如何为lmer模型结果绘制预测值和标准误差的图表？

Question

我有一个移植实验，涉及四个地点和四种基质（每个地点取自不同基质）。我已经确定了每个种群在每个地点和基质组合中的存活情况。该实验已被复制三次。

我创建了一个lmm，如下所示：

Survival.model <- lmer(Survival ~ Location + Substrate + Location:Substrate + (1|Replicate), data=Transplant.Survival,, REML = TRUE)

我想使用predict命令提取预测结果，例如：

Survival.pred <- predict(Survival.model)

然后提取标准误差，以便我可以将它们与预测一起绘制，生成类似以下图表的内容：

我知道如何使用标准glm来完成这个任务（这也是我创建示例图的方式），但我不确定是否可以或应该使用lmm。

作为线性混合模型的新用户，我是否缺少了一些基础知识？

我在Stack Overflow上找到了post，但并没有得到帮助。

根据RHertel的评论，也许我应该这样表达问题：如何绘制我的lmer模型结果的模型估计值和置信区间，以便我可以获得与上面创建的类似的图形？

示例数据：

Transplant.Survival <- structure(list(Location = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("Steninge", "Molle", 
"Kampinge", "Kaseberga"), class = "factor"), Substrate = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 3L, 3L, 4L, 4L, 4L, 1L, 1L, 2L, 2L, 2L, 3L, 
3L, 3L, 4L, 4L, 4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 
4L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L), .Label = c("Steninge", 
"Molle", "Kampinge", "Kaseberga"), class = "factor"), Replicate = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1", 
"2", "3"), class = "factor"), Survival = c(1, 1, 1, 0.633333333333333, 
0.966666666666667, 0.5, 0.3, 0.233333333333333, 0.433333333333333, 
0.966666666666667, 0.866666666666667, 0.5, 0.6, 0.266666666666667, 
0.733333333333333, 0.6, 0.3, 0.5, 0.3, 0.633333333333333, 0.9, 
0.266666666666667, 0.633333333333333, 0.7, 0.633333333333333, 
0.833333333333333, 0.9, 0.6, 0.166666666666667, 0.333333333333333, 
0.433333333333333, 0.6, 0.9, 0.6, 0.133333333333333, 0.566666666666667, 
0.633333333333333, 0.633333333333333, 0.766666666666667, 0.766666666666667, 
0.0333333333333333, 0.733333333333333, 0.3, 1.03333333333333, 
0.6, 1)), .Names = c("Location", "Substrate", "Replicate", "Survival"
), class = "data.frame", row.names = c(NA, -46L))

- Keith W. Larson

---

RHertel：非常抱歉问这个愚蠢的问题，但对于这样一个小型研究，将数据分成训练集和测试集是否有意义？如果这样做的话，用置信区间绘制估计值是否更合适？ - Keith W. Larson

阅读r-sig-mixed-models FAQ。其中包含一个例子。此外，还有我的这个答案，它指向了一个讨论该问题更深层次概念问题的邮件列表线程。 - Roland

@Roland：感谢您的回复。我已经阅读了这些页面并理解了方法，但是由于数据集太小，将其分成训练和测试数据集似乎有问题。我需要重新考虑如何在图表中呈现模型结果或整个建模方法！ - Keith W. Larson

问题在于SE（标准误差）对于MLM（多层次线性模型）并不是标准的，这就是为什么lme4没有为它们实现函数的原因。您只需要前端标准误差吗？后端标准误差？还是组合标准误差？此外，您的方差矩阵结构也很重要。 - alexwhitworth

@KeithLarson，MLM FAQ中已经提供了一个示例，您可以参考。任何标准文本（Goldstein、Pinheiro和Bates等）都会有详细说明。 - alexwhitworth

Answer 1

编辑：修复函数/图中的错误。

如果您想使用置信区间绘制估计值，可以看一下sjPlot包中的sjp.lmer函数。在这里可以查看各种绘图类型的示例此处。

此外，arm包提供了用于计算标准误差的函数（arm::se.fixef和arm::se.ranef）。

sjp.setTheme("forestgrey") # plot theme
sjp.lmer(Survival.model, type = "fe")

会给出以下图表：

- Daniel

---

刚刚发现了“type =“ fe””的x轴标签错误，会进行修复。 - Daniel

好的，sjPlot包已在CRAN上更新，同时修复了错误。 - Daniel