首页 文章

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

提问于
浏览
1

我可以得到这样的混合模型的预测值:

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

2 回答

  • 0

    predictInterval 中将 which 设置为 "fixed" 就足够了,但事实并非如此 . 所以,它看起来像一个bug . 但是,如果我们为分组变量提供任何值,那么除了这个参数外,一切都有效 .

    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.
    
  • 0

    您也可以使用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
    

相关问题