【发布时间】:2017-06-21 16:47:48
【问题描述】:
library(nlme)
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
> summary(fm1)
Nonlinear mixed-effects model fit by maximum likelihood
Model: height ~ SSasymp(age, Asym, R0, lrc)
Data: Loblolly
AIC BIC logLik
239.4856 251.6397 -114.7428
Random effects:
Formula: Asym ~ 1 | Seed
Asym Residual
StdDev: 3.650642 0.7188625
Fixed effects: Asym + R0 + lrc ~ 1
Value Std.Error DF t-value p-value
Asym 101.44960 2.4616951 68 41.21128 0
R0 -8.62733 0.3179505 68 -27.13420 0
lrc -3.23375 0.0342702 68 -94.36052 0
Correlation:
Asym R0
R0 0.704
lrc -0.908 -0.827
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.23601930 -0.62380854 0.05917466 0.65727206 1.95794425
Number of Observations: 84
Number of Groups: 14
我有兴趣从 NLME 拟合的摘要输出中提取信息。
我想提取
- 随机效应的 StdDev(即 Asym 的 StdDev,这 = 3.65)对于这个我试过
fm1$apVar但没有运气。 - 固定效应的参数估计值(即 Asym = 101.44960,R0 = -8.62733 等),可通过
fixef(fm1)提取 - 固定效果的Std.Error(即2.46、0.317、0.034)。对于这个我已经尝试过
sqrt(diag(fm1$varFix)),但这些值与固定效果下 Std.Error 列下的值不完全匹配? - logLikelihood(即-114.7428,可以使用
fm1$logLik提取) - 残差(即 0.7188625,可以使用
fm1$Residuals提取)
我的最终目标是拟合多个模型并将它们各自的汇总估计存储到一个有组织的data.frame。
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
fm2 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
data = Loblolly,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1,
start = c(Asym = 103, R0 = -5.4, lrc = -3.3))
summary(fm1)
summary(fm2)
mylist = list(NULL, summary(fm1), NULL, summary(fm2), NULL, NULL)
假设我的列表对象看起来像mylist。现在我想创建一个data.frame,看起来像:
model FixedAsym FixedAsymStdError FixedR0 ... Residual
1 101.44960 2.4616951 -8.62733 0.7188625
2 101.44934 2.4616788 -8.62736 ... 0.7188625
为了创建这个 data.frame(行数对应于我在 mylist 中的模型摘要数量),我需要从模型摘要输出中系统地提取这些值(编号 1-5)。
【问题讨论】:
-
创建一个函数,该函数返回您想要为单个模型返回的所有统计信息的命名向量,然后使用
sapply在模型列表上运行它。 -
@lmo 听起来不错。你知道我如何提取标准。固定效应的误差?
-
您可以使用 broom 包中的
tidy和glance来获得所需的大部分模型估计值。tidy的替代方法是简单地将tTable从摘要输出中提取出来,然后从那里提取您想要的(例如 SE);例如,summary(fm1)$tTable. -
对数似然和残差都有提取函数:
logLik&residuals。您可以匹配固定效果 st。如果您不调整summary(fm1, adjustSigma = FALSE)或使用method="REML"拟合您的模型,则sqrt(diag(vcov(fm1)))(或sqrt(diag(fm1$varFix),但最好使用提取器函数)返回到摘要输出的错误 - 请参阅?nlme:::summary.lme