【问题标题】：How to fit confidence intervals using predict function for glmmTMB如何使用 glmmTMB 的预测函数拟合置信区间
【发布时间】：2020-12-03 13:01:06
【问题描述】：

我正在使用 glmmTMB 包运行混合模型，并使用 predict 函数使用以下代码计算预测均值：

运行模型

model_1 <- glmmTMB(Step.rate ~ Treatment*Week + 
    (1|Treatment.Group/Lamb.ID) +  (1|Plot),
     data = data.df, family = nbinom1)

创建新数据框

new.dat <- data.frame(Treatment = data.df$Treatment,
                      Week = data.df$Week, Plot = data.df$Plot, 
                      Treatment.Group = data.df$Treatment.Group,
                      Lamb.ID = data.df$Lamb.ID)

预测平均值

new.dat$prediction <- predict(model_1, new.data = new.dat, 
       type = "response", re.form = NA)

这段代码运行良好，但是当我添加 interval = "confidence" 来计算置信区间时，它似乎不起作用。 R 忽略了代码的最后一部分，只计算预测的均值。

new.dat$prediction <- predict(model_1, new.data = new.dat, 
     type = "response", re.form = NA, intervals = "confidence")

为什么intervals = "confidence" 不起作用？这可能是与 glmmTMB 软件包相关的问题吗？

【问题讨论】：

一些可能有帮助的事情：要查看predict() 的 glmmTMB 特定版本的帮助页面，请转到 ?predict.glmmTMB。你会看到那里没有intervals 参数。您可以“手动”制作近似 CI。 Bolker 的 GLMM 常见问题解答显示了一个示例 here。我还使用了 ggeffects 包中的ggpredict() 来获得近似 CI 和预测，尽管我承认我在尝试定义我想要的确切数据集时会有些沮丧得到预测。
@aosmith，这可能值得回答（特别是如果您澄清通用方法（?predict）、最广泛使用的方法（?predict.lm）和 OP 方法之间的区别正在使用 (?predict.glmmTMB) ...
对于?predict.glmmTMB，您可以选择se.fit=TRUE，这使您构建置信区间的工作变得更加容易。

标签： r predict confidence-interval glmmtmb

【解决方案1】：

有一些包可以为您工作，例如 emmeans 或 ggeffects，或 effects 包（可能还有更多包):

library(ggeffects)
library(glmmTMB)
library(emmeans)
data("Salamanders")
m <- glmmTMB(count ~ spp * mined + sample + (1 | site), family = nbinom1, data = Salamanders)

emmeans(m, c("spp", "mined"), type = "response")
#>  spp   mined response     SE  df lower.CL upper.CL
#>  GP    yes     0.0368 0.0373 627  0.00504    0.269
#>  PR    yes     0.1099 0.0661 627  0.03368    0.358
#>  DM    yes     0.3842 0.1397 627  0.18808    0.785
#>  EC-A  yes     0.1099 0.0660 627  0.03377    0.357
#>  EC-L  yes     0.3238 0.1222 627  0.15437    0.679
#>  DES-L yes     0.4910 0.1641 627  0.25468    0.947
#>  DF    yes     0.5561 0.1764 627  0.29822    1.037
#>  GP    no      2.2686 0.4577 627  1.52646    3.372
#>  PR    no      0.4582 0.1515 627  0.23940    0.877
#>  DM    no      2.4201 0.4835 627  1.63472    3.583
#>  EC-A  no      0.8931 0.2373 627  0.53005    1.505
#>  EC-L  no      3.2017 0.6084 627  2.20451    4.650
#>  DES-L no      3.4921 0.6517 627  2.42061    5.038
#>  DF    no      1.8495 0.3948 627  1.21623    2.813
#> 
#> Confidence level used: 0.95 
#> Intervals are back-transformed from the log scale
ggpredict(m, c("spp", "mined"))
#> 
#> # Predicted counts of count
#> # x = spp
#> 
#> # mined = yes
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      0.04 | 1.01 | [0.01, 0.27]
#> PR   |      0.11 | 0.60 | [0.03, 0.36]
#> DM   |      0.38 | 0.36 | [0.19, 0.78]
#> EC-A |      0.11 | 0.60 | [0.03, 0.36]
#> EC-L |      0.32 | 0.38 | [0.15, 0.68]
#> DF   |      0.56 | 0.32 | [0.30, 1.04]
#> 
#> # mined = no
#> 
#> x    | Predicted |   SE |       95% CI
#> --------------------------------------
#> GP   |      2.27 | 0.20 | [1.53, 3.37]
#> PR   |      0.46 | 0.33 | [0.24, 0.88]
#> DM   |      2.42 | 0.20 | [1.64, 3.58]
#> EC-A |      0.89 | 0.27 | [0.53, 1.50]
#> EC-L |      3.20 | 0.19 | [2.21, 4.65]
#> DF   |      1.85 | 0.21 | [1.22, 2.81]
#> 
#> Adjusted for:
#> * sample = 2.50
#> *   site = NA (population-level)
#> Standard errors are on the link-scale (untransformed).

