I recommend sending questions like this to the R-SIG-Mixed-Models at R-project.org mailing list which I have cc:'d on this reply. Many of those who read that list can reply faster than I am able to. On Thu, Sep 11, 2014 at 9:06 AM, Hung-Chih Ku <
Hung-Chih.Ku at utsouthwestern.edu> wrote:
Dear Dr. Bates,
My name is Hung-Chih Ku and I am a postdoctoral fellow at UT
Southwestern Medical Center. I am trying to estimate variance components
from a linear mixed model Y=Xa+Zb+e where Z is an nxp unstructured matrix.
The following is an example of R code.
n <- 50
p <- 5
X <- runif(n)
Z <- NULL
for (i in 1:p) {
Z <- cbind(Z, rbinom(n, 2, .3))
}
a <- 0.5
b <- rnorm(p)
e<- rnorm(n)
Y <- X*a + Z*b + e
library(nlme)
group=rep(1,n)
fit <- lme(Y~X, random=list(group=pdIdent(~-1+Z)))
VarCorr(fit)
However, when n and p increase, let's say n=1000 and p=500, the lme will
be taking a long time to estimate two variances (variances of b and e). Is
there a way to speed up? Thank you for your time and I am looking forward
to hearing from you.
When n=1000 and p=500 you are trying to estimate 500,000 random effects coefficients from 1000 observations which seems a bit optimistic You may want to reexamine your model specification. It doesn't make sense to me.