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 Liu
------------------------------
Peng Liu |
Division of Statistics |
Northern Illinois University |
De Kalb, IL 60115, USA |
E-mail: pliu at math.niu.edu |
------------------------------
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 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
--- Renaud Lancelot <renaud.lancelot at cirad.fr> wrote:> Hi Peng, > > Peng wrote: > [snip] > > 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>library(nlme) > R>data(Orthodont) > R>fm1 <- lme(distance ~ age, data = Orthodont) # > random is ~ age > R> > R>lapply(pdMatrix(fm1$modelStruct$reStruct), "*", > fm1$sigma^2) > $Subject > (Intercept) age > (Intercept) 5.414722 -0.32102403 > age -0.321024 0.05126664 > > Comment: pdMatrix returns the list of var-cov > matrices of the random > effects (one matrix for each level of the possibly > multilevel > structure). For some computational reason, these > matrices are scaled by > the residuals variance: you have to multiply them by > the residuals variance. > > If you really want a matrix (the above is a list > with one component), > you can do: > > R>lapply(pdMatrix(fm1$modelStruct$reStruct), "*", > fm1$sigma^2)$Subject > (Intercept) age > (Intercept) 5.414722 -0.32102403 > age -0.321024 0.05126664Hi, Renaud, I think reStruct gives the structure of random effects, not the within group errors. And I happened to want the structure for the later one, especially when the within group errors are differently structured between groups. Regards, Peng ------------------------------ Peng Liu | Division of Statistics | Northern Illinois University | De Kalb, IL 60115, USA | E-mail: pliu at math.niu.edu | ------------------------------> > > Best, > > Renaud > > -- > Dr Renaud Lancelot, vétérinaire > CIRAD, Département Elevage et Médecine Vétérinaire > (CIRAD-Emvt) > Programme Productions Animales > http://www.cirad.fr/fr/pg_recherche/page.php?id=14 > > ISRA-LNERV tel +221 832 49 > 02 > BP 2057 Dakar-Hann fax +221 821 18 > 79 (CIRAD) > Senegal e-mail > renaud.lancelot at cirad.fr >