R devel - Nov 2001 - lme: how to extract the variance components?

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?

      Thank you very much!

     

Jean

"Jean Wu" <darmy16@hotmail.com> writes:
>       Dear all, 
>       Here is the question:
>       For example, using the "petrol" data offered with R.
I think you mean the "petrol" data from the MASS library for 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)".
You can get a character representation of the number from
VarCorr
> VarCorr(pet3.lme)
No = pdLogChol(1) 
            Variance StdDev  
(Intercept) 2.088107 1.445028
Residual    3.504931 1.872146
> VarCorr(pet3.lme)["(Intercept)", "Variance"]
[1] "2.088107"
> as.numeric(VarCorr(pet3.lme)["(Intercept)",
"Variance"])
[1] 2.088107

A better but more complicated approach is to extract the reStruct
component of the modelStruct component of the fitted model and convert
convert to a matrix.
> str(as.matrix(pet3.lme$modelStruct$reStruct))
List of 1
 $ No: num 0.596
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr "(Intercept)"
  .. ..$ : chr "(Intercept)"

The result is actually a list of matrices because we must allow for
cases of nested grouping factors.  The first (and only, in this case)
element of the list is the relative variance-covariance matrix of the
random effects for the first grouping factor.  To convert it to the
variance-covariance matrix you multiply by sigma^2.
> as.matrix(pet3.lme$modelStruct$reStruct)[[1]] * pet3.lme$sigma^2
            (Intercept)
(Intercept)    2.088107

That number prints the same as the previous value derived from VarCorr
but this extraction retains the full numerical precision because it
does not convert to character and back to numeric.
>       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?
That's from a very old version of the NLME package that uses different
computational methods from those used in NLME 3.0 and higher.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._