[R-sig-ME] lme4::mcmcsamp + coda::HPDinterval

[Sundar Dorai-Raj]

Hi, lmer+coda Users,

(reproducible code at end)

I have a question regarding confidence intervals on the random effects. 
The data I'm working with is a highly unbalanced, nested design with two 
random effects (say, A/B), where the B variance is expected to be larger 
than A variance. I need to compute confidence bounds on the standard 
deviations of each effect (A and A:B). To do this, I use lme4::mcmcsamp 
with coda::HPDinterval. As in,

f <- lmer(...)
m <- mcmcsamp(f, n = 1000)
(s <- sqrt(exp(HPDinterval(m))))

However, one of the point estimates falls below the lower bound of the 
confidence interval. My guess is that this is due to the correlation 
between A and A:B (due to imbalance? relative magnitude of A vs. A:B?) 
leading to an unstable model fit. Are the point estimates completely 
untrustworthy? In this case, should I simply remove the offending random 
effect and refit? Is there a reference that describes situations such as 
these (point estimates outside the C.I.)?

Also, this may be an example where the nlme:::intervals.lme function 
produces intervals that are nonsensical (at least to me) if you change 
sd.A below to, say, 0.5.

Thanks,

--sundar

<code>
set.seed(2)
## simulate data
A <- factor(rep(1:10, each = 4))
n <- length(A)
B <- factor(rep(1:2, 20))
sd.A <- 5
sd.B <- 10
sd.e <- 0.5
rA <- rnorm(nlevels(A), 0, sd.A)
rB <- rnorm(nlevels(A:B), 0, sd.B)
e <- rnorm(n, 0, sd.e)
y <- 0.5 + rA[A] + rB[A:B] + e
## create unbalance
out <- seq(length(y)) %in% sample(length(y), 20)

## fit model
library(lme4)
library(coda)
f <- lmer2(y ~ (1 | A) + (1 | A:B), subset = !out)
m <- mcmcsamp(f, 1000)
v <- VarCorr(f)
s.lmer <- cbind(sqrt(exp(HPDinterval(m[, -1]))),
                 est. = c(attr(v, "sc"), sqrt(sapply(v, as.matrix))),
                 true = c(sd.e, sd.B, sd.A))
s.lmer <- s.lmer[, c("lower", "est.", "upper", "true")]
rownames(s.lmer) <- c("sigma", "A:B", "A")
print(zapsmall(s.lmer), digits = 4)
</code>