How trustworthy are the confidence intervals for lmer objects through effects package?

Question

Effects package provides a very fast and convenient way for plotting linear mixed effect model results obtained through lme4 package. The effect function calculates confidence intervals (CIs) very quickly, but how trustworthy are these confidence intervals?

For example:

library(lme4)
library(effects)
library(ggplot)

data(Pastes)

fm1  <- lmer(strength ~ batch + (1 | cask), Pastes)
effs <- as.data.frame(effect(c("batch"), fm1))
ggplot(effs, aes(x = batch, y = fit, ymin = lower, ymax = upper)) + 
  geom_rect(xmax = Inf, xmin = -Inf, ymin = effs[effs$batch == "A", "lower"],
        ymax = effs[effs$batch == "A", "upper"], alpha = 0.5, fill = "grey") +
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

According to CIs calculated using effects package, batch "E" does not overlap with batch "A".

If I try the same using confint.merMod function and the default method:

a <- fixef(fm1)
b <- confint(fm1)
# Computing profile confidence intervals ...
# There were 26 warnings (use warnings() to see them)

b <- data.frame(b)
b <- b[-1:-2,]

b1 <- b[[1]]
b2 <- b[[2]]

dt <- data.frame(fit   = c(a[1],  a[1] + a[2:length(a)]), 
                 lower = c(b1[1],  b1[1] + b1[2:length(b1)]), 
                 upper = c(b2[1],  b2[1] + b2[2:length(b2)]) )
dt$batch <- LETTERS[1:nrow(dt)]

ggplot(dt, aes(x = batch, y = fit, ymin = lower, ymax = upper)) +
  geom_rect(xmax = Inf, xmin = -Inf, ymin = dt[dt$batch == "A", "lower"], 
        ymax = dt[dt$batch == "A", "upper"], alpha = 0.5, fill = "grey") + 
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

I see that all of the CIs overlap. I also get warnings indicating that the function failed to calculate trustworthy CIs. This example, and my actual dataset, makes me to suspect that effects package takes shortcuts in CI calculation that might not entirely be approved by statisticians. How trustworthy are the CIs returned by effect function from effects package for lmer objects?

What have I tried: Looking into the source code, I noticed that effect function relies on Effect.merMod function, which in turn directs to Effect.mer function, which looks like this:

effects:::Effect.mer
function (focal.predictors, mod, ...) 
{
    result <- Effect(focal.predictors, mer.to.glm(mod), ...)
    result$formula <- as.formula(formula(mod))
    result
}
<environment: namespace:effects>

mer.to.glm function seems to calculate Variance-Covariate Matrix from the lmerobject:

effects:::mer.to.glm

function (mod) 
{
...
mod2$vcov <- as.matrix(vcov(mod))
...
mod2
}

This, in turn, is probably used in Effect.default function to calculate CIs (I might have misunderstood this part):

effects:::Effect.default
...
     z <- qnorm(1 - (1 - confidence.level)/2)
        V <- vcov.(mod)
        eff.vcov <- mod.matrix %*% V %*% t(mod.matrix)
        rownames(eff.vcov) <- colnames(eff.vcov) <- NULL
        var <- diag(eff.vcov)
        result$vcov <- eff.vcov
        result$se <- sqrt(var)
        result$lower <- effect - z * result$se
        result$upper <- effect + z * result$se
...

I do not know enough about LMMs to judge whether this is a right approach, but considering the discussion around confidence interval calculation for LMMs, this approach appears suspiciously simple.

Answer 1

All of the results are essentially the same (for this particular example). Some theoretical differences are:

• as @rvl points out, your reconstruction of CIs without taking account of covariance among parameters is just wrong (sorry)
• confidence intervals for parameters can be based on Wald confidence intervals (assuming a quadratic log-likelihood surface): lsmeans, effects, confint(.,method="Wald"); except for lsmeans, these methods ignore finite-size effects ("degrees of freedom"), but in this case it barely makes any difference (df=40 is practically indistinguishable from infinite df)
• ... or on profile confidence intervals (the default method; ignores finite-size effects but allows for non-quadratic surfaces)
• ... or on parametric bootstrapping (the gold standard -- assumes the model is correct [responses are Normal, random effects are Normally distributed, data are conditionally independent, etc.], but otherwise makes few assumptions)

