如何使用 glmmTMB 的预测函数来拟合置信区间
How to fit confidence intervals using predict function for glmmTMB
我运行正在使用包 glmmTMB 构建一个混合模型，并使用预测函数使用以下代码计算预测均值：
运行 型号
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)
这段代码工作正常，但是当我添加 intervals = "confidence" 来计算置信区间时，它似乎不起作用。 R忽略代码的最后部分，只计算预测均值。
new.dat$prediction <- predict(model_1, new.data = new.dat,
type = "response", re.form = NA, intervals = "confidence")
为什么 intervals = "confidence" 不起作用？这可能是与包 glmmTMB 相关的问题吗？
您可以使用参数 se.fit = TRUE 获取预测值的标准误差，然后使用这些来计算置信区间。
https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/predict.glmmTMB
我认为另一个答案为您提供了一个 work-around，用于使用 se.fit 参数为 glmmTMB 对象获取 CI。但是，为不同的对象类型（由对象的 class 定义的函数提供特定版本的问题在过去曾让我有些悲伤，所以它可能是值得在这里展开。
无需过多详细说明，R 中在许多对象类型中通用的函数都有通用版本。例如，如果您访问 ?predict 的文档，您将看到该函数的通用版本的帮助页面。在那里你会看到一些关于函数通常如何工作的一般性陈述，但几乎没有对特定参数的解释，因为可用的参数取决于你正在使用的对象的类型。通用 predict():
帮助页面的描述
predict is a generic function for predictions from the results of
various model fitting functions. The function invokes particular
methods which depend on the class of the first argument.
特定模型拟合函数可以具有特定版本的predict()，用于生成的模型对象。例如，对于符合 lm() 的模型，有一个特定的 predict()。从 lm() 返回的对象属于 class lm。您可以在 ?predict.lm 查看 lm 对象函数版本的文档。此函数包含一个 intervals 参数，用于计算置信区间和预测区间。虽然我们中的许多人从 lm 对象开始学习 intervals，但事实证明许多（大多数？）其他 predict() 函数没有此选项。
进入您正在使用的特定 predict() 函数的帮助页面的关键是了解您正在使用的模型拟合函数返回的对象的 class。例如，适合 glmmTMB() 的模型属于 class glmmTMB，因此您可以转到 ?predict.glmmTMB。适合 lme4::lmer() 的模型属于 class merMod 所以你可以去 ?predict.merMod. 如果你不知道 class model-fitting 函数 returns，看起来你经常可以在 Value 部分下的文档中找到信息。至少 lm() 和 lmer() 是这样。
最后，如果您需要知道某个 class 对象是否具有与之关联的通用函数的特定版本，您可以查看该 class 可用的方法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
有些软件包可以为您完成这项工作，例如 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-09-14 创建
我运行正在使用包 glmmTMB 构建一个混合模型，并使用预测函数使用以下代码计算预测均值：
运行 型号
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)
这段代码工作正常，但是当我添加 intervals = "confidence" 来计算置信区间时，它似乎不起作用。 R忽略代码的最后部分，只计算预测均值。
new.dat$prediction <- predict(model_1, new.data = new.dat,
type = "response", re.form = NA, intervals = "confidence")
为什么 intervals = "confidence" 不起作用？这可能是与包 glmmTMB 相关的问题吗？
您可以使用参数 se.fit = TRUE 获取预测值的标准误差，然后使用这些来计算置信区间。
https://www.rdocumentation.org/packages/glmmTMB/versions/1.0.2.1/topics/predict.glmmTMB
我认为另一个答案为您提供了一个 work-around，用于使用 se.fit 参数为 glmmTMB 对象获取 CI。但是，为不同的对象类型（由对象的 class 定义的函数提供特定版本的问题在过去曾让我有些悲伤，所以它可能是值得在这里展开。
无需过多详细说明，R 中在许多对象类型中通用的函数都有通用版本。例如，如果您访问 ?predict 的文档，您将看到该函数的通用版本的帮助页面。在那里你会看到一些关于函数通常如何工作的一般性陈述，但几乎没有对特定参数的解释，因为可用的参数取决于你正在使用的对象的类型。通用 predict():
predict is a generic function for predictions from the results of various model fitting functions. The function invokes particular methods which depend on the class of the first argument.
特定模型拟合函数可以具有特定版本的predict()，用于生成的模型对象。例如，对于符合 lm() 的模型，有一个特定的 predict()。从 lm() 返回的对象属于 class lm。您可以在 ?predict.lm 查看 lm 对象函数版本的文档。此函数包含一个 intervals 参数，用于计算置信区间和预测区间。虽然我们中的许多人从 lm 对象开始学习 intervals，但事实证明许多（大多数？）其他 predict() 函数没有此选项。
进入您正在使用的特定 predict() 函数的帮助页面的关键是了解您正在使用的模型拟合函数返回的对象的 class。例如，适合 glmmTMB() 的模型属于 class glmmTMB，因此您可以转到 ?predict.glmmTMB。适合 lme4::lmer() 的模型属于 class merMod 所以你可以去 ?predict.merMod. 如果你不知道 class model-fitting 函数 returns，看起来你经常可以在 Value 部分下的文档中找到信息。至少 lm() 和 lmer() 是这样。
最后，如果您需要知道某个 class 对象是否具有与之关联的通用函数的特定版本，您可以查看该 class 可用的方法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
有些软件包可以为您完成这项工作，例如 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-09-14 创建