[R] How to extract the variance componets from lme
joris.dewolf at cropdesign.com
joris.dewolf at cropdesign.com
Tue Jan 2 15:12:04 CET 2007
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 guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list