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

【问题标题】：Extracting coefficients and their standard error for each unit in an lme model fit提取 lme 模型拟合中每个单元的系数及其标准误差
【发布时间】：2014-11-29 16:38:14
【问题描述】：

如何在线性混合模型中提取系数（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

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

【问题讨论】：

标签： r mixed-models nlme coefficients

【解决方案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

【讨论】：

优秀的@Ben！这就是我要找的……我会去那个网站看看。
两个问题：1）每个分组级别的 SE 是否始终相同，或者是否存在每个分组级别的 SE 值不同的模型？（换句话说：只返回一个 SE 值作为截距和一个斜率是否足够）2）这是否也适用于glmer？

【解决方案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 等）。但是，我不确定应该如何调整每个图的标准误差。

【讨论】：