如何使用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) 

这段代码可以正常工作,但是当我添加interval =“ confidence”来计算置信区间时,它似乎不起作用。 R忽略了代码的最后一部分,只计算了预测的均值。

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

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

解决方法

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

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

,

我认为另一个答案为您提供了一种解决方法,可以使用se.fit参数获取glmmTMB对象的CI。但是,针对不同的对象类型(由对象的 class 定义)具有特定版本的功能的问题在过去让我有些痛苦,因此可能值得在此进行扩展。 >

R中的函数在许多对象类型中都是通用的,因此无需过多详细说明。例如,如果您获得?predict的文档,则将看到该函数的通用版本的帮助页面。在那里,您将看到一些有关该函数正常工作方式的一般说明,但对特定参数的解释很少甚至没有解释,因为可用的参数取决于您使用的对象的类型。通用predict()的帮助页面中的描述:

predict是用于根据以下结果进行预测的通用函数 各种模型拟合功能。该函数调用特定的 取决于第一个参数的类的方法。

特定的模型拟合函数可以具有与生成的模型对象一起使用的特定版本的predict()。例如,对于符合predict()的模型,有一个特定的lm()。从lm()返回的对象属于 lm 类。您可以在?predict.lm上查看lm对象功能版本的文档。该函数包含一个intervals自变量,用于计算置信度和预测间隔。虽然我们很多人都是从lm对象开始的,所以他们学习intervals,但事实证明,很多(大多数?)其他predict()函数没有此选项。

进入正在使用的特定predict()函数的帮助页面的关键是了解所使用的模型拟合函数返回的对象的类。例如,符合glmmTMB()的模型属于 glmmTMB 类,因此您可以转到?predict.glmmTMB。符合lme4::lmer()的模型属于 merMod 类,因此您可以转到?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      
,

有些软件包可以为您完成工作,例如 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,"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创建

相关问答

依赖报错 idea导入项目后依赖报错,解决方案:https://blog....
错误1:代码生成器依赖和mybatis依赖冲突 启动项目时报错如下...
错误1:gradle项目控制台输出为乱码 # 解决方案:https://bl...
错误还原:在查询的过程中,传入的workType为0时,该条件不起...
报错如下,gcc版本太低 ^ server.c:5346:31: 错误:‘struct...