如何基于glmmTMB截断nbinom2模型的连接尺度进行预测？

我在用glmmTMB拟合零截断负二项2（truncated_nbinom2）族模型，想先在连接尺度上预测并计算95%置信区间，再转换到响应尺度，但不知道怎么针对这个截断族实现操作。

以下是示例代码：

library(dplyr)
library(glmmTMB)

set.seed(1)
df <- data.frame(Group = rep(c("a", "b"), each = 20))

## 生成零截断负二项数据的粗糙方法
tnb <- rep(0, 40)
z <- (tnb==0)
while(any(z)) {
    tnb[z] <- rnbinom(sum(z), mu = 1, size = 1)
    z <- (tnb==0)
}
df$Nnb <- tnb
## 分组统计均值
df %>% group_by(Group) %>% summarize(across(starts_with("N"), mean))

m <- glmmTMB(Nnb ~ Group, data = df, family = "truncated_nbinom2")
df %>% group_by(Group) %>% summarize(mean(Nnb))
predict(m, newdata = data.frame(Group = c("a", "b")), type = "response")

对于非截断模型，我可以通过连接尺度的条件值乘以非零概率得到响应尺度结果，但零截断模型这么做的结果和predict(..., type="response")的输出不匹配。请问怎么针对truncated_nbinom2族实现类似下面的操作？

mu <- predict(m, newdata = data.frame(Group = c("a", "b")), type = "link", se.fit = TRUE)
zi <- predict(m, newdata = data.frame(Group = c("a", "b")), type = "zlink", se.fit = TRUE)

Pred.response <- exp(mu$fit)*(1 - plogis(zi$fit))
Lower.response <- exp(mu$fit - 1.96*mu$se.fit) * (1 - plogis(zi$fit - 1.96*zi$se.fit))
Upper.response <- exp(mu$fit + 1.96*mu$se.fit) * (1 -plogis(zi$fit + 1.96*zi$se.fit))

Pred.response

零截断负二项2模型的响应尺度均值不能直接用非截断模型的逻辑计算，必须遵循零截断分布的均值转换规则，核心是用原始分布均值除以（1-原始分布取0的概率）。

核心逻辑 #

零截断负二项分布的均值公式为：
$$\text{截断均值} = \frac{\text{原始负二项均值}}{1 - P(Y=0)}$$
其中：

• $\text{原始负二项均值}$：连接尺度预测值转换后的值（即exp(连接尺度拟合值)）
• $P(Y=0)$：原始（未截断）负二项分布取0的概率

具体实现代码 #

library(dplyr)
library(glmmTMB)

# 复用示例代码拟合模型
set.seed(1)
df <- data.frame(Group = rep(c("a", "b"), each = 20))
tnb <- rep(0, 40)
z <- (tnb==0)
while(any(z)) {
    tnb[z] <- rnbinom(sum(z), mu = 1, size = 1)
    z <- (tnb==0)
}
df$Nnb <- tnb
m <- glmmTMB(Nnb ~ Group, data = df, family = "truncated_nbinom2")

# 1. 获取连接尺度的预测值与标准误
mu_link <- predict(m, newdata = data.frame(Group = c("a", "b")), type = "link", se.fit = TRUE)
# 提取负二项模型的离散参数theta
theta <- m$fit$par["theta"]

# 2. 计算原始负二项分布的零概率（分别对应拟合值、置信上下限）
p0_fit <- dnbinom(0, mu = exp(mu_link$fit), size = theta)
p0_lower <- dnbinom(0, mu = exp(mu_link$fit - 1.96*mu_link$se.fit), size = theta)
p0_upper <- dnbinom(0, mu = exp(mu_link$fit + 1.96*mu_link$se.fit), size = theta)

# 3. 转换到响应尺度并计算95%置信区间
Pred.response <- exp(mu_link$fit) / (1 - p0_fit)
Lower.response <- exp(mu_link$fit - 1.96*mu_link$se.fit) / (1 - p0_lower)
Upper.response <- exp(mu_link$fit + 1.96*mu_link$se.fit) / (1 - p0_upper)

# 验证：与直接predict(type="response")结果一致
all.equal(Pred.response, predict(m, newdata = data.frame(Group = c("a", "b")), type = "response"))
# 输出结果
data.frame(Group = c("a", "b"),
           Predicted = Pred.response,
           Lower_95CI = Lower.response,
           Upper_95CI = Upper.response)

关键注意点 #

• truncated_nbinom2是零截断模型，不是零膨胀模型，因此没有零膨胀参数，不需要处理zlink相关的预测结果。
• 置信区间的计算必须在连接尺度完成上下限推导后，再转换到响应尺度应用截断公式，不能直接对响应尺度均值加减标准误（转换是非线性的，直接操作会导致偏差）。

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