I think all of these approaches are reasonable (some are more approximate than others), but in this case it barely makes any difference which one you use. If you're concerned, try out several contrasting methods on your data, or on simulated data that resemble your own, and see what happens ...

(PS: I wouldn't put too much weight on the fact that the confidence intervals of A and E don't overlap. You'd have to do a proper pairwise comparison procedure to make reliable inferences about the differences between this particular pair of estimates ...)

95% CIs:

Comparison code:

library(lme4)
fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
c0 <- confint(fm2,method="Wald")
c1 <- confint(fm2)
c2 <- confint(fm2,method="boot")
library(effects)
library(lsmeans)
c3 <- with(effect("batch",fm2),cbind(lower,upper))
c4 <- with(summary(lsmeans(fm2,spec="batch")),cbind(lower.CL,upper.CL))
tmpf <- function(method,val) {
    data.frame(method=method,
               v=LETTERS[1:10],
               setNames(as.data.frame(tail(val,10)),
                        c("lwr","upr")))
}
library(ggplot2); theme_set(theme_bw())
allCI <- rbind(tmpf("lme4_wald",c0),
      tmpf("lme4_prof",c1),
      tmpf("lme4_boot",c2),
      tmpf("effects",c3),
               tmpf("lsmeans",c4))
ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
    geom_linerange(position=position_dodge(width=0.8))

ggsave("pastes_confint.png",width=10)

Answer 2

It looks like what you have done in the second method is to have computed confidence intervals for the regression coefficients, then transformed those to obtain CIs for the predictions. This ignores the covariances between the regression coefficients.

Try fitting the model without an intercept, so that the batch effects will actually be the predictions, and confint will return the intervals you need.

Addendum 1

I did exactly what I suggested above:

> fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
> confint(fm2)
Computing profile confidence intervals ...
           2.5 %    97.5 %
.sig01  0.000000  1.637468
.sigma  2.086385  3.007380
batchA 60.234772 64.298581
batchB 57.268105 61.331915
batchC 60.018105 64.081915
batchD 57.668105 61.731915
batchE 53.868105 57.931915
batchF 59.001439 63.065248
batchG 57.868105 61.931915
batchH 61.084772 65.148581
batchI 56.651439 60.715248
batchJ 56.551439 60.615248

These intervals seem to jibe with the results from effects.

Addendum 2

Another alternative is the lsmeans package. It obtains degrees of freedom and an adjusted covariance matrix from the pbkrtest package.

> library("lsmeans")
> lsmeans(fm1, "batch")
Loading required namespace: pbkrtest
 batch   lsmean       SE    df lower.CL upper.CL
 A     62.26667 1.125709 40.45 59.99232 64.54101
 B     59.30000 1.125709 40.45 57.02565 61.57435
 C     62.05000 1.125709 40.45 59.77565 64.32435
 D     59.70000 1.125709 40.45 57.42565 61.97435
 E     55.90000 1.125709 40.45 53.62565 58.17435
 F     61.03333 1.125709 40.45 58.75899 63.30768
 G     59.90000 1.125709 40.45 57.62565 62.17435
 H     63.11667 1.125709 40.45 60.84232 65.39101
 I     58.68333 1.125709 40.45 56.40899 60.95768
 J     58.58333 1.125709 40.45 56.30899 60.85768

Confidence level used: 0.95 

These are even more in line with the effect results: the standard errors are identical, but effect uses different d.f. The confint results in Addendum 1 are even narrower than the asymptotic ones based on using $\pm1.96\times\mbox{se}$. So now I think those are not very trustworthy.

Results from effect and lsmeans are similar, but with an unbalanced multi-factor situation, lsmeans by default averages over unused factors with equal weights, whereas effect weights by the observed frequencies (available as an option in lsmeans).