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

Question

在类似于以下的线性混合模型中，我如何提取每个实验单位（plot）的系数（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

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

- Juanchi

Answer 1

首先要说明的是，这实际上是一个不太简单的理论问题：在r-sig-mixed-models上有一个相当长的线程，其中涉及了一些技术细节；你应该看一下，尽管可能会有点吓人。基本问题在于，每个组的估计系数值是固定效应参数和该组的BLUP/条件模式之和，它们是不同类别的对象（一个是参数，一个是随机变量的条件均值），这就产生了一些技术难题。

其次，（不幸的是）我不知道在lme中有什么简单的方法来解决这个问题，所以我的答案使用了lmer（从lme4包中获得）。

如果你愿意做最简单的事情，并忽略固定效应参数和BLUP之间（可能存在的）不良协方差，你可以使用下面的代码。

两种选择：（1）使用贝叶斯分层方法（例如MCMCglmm包）拟合模型，并计算每个级别的后验预测标准差；（2）使用参数自助法计算BLUPs/条件模式，然后取自助分布的标准差。请记住，如往常一样，此建议不保证正确性。

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

- Ben Bolker

---

两个问题：1）每个分组水平的标准误是否总是相同的，或者可能存在标准误值因分组水平而异的模型？换句话说，仅返回拦截项和斜率的一个标准误值是否足够？2）这也适用于 glmer 吗？ - Daniel

Answer 2

使用nlme，您可以部分地完成以下操作...

您可以使用以下方法提取summary()的组件：

summary(fitL1)$tTable[,1] #fixed-effect parameter estimates
summary(fitL1)$tTable[,2] #fixed-effect parameter standard errors

你可以通过行进一步分割这些数据集：

summary(fitL1)$tTable[1,1] #the first fixed-effect parameter estimate
summary(fitL1)$tTable[1,2] #the first fixed-effect parameter standard error

例如，使用以下代码从结果中提取单个参数或标准误差并将其组合成数据框：

df<-data.frame(cbind(summary(fitL1)$tTable[1,1], summary(fitL1)$tTable[1,2]))
names(df)<-c("Estimate","SE")
df

要为每个图形调整这些参数（我想是随机效应），您可以使用以下方法提取随机系数：

fitL1$coefficients$random

将它们添加到参数估计值（B0（截距），B1等）中。然而，我不确定每个图的标准误差应该如何调整。

- AKR