[R-sig-ME] negatively biased variance estimates of random factors with few levels

[kalman.toth]

Dear List Members,

On the glmmFAQ webpage Ben Bolker writes:

“Treating factors with small numbers of levels as random will in the best case lead to very small and/or imprecise estimates of random effects; in the worst case it will lead to various numerical difficulties such as lack of convergence, zero variance estimates, etc.. (A small [simulation exercise](https://rpubs.com/bbolker/4187) shows that at least the estimates of the standard deviation are downwardly biased in this case; it’s not clear whether/how this bias would affect the point estimates of fixed effects or their estimated confidence intervals.) In the classical method-of-moments approach these problems may not arise (because the sums of squares are always well defined as long as there are at least two units), but the underlying problems of lack of power are there nevertheless. “

I am struggling to understand the reasons behind the potential bias in REML estimates when dealing with a small number of factor levels. To explore this issue, I modified Ben Bolker’s simulation as follows:

- set the sd1/sd2 ratio to 5/1 (to reduce the number of zero variance estimates and eliminate the possibility of causing bias)

- compared the REML estimates with ANOVA because 1) the ANOVA variance estimates are expected to be unbiased and 2) for these datasets the two estimates should be the same

- extended the simulation to sample sizes 3, 5, 10 and 20

My conclusions based on the extended simulation are:

- Consistency of REML and ANOVA estimates: the estimates are indeed the same

- Negative bias: Both methods give negatively biased estimates!?

- Impact of sample size: the magnitude of bias decreases as the sample size increases (11, 7, 4, and 1% for sample sizes 3, 5, 10 and 20, respectively) but sample sizes as large as 20 or 50 might be still negatively biased

I am particularly interested in addressing the treatment of factors conceptually viewed as random effects but having only a few factor levels. This scenario is common in the datasets I work with. I understand that there can be issue with the power and the SE can be large if the sample size is small. But why would both ANOVA and REML estimates be biased? And what could be done about the bias apart from treating the factor as fixed? (I prefer to treat them as fixed if they really have few levels e.g. 3-4 but not with 8-10 levels)

Regards,

Kalman

Please, find here the modified simulation:

nrep <- 5

simfun <- function(n1 = 20, n2 = nrep, sd1 = 1, sd2 = 0.2) {
d <- expand.grid(f1 = factor(seq(n1)), f2 = factor(seq(n2)))
u1 <- rnorm(n1, mean = 0, sd = sd1)
d$y <- rnorm(n1 * n2, mean = u1, sd = sd2)
d
}
require(lme4)
fitfun <- function(d = simfun()) {
lme <- sqrt(unlist(VarCorr(lmer(y ~ (1 | f1), data = d))))
MSb <- summary(aov(y ~ f1, data = d))[[1]]['Mean Sq']['f1',]
MSw <- summary(aov(y ~ f1, data = d))[[1]]['Mean Sq']['Residuals',]
an <- sqrt((MSb-MSw)/nrep)
list(lme = lme, anova = an)
}

set.seed(101)
nsim = 500
sd_dist0 <- replicate(nsim, fitfun())
sd_dist1 <- replicate(nsim, fitfun(simfun(n1 = 10)))
sd_dist2 <- replicate(nsim, fitfun(simfun(n1 = 5)))
sd_dist3 <- replicate(nsim, fitfun(simfun(n1 = 3)))
sd_List <- list(lme.20 = unlist(sd_dist0[1,]), lme.10 = unlist(sd_dist1[1,]), lme.5 = unlist(sd_dist2[1,]), lme.3 = unlist(sd_dist3[1,]),
anova.20 = unlist(sd_dist0[2,]), anova.10 = unlist(sd_dist1[2,]), anova.5 = unlist(sd_dist2[2,]), anova.3 = unlist(sd_dist3[2,]))

#sd_List <- list(n1.5 = sd_dist1, n1.3 = sd_dist2)

plotfun <- function(x, title) {
par(las = 1, bty = "l")
hist(x, breaks = 50, col = "gray", main = title, xlab = "est. sd", freq = FALSE)
}

par(mfrow = c(2, 4))
invisible(mapply(plotfun, sd_List, names(sd_List)))

sapply(sd_List, function(x) mean(x == 0 | is.na(x)))

sfun <- function(x) {
r <- list(mean = mean(x, na.rm=T), mean2 = mean(x[which(x > 1e-5 )]), sem = sd(x, na.rm=T)/sqrt(length(x)))
r <- with(r, c(r, list(lwr = mean - 2 * sem, upr = mean + 2 * sem)))
unlist(r)
}
print(s_tab <- sapply(sd_List, sfun), digits = 3)

bias_pct <- round((s_tab["mean", ] - 1) * 100)
bias_pct2 <- round((s_tab["mean2", ] - 1) * 100)
bias_pct
bias_pct2 # bias without zero SD estimates
	[[alternative HTML version deleted]]