BXC (Bendix Carstensen)
2007-Jan-02 11:50 UTC
[R] How to extract the variance componets from lme
Here is a piece of code fitting a model to a (part) of a dataset, just
for
illustration. I can extract the random interaction and the residual
variance
in group meth==1 using VarCorr, but how do I get the other residual
variance?
Is there any way to get the other variances in numerical form directly -
it
seems a litte contraintuitive to use "as.numeric" when extracting
estimates, it's
a bit like good old days writing SAS-programs that reads the SAS output
files...
Bendix
------------------------------------------------------------------------
-------
library( nlme )
dfr <-
structure(list(meth = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), .Label = c("CO",
"pulse"), class = "factor"), item = c(1, 1, 1, 2, 2, 2, 3,
3,
3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4), repl = c(1,
2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1,
2, 3), y = c(78, 76.4, 77.2, 68.7, 67.6, 68.3, 82.9, 80.1, 80.7,
62.3, 65.8, 67.5, 71, 72, 73, 68, 67, 68, 82, 77, 77, 43, 69,
77)), .Names = c("meth", "item", "repl",
"y"), row.names = c("1",
"2", "3", "4", "5", "6",
"7", "8", "9", "10", "11",
"12", "184",
"185", "186", "187", "188",
"189", "190", "191", "192",
"193",
"194", "195"), class = "data.frame")
m1 <-
lme( y ~ factor( meth ) + factor( item ),
random = ~1 | MI,
weights = varIdent( form = ~1 | meth ),
method ="REML",
data = cbind( dfr, MI=interaction( dfr$meth, dfr$item ) ) )
m1
# The MI std and the residual std for meth==1
as.numeric(VarCorr(m1)[,2])
______________________________________________
Bendix Carstensen
Senior Statistician
Steno Diabetes Center
Niels Steensens Vej 2-4
DK-2820 Gentofte
Denmark
+45 44 43 87 38 (direct)
+45 30 75 87 38 (mobile)
+45 44 43 73 13 (fax)
bxc at steno.dk http://www.biostat.ku.dk/~bxc
joris.dewolf at cropdesign.com
2007-Jan-02 14:12 UTC
[R] How to extract the variance componets from lme
I advice you strongly to use VarCorr(), but if you insist tmp <- as.matrix(m1$modelStruct$reStruct$MI) c(sqrt(diag(tmp)), Residual = 1) * m1$sigma Joris r-help-bounces at stat.math.ethz.ch wrote on 02/01/2007 12:50:05:> Here is a piece of code fitting a model to a (part) of a dataset, just > for > illustration. I can extract the random interaction and the residual > variance > in group meth==1 using VarCorr, but how do I get the other residual > variance? > > Is there any way to get the other variances in numerical form directly - > it > seems a litte contraintuitive to use "as.numeric" when extracting > estimates, it's > a bit like good old days writing SAS-programs that reads the SAS output > files... > > Bendix > ------------------------------------------------------------------------ > ------- > > library( nlme ) > > dfr <- > structure(list(meth = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, > 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), .Label = c("CO", > "pulse"), class = "factor"), item = c(1, 1, 1, 2, 2, 2, 3, 3, > 3, 4, 4, 4, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4), repl = c(1, > 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, > 2, 3), y = c(78, 76.4, 77.2, 68.7, 67.6, 68.3, 82.9, 80.1, 80.7, > 62.3, 65.8, 67.5, 71, 72, 73, 68, 67, 68, 82, 77, 77, 43, 69, > 77)), .Names = c("meth", "item", "repl", "y"), row.names = c("1", > "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "184", > "185", "186", "187", "188", "189", "190", "191", "192", "193", > "194", "195"), class = "data.frame") > > m1 <- > lme( y ~ factor( meth ) + factor( item ), > random = ~1 | MI, > weights = varIdent( form = ~1 | meth ), > method ="REML", > data = cbind( dfr, MI=interaction( dfr$meth, dfr$item ) ) ) > > m1 > > # The MI std and the residual std for meth==1 > as.numeric(VarCorr(m1)[,2]) > > ______________________________________________ > > Bendix Carstensen > Senior Statistician > Steno Diabetes Center > Niels Steensens Vej 2-4 > DK-2820 Gentofte > Denmark > +45 44 43 87 38 (direct) > +45 30 75 87 38 (mobile) > +45 44 43 73 13 (fax) > bxc at steno.dk http://www.biostat.ku.dk/~bxc > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code.