Lmer 输出的 Ggplot 可视化

Question

我正在尝试使用 ggplot 可视化 lmer 输出，但无法得到我想要的结果。

代码如下：

model <- lmer(Fluency ~ Roundnr * relevel(as.factor(Condition), ref='C3') * block + (1|Participant), data=data_aggregate, REML=FALSE)

ggplot(data_aggregate,aes(Roundnr, Fluency, group=interaction(Participant, Condition), col=Condition)) + 
  facet_grid(~block, scales="free", space="free_x") +
  geom_line(aes(y=predict(model), lty=Condition), linewidth=0.8) +
  geom_point(alpha = 0.3) + 
  geom_hline(yintercept=0, linetype="dashed")

结果如下所示：

（代码改编自这个问题的答案：在ggplot中绘制混合效应模型）

但是，我想为每个条件显示 1 行，而不是为每个参与者显示一行。我无法真正分享完整的数据，但这是大约 60 名参与者的数据，他们每人执行一项任务 8 次（A 中 4 次，B 中 4 次）。更改它使得 group=Condition 不起作用，它会给出以下奇怪的结果：

对于如何实现这一目标有什么建议吗？

Answer 1

如果您想在人口级别绘制预测，最好在人口级别自行生成预测（见下文）：这也可以使用各种包

emmeans

、

ggeffects

，也许

sjPlot

更自动地实现。 ..

这个示例与您的示例不太匹配，但应该很容易适应。

library(lme4)
model <- lmer(r ~ Myc * N * TRTYEAR + (1|site), data=tempEf)
tempEf$fit <- predict(model)   #Add model fits to dataframe


library(ggplot2)
gg0 <- ggplot(tempEf,aes(TRTYEAR, r, colour = Myc)) +
    facet_grid(Myc~N, scales = "free") +
    geom_point(alpha = 0.3) + 
    geom_hline(yintercept=0, linetype="dashed") +
    geom_line(aes(y=fit, group = interaction(site,Myc)), alpha = 0.5)


## generate predictions at the population level
pp <- with(tempEf,
   ## n=51 is overkill for this example but will be useful if you have
   ## something nonlinear (confidence intervals, GLMM responses, etc.)
           expand.grid(TRTYEAR = seq(min(TRTYEAR), max(TRTYEAR), length.out = 51),
                       Myc = unique(Myc),
                       N = unique(N)))

pp$r <- predict(model, newdata = pp, re.form = NA)
gg0 + geom_line(data = pp, lwd = 2)

模拟数据

## https://stackoverflow.com/questions/31075407/plot-mixed-effects-model-in-ggplot
set.seed(101)
tempEf <- data.frame(
  N = rep(c("Nlow", "Nhigh"), each=300),
  Myc = rep(c("AM", "ECM"), each=150, times=2),
  TRTYEAR = runif(600, 1, 15),
  site = rep(c("A","B","C","D","E"), each=10, times=12)   #5 sites
)
tempEf$r <- simulate( ~ Myc * N * TRTYEAR + (1|site),
                     newdata = tempEf,
                     newparams = list(beta = c(2, -2, 1, 0.1, 0.2, -1, 0.5, 0),
                                      theta = 1,
                                      sigma = 1))[[1]]