A few more points to add to @Nate's answer.
The idea that the bootstrap function from Hmisc (which is what mean_cl_boot uses) is wrong because it doesn't take the grouping structure into account is basically correct.
I modified your fitting function slightly to make it more convenient to look at the confidence intervals for species A (suppressing the intercept by including -1 in the formula. I also tried it with and without lmerTest, for the purpose of making some comparisons discussed in more detail below.
library(lme4)
mod0 <- lmer(oviposition.index ~ species-1 + (1|month/plot), df)
library(lmerTest)
mod1 <- as(mod0, "lmerModLmerTest")
library(broom.mixed)
f <- function(m, mod = mod0, ...) {
tt <- tidy(mod, conf.int = TRUE, effects = "fixed", conf.method = m, ...)
as.data.frame(tt)[1, c("estimate", "conf.low", "conf.high")]
}
ctab <- rbind(
hmboot = Hmisc::smean.cl.boot(oviposition.index[1:10]),
hmwald = Hmisc::smean.cl.normal(oviposition.index[1:10]),
wald = f("Wald"),
wald_t_satt = f("Wald", mod1),
wald_t_kr = f("Wald", mod1, ddf.method = "Kenward-Roger"),
profile = f("profile"),
pboot = f("boot")
)
print(ctab,digits =3)
- a Wald test is based on the estimated curvature of the likelihood surface at the ML estimate. It's usually fastest but least accurate; it always gives symmetric CIs. It can be based on the assumption of a Normal sampling distribution of the estimate or based on a t-distribution; if the latter, then you need to specify some method of approximating the 'degrees of freedom' parameter of the t distribution.
- profile likelihood is based on measuring the whole likelihood surface. It's more reliable (and slower) than Wald, but doesn't take small sample sizes into account
- parametric bootstrap is the most reliable, but slowest method. It is based on simulating new data sets from the model.
Conclusions here are that the methods all give approximately the same estimates for the CI. The naive bootstrap (as you've used above) gives the (slightly) narrowest CIs, and the Wald estimate with Kenward-Roger degrees of freedom gives the widest (probably overconservative, as the parametric bootstrap (pboot) probably gives the best answer). (The Satterthwaite ddf approximation completely breaks down in this example.)
estimate conf.low conf.high
hmboot 1.13 0.4397 1.68 ## naive bootstrap
hmwald 1.13 0.4005 1.86 ## naive Wald (t-distrib)
wald_lmer 1.13 0.4082 1.85 ## mixed-model Wald (Z-distrib)
wald_t_satt 1.13 NaN NaN ## mixed-model Wald (Satterthwaite)
wald_t_kr 1.13 0.0586 2.20 ## mixed-model Wald (Kenward-Roger)
profile 1.13 0.3600 1.90 ## likelihood profile CI
pboot 1.13 0.4111 1.82 ## parametric bootstrap CI
If we get a little fancier (code below) we can get CIs for both groups:
library(Hmisc)
f <- function(m, mod = mod0, w = 1:2, ...) {
tt <- tidy(mod, conf.int = TRUE, effects = "fixed", conf.method = m, ...)
tt[1:2, c("term","estimate", "conf.low", "conf.high")]
}
h <- function(sfun) {
tab <- do.call(rbind, lapply(split(df, species),
function(d) sfun(d$oviposition.index)))
tab <- data.frame(term = paste0("species", c("A","B")),
setNames(as.data.frame(tab), c("estimate", "conf.low", "conf.high")))
return(tab)
}
h(smean.cl.normal)
tab2 <- dplyr::bind_rows(list(
hmisc_boot = h(smean.cl.boot),
hmisc_normal = h(smean.cl.normal),
wald_lmer = f("Wald"),
wald_t_satt = f("Wald", mod1),
wald_t_kr = f("Wald", mod1, ddf.method = "Kenward-Roger"),
profile = f("profile"),
boot = f("boot")),
.id = "method")
tab2$method <- factor(tab2$method, levels = unique(tab2$method))
ggplot(tab2, aes(x=term, y = estimate, colour = method)) +
geom_pointrange(aes(ymin=conf.low, ymax = conf.high), position = position_dodge(width=0.25)) +
geom_hline(yintercept = 1, lty = 2)
