6
votes
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.

1
Create a function that returns a named vector of all the stats you want to return for a single model and then use sapply to run it over a list of the models. - lmo
@lmo That sounds good. Do you have any idea how I can extract the std. errors of the fixed effects? - Adrian
You can use tidy and glance from package broom to get most of of the model estimates you want. An alternative to tidy is to simply pull the tTable out of the summary output and extract what you want (like SE) from there; e.g., summary(fm1)$tTable . - aosmith
The log-likelihood and residuals both have extractor functions: logLik & residuals. You can match the fixed effects st. errors returned by sqrt(diag(vcov(fm1))) (or sqrt(diag(fm1$varFix), but again better to use extractor functions) to the summary output if you dont adjust them summary(fm1, adjustSigma = FALSE) or fit your model using method="REML" - see ?nlme:::summary.lme - user20650

1 Answers

0
votes

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.