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

问题描述 投票:7回答:2

如何在线性混合模型中提取每个实验单位的系数(b0和b1)及其各自的标准误差(图)。

更好地拟合线性模型

对于相同的数据集(df)和拟合模型(fitL1):我怎么能得到一个像这样的数据框......

   plot    b0      b0_se   b1    b1_se 
    1    2898.69   53.85   -7.5  4.3

   ...    ...       ...     ...   ...
r mixed-models nlme coefficients
2个回答
8
投票

首先要说的是,这其实是一个非平凡的理论问题:有一个相当的 关于R-Sig-混合模型的长线。 这本书详细介绍了一些技术细节,你一定要看看,尽管它有点吓人。基本的问题是,每组的估计系数值是该组的固定效应参数和BLUPconditional模式之和,它们是不同类别的对象(一个是参数,一个是随机变量的条件均值),这就造成了一些技术上的困难。

第二点是,(很遗憾)我不知道有什么简单的方法可以在 lme所以我的答案使用 lmer (来自 lme4 包)。)

如果 您可以用最简单的方法忽略固定效应参数和BLUPs之间的(可能是定义不清的)协方差,您可以使用下面的代码。

两种选择是:(1)用贝叶斯分层方法来拟合你的模型(如 MCMCglmm 包),并计算每个层次的后验预测的标准差(2)使用参数引导计算BLUPsconditional模式,然后取引导分布的标准差。

请记住,和往常一样,这个建议是没有保证的。

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

4
投票

为了让你使用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等)。但是,我不知道应该如何调整每个小区的标准误差。

© www.soinside.com 2019 - 2024. All rights reserved.