[R] nested analysis with lme - odd result?

[Gary Allison]

Hello!

When I simulate variance at only a single level in a nested analysis 
using lme (all levels are random effects), the results confuse me. 
Instead of lme reporting high variance in only that simulated level, 
substantial variance (>10% of simulated level) often appears in other 
levels -- in some configurations, as often as 50% of the time. Usually 
this spurious variance shows up in levels of nesting above the level 
with high simulated variance.

Am I doing something wrong or should I expect such 'over-detection'? I 
was expecting near zero variance in those other levels except say, 5% of 
the time.

Here's some example code with variance simulated in level 2 and spurious 
variance appearing frequently in level 1:

library(nlme)
# 4 level nesting
simVar <- c(0,1,0,.1) # first is level one, last becomes error
nAtLevel <- c(5,5,5,5)  # number of replicates at each level

F1V <- F2V <- F3V  <- residV <- numeric(0)
for (rep in 1:100){
   F1f <- F2f <- F3f <- value <- numeric(0)
   mn <- numeric(length(nAtLevel))
   # data generator
   for (F1 in 1:nAtLevel[1]) {
     mn[1] <- rnorm(1,sd=simVar[1]) # set mean for level 1
     for (F2 in 1:nAtLevel[2]) {
       mn[2] <- rnorm(1,sd=simVar[2]) # set mean for level 2
       for (F3 in 1:nAtLevel[3]) {
         mn[3] <- rnorm(1,sd=simVar[3]) # set mean for level 3
         for (F4 in 1:nAtLevel[4]) {
           mn[4] <- rnorm(1,sd=simVar[4]) #set mean for lowest level
           value <- c(value, sum(mn))
           F1f <- c(F1f,F1)
           F2f <- c(F2f,F2)
           F3f <- c(F3f,F3)
         }
       }
     }
   }

   y.lme <- lme(value ~ 1,random = ~ 1 | as.factor(F1f)/
                                         as.factor(F2f)/
                                         as.factor(F3f))
   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)
}
var.df <- data.frame(F1V,F2V,F3V,residV)
par(mfrow=c(2,2))
hist(var.df$F1V,main='Level 1 StdDev')
hist(var.df$F2V,main='Level 2 StdDev')
hist(var.df$F3V,main='Level 3 StdDev')
hist(var.df$residV,main='Residual StdDev')
txt <- 'Simulated stddev:'
for (i in 1:length(simVar)) {
   txt <- paste(txt,' level ',i,'=',simVar[i],',',sep='')}
mtext(txt, outer=T,line=-1,side=3)


Summary plots for several sample sizes is at:
http://david.science.oregonstate.edu/~allisong/R/sampSize.pdf

Thanks for any suggestions!
Gary

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