[R] lme: extract variance estimate
spencer.graves at pdf.com
Wed Jul 7 22:02:11 CEST 2004
I just tried it in both lme4 and nlme. I got it in nlme but not
lme4. I'm sure lme4 is better in many ways, but I could not figure out
how to get what you want in lme4.
Specifically, I tried the following:
DF <- data.frame(x=rep(letters[1:2], 2), y=rep(1:2, 2)+0.01*(1:4))
fit <- lme(y~1, random=~1|x, data=DF)
VC <- VarCorr(fit)
When I did this in library(nlme), I got the following:
x = pdLogChol(1)
(Intercept) 0.5101563004 0.71425227
Residual 0.0001999596 0.01414071
However, when I did it in library(lme4), I got something different:
Groups Name Variance Std.Dev.
x (Intercept) 0.50995 0.71411
Residual 2e-04 0.014142
Error in VC[1, 2] : incorrect number of dimensions
I was concerned by the differences in the estimates in this case,
so I ported this problem to S-Plus 6.2. There I got the "lme4" answers
from both "lme" and "varcomp".
hope this helps. spencer graves
Stephen Ellner wrote:
>Spencer Graves wrote:
>> Have you considered "VarCorr"? I've used it with "lme", and the
>>documentation in package lme4 suggests it should work with GLMM, which
>>might also do what you want from glmmPQL.
>Thanks for the pointer (I was not aware that nlme and lme4 had different
>versions of lme), but I'm still stuck at the same place using lme4:> VarCorr(fit)
> Groups Name Variance Std.Dev.
> yeart (Intercept) 0.040896 0.20223
> Residual 0.091125 0.30187
>The number I need to extract and store is the .20223, but all the
>components I can find in VarCorr(fit) are something else.
> "scale" "reSumry" "useScale"
>>u at scale
>>u at reSumry
>An object of class "corrmatrix"
>>u at useScale
>In glmmML the estimate is returned as the $sigma component
>of the model, but I also need the same info from 'family=gaussian'
>Stephen P. Ellner (spe2 at cornell.edu)
>Department of Ecology and Evolutionary Biology
>Corson Hall, Cornell University, Ithaca NY 14853-2701
>Phone (607) 254-4221 FAX (607) 255-8088
>R-help at stat.math.ethz.ch mailing list
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help