混合模型中估计值的置信区间

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

在predictInterval中将which设置为"fixed"应该就足够了，但事实并非如此。所以，它看起来像一个错误。但是，如果我们为分组变量提供任何值，则连同此参数，一切正常。

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.

Answer 2

您也可以使用 ggeffects-package (examples, for instance, in 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

混合模型中估计值的置信区间

confidence intervals of estimates in mixed models

r

confidence-interval

mixed-models