For a Monte Carlo study I need to extract from an lme model the estimated standard deviation of a random effect and store it in a vector. If I do a print() or summary() on the model, the number I need is displayed in the Console [it's the 0.1590195 in the output below]>print(fit) >Linear mixed-effects model fit by maximum likelihood > Data: datag2 > Log-likelihood: -145.0028 > Fixed: lst1 ~ lst >(Intercept) lst > 1.1080629 0.7582595 > >Random effects: > Formula: ~1 | yeart > (Intercept) Residual >StdDev: 0.1590195 0.3313497However I have not been able to find that number anywhere in the model object returned by lme. Can anyone tell me how to pull it or compute it from a model returned by lme? Ditto for models with random intercept fitted by glmmPQL with family=binomial. Details: The relevant data are the first 3 colums from the data frame extracted below. Yeart is a factor (year of study) ranging from 1 to 16. The model is linear regression of lst1 on lst with constant slope, random intercept varying across years. I'm using R 1.9.1 in WinXP and the current version of nlme from CRAN. Code: datag2=groupedData(lst1~1|yeart,data=cdata); mixed.grow=lme(fixed=lst1~lst,random=~1, method="ML",data=datag2); mixed.surv=glmmPQL(fixed=surv~lst, random=~1|yeart, family=binomial, data=cdata,niter=50); Some of the data: yeart lst lst1 surv flow 1 1 3.158088 3.331216 1 0 2 1 2.472618 2.486410 1 0 3 1 3.582950 3.417807 1 0 4 1 3.554819 3.377117 1 0 5 1 2.830049 2.992426 1 0 6 1 2.779616 3.240449 1 0 7 1 2.580484 2.930154 1 0 Thanks in advance, Steve 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
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.
hope this helps. spencer graves
Stephen Ellner wrote:
>For a Monte Carlo study I need to extract from an lme model
>the estimated standard deviation of a random effect
>and store it in a vector. If I do a print() or summary()
>on the model, the number I need is displayed in the Console
>[it's the 0.1590195 in the output below]
>
>
>
>>print(fit)
>>Linear mixed-effects model fit by maximum likelihood
>> Data: datag2
>> Log-likelihood: -145.0028
>> Fixed: lst1 ~ lst
>>(Intercept) lst
>> 1.1080629 0.7582595
>>
>>Random effects:
>>Formula: ~1 | yeart
>> (Intercept) Residual
>>StdDev: 0.1590195 0.3313497
>>
>>
>
>However I have not been able to find that number anywhere in the
>model object returned by lme. Can anyone tell me how to pull it
>or compute it from a model returned by lme? Ditto for models with
>random intercept fitted by glmmPQL with family=binomial.
>
>Details:
>The relevant data are the first 3 colums from the
>data frame extracted below. Yeart is a factor (year of
>study) ranging from 1 to 16. The model is linear regression of
>lst1 on lst with constant slope, random intercept varying
>across years. I'm using R 1.9.1 in WinXP and the current
>version of nlme from CRAN.
>
>Code:
>datag2=groupedData(lst1~1|yeart,data=cdata);
>mixed.grow=lme(fixed=lst1~lst,random=~1, method="ML",data=datag2);
>
>mixed.surv=glmmPQL(fixed=surv~lst, random=~1|yeart, family=binomial,
>data=cdata,niter=50);
>
>Some of the data:
> yeart lst lst1 surv flow
>1 1 3.158088 3.331216 1 0
>2 1 2.472618 2.486410 1 0
>3 1 3.582950 3.417807 1 0
>4 1 3.554819 3.377117 1 0
>5 1 2.830049 2.992426 1 0
>6 1 2.779616 3.240449 1 0
>7 1 2.580484 2.930154 1 0
>
>Thanks in advance,
>Steve
>
>
>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
>
>
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