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
- StdDev of the random effects (i.e. StdDev of Asym, which = 3.65) For this one I've tried
fm1$apVarbut no luck. - Parameter estimates of the fixed effects (i.e. Asym = 101.44960, R0 = -8.62733, etc), which can be extracted via
fixef(fm1) - 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? - logLikelihood (i.e. -114.7428, which can be extracted using
fm1$logLik) - 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.
sapplyto run it over a list of the models. - lmotidyandglancefrom package broom to get most of of the model estimates you want. An alternative totidyis to simply pull thetTableout of the summary output and extract what you want (like SE) from there; e.g.,summary(fm1)$tTable. - aosmithlogLik&residuals. You can match the fixed effects st. errors returned bysqrt(diag(vcov(fm1)))(orsqrt(diag(fm1$varFix), but again better to use extractor functions) to the summary output if you dont adjust themsummary(fm1, adjustSigma = FALSE)or fit your model usingmethod="REML"- see?nlme:::summary.lme- user20650