在R中拟合零膨胀右截断负二项分布及拟合优度检验咨询

当然可以！在R里拟合零膨胀右截断负二项分布完全没问题，我会一步步带你搞定模型拟合和拟合优度检验的步骤，都是实际工作中常用的方法~

我们可以用glmmTMB包来实现这个需求——它支持零膨胀模型、负二项分布，还能直接指定截断参数，非常方便。

首先安装并加载所需的包：

# 安装包（如果还没安装的话）
install.packages(c("glmmTMB", "DHARMa", "pscl"))
# 加载包
library(glmmTMB)
library(DHARMa)
library(pscl)

假设你的数据框是dat，响应变量是y，还有协变量x1、x2（你可以替换成自己的变量）。拟合模型的代码如下：

# 拟合零膨胀右截断负二项分布
model <- glmmTMB(
  y ~ x1 + x2,  # 计数部分的协变量（解释y的非零计数变化）
  ziformula = ~ x1 + x2,  # 零膨胀部分的协变量（解释额外零的产生）
  family = nbinom2,  # 选择负二项分布（nbinom2是方差随均值平方变化的参数化）
  data = dat,
  truncation = upper(2000)  # 指定右截断点为2000
)

关键参数解释： #

• family = nbinom2：如果你的数据方差和均值线性相关，可以换成nbinom1，后续可以通过分散性检验确认哪种更合适。
• ziformula：如果零膨胀部分不需要协变量，可以写成~1，表示只有截距项。
• truncation = upper(2000)：告诉模型我们的响应变量在2000处右截断——也就是所有真实值大于2000的观测都被记录为2000（符合你说的物理上限）。

拟合完成后，查看模型结果：

summary(model)

接下来我们用几种方法来检验模型的拟合效果，从直观诊断到正式检验都有：

2.1 基于模拟残差的直观诊断（DHARMa包） #

DHARMa生成的模拟残差能帮我们直观判断模型是否捕捉到了数据的分布特征，是我最常用的诊断工具：

# 生成1000次模拟的残差
sim_res <- simulateResiduals(model = model, n = 1000)
# 绘制诊断图
plot(sim_res)

这个图包含四个子图：

• QQ图：如果模型拟合得好，点应该沿着直线分布；
• 残差 vs 预测值：如果没有明显趋势，说明模型的均值预测没问题；
• 残差分布：应该接近均匀分布；
• 离群点检验：标记出可能的异常值。

还可以做正式的统计检验：

# 检验残差是否服从均匀分布（拟合优度核心检验）
testUniformity(sim_res)
# 检验是否存在过度/不足分散
testDispersion(sim_res)
# 检验零膨胀是否有必要（如果不确定的话）
testZeroInflation(sim_res)

2.2 皮尔逊卡方检验 #

我们可以把观测值分组，对比观测频数和模型预测的期望频数，用卡方检验来量化拟合程度：

# 定义一个函数，计算零膨胀右截断负二项分布的概率
zi_trunc_nbinom_prob <- function(y, mu, theta, zi_prob, upper_trunc) {
  # 计算非零膨胀部分的负二项概率
  nb_prob <- dnbinom(y, mu = mu, size = theta)
  # 计算截断因子（真实值<=截断点的概率）
  trunc_factor <- pnbinom(upper_trunc, mu = mu, size = theta)
  # 截断后的非零概率
  trunc_nb_prob <- nb_prob / trunc_factor
  # 零膨胀部分的概率：y=0时是结构零概率+截断后的零概率；否则是截断后的非零概率
  ifelse(y == 0,
         zi_prob + (1 - zi_prob) * trunc_nb_prob[y == 0],
         (1 - zi_prob) * trunc_nb_prob)
}

# 对每个观测计算预测概率
dat$pred_prob <- mapply(
  zi_trunc_nbinom_prob,
  y = dat$y,
  mu = predict(model, type = "response"),
  theta = model$fit$par["theta"],
  zi_prob = plogis(predict(model, type = "zprob")),
  upper_trunc = 2000
)

# 分组（根据数据分布调整断点，确保每组期望频数>=5）
y_groups <- cut(dat$y, breaks = c(-1, 0, 100, 500, 1000, 1500, 2000))
# 计算观测频数
obs_counts <- table(y_groups)
# 计算期望频数
exp_counts <- tapply(dat$pred_prob, y_groups, sum)

# 做皮尔逊卡方检验
chisq.test(obs_counts, p = exp_counts / sum(exp_counts))

注意： #

如果分组后有些组的期望频数太小（<5），卡方检验的结果会不可靠，这时候可以调整分组断点，或者优先用DHARMa的模拟检验。

2.3 模型对比（AIC准则） #

如果有备选模型（比如不考虑截断的零膨胀负二项、不考虑零膨胀的截断负二项），可以通过AIC值来比较——AIC越小，模型拟合越好：

# 拟合普通零膨胀负二项（不考虑截断）
model_no_trunc <- zeroinfl(y ~ x1 + x2 | x1 + x2, data = dat, dist = "negbin")
# 拟合截断负二项（不考虑零膨胀）
model_no_zi <- glmmTMB(y ~ x1 + x2, family = nbinom2, data = dat, truncation = upper(2000))

# 比较三个模型的AIC
AIC(model, model_no_trunc, model_no_zi)

如果我们的目标模型（零膨胀+截断）AIC最小，说明加入这两个成分是必要的。

• 先检查数据的零比例：用table(dat$y)看零的占比，如果很高，零膨胀模型才是必要的；
• 模型收敛问题：如果模型报错不收敛，可以尝试简化协变量，或者给glmmTMB加start参数指定初始值；
• 截断点确认：如果你的数据中没有大于2000的值（因为物理上限），这个截断设置完全适用；如果是大于2000的值被删除了（删失），那需要用生存分析的方法处理，但你描述的情况更符合截断。

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