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