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

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）...你怎么看认为？

Answer 1

不确定何时添加，但现在 lme4 中实现了confint()。 例如，以下示例有效：

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

Answer 2

有两个包，lmerTest 和 emmeans，可以计算

lmer

和

glmer

输出的 95% 置信限。也许你可以看看那些？而且 coefplot2，我认为也可以做到（尽管正如 Ben 在下面指出的那样，以一种不太复杂的方式，来自 Wald 统计数据的标准误差，而不是

 中使用的 Kenward-Roger 和/或 Satterthwaite df 近似） lmerTest

和

emmeans

）...

Answer 3

我要在这里补充一点。如果

m

是拟合的

(g)lmer

模型（其中大部分也适用于

lme

）：

• fixef(m)

  是从混合模型中提取系数的规范方法（该约定从

  nlme

  开始，一直延续到

  lme4

  ）
• 您可以通过

  coef(summary(m))

  获得完整的系数表；如果您在拟合模型之前加载了

  lmerTest

  ，或者在拟合后通过

  lmerTest

  转换模型（然后加载

  coef(summary(as(m,"merModLmerTest")))

  ），则系数表将包含 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) CI
• 添加

  conf.method="profile"

  （与

  conf.int=TRUE

  一起）给出似然剖面 CI

您还可以通过参数引导程序 (

method="boot"

) 获得置信区间，这在某些情况下速度要慢得多，但更准确。

Answer 4

假设固定效应的正态近似（限制也可以做到），我们可以通过

获得 95% 的置信区间

估计 + 1.96*标准误差。

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

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

Answer 5

我建议你使用好的旧lme（在nlme包中）。它有限制，如果你需要对比限制，有一系列选择（gmodels中的esimable，contrasts中的contrast，multcomp中的glht）。

为什么 lmer 中缺少 p 值和限制：请参阅 http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html .

Answer 6

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

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)

Answer 7

我建议使用

tab_model()

包中的

sjPlot

函数作为替代方案。干净且可读的输出可供降价使用。参考这里和示例这里。

对于那些更具视觉倾向的人来说，同一个包装中的

plot_model()

也可能会派上用场。

替代解决方案是通过

parameters

package 使用

model_parameters()

function。