Search code examples
rarmmultilevel-analysis

Why are the coinfidence intervals predicted by arm::sim vs merTools::predictInterval different?

I'm comparing the (CI)s produced by arm's sim() function and predictInterval() from merTools. I'm using the sleepstudy dataset from lme4 as an example. I am expecting the same result from the two methods but that is not the case. What is the fundamental difference between the two methods what I am missing?

The code is the following:

importing test data

sleepstudy <- as_tibble(sleepstudy) %>%
  mutate(id = rep(1:18, each = 10)) %>%
  dplyr::select(id, Days, Reaction) %>%
  filter(id <= 16)

the multilevel model from lme4

lmerfit <- lmer(Reaction ~ Days + (Days | id), data = sleepstudy)

generating prediction

This is to compare the median values generated later by sim and preditInterval.

sleepstudy$predicted <- predict(lmerfit, newdata=sleepstudy, allow.new.levels=T)

CIs using arm: Individual level

sims <- sim(lmerfit, n.sims = 1000)
yhat <- fitted(sims, lmerfit)
sleepstudy$lower <- apply(yhat, 1, quantile, prob=0.025)
sleepstudy$median <- apply(yhat, 1, quantile, prob=0.5)
sleepstudy$upper <- apply(yhat, 1, quantile, prob=0.975)

CIs using merTols

preds <- predictInterval(lmerfit, 
newdata = sleepstudy, 
n.sims = 1000, 
include.resid.var=FALSE, 
level=0.95, 
stat="median")
sleepstudy <- cbind(sleepstudy, preds)

As an example I am plotting the first data together with the two different CI prediction. Black points are the data. Red points are the predicted values from lmerfit. Black line and black dashed lines are the median and 95% CIs from arm::sim respectively. Red line and dashed lines are the median and 95% CIs from merTools::predictInterval respectively.

The predicted values and simulated median values are identical, but the CIs are rather different. What could be the reason? Which one is accurate?

ggplot(data =  filter(sleepstudy, id == 1), aes(x=Days, y=Reaction)) +
  geom_point() +
  geom_point(aes(y=predicted), col = "red") +
  geom_line(aes(y=median), col ="black" ) +
  geom_line(aes(y=lower), col ="black", lty = 2) +
  geom_line(aes(y=upper), col ="black", lty = 2) +
  geom_line(aes(y=fit), col = "red") +
  geom_line(aes(y=lwr), col = "red", lty = 2) +
  geom_line(aes(y=upr), col = "red", lty = 2)
Solution

  • The merTools CRAN page goes into this (https://cran.r-project.org/web/packages/merTools/vignettes/Using_predictInterval.html), making a direct comparison between sim and predictInterval. Basically, my understanding is that sim ignores uncertainty about the random intercepts, using the mode as a point estimate. The intervals for predictInterval are wider because they account for this additional uncertainty, and are therefore probably more realistic.