доверительные интервалы оценок в смешанных моделях

Question

Я могу получить прогнозируемые значения для смешанной модели, например:

mod <- lmer(sales1 ~ price1 + (1|store), oranges)
X <- with(oranges, expand.grid(price1=c(30,50,70)))
X$pred <- predict(mod, newdata=X, re.form=NA)

> X
      price1      pred
    1     30 23.843916
    2     50 11.001901
    3     70 -1.840114

, но как я могу получить нижний и верхний доверительные интервалы этих трех оценок?

Я установил merTools пакет и попытался

predictInterval(mod, newdata = X, n.sims = 999) 

, но получил ошибку

Error in eval(predvars, data, env) : object 'store' not found

Answer 1

Вы также можете использовать ggeffects-package (например, в this package-vignette ), что экономит ваше время, поскольку вам не нужно создаватьфрейм данных для newdata:

library(ggeffects)
library(lme4)
#> Loading required package: Matrix
data("sleepstudy")
m <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
ggpredict(m, "Days")
#> 
#> # Predicted values of Reaction 
#> # x = Days 
#> 
#>  x predicted std.error conf.low conf.high
#>  0   251.405     6.825  238.029   264.781
#>  1   261.872     6.787  248.570   275.174
#>  2   272.340     7.094  258.435   286.244
#>  3   282.807     7.705  267.705   297.909
#>  5   303.742     9.581  284.963   322.520
#>  6   314.209    10.732  293.174   335.244
#>  7   324.676    11.973  301.210   348.142
#>  9   345.611    14.629  316.939   374.283
#> 
#> Adjusted for:
#> * Subject = 308

# example solution for the case mentioned
# in the comments...
r <- c(2,4,6)
s <- paste0("Days [", toString(sprintf("%s", r)), "]", collapse = "")

ggpredict(m, s)
#> 
#> # Predicted values of Reaction 
#> # x = Days 
#> 
#>  x predicted std.error conf.low conf.high
#>  2   272.340     7.094  258.435   286.244
#>  4   293.274     8.556  276.506   310.043
#>  6   314.209    10.732  293.174   335.244
#> 
#> Adjusted for:
#> * Subject = 308

Answer 2

Установка which в "fixed" в predictInterval должна быть достаточной, но это не так.Итак, это похоже на ошибку.Однако, наряду с этим параметром, если мы предоставим какое-либо значение для переменной группировки, все будет работать.

library(lme4)
library(merTools)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
X1 <- data.frame(Reaction = 250, Days = 4, Subject = 309)
predictInterval(fm1, newdata = X1, n.sims = 999, seed = 1)
#        fit      upr      lwr
# 1 216.8374 256.8839 181.1969
X1 <- data.frame(Reaction = 250, Days = 4, Subject = 310)
predictInterval(fm1, newdata = X1, n.sims = 999, seed = 1)
#       fit      upr      lwr
# 1 230.959 271.0055 195.3185

Как и ожидалось, разные субъекты дают разные прогнозы.Однако установка which на "fixed" помогает:

X1 <- data.frame(Reaction = 250, Days = 4, Subject = 309)
predictInterval(fm1, newdata = X1, n.sims = 999, seed = 1, which = "fixed")
#        fit      upr      lwr
# 1 291.9062 328.5429 256.2472
X1 <- data.frame(Reaction = 250, Days = 4, Subject = 310)
predictInterval(fm1, newdata = X1, n.sims = 999, seed = 1, which = "fixed")
#        fit      upr      lwr
# 1 291.9062 328.5429 256.2472

Значение группировки даже не должно быть значимым, так как в конечном итоге оно игнорируется:

X1 <- data.frame(Reaction = 250, Days = 4, Subject = -1)
predictInterval(fm1, newdata = X1, n.sims = 999, seed = 1, which = "fixed")
#        fit      upr      lwr
# 1 291.9062 328.5429 256.2472
# Warning message:
#      The following levels of Subject from newdata 
#  -- -1 -- are not in the model data. 
#      Currently, predictions for these values are based only on the 
#  fixed coefficients and the observation-level error.