> 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
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
