Rense Nieuwenhuis
2007-Jul-31 09:22 UTC
[R] Extracting random parameters from summary lme and lmer
LS,
I'm estimating multilevel regression models, using the lme-function
from the nlme-package. Let's say that I estimated a model and stored
it inside the object named 'model'. The summary of that model is
shown below:
Using summary(model)$tTable , I receive the following output:
> summary(model)$tTable
Value Std.Error DF t-value p-value
(Intercept) 0.23268607 0.09350662 3990 2.4884449 1.287080e-02
sexM -0.15338225 0.03169762 3990 -4.8389206 1.354802e-06
standLRT 0.38593558 0.01677195 3990 23.0107762 4.005182e-110
vrmid 50% 0.07606394 0.09389376 61 0.8101064 4.210281e-01
vrtop 25% 0.24561327 0.10483374 61 2.3428838 2.241317e-02
intakemid 50% -0.41469716 0.03177240 3990 -13.0521199 3.698344e-38
intaketop 25% -0.75920783 0.05357980 3990 -14.1696648 1.666780e-44
typeSngl 0.15680532 0.07173835 61 2.1857949 3.267903e-02
All looks fine to me. The output above is simply a section from the
full summary shown below. Now, I want to extract from the summary (or
the full model) the part stating the random parameters. More
specifically, I want to extract from the summary the following:
(Intercept) 0.2869401 (Intr)
typeSngl 0.2791040 -0.617
Residual 0.7302233
How could this be done, and how can the same be done using the lmer-
function from the lme4-package?
Thanks for the effort,
Rense Nieuwenhuis
Linear mixed-effects model fit by REML
Data: Exam
AIC BIC logLik
9158.56 9234.241 -4567.28
Random effects:
Formula: ~type | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.2869401 (Intr)
typeSngl 0.2791040 -0.617
Residual 0.7302233
Fixed effects: normexam ~ sex + standLRT + vr + intake + type
Value Std.Error DF t-value p-value
(Intercept) 0.2326861 0.09350662 3990 2.488445 0.0129
sexM -0.1533822 0.03169762 3990 -4.838921 0.0000
standLRT 0.3859356 0.01677195 3990 23.010776 0.0000
vrmid 50% 0.0760639 0.09389376 61 0.810106 0.4210
vrtop 25% 0.2456133 0.10483374 61 2.342884 0.0224
intakemid 50% -0.4146972 0.03177240 3990 -13.052120 0.0000
intaketop 25% -0.7592078 0.05357980 3990 -14.169665 0.0000
typeSngl 0.1568053 0.07173835 61 2.185795 0.0327
Correlation:
(Intr) sexM stnLRT vrm50% vrt25% int50% int25%
sexM -0.201
standLRT -0.125 0.028
vrmid 50% -0.742 0.028 -0.035
vrtop 25% -0.652 0.051 -0.065 0.649
intakemid 50% -0.246 -0.011 0.541 -0.002 0.007
intaketop 25% -0.218 -0.018 0.676 0.014 0.013 0.660
typeSngl -0.421 0.080 0.007 0.033 -0.027 -0.001 0.001
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.59074329 -0.63776965 0.03829878 0.67303837 3.33952680
Number of Observations: 4059
Number of Groups: 65
[[alternative HTML version deleted]]
Chuck Cleland
2007-Jul-31 10:58 UTC
[R] Extracting random parameters from summary lme and lmer
Rense Nieuwenhuis wrote:> LS, > > I'm estimating multilevel regression models, using the lme-function > from the nlme-package. Let's say that I estimated a model and stored > it inside the object named 'model'. The summary of that model is > shown below: > > Using summary(model)$tTable , I receive the following output: > > > summary(model)$tTable > Value Std.Error DF t-value p-value > (Intercept) 0.23268607 0.09350662 3990 2.4884449 1.287080e-02 > sexM -0.15338225 0.03169762 3990 -4.8389206 1.354802e-06 > standLRT 0.38593558 0.01677195 3990 23.0107762 4.005182e-110 > vrmid 50% 0.07606394 0.09389376 61 0.8101064 4.210281e-01 > vrtop 25% 0.24561327 0.10483374 61 2.3428838 2.241317e-02 > intakemid 50% -0.41469716 0.03177240 3990 -13.0521199 3.698344e-38 > intaketop 25% -0.75920783 0.05357980 3990 -14.1696648 1.666780e-44 > typeSngl 0.15680532 0.07173835 61 2.1857949 3.267903e-02 > > > All looks fine to me. The output above is simply a section from the > full summary shown below. Now, I want to extract from the summary (or > the full model) the part stating the random parameters. More > specifically, I want to extract from the summary the following: > > (Intercept) 0.2869401 (Intr) > typeSngl 0.2791040 -0.617 > Residual 0.7302233 > > How could this be done, and how can the same be done using the lmer- > function from the lme4-package??VarCorr> library(nlme)> fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age> summary(fm1)Linear mixed-effects model fit by REML Data: Orthodont AIC BIC logLik 454.6367 470.6173 -221.3183 Random effects: Formula: ~age | Subject Structure: General positive-definite StdDev Corr (Intercept) 2.3270339 (Intr) age 0.2264276 -0.609 Residual 1.3100399 Fixed effects: distance ~ age Value Std.Error DF t-value p-value (Intercept) 16.761111 0.7752461 80 21.620375 0 age 0.660185 0.0712533 80 9.265334 0 Correlation: (Intr) age -0.848 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -3.223106017 -0.493760867 0.007316632 0.472151090 3.916032743 Number of Observations: 108 Number of Groups: 27> VarCorr(fm1)Subject = pdSymm(age) Variance StdDev Corr (Intercept) 5.41508688 2.3270339 (Intr) age 0.05126947 0.2264276 -0.609 Residual 1.71620459 1.3100399> library(lme4)> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))Linear mixed-effects model fit by REML Formula: Reaction ~ Days + (Days | Subject) Data: sleepstudy AIC BIC logLik MLdeviance REMLdeviance 1754 1770 -871.8 1752 1744 Random effects: Groups Name Variance Std.Dev. Corr Subject (Intercept) 610.835 24.7151 Days 35.056 5.9208 0.067 Residual 655.066 25.5943 number of obs: 180, groups: Subject, 18 Fixed effects: Estimate Std. Error t value (Intercept) 251.405 6.820 36.86 Days 10.467 1.546 6.77 Correlation of Fixed Effects: (Intr) Days -0.137> VarCorr(fm1)$Subject 2 x 2 Matrix of class "dpoMatrix" (Intercept) Days (Intercept) 610.834546 9.738707 Days 9.738707 35.056337 attr(,"sc") scale 25.59426> Thanks for the effort, > > Rense Nieuwenhuis > > Linear mixed-effects model fit by REML > Data: Exam > AIC BIC logLik > 9158.56 9234.241 -4567.28 > > Random effects: > Formula: ~type | school > Structure: General positive-definite, Log-Cholesky parametrization > StdDev Corr > (Intercept) 0.2869401 (Intr) > typeSngl 0.2791040 -0.617 > Residual 0.7302233 > > Fixed effects: normexam ~ sex + standLRT + vr + intake + type > Value Std.Error DF t-value p-value > (Intercept) 0.2326861 0.09350662 3990 2.488445 0.0129 > sexM -0.1533822 0.03169762 3990 -4.838921 0.0000 > standLRT 0.3859356 0.01677195 3990 23.010776 0.0000 > vrmid 50% 0.0760639 0.09389376 61 0.810106 0.4210 > vrtop 25% 0.2456133 0.10483374 61 2.342884 0.0224 > intakemid 50% -0.4146972 0.03177240 3990 -13.052120 0.0000 > intaketop 25% -0.7592078 0.05357980 3990 -14.169665 0.0000 > typeSngl 0.1568053 0.07173835 61 2.185795 0.0327 > Correlation: > (Intr) sexM stnLRT vrm50% vrt25% int50% int25% > sexM -0.201 > standLRT -0.125 0.028 > vrmid 50% -0.742 0.028 -0.035 > vrtop 25% -0.652 0.051 -0.065 0.649 > intakemid 50% -0.246 -0.011 0.541 -0.002 0.007 > intaketop 25% -0.218 -0.018 0.676 0.014 0.013 0.660 > typeSngl -0.421 0.080 0.007 0.033 -0.027 -0.001 0.001 > > Standardized Within-Group Residuals: > Min Q1 Med Q3 Max > -3.59074329 -0.63776965 0.03829878 0.67303837 3.33952680 > > Number of Observations: 4059 > Number of Groups: 65 > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- Chuck Cleland, Ph.D. NDRI, Inc. 71 West 23rd Street, 8th floor New York, NY 10010 tel: (212) 845-4495 (Tu, Th) tel: (732) 512-0171 (M, W, F) fax: (917) 438-0894