【问题标题】:R: how to extract information from summary model fitR:如何从摘要模型拟合中提取信息
【发布时间】: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 拟合的摘要输出中提取信息。

我想提取

  1. 随机效应的 StdDev(即 Asym 的 StdDev,这 = 3.65)对于这个我试过 fm1$apVar 但没有运气。
  2. 固定效应的参数估计值(即 Asym = 101.44960,R0 = -8.62733 等),可通过fixef(fm1) 提取
  3. 固定效果的Std.Error(即2.46、0.317、0.034)。对于这个我已经尝试过sqrt(diag(fm1$varFix)),但这些值与固定效果下 Std.Error 列下的值不完全匹配?
  4. logLikelihood(即-114.7428,可以使用fm1$logLik提取)
  5. 残差(即 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 包中的 tidyglance 来获得所需的大部分模型估计值。 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

标签: r nlme


【解决方案1】:

这里还有一些...

as.numeric(VarCorr(fm1)[,2])
# [1] 3.6506418 0.7188625

summary(fm1)$tTable[,2]
#       Asym         R0        lrc 
# 2.46169512 0.31795045 0.03427017 

# looks like you don't need this one anymore, but here's a way of getting it
summary(fm1)$corFixed
#            Asym         R0        lrc
# Asym  1.0000000  0.7039498 -0.9077793
# R0    0.7039498  1.0000000 -0.8271022
# lrc  -0.9077793 -0.8271022  1.0000000

抱歉,这不是一个完整的答案 - 创建一个像您描述的那样的汇总表可能很困难,因为每个潜在行的结构都会不同,并且取决于哪些变量被包含为固定变量和随机效应。

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-04-19
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多