Я сравниваю доверительный интервал (CI) с, полученный с помощью функции sim() руки и predictInterval() из merTools.
Я использую sleepstudy набор данных из lme4 в качестве примера.
Я ожидаю одинакового результата от двух методов, но это не так. В чем принципиальная разница между двумя методами, которые мне не хватает?
Код следующий:
импорт данных испытаний
sleepstudy <- as_tibble(sleepstudy) %>%
mutate(id = rep(1:18, each = 10)) %>%
dplyr::select(id, Days, Reaction) %>%
filter(id <= 16)
многоуровневая модель от lme4
lmerfit <- lmer(Reaction ~ Days + (Days | id), data = sleepstudy)
генерация прогноза
Это для сравнения медианных значений, сгенерированных позже sim и preditInterval .
sleepstudy$predicted <- predict(lmerfit, newdata=sleepstudy, allow.new.levels=T)
ДИ, использующие руку: индивидуальный уровень
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)
CI, использующие merTols
preds <- predictInterval(lmerfit,
newdata = sleepstudy,
n.sims = 1000,
include.resid.var=FALSE,
level=0.95,
stat="median")
sleepstudy <- cbind(sleepstudy, preds)
В качестве примера я строю первые данные вместе с двумя различными прогнозами CI. Черные точки - это данные. Красные точки - это прогнозные значения от lmerfit.
Черная линия и черные пунктирные линии - это медиана и 95% ДИ от arm::sim соответственно.
Красная линия и пунктирная линия - это медиана и 95% ДИ от merTools::predictInterval соответственно.
Прогнозируемые значения и моделируемые медианные значения идентичны, но CI довольно различны. Что может быть причиной? Какой из них точный?
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)