如何在混合效应模型中获取系数及其置信区间？

Question

在 lm 和 glm 模型中，我使用函数 coef 和 confint 来达到目标:

m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
coef(m)
confint(m)

我现在将随机效应添加到模型中，使用lme4包中的lmer函数进行混合效应建模。但是，coef和confint函数不再起作用！

> mix1 = lmer(resp ~ 0 + var1 + var1:var2 + (1|var3)) 
                                      # var1, var3 categorical, var2 continuous
> coef(mix1)
Error in coef(mix1) : unable to align random and fixed effects
> confint(mix1)
Error: $ operator not defined for this S4 class

我尝试过谷歌和文档，但没有结果。请指点我正确的方向。

编辑：我也在思考这个问题是否更适合https://stats.stackexchange.com/，但我认为它比统计学更具技术性，所以我得出结论最适合在这里提问（SO）......你怎么看？

- Tomas

---

为了帮助你入门，直到像@BenBolker这样的专家出现：?lmer列出了方法fixef和ranef，除了coef。由于您的错误显示它在尝试组合两者时遇到了问题，所以问题很可能是您的模型规范在某种程度上是“不寻常的”。 - joran

谢谢@joran。我的模型规范可能不太寻常，因为我想省略截距 - 我想这样做，因为否则系数是无意义的。var1是分类变量，我想为每个类别设置“组特定截距”。如果我允许截距（从公式中删除0 +），coef可以运行但不会给出我期望的结果。fixef非常好用，谢谢！然而，confint根本不起作用。 - Tomas

我会直接从S4对象中提取您需要的数据--请参阅此帖子的答案：https://dev59.com/jWoy5IYBdhLWcg3wdt8L - baha-kev

谢谢@baha-kev，但是你确定置信区间在这个对象里吗？我觉得不是... - Tomas

不是直接的，但你只需要将标准误差乘以+/-1.96并加到系数估计值上，就可以得到一个95%置信区间。例如，S4对象具有标准误差和系数估计值。 - baha-kev

我正在修复lme4的r-forge版本（lme4.0，目前稳定的分支对应于CRAN-lme4）中coef中的bug(let)问题，以及开发分支中的lme4。confint是一个更大的问题，正如已经讨论过的那样，尽管lme4的开发分支可以计算轮廓置信区间... - Ben Bolker

Answer 1

不确定是何时添加的，但现在lme4中已经实现了confint()。例如，以下示例可正常工作：

library(lme4)
m = lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
confint(m)

- jciloa

Answer 2

有两个新的包lmerTest和lsmeans可以计算lmer和glmer输出的95％置信区间。也许你可以研究一下这些？而coefplot2，我认为也可以（尽管正如下面的Ben指出的那样，采用的是Wald统计量的标准误差，而不是lmerTest和lsmeans中使用的Kenward-Roger和/或Satterthwaite df逼近的方法）......只是遗憾的是，在包lsmeans中仍然没有内置的绘图功能（就像在包effects()中一样，它还返回lmer和glmer对象的95％置信区间，但是通过重新拟合一个没有任何随机因素的模型来实现，显然是不正确的）。

- Tom Wenseleers

---

coefplot2非常天真地计算了1.96倍的Wald标准误差--它没有解决置信区间有限大小修正的非常重要的问题。 - Ben Bolker

请查看此帖子 http://stats.stackexchange.com/questions/117641/how-trustworthy-are-the-confidence-intervals-for-lmer-objects-through-effects-pa 以获取更详细的答案。 - Tom Wenseleers

lmerTest现在在JoSS中得到了很好的描述 https://www.jstatsoft.org/article/view/v082i13 - radek

请注意，这些评论现在很过时。使用emmeans或lmerTest是正确的方法，并且现在有绘图方法可用。 - Axeman

Answer 3

我来补充一下。如果m是已拟合的(g)lmer模型（大部分也适用于lme）：

• fixef(m)是从混合模型中提取系数的规范方法（此约定始于nlme并延续到lme4）
• 您可以使用coef(summary(m))获取完整的系数表；如果在拟合模型之前加载了lmerTest，或者在拟合模型之后转换模型（然后加载lmerTest），则系数表将包括p值。 （系数表是矩阵；您可以通过例如ctab [，“Estimate”]、ctab[，“Pr(>|t|)”]提取列，或者将矩阵转换为数据框并使用$-索引。）
• 如上所述，您可以通过confint(m)获取似然轮廓置信区间；这可能需要计算较多的时间。 如果您使用confint(m, method="Wald")，则会获得标准的+/- 1.96SE置信区间。 (lme使用intervals(m)而不是confint()。)

