我正在尝试计算零通胀的 glmmTMB 模型预测的置信区间。
我浏览了 github 上发布的一些问题以及描述 glmmTMB 的原始论文。然而,glmmTMB 似乎发生了轻微的变化,我对正确的使用方法感到非常困惑。
这是我到目前为止所看到的:
-https://github.com/glmmTMB/glmmTMB/issues/324 https://github.com/glmmTMB/glmmTMB/issues/324
- https://github.com/glmmTMB/glmmTMB/issues/378 https://github.com/glmmTMB/glmmTMB/issues/378
- https://www.biorxiv.org/content/biorxiv/suppl/2017/05/01/132753.DC1/132753-2.pdf https://www.biorxiv.org/content/biorxiv/suppl/2017/05/01/132753.DC1/132753-2.pdf
目前,我只是在响应尺度上进行预测和 se(因为这也应该包括 ZI 项),转换回对数尺度,计算 0.95 置信区间并转换回响应尺度。 (参见下面的代码)。
library(glmmTMB)
data("Salamanders")
fit = glmmTMB(count~spp + cover + mined +(1|site),
ziformula=~spp + mined,
dispformula = ~DOY,
data = Salamanders,
family=nbinom2)
newdata <- expand.grid(
site = "new",
spp = unique(Salamanders$spp),
mined = factor(c("yes", "no"), levels = levels(Salamanders$mined)),
cover = seq(from = min(Salamanders$cover),
to = max(Salamanders$cover),
length = 25
),
DOY = mean(Salamanders$DOY)
)
preds <- predict(fit, newdata, se=T, allow.new.levels = T, type='response')
newdata$pred=preds$fit
newdata$se = preds$se.fit
newdata$ulimit = exp(log(newdata$pred)+(qnorm(0.975)*log(newdata$se)))
newdata$llimit = exp(log(newdata$pred)-(qnorm(0.975)*log(newdata$se)))
library(ggplot2)
ggplot(data=newdata, aes(x=cover, y = pred))+
geom_ribbon(aes(ymin=llimit, ymax=ulimit, fill=mined), alpha = 0.25)+
geom_line(aes(color = mined), size=1)+
facet_wrap(~spp)
由于文档显示了执行此操作的其他方法,因此我想保证我正确计算了上限和下限。