[R] cross-classified random factors in lme without blocking - success

[Gordon Smyth]

Sundar's advice seems to do the trick. Here is a small simulation of a 
cross-classified random model with extraction of the fitted variance 
components from lme:

 > library(nlme)
 > set.seed(18112003)
 > na <- 20
 > nb <- 20
 > sigma.a <- 2
 > sigma.b <- 3
 > sigma.res <- 1
 > mu <- 0.5
 > n <- na*nb
 > a <- gl(na,1,n)
 > b <- gl(nb,na,n)
 > u <- gl(1,1,n)
 > ya <- rnorm(na,mean=0,sd=sigma.a)
 > yb <- rnorm(nb,mean=0,sd=sigma.b)
 > y <- model.matrix(~a-1) %*% ya + model.matrix(~b-1) %*% yb + 
rnorm(n,mean=mu,sd=sigma.res)
 > fm <- lme(y~1,random=list(u=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1)))))
 > v <- sqrt(diag(as.matrix(fm$modelStruct$reStruct$u))[c(1,na+1)])
 > v <- fm$sigma * c(v, 1)
 > names(v) <- c("Sigma.a","Sigma.b","Sigma.res")
 > print(v)
   Sigma.a   Sigma.b Sigma.res
  2.450700  3.650427  1.046689

Gordon

At 05:29 AM 16/11/2003, Sundar Dorai-Raj wrote:
>Gordon Smyth wrote:
>
>>Many thanks for help from Peter Dalgaard, Douglas Bates and Bill 
>>Venables. As a result of their help, here is a working example of using 
>>lme to fit an additive random effects model. The model here is 
>>effectively y~a+b with a and b random:
>>y <- rnorm(12)
>>a <- gl(4,1,12)
>>b <- gl(3,4,12)
>>u <- gl(1,1,12)
>>library(nlme)
>>fm <- lme(y~1,random=list(u=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1)))))
>>summary(fm)
>>I still can't see though how to extract the three variance components 
>>(for a and b and for the residual) from the fitted object 'fm'.
>>Thanks
>>Gordon
>
>Would this work for you? It still needs some parsing to extract sigma_a 
>and sigma_b, but that should be simple.
>
>v <- as.matrix(fm$modelStruct$reStruct$u)
>c(sqrt(diag(v)), Residual = 1) * fm$sigma
>
>Regards,
>Sundar