[Rd] lme: how to extract the variance components?

Peter Dalgaard BSA p.dalgaard@biostat.ku.dk
14 Nov 2001 19:19:39 +0100


"Jean Wu" <darmy16@hotmail.com> writes:

>       Dear all, 
>       Here is the question:
>       For example, using the "petrol" data offered with R.
>       pet3.lme<-lme(Y~SG+VP+V10+EP,random=~1|No,data=petrol)
>       pet3.lme$sigma gives the residual StdDev.
>       But I can't figure out how to extract the "(intercept) StdDev", 
>       although it is in the print out if I do "summary(pet3.lme)".
> 
>       In S-plus3.4 , $var.ran is the estimate of the variance of the random 
>       effects between clusters. what is the equivalent thing to do in R?

It's been a while, but I have some old code saying


lme.obj <- lme(log(Height)~ns(sqrt(Age),knots=c(0.25,.5,1,5),
  Boundary.knots=c(0,10)), random=~sqrt(Age)|ID,
  correlation=corExp(form=~sqrt(Age),nugget=F))
Age.new <- seq(0,10,0.01)  
C.mat  <- lme.obj$sigma^2*as.matrix(lme.obj$modelStruct$reStruct$ID)
SD <- sqrt(sapply(Age.new,
  function(a){x<-c(1,sqrt(a)); t(x) %*% C.mat %*% x})
  +lme.obj$sigma^2)

so I would expect that you want 

pet3.lme$sigma^2 * as.matrix(pet3.lme$modelStruct$reStruct$No)


-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._