在ggplot中绘制混合效应模型

Question

我是混合效果模型的新手，需要您的帮助。我在 ggplot 中绘制了下图：

ggplot(tempEf,aes(TRTYEAR,CO2effect,group=Myc,col=Myc)) + 
  facet_grid(~N) +
  geom_smooth(method="lm",se=T,size=1) +
  geom_point(alpha = 0.3) + 
  geom_hline(yintercept=0, linetype="dashed") +
  theme_bw()

但是，我想在 geom_smooth 中表示混合效应模型而不是 lm，因此我可以将 SITE 作为随机效应包括在内。

模型如下：

library(lme4)
tempEf$TRTYEAR <- as.numeric(tempEf$TRTYEAR)
mod <- lmer(r ~ Myc * N * TRTYEAR + (1|SITE), data=tempEf)

我已经包括了 TRTYEAR（治疗年份），因为我也对效果的模式感兴趣，对于某些组来说，效果可能会随着时间的推移而增加或减少。

下一步是我迄今为止从模型中提取绘图变量的最佳尝试，但只提取了 TRTYEAR= 5、10 和 15 的值。

library(effects)
ef <- effect("Myc:N:TRTYEAR", mod)
x <- as.data.frame(ef)
> x
   Myc     N TRTYEAR        fit         se       lower     upper
1   AM  Nlow       5 0.04100963 0.04049789 -0.03854476 0.1205640
2  ECM  Nlow       5 0.41727928 0.07342289  0.27304676 0.5615118
3   AM Nhigh       5 0.20562700 0.04060572  0.12586080 0.2853932
4  ECM Nhigh       5 0.24754017 0.27647151 -0.29556267 0.7906430
5   AM  Nlow      10 0.08913042 0.03751783  0.01543008 0.1628307
6  ECM  Nlow      10 0.42211957 0.15631679  0.11504963 0.7291895
7   AM Nhigh      10 0.30411129 0.03615213  0.23309376 0.3751288
8  ECM Nhigh      10 0.29540744 0.76966410 -1.21652689 1.8073418
9   AM  Nlow      15 0.13725120 0.06325159  0.01299927 0.2615031
10 ECM  Nlow      15 0.42695986 0.27301163 -0.10934636 0.9632661
11  AM Nhigh      15 0.40259559 0.05990085  0.28492587 0.5202653
12 ECM Nhigh      15 0.34327471 1.29676632 -2.20410343 2.8906529

欢迎提出完全不同的方法来表示此分析的建议。我认为这个问题更适合 Whosebug，因为它是关于 R 中的技术细节而不是背后的统计数据。谢谢

Answer 1

您可以用多种不同的方式表示您的模型。最简单的方法是使用不同的绘图工具（颜色、形状、线型、面）根据各种参数绘制数据，这就是您对示例所做的，除了随机效应 site。还可以绘制模型残差以传达结果。就像@MrFlick 评论的那样，这取决于您想传达的内容。如果您想在估计值周围添加 confidence/prediction 个区间，则必须更深入地挖掘并考虑更大的统计问题 (example1, example2)。

这里有一个例子让你更进一步。
此外，在您的评论中，您说您没有提供可重现的示例，因为数据不属于您。这并不意味着您不能提供虚构数据的示例。请在以后的帖子中考虑这一点，以便您可以更快地得到答案。

#Make up data:
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
  )

# Make up some response data
tempEf$r <- 2*tempEf$TRTYEAR +                   
            -8*as.numeric(tempEf$Myc=="ECM") +
            4*as.numeric(tempEf$N=="Nlow") +
            0.1*tempEf$TRTYEAR * as.numeric(tempEf$N=="Nlow") +
            0.2*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM") +
           -11*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
            0.5*tempEf$TRTYEAR*as.numeric(tempEf$Myc=="ECM")*as.numeric(tempEf$N=="Nlow")+ 
           as.numeric(tempEf$site) +  #Random intercepts; intercepts will increase by 1
           tempEf$TRTYEAR/10*rnorm(600, mean=0, sd=2)    #Add some noise

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

顺便说一句，与上面的系数相比，该模型很好地拟合了数据：

model

#Linear mixed model fit by REML ['lmerMod']
#Formula: r ~ Myc * N * TRTYEAR + (1 | site)
#   Data: tempEf
#REML criterion at convergence: 2461.705
#Random effects:
# Groups   Name        Std.Dev.
# site     (Intercept) 1.684   
# Residual             1.825   
#Number of obs: 600, groups:  site, 5
#Fixed Effects:
#         (Intercept)                MycECM                 NNlow               
#             3.03411              -7.92755               4.30380               
#             TRTYEAR          MycECM:NNlow        MycECM:TRTYEAR  
#             1.98889             -11.64218               0.18589  
#       NNlow:TRTYEAR  MycECM:NNlow:TRTYEAR  
#             0.07781               0.60224

调整您的示例以显示覆盖在数据上的模型输出

   library(ggplot2)
    ggplot(tempEf,aes(TRTYEAR, r, group=interaction(site, Myc), col=site, shape=Myc )) + 
      facet_grid(~N) +
      geom_line(aes(y=fit, lty=Myc), size=0.8) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

请注意，我所做的只是将您的颜色从 Myc 更改为 site，并将线型更改为 Myc.

我希望这个例子能提供一些关于如何可视化混合效果模型的想法。

在ggplot中绘制混合效应模型

plot mixed effects model in ggplot

r

ggplot2

lmer