^{由reprex package (v0.3.0) 于 2020 年 9 月 14 日创建}

【讨论】：

【解决方案2】：

我认为另一个答案为您提供了一种解决方法，可以使用 se.fit 参数为 glmmTMB 对象获取 CI。但是对于不同的对象类型（由对象的 class 定义）具有特定版本的函数的问题在过去让我感到有些悲伤，因此可能值得在这里扩展。

无需过多详细介绍，R 中许多对象类型通用的函数都有通用版本。例如，如果您访问了 ?predict 的文档，您将看到该函数的通用版本的帮助页面。在那里，您将看到一些关于函数通常如何工作的一般性陈述，但几乎没有对特定参数的解释，因为可用的参数取决于您正在使用的对象的类型。通用predict() 帮助页面中的描述：

predict 是一个通用函数，用于根据各种模型拟合功能。该函数调用特定的依赖于第一个参数的类的方法。

特定的模型拟合函数可以有特定版本的predict() 与生成的模型对象一起使用。例如，对于与 lm() 匹配的模型，有一个特定的 predict()。从lm() 返回的对象属于 lm 类。您可以在?predict.lm 上查看有关 lm 对象的函数版本的文档。正是这个函数包含一个intervals 参数，用于计算置信区间和预测区间。虽然我们中的许多人从 lm 对象开始学习 intervals，但事实证明许多（大多数？）其他 predict() 函数没有这个选项。

访问您正在使用的特定predict() 函数的帮助页面的关键是了解您正在使用的模型拟合函数返回的对象的类。例如，符合glmmTMB() 的模型属于glmmTMB 类，因此您可以转到?predict.glmmTMB。适合 lme4::lmer() 的模型属于 merMod 类，因此您可以转到 ?predict.merMod. 如果您不知道模型拟合函数返回的类，看起来您经常可以在Value 部分下的文档。至少 lm() 和 lmer() 是这样。

最后，如果您需要知道某个对象类是否具有与其关联的通用函数的特定版本，您可以使用methods() 函数查看该类可用的方法。 lm 示例：

methods(class = "lm")
 [1] add1           alias          anova          case.names     coerce        
 [6] confint        cooks.distance deviance       dfbeta         dfbetas       
[11] drop1          dummy.coef     effects        extractAIC     family        
[16] formula        hatvalues      influence      initialize     kappa         
[21] labels         logLik         model.frame    model.matrix   nobs          
[26] plot           predict        print          proj           qqnorm        
[31] qr             residuals      rstandard      rstudent       show          
[36] simulate       slotsFromS3    summary        variable.names vcov

【讨论】：

【解决方案3】：

您可以使用参数se.fit = TRUE 来获取预测值的标准误，然后使用这些来计算置信区间。

https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/predict.glmmTMB

【讨论】：

R 区分大小写，它是TRUE，而不是True。
谢谢，这真的很有帮助！