如果您喜欢使用broom.mixed：

• tidy(m,effects="fixed")会给出一个包含估计值、标准误等的表格。
• tidy(as(m,"merModLmerTest"), effects="fixed") （或者一开始就用lmerTest进行拟合）包括p值
• 添加conf.int=TRUE会给出（Wald）置信区间
• 添加conf.method="profile"（以及conf.int=TRUE）会给出似然轮廓置信区间

您还可以通过参数Bootstrap（method="boot"）来获得置信区间，这在某些情况下会更加准确，但速度较慢。

- Ben Bolker

---

嗨Ben，谢谢！我有点困惑，独立的点“.”是什么意思？如果它只是一个模型变量名，为什么不使用例如“m”？ :-) - Tomas

我可以使用 m。有时我会使用 . 作为占位符。 - Ben Bolker

Answer 4

假设对固定效应采用正常逼近（confint 也是这样做的），我们可以通过 估计值 + 1.96 * 标准误差 来获得95%的置信区间。

以下不适用于方差分量/随机效应。

library("lme4")
mylm <- lmer(Reaction ~ Days + (Days|Subject),  data =sleepstudy)

# standard error of coefficient

days_se <- sqrt(diag(vcov(mylm)))[2]

# estimated coefficient

days_coef <- fixef(mylm)[2]

upperCI <-  days_coef + 1.96*days_se
lowerCI <-  days_coef  - 1.96*days_se

- julieth

---

你好julieth，很不错的想法，然而真正的置信区间（由confint计算）与这些区间存在差异...。也许t分布会给出与confint相同的结果（不确定），但在这种情况下我不知道应该使用哪个df。 - Tomas

换句话说，这就是为什么我更喜欢使用像confint等函数来为我完成所有这些操作的原因...（特别是如果我不确定系数的正态分布）。 - Tomas

t分布渐近于正态分布，而在许多多层次设计中误差项的自由度非常高，因此在该点处误差分布是正态的。因此，如果您有一个具有大量自由度的设计，则这是一个完全合理的置信区间估计。 - John

Answer 5

我建议您使用已有的lme软件包（在nlme中）。它有confint函数，如果您需要对比度的置信区间，可选择一系列方案（在gmodels中为estimable，在contrasts中为contrast，在multcomp中为glht）。

为什么lmer中没有p值和置信区间：请参见http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html 。

- Dieter Menne

---

谢谢Dieter，我会尝试使用旧版本的软件包。缺少p值并且无法立即确定显著性也让我感到担忧！如果我能够获得置信区间，那么我只需要查看是否包含零，就可以得出显著性结论了！祝好！ - Tomas

我忘了提到multcomp包中的confint(glht...)可以为lmer提供渐近置信区间。Douglas Bates的警告仍然适用，但他大胆地将lmer / gaussian中的p值省略出去无疑激起了轩然大波。 - Dieter Menne

Dieter，你说的“confint(glht”是什么意思？multcomp包中没有confint函数... - Tomas

使用 intervals(mix1) 将会给出类似于 @julieth 的下面回答中的渐近置信区间；intervals(mix1)$fixed 提取固定效应区间。这些基于正态近似，而不是 t 分布或任何更奇特的东西... - Ben Bolker

你对于默认表是正确的；我想的是非默认对比。 - Dieter Menne

Answer 6

我建议使用sjPlot包中的tab_model()函数作为替代方案。它可以生成干净易读的输出，适用于markdown。参考此处和示例此处。

对于那些更注重视觉效果的人，同一包中的plot_model()也可能很有用。

另一种选择是使用parameters包，使用model_parameters()函数function。 package

- radek

---

或者 broom.mixed::tidy() - Ben Bolker

在我看来，@Ben Bolker的回答应该是最高评分的答案 :) - Extrapolator

Answer 7

要找到系数，您可以简单地使用lme4的summary函数

m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
m_summary <- summary(m)

要求所有系数：

m_summary$coefficient

如果您想要置信区间，请将标准误差乘以1.96：

CI <- m_summary$coefficient[,"Std. Error"]*1.96
print(CI)

- denis

---

这里的1.96因子是用于95%的置信区间。 - denis