[R] lme: extract variance estimate
Spencer Graves
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)
VC
VC[1,2]
When I did this in library(nlme), I got the following:
> VC
x = pdLogChol(1)
Variance StdDev
(Intercept) 0.5101563004 0.71425227
Residual 0.0001999596 0.01414071
> VC[1,2]
[1] "0.71425227"
However, when I did it in library(lme4), I got something different:
> VC
Groups Name Variance Std.Dev.
x (Intercept) 0.50995 0.71411
Residual 2e-04 0.014142
> VC[1,2]
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.
>
>u=VarCorr(fit); slotNames(u)
>[1] "scale" "reSumry" "useScale"
>
>
>>u at scale
>>
>>
>[1] 0.3018693
>
>
>>u at reSumry
>>
>>
>$yeart
>An object of class "corrmatrix"
> (Intercept)
>(Intercept) 1
>Slot "stdDev":
>(Intercept)
> 0.6699156
>
>
>>u at useScale
>>
>>
>[1] TRUE
>
>In glmmML the estimate is returned as the $sigma component
>of the model, but I also need the same info from 'family=gaussian'
>models.
>
>
>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
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>
More information about the R-help
mailing list