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
Running lme on your data set results exactly in what you expect - or do
you expect something different?
Pascal
> L1<-factor(F1f)
> L2<-factor(F2f)
> L3<-factor(F3f)
> lme(value ~ 1,random = ~ 1 | L1/L2/L3)
Linear mixed-effects model fit by REML
Data: NULL
Log-restricted-likelihood: 438.9476
Fixed: value ~ 1
(Intercept)
0.2955631
Random effects:
Formula: ~1 | L1
(Intercept)
StdDev: 0.02472988 <== level F1 which is 0
Formula: ~1 | L2 %in% L1
(Intercept)
StdDev: 1.140782 <== level F2 which is 1
Formula: ~1 | L3 %in% L2 %in% L1
(Intercept) Residual
StdDev: 0.0005512791 0.1020479 <== F3 which is 0 , and F4 which is 0.1
Number of Observations: 625
Number of Groups:
L1 L2 %in% L1 L3 %in% L2 %in% L1
5 25 125
On Tue, 16 Dec 2003, Gary Allison wrote:> 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. >While I haven't seen any literature on this precise problem I would not be particularly surprised by it. Your model says that there is a lot of variability at level two but that it all averages out to zero within a level 1 unit. It would not be particularly surprising if some of the level 2 variability leaked up to level one, since it is only being averaged over 5 level 2 units per level 1 unit. Given sufficient data this will eventually stop happening, but you don't have very many level 1 replicates either. I think that if you really need to estimate a small variance component with large variances nested within it, you need a lot more data or a prior. -thomas