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

Question

在lm和glm模型中,我使用函数coef并confint实现目标:

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

现在我添加了随机效果到模型 - 使用lmerlme4包的函数使用混合效果模型.但是,功能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

有两个新的包,lmerTest和lsmeans,可以计算95%的置信限lmer和glmer输出.也许你可以看看那些？而coefplot2,我认为也可以这样做(尽管Ben在下面指出,以不太复杂的方式,从Wald统计数据的标准误差,而不是Kenward-Roger和/或Satterthwaite df近似用于lmerTest和lsmeans). ..只是一个耻辱,还有在包中没有内置的绘图设备lsmeans(因为有封装effects(),其中顺便说一句也返回上95%置信区间lmer和glmer对象,但通过改装的模型没有任何随机的因素确实如此,这显然是不正确).

Answer 2

我要在这里补充一点。如果m是一个拟合(g)lmer模型（其中大多数lme也适用）：

• fixef(m)是从混合模型中提取系数的规范方法（该约定始于nlme并延续到lme4）
• 您可以获得完整的系数表coef(summary(m))；如果您在拟合模型之前已加载lmerTest，或通过 拟合后转换模型（然后加载lmerTest）coef(summary(as(m,"merModLmerTest")))，则系数表将包含 p 值。（系数表是一个矩阵；您可以通过例如ctab[,"Estimate"]、来提取列ctab[,"Pr(>|t|)"]，或者将矩阵转换为数据框并使用$-indexing。）
• 如上所述，您可以通过以下方式获得可能性剖面置信区间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 3

我建议你使用好的老lme(在包装nlme).它有一些特权,如果你需要对比,那就有一系列的选择(在gmodels中可以估算,对比中的对比,在multcomp中的glht).

为什么在lmer中不存在p值和confint:请参阅http://finzi.psych.upenn.edu/R/Rhelp02a/archive/76742.html.

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

不确定它何时添加,但现在confint()在lme4中实现.例如,以下示例有效:

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