Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: how to extract information from summary model fit

Tags:

r

nlme

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 

I'm interested in extracting information from the summary output of an NLME fit.

I would like to extract the

  1. StdDev of the random effects (i.e. StdDev of Asym, which = 3.65) For this one I've tried fm1$apVar but no luck.
  2. Parameter estimates of the fixed effects (i.e. Asym = 101.44960, R0 = -8.62733, etc), which can be extracted via fixef(fm1)
  3. Std.Error of the fixed effects (i.e. 2.46, 0.317, 0.034). For this one I've tried sqrt(diag(fm1$varFix)) but those values don't exactly match the ones under Std.Error column under fixed effects?
  4. logLikelihood (i.e. -114.7428, which can be extracted using fm1$logLik)
  5. Residuals (i.e. 0.7188625, which can be extracted using fm1$Residuals)

My ultimate goal is to fit multiple models and store their respective summary estimates into an organized 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)

Suppose my list object looks like mylist. Now I would like to create a data.frame that looks like:

model    FixedAsym    FixedAsymStdError   FixedR0      ...     Residual
 1       101.44960        2.4616951       -8.62733            0.7188625
 2       101.44934        2.4616788       -8.62736     ...    0.7188625

And to create this data.frame (the number of rows corresponds to how many model summaries I have in mylist) I would need to systematically extract those values (numbered 1-5) from the model summary output.

like image 527
Adrian Avatar asked Jun 21 '17 16:06

Adrian


1 Answers

Here are a few more pieces...

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

Apologies that this isn't a complete answer -- It may prove difficult to create a summary table like you're describing, since the structure of each potential row will be different, and will depend on which variables are included as fixed and random effects.

like image 76
Matt Tyers Avatar answered Nov 29 '22 20:11

Matt Tyers