r - 在 lme 模型拟合中提取每个单元的系数及其标准误差

Question

我如何在线性混合模型中提取系数(b0 和 b1)及其各自的标准误差(图)，例如:

Better fits for a linear model

使用相同的数据集(df)，以及拟合模型(fitL1):我怎么能得到这样的数据框......

   plot    b0      b0_se   b1    b1_se 
    1    2898.69   53.85   -7.5  4.3

   ...    ...       ...     ...   ...

Answer 1

第一个评论是，这实际上是一个不平凡的理论问题:有一个相当long thread on r-sig-mixed-models这涉及到一些技术细节；你绝对应该看看，即使它有点吓人。基本问题是每个组的估计系数值是该组的固定效应参数和 BLUP/条件模式的总和，它们是不同类别的对象(一个是参数，一个是条件均值)随机变量)，这造成了一些技术困难。

第二点是(不幸的是)我不知道在 lme 中有什么简单的方法可以做到这一点。 ，所以我的回答使用 lmer (来自 lme4 包)。

如 您可以轻松地做最简单的事情并忽略固定效果参数和 BLUP 之间的(可能定义不明确的)协方差，您可以使用下面的代码。

两种选择是 (1) 使用贝叶斯分层方法(例如 MCMCglmm 包)拟合您的模型并计算每个级别的后验预测的标准偏差 (2) 使用参数自举来计算 BLUP/条件模式，然后取自举分布的标准偏差。

请记住，像往常一样，此建议不提供任何保证。

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
cc <- coef(fm1)$Subject
## variances of fixed effects
fixed.vars <- diag(vcov(fm1))
## extract variances of conditional modes
r1 <- ranef(fm1,condVar=TRUE)
cmode.vars <- t(apply(cv <- attr(r1[[1]],"postVar"),3,diag))
seVals <- sqrt(sweep(cmode.vars,2,fixed.vars,"+"))
res <- cbind(cc,seVals)
res2 <- setNames(res[,c(1,3,2,4)],
                 c("int","int_se","slope","slope_se"))
##          int   int_se     slope slope_se
## 308 253.6637 13.86649 19.666258   2.7752
## 309 211.0065 13.86649  1.847583   2.7752
## 310 212.4449 13.86649  5.018406   2.7752
## 330 275.0956 13.86649  5.652955   2.7752
## 331 273.6653 13.86649  7.397391   2.7752
## 332 260.4446 13.86649 10.195115   2.7752