[R] variance estimates in lme biased?

Gary Allison allison.100 at osu.edu
Wed Dec 17 02:35:37 CET 2003


Hi all,
I didn't get a response to my post of this issue a week ago, so I've 
tried to clarify:

When I use lme to analyze a model of nested random effects, the variance 
estimates of levels higher in the hierarchy appear to have much more 
variance than they should.

In the example below with 4 levels, I simulate variance in level 2 
(sd=1.0) and level 4 (sd=0.1), but levels 1 and 3 do not vary.  Although 
I expected to see a small amount of variability in the lme estimates of 
levels 1 and 3, I am confused by what happens in the level 1 estimates: 
  more than 10% of the runs produce stdDev of >0.4 and in many cases the 
level 1 estimate is close to level 2's. Nothing like that happens in 
level 3.

I use R 1.8.1 & Win2000.  Since I last posted, I have also seen these 
symptoms in SAS when I use PROC VARCOMP, although I haven't explored it 
much there.

library(nlme)
F1V <- F2V <- F3V  <- residV <- numeric(0)
for (rep in 1:500){
   F1f <- F2f <- F3f <- value <- numeric(0)
   # data generator
   for (F1 in 1:5) {
     for (F2 in 1:5) {
       lev2 <- rnorm(1,sd=1.0)
       for (F3 in 1:5) {
         for (F4 in 1:5) {
           lev4 <- rnorm(1,sd=0.1)
           value <- c(value, lev2+lev4)
           F1f <- c(F1f,F1);F2f <- c(F2f,F2);F3f <- c(F3f,F3)
         }
       }
     }
   }
   L1 <- as.factor(F1f); L2 <- as.factor(F2f); L3 <- as.factor(F3f)
   y.lme <- lme(value ~ 1,random = ~ 1 | L1/L2/L3)
   v <- as.numeric(VarCorr(y.lme)[,2])
   v <- as.numeric(na.omit(v))
   F1V <- c(F1V,v[1]);  F2V <- c(F2V,v[2]);  F3V <- c(F3V,v[3])
   residV <- c(residV,y.lme$sigma)
}

Because this code can take several minutes to run (probably because of 
my heavy use of for loops), I've posted a set of results at:
http://david.science.oregonstate.edu/~allisong/R/nestedVar.pdf
   (see last page for summary)

If anyone can help me make sense of this, point out an error I'm making 
or point me to some literature, I would greatly appreciate it!

Thanks,
Gary

--
Gary Allison
Evolution, Ecology and Organismal Biology
Ohio State University




More information about the R-help mailing list