如何为glmmTMB拟合的广义泊松模型预测概率？

问题 #

使用pscl::predprob()函数获取glmmTMB拟合的广义泊松模型（model.gp）的预测概率时会报错，因为predprob()仅支持pscl包原生的模型对象（如zeroinfl、hurdle），不兼容glmmTMB模型。

解决方案 #

手动计算广义泊松模型的预测概率，步骤如下：

• 确定需要计算概率的计数范围（与其他模型保持一致，0到观测数据的最大值）
• 对每个观测，基于模型的线性预测值，用广义泊松的概率质量函数计算每个计数的概率
• 对所有观测的概率取均值，得到与predprob()输出一致的结果

修改后的完整代码 #

# MODEL COMPARISONS
library(pscl)
library(glmmTMB)
library(dplyr)
library(tidyr)
library(ggplot2)

# 原泊松和负二项模型的预测概率（保留原代码）
po.p <- predprob(model.p) %>% colMeans
po.nb <- predprob(model.nb) %>% colMeans

# 手动计算glmmTMB广义泊松模型的预测概率
count_range <- 0:max(na.omit(bo$Grav))
# 获取模型的线性预测值（link尺度）
lin_pred <- predict(model.gp, type = "link")
# 对每个观测，计算所有count的概率
gp_probs <- sapply(count_range, function(k) {
  dgenpois(k, mu = exp(lin_pred), theta = model.gp$fit$par["theta"])
})
# 取列均值，得到每个count的平均预测概率
po.gp <- colMeans(gp_probs)

# 后续数据合并和绘图代码（保留原结构）
df <- data.frame(x = count_range, Poisson = po.p, 
                 NegBin = po.nb, Gen_Poisson = po.gp)

obs <- table(bo$Grav) %>% prop.table() %>% data.frame # Observed
names(obs) <- c("x", 'Observed')

p1 <- predict(linear) %>% round() %>% table %>% prop.table %>% data.frame # for OLS
names(p1) <- c('x', 'OLS')

tmp <- merge(p1, obs, by = 'x', all = T)
tmp$x <- as.numeric(as.character(tmp$x))

comb <- merge(tmp, df, by = 'x', all = T)
comb[is.na(comb)] <- 0

comb2 <- comb[1:11, ] # just for the first 11 results, including zero

mm <- melt(comb2, id.vars = 'x', value.name = 'prob', variable.name = 'Model')
mm <- filter(mm, Model != "OLS") # can include the linear model too if you want
# the SAS note does not, so I am not including it

ggplot(mm, aes(x = x, y = prob, group = Model, col = Model)) +
  geom_line(aes(lty = Model), lwd = 1) +
  theme_bw() +
  labs(x = "Number of pregnancies", y = 'Probability',
       title = "Models for number of pregnancies") +
  scale_color_manual(values = c('black', 'blue', 'red', 'green')) +
  scale_linetype_manual(values = c('solid', 'dotted', 'dotted', 'dotted')) +
  theme(legend.position=c(.2, .80), axis.title.y = element_text(angle = 0))

注意事项 #

• 确保已加载glmmTMB包，因为dgenpois()是该包提供的广义泊松概率函数
• theta参数直接从模型对象中提取，对应广义泊松的离散参数
• 计算逻辑与predprob()一致：先得到每个观测的各计数概率，再对所有观测取均值，确保结果可与其他模型直接对比

内容的提问来源于stack exchange，提问作者Katlego Motsatsing