[R] How to extract R{i} from lme object?
Peng
l12345p2000 at yahoo.com
Tue Mar 4 23:11:54 CET 2003
--- 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
More information about the R-help
mailing list