Re: [R] lme: extract variance estimate

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Thu 08 Jul 2004 - 06:02:11 EST

      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@scale
>>
>>
>[1] 0.3018693
>
>
>>u@reSumry
>>
>>
>$yeart
>An object of class "corrmatrix"
> (Intercept)
>(Intercept) 1
>Slot "stdDev":
>(Intercept)
> 0.6699156
>
>
>>u@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@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@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
>
>



R-help@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 Received on Thu Jul 08 06:06:44 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 02:34:54 EST