[R] lme: extract variance estimate

From: Stephen Ellner <spe2_at_cornell.edu>
Date: Wed 07 Jul 2004 - 06:26:52 EST


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@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 Received on Wed Jul 07 06:31:51 2004

This archive was generated by hypermail 2.1.8 : Fri 18 Mar 2005 - 09:45:59 EST