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

______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

Reply via email to