[R] How to extract the variance componets from lme

BXC (Bendix Carstensen) bxc at steno.dk
Tue Jan 2 12:50:05 CET 2007


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



More information about the R-help mailing list