Why are there large discrepancies between Wald and bootstrapped confidence intervals for parameters of a lmer model in R?

Question

I am dealing with a multilevel model with gaussian error distribution that has ~21,000 observations and 5000 clusters. The model is of the simple form:

lmer(y ~ a + b + a:b + (1|z), weights=b)

and has weights applied in proportion to the variable b.

The behavior I'm encountering occurs when calculating confidence intervals of the fixed effects parameters with the following:

confint(mod, method="Wald")
confint(mod, method="profile") 
confint(mod1, method="boot", nsim=1000, parm="beta_")

The results from bootstrapping give confidence intervals that are ~3 times wider than the Wald results. The profile results throw a number of warnings such as:

1: In profile.merMod(object, which = parm, signames = oldNames, ...) : non-monotonic profile for (Intercept)

6: In confint.thpr(pp, level = level, zeta = zeta) : bad spline fit for (Intercept): falling back to linear interpolation

I have searched through many old threads that compare these methods, and I do expect the results from these methods to be different. However previous posts (and my own experience) suggests that usually these methods produce results that are not quite so far apart. I am looking for some intuition as to why this might occur and whether I am right to think that the bootstrapped results are probably more realistic (I suppose this means I should be skeptical of the SEs reported in summary() as well?).

I apologize for not providing a reproducible example, as I am not able to share the original data and my attempts to simulate the same issue only lead to situations in which the Wald and bootstrapped CIs were similar.

EDIT: Plot of profile:

Edit #2: In response to Ben Bolker's comment below, the following code (stolen from elsewhere on the internet) simulates data for a mixed model and demonstrates that confint(..., method="profile") fails when weights are included as well as some differences in CIs from the different methods.

library(mvtnorm)
set.seed(2345)

N <- 150
unit.df <- data.frame(unit = c(1:N), a = rnorm(N))
unit.df <-  within(unit.df, {
  E.alpha.given.a <-  1 - 0.15 * a
  E.beta.given.a <-  3 + 0.3 * a
})

q = 0.2
r = 0.9
s = 0.5
cov.matrix <- matrix(c(q^2, r * q * s, r * q * s, s^2), nrow = 2,
                     byrow = TRUE)
random.effects <- rmvnorm(N, mean = c(0, 0), sigma = cov.matrix)
unit.df$alpha <- unit.df$E.alpha.given.a + random.effects[, 1]
unit.df$beta <- unit.df$E.beta.given.a + random.effects[, 2]

J <- 300
M = J * N  #Total number of observations
x.grid = seq(-4, 4, by = 8/J)[0:30]
within.unit.df <-  data.frame(unit = sort(rep(c(1:N), J)), j = rep(c(1:J),N), x =rep(x.grid, N))
flat.df = merge(unit.df, within.unit.df)

flat.df <-  within(flat.df, y <-  alpha + x * beta + 0.75 * rnorm(n = M))
simple.df <-  flat.df[, c("unit", "a", "x", "y")]
simple.df$wht <- rpois(n = dim(simple.df)[1], lambda = 5)+1

my.lmer <-  lmer(y ~ x + a + x * a + wht + (1 | unit), data = simple.df, weights = wht)
summary(my.lmer)

confint(my.lmer, method="Wald")
confint(my.lmer, method="profile")
confint(my.lmer, method="boot", nsim=100)

Answer 1

Basically, the Wald statistic is not good and you shouldn't trust it for mixed models. It uses a much cruder approximation to the actual likelihood than you get with the profile and boot.ci methods. If R (and SAS and JMP and...) would have been written today, they would not have bothered implementing Wald stats. That's why the summary.merMod method intentionally omits $p$-values from the fixed effect coefficient output. The computational intensity of profile/bootstrap is at most on the scale of minutes by today's standards, but in the olden days, it would take weeks. So, the analyst was expected to do massive amounts of testing and variable transformation methods so that Wald stat might have good-ish properties.

EDIT: below is a snippet of a conversation between me, David Dahl, and Douglas Bates back in 2010 when I tried to suggest using the Wald $p$-values for xtable.

A user of your lme4 package would like to use xtable on mer objects from lme4. That means defining a function "xtable.mer". He suggests the implementation below. I regrettably am not very familar with lme4. Do you have any suggestions?

I appreciate Adam's suggestion and his providing an implementation. Regrettably, I think that the implementation would be controversial, to say the least, and I would prefer not to be the recipient of the fallout. There is a long-standing issue with lme4 regarding p-values on tests of the fixed-effects parameters. For linear mixed models there is a widespread belief that you can calculate a t-statistic (what is labelled here as a "z value") and convert it to a p-value by the simple expedient of determining an approximate number of degrees of freedom. In fact, SAS PROC MIXED offers several (6, I believe) different, and incompatible, ways of determining such degrees of freedom and the corresponding p-values. The fact that these give different answers doesn't deter people from regarding such approximations as "absolute truth".

In reality the distribution of such a statistic is not a Student's T. It is much more complicated than that and I advocate other ways of calculating confidence intervals or testing hypotheses. In the case of a generalized linear mixed model I do calculate a p-value from the standard normal distribution, not because the approximation is better for GLMMs than for LMMs but because it is worse.

I am writing a book for Springer on lme4 (chapter drafts are available at http://lme4.R-forge.R-project.org/book/) where I describe using likelihood ratio tests for hypothesis tests and techniques based on profiling the LRT statistic to produce confidence intervals on parameters. The examples in that book are based on the development version of the package which uses a different representation of the model. The implementation is not complete, which is why I haven't released it as lme4, but right now I need to concentrate on the writing because the book is going to be used in a seminar which starts next week.