[R] How to extract R{i} from lme object?

[Peng]

--- Sundar Dorai-Raj <sundar.dorai-raj at pdf.com> wrote:
> 
> 
> Peng wrote:
> > Hi, lme() users,
> > 
> > Can some one tell me how to do this.
> > I model Orthodont with the same G for random
> > variables, but different R{i}'s for boys and
> girls, so
> > that I can get sigma1_square_hat for boys and
> > sigma2_square_hat for girls.
> > 
> > The model is Y{i}=X{i}beta + Z{i}b + e{i}
> > b ~ iid N(0,G) and e{i} ~ iid N(0,R{i}) i=1,2
> > orth.lme <- lme(distance ~ Sex * age,
> data=Orthodont,
> > random=~age|Subject,
> weights=varIdent(form=~1|Sex),
> > method="ML")
> > 
> > I can see the numbers I need from summary(), but
> how
> > can I extract them? I tried several functions in
> nlme,
> > but I cannot find a correct one.
> > 
> > Peng
> > 
> 
> Peng,
> 
> Is this what you need?
> 
> R> data(Orthodont)
> R> orth.lme <- lme(distance ~ Sex * age,
> +                  data=Orthodont,
> +                  random=~age|Subject,
> +                  weights=varIdent(form=~1|Sex),
> +                  method="ML")
> R> orth.lme$modelStruct$varStruct
> Variance function structure of class varIdent
> representing
>       Male    Female
> 1.0000000 0.4112708

Sundar, 

Thanks a lot!
In the above model, I assume the error term for boy
has iid N(0, sigma1), and that for girl has iid
N(0,sigma2). I hope that I wrote the correct lme(). If
so, what I want is:
(1/unique(attributes(orth.lme$modelStruct$varStruct)$weights)*orth.lme$sigma)^2
[1] 2.6292575 0.4447223
, where 2.629 for boys, and 0.445 for girls.

I am still wondering whether there is a function like
getVarCov(), so that I can get R var-cov matrix
directly. I looked through the function list of nlme,
but I cannot find.

Peng

> R> rcov.unscaled <-
> as.matrix(orth.lme$modelStruct$reStruct[[1]])
> R> rcov.scaled <- rcov.unscaled * orth.lme$sigma^2
> R> Stdev <- sqrt(diag(rcov.scaled))
> R> Stdev
> (Intercept)         age
>    1.7914848   0.1408001
> R>
> 
> If there's nesting, look for more elements in
> $reStruct.
> 
> Regards,
> Sundar
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help