0

我正在使用 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")

为什么间隔=“信心”不起作用？这可能是与 glmmTMB 软件包相关的问题吗？

4

3 回答 3

3

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

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

于 2020-08-13T16:39:17.753 回答
0

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

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

predict 是一个通用函数，用于根据各种模型拟合函数的结果进行预测。该函数调用依赖于第一个参数的类的特定方法。

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

进入您正在使用的特定函数的帮助页面的关键predict()是了解您正在使用的模型拟合函数返回的对象的类。例如，适合的模型属于glmmTMBglmmTMB()类，因此您可以转到. 适合的模型属于merMod类，因此您可以转到 如果您不知道模型拟合函数返回的类，看起来您通常可以在“值”部分下的文档中找到信息。至少对于和是这样。?predict.glmmTMBlme4::lmer()?predict.merMod.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
于 2020-08-13T17:21:27.930 回答
0

有一些包可以为你工作，比如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 包（v0.3.0）于 2020 年 9 月 14 日创建

于 2020-09-14T09:39:26.187 回答