[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