The V Lab
适用对象公共卫生、流行病学、医学与健康数据科学学习者
学习时长约 150–210 分钟
先修要求描述统计、线性/逻辑回归基础与基本 R 语法

本教程的数据与目标 所有数据均由固定随机种子模拟：2,400 名成人在不同人年中发生的呼吸系统急诊（ED）就诊次数。部分人不在参与医疗机构的资料捕获范围内，因此必为零；其他人即使在范围内，也可能只是观察期内恰好没有就诊。这个刻意设置让我们区分结构零和抽样零，但绝不把模拟中的机制直接当成真实世界结论。

先问研究设计，再选模型 Poisson、负二项、hurdle 和 ZIP 都是结局分布模型；它们不会自动处理混杂、重复测量、选择进入队列、缺失或错误的随访起点。若“预防项目”不是随机分配，本页调整后的 IRR 仍是条件关联，不能仅凭模型名称解读为因果效应。

如何使用本教程

建议沿着以下链条学习：

计数问题与分母 → Poisson 模型 → 率比和绝对预测 → 离散度 → 负二项/稳健推断 → 零的机制 → ZIP 估计与预测 → 频数校准 → 报告与研究边界

学习目标

完成本教程后，你应能够：

  • 区分事件次数、事件率、二分类风险和随访暴露时间；
  • 用 offset(log(person_years)) 建立以人年为分母的 Poisson 率模型；
  • 把回归系数转换为 IRR、置信区间，并报告有政策意义的绝对预测；
  • 用 Pearson 离散度评估 Poisson 的均值—方差约束，理解过度离散的常见来源；
  • 在 glm() Poisson、quasi-Poisson 和 MASS::glm.nb() 之间作有目的的选择；
  • 解释“零很多”为什么不是 ZIP 的充分证据；
  • 写出 ZIP 的双部分似然，并使用 base R 的 optim() 实际估计一个 ZIP 模型；
  • 分别解释在风险组的条件均值 μ\mu、结构零概率 π\pi、边际均值及零概率；
  • 用观测—预测频数、零比例和分组校准检查模型，而不迷信单一 AIC；
  • 说明 hurdle、聚类模型和因果设计与普通 ZIP 的不同职责；
  • 写出可复核的分析流程、结果段落与敏感性分析清单。

1 计数、率、风险：先把问题说清楚

1.1 同一个“零”可回答不同的问题

令 YiY_i 为第 ii 人在观察窗内的呼吸系统 ED 就诊次数，tit_i 为实际观察到的人年。

量 例子 合适的分母/模型 常见误解
次数 一年内 ED 就诊 0、1、2 次 计数模型 把 2 次当成“有病”二元变量
率 每人年 ED 就诊次数 offset(log(t)) 的 Poisson/负二项 忽略 0.3 年与 3 年随访不同
风险 一年内是否至少一次就诊 二项/逻辑或风险模型 将“至少一次”与次数混用
强度/发生率 瞬时事件过程 生存/复发事件方法 把删失当成完整暴露

本页的主要 estimand 是条件发生率比：在相同年龄、吸烟、疾病严重度与居住地条件下，预防项目组的单位人年期望 ED 就诊率与常规服务组之比。它不是“至少一次就诊”的风险比，也不是个体层面的必然变化。

1.2 为什么 offset 不是一个普通协变量

最常用的 Poisson 率模型是：

Yi∣Xi∼Poisson⁡(μi),log⁡(μi)=log⁡(ti)+XiTβ. Y_i\mid X_i \sim \operatorname{Poisson}(\mu_i),\qquad \log(\mu_i)=\log(t_i)+X_i^T\beta.

因此 E(Yi∣Xi)/ti=exp⁡(XiTβ)E(Y_i\mid X_i)/t_i=\exp(X_i^T\beta) 是单位人年的期望率。offset(log(person_years)) 的系数固定为 1；若把 log(person_years) 当普通自变量，其系数会被数据估计，回答的是另一问题，通常没有研究设计上的理由。

overview <- within(ed_data, {
  rate_per_py <- ed_visits / person_years
})
summary_table <- data.frame(
  指标 = c("样本量", "总人年", "总 ED 就诊", "观察到的零比例", "粗发生率（每人年）"),
  数值 = c(
    nrow(overview), sum(overview$person_years), sum(overview$ed_visits),
    mean(overview$ed_visits == 0), sum(overview$ed_visits) / sum(overview$person_years)
  )
)
knitr::kable(summary_table, digits = 3, caption = "模拟队列的基本规模与未调整发生率")
模拟队列的基本规模与未调整发生率
指标 数值
样本量 2400.000
总人年 2926.229
总 ED 就诊 1717.000
观察到的零比例 0.605
粗发生率（每人年） 0.587 
by_program <- aggregate(
  cbind(ed_visits, person_years) ~ program, data = ed_data, FUN = sum
)
by_program$粗发生率 <- by_program$ed_visits / by_program$person_years
knitr::kable(by_program, digits = 3, caption = "按项目分组的粗次数、人年与发生率")
按项目分组的粗次数、人年与发生率
program ed_visits person_years 粗发生率
常规服务 941 1402 0.671
预防项目 776 1525 0.509

分母审计 在拟合前核对：人年是否从同一个时间零点累计？死亡、迁出、失访和行政截止是否正确截断？结局发生是否会缩短可观察时间？若答案是否定或不确定，先修复队列定义；模型无法弥补错误分母。

2 Poisson 回归：从系数到可解释结果

2.1 拟合一个率模型

poisson_fit <- glm(
  ed_visits ~ program + smoking + severity + rural + age10 +
    offset(log(person_years)),
  family = poisson(link = "log"), data = ed_data
)
summary(poisson_fit)
## 
## Call:
## glm(formula = ed_visits ~ program + smoking + severity + rural + 
##     age10 + offset(log(person_years)), family = poisson(link = "log"), 
##     data = ed_data)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -0.5881     0.0411  -14.30  < 2e-16 ***
## program预防项目  -0.2538     0.0485   -5.23  1.7e-07 ***
## smoking吸烟       0.3709     0.0521    7.12  1.1e-12 ***
## severity          0.3675     0.0243   15.12  < 2e-16 ***
## rural农村        -0.1756     0.0521   -3.37  0.00076 ***
## age10             0.1243     0.0208    5.97  2.3e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3492.4  on 2399  degrees of freedom
## Residual deviance: 3110.8  on 2394  degrees of freedom
## AIC: 5400
## 
## Number of Fisher Scoring iterations: 6
poisson_coef <- summary(poisson_fit)$coefficients
poisson_irr <- data.frame(
  变量 = rownames(poisson_coef),
  IRR = exp(poisson_coef[, "Estimate"]),
  下限 = exp(poisson_coef[, "Estimate"] - 1.96 * poisson_coef[, "Std. Error"]),
  上限 = exp(poisson_coef[, "Estimate"] + 1.96 * poisson_coef[, "Std. Error"]),
  P值 = poisson_coef[, "Pr(>|z|)"]
)
rownames(poisson_irr) <- NULL
poisson_effects <- poisson_irr[poisson_irr[["变量"]] != "(Intercept)", ]
knitr::kable(poisson_effects, digits = 3, caption = "Poisson 协变量率比（IRR）及 95% 置信区间")
Poisson 协变量率比（IRR）及 95% 置信区间
变量 IRR 下限 上限 P值
2 program预防项目 0.776 0.705 0.853 0.000
3 smoking吸烟 1.449 1.308 1.605 0.000
4 severity 1.444 1.377 1.515 0.000
5 rural农村 0.839 0.757 0.929 0.001
6 age10 1.132 1.087 1.180 0.000

截距未列入 IRR 表；其指数化值是所有参考水平、连续协变量为 0 时的基线率，而不是协变量效应比。

在其他协变量相同且人年相同时，program预防项目 的 exp⁡(β)\exp(\beta) 是预防项目相对常规服务的条件 IRR。例如 IRR 为 0.80，表示模型中的期望发生率低 20%，而不是“每个人少 0.20 次”，也不必然是因果效果。

2.2 IRR 必须回到绝对尺度

相对量需要参照风险（基线率）才容易用于资源规划。以下先预测一个明确画像在 1 人年内的期望次数，再把该率线性换算为每 100 人年的期望次数。这里的普通 Poisson 预测假设所有人都来自同一个 Poisson 过程；随后 ZIP 会放宽这个假设。

profile <- data.frame(
  program = factor(c("常规服务", "预防项目"), levels = levels(ed_data$program)),
  smoking = factor("不吸烟", levels = levels(ed_data$smoking)),
  severity = 0,
  rural = factor("城市", levels = levels(ed_data$rural)),
  age10 = 0, # 55 岁
  access_barrier = factor("无明显障碍", levels = levels(ed_data$access_barrier)),
  person_years = 1
)
profile$预测次数_每人年 <- predict(poisson_fit, newdata = profile, type = "response")
profile$预测次数_每100人年 <- 100 * profile$预测次数_每人年
profile$项目 <- profile$program
knitr::kable(profile[, c("项目", "预测次数_每人年", "预测次数_每100人年")], digits = 2,
             caption = "指定画像下 Poisson 的绝对预测")
指定画像下 Poisson 的绝对预测
项目 预测次数_每人年 预测次数_每100人年
常规服务 0.56 55.5
预防项目 0.43 43.1 
set_count_plot_font()
barplot(
  height = profile$预测次数_每100人年,
  names.arg = profile$项目, col = c(palette_count["gray"], palette_count["teal"]),
  ylab = "预测 ED 就诊次数（每 100 人年）", ylim = c(0, max(profile$预测次数_每100人年) * 1.2),
  main = "绝对预测须说明人群画像与暴露时间"
)
柱状图比较常规服务与预防项目在同一参考画像下的预测 ED 就诊次数，纵轴为每一百人年次数。

指定画像中两种项目的 Poisson 预测 ED 就诊次数（每 100 人年）

2.3 置信区间与模型尺度

predict(..., type = "response") 给出均值预测，不自动回答所有不确定性问题。对线性预测子 η=XTβ̂+log⁡(t)\eta=X^T\hat\beta+\log(t)，可用 predict(..., type = "link", se.fit = TRUE) 在 log 均值尺度构造区间，再指数变换。它是平均次数的模型不确定性区间，不是单个新人的 prediction interval；后者还包含随机事件变异和参数不确定性。

link_pred <- predict(poisson_fit, newdata = profile, type = "link", se.fit = TRUE)
profile$均值下限 <- exp(link_pred$fit - 1.96 * link_pred$se.fit)
profile$均值上限 <- exp(link_pred$fit + 1.96 * link_pred$se.fit)
profile$每100人年均值下限 <- 100 * profile$均值下限
profile$每100人年均值上限 <- 100 * profile$均值上限
knitr::kable(profile[, c("项目", "预测次数_每人年", "均值下限", "均值上限",
                         "预测次数_每100人年", "每100人年均值下限", "每100人年均值上限")],
             digits = 2, caption = "Poisson 预测均值的近似 95% 置信区间")
Poisson 预测均值的近似 95% 置信区间
项目 预测次数_每人年 均值下限 均值上限 预测次数_每100人年 每100人年均值下限 每100人年均值上限
常规服务 0.56 0.51 0.60 55.5 51.2 60.2
预防项目 0.43 0.40 0.47 43.1 39.6 46.9

3 Poisson 假设与过度离散

3.1 均值等于方差是可检验的工作假设

Poisson 要求条件方差等于条件均值：Var⁡(Yi∣Xi)=μi\operatorname{Var}(Y_i\mid X_i)=\mu_i。Pearson 离散度是常用快速诊断：

ϕ̂P=∑i(yi−μ̂i)2/μ̂in−p. \hat\phi_P=\frac{\sum_i (y_i-\hat\mu_i)^2/\hat\mu_i}{n-p}.

约为 1 与 Poisson 相容；明显大于 1 说明残差比假设更分散。它不是机械阈值：大样本下轻微偏离也会显著，且遗漏非线性、相关性或零机制都会影响它。

pearson_resid <- residuals(poisson_fit, type = "pearson")
phi_pearson <- sum(pearson_resid^2) / df.residual(poisson_fit)
dispersion_table <- data.frame(
  指标 = c("Pearson 卡方", "残差自由度", "Pearson 离散度"),
  数值 = c(sum(pearson_resid^2), df.residual(poisson_fit), phi_pearson)
)
knitr::kable(dispersion_table, digits = 2, caption = "Poisson 的 Pearson 离散度诊断")
Poisson 的 Pearson 离散度诊断
指标 数值
Pearson 卡方 3262.22
残差自由度 2394.00
Pearson 离散度 1.36

常见的过度离散来源包括：未测量的易感性差异、遗漏交互或非线性、同一诊所/家庭的聚类、重复事件的依赖性、混合人群，以及零过多。相反，不能因为 ϕ̂P>1\hat\phi_P>1 就断言“必须 ZIP”；先检查数据质量、均值结构和聚类设计。

set_count_plot_font()
plot(
  fitted(poisson_fit), pearson_resid, pch = 16, cex = 0.55,
  col = grDevices::adjustcolor(palette_count["navy"], alpha.f = 0.35),
  xlab = "Poisson 拟合的期望次数", ylab = "Pearson 残差",
  main = "残差图不是正式的分布检验，但能暴露模式"
)
abline(h = 0, lty = 2, col = palette_count["vermillion"])
散点图横轴为 Poisson 拟合的期望 ED 次数，纵轴为 Pearson 残差，并有零水平参考线。

Poisson Pearson 残差与拟合均值：用于发现均值结构或离散度问题

3.2 quasi-Poisson 与负二项：解决的不是同一件事

quasi-Poisson 保留 E(Y∣X)=μE(Y\mid X)=\mu 和 log 链接，但把方差设为 ϕμ\phi\mu，主要调整标准误；它没有完整似然，因此不能与 Poisson 用 AIC 比较。负二项常设：

Var⁡(Yi∣Xi)=μi+μi2/θ, \operatorname{Var}(Y_i\mid X_i)=\mu_i+\mu_i^2/\theta,

其中 θ\theta 越小，额外异质性越大；它改变了分布与预测尾部，能用似然型 AIC。两者都不是对结构零的直接解释。

quasi_fit <- update(poisson_fit, family = quasipoisson(link = "log"))
quasi_program <- summary(quasi_fit)[["coefficients"]]["program预防项目", ]

nb_available <- requireNamespace("MASS", quietly = TRUE)
if (nb_available) {
  nb_fit <- MASS::glm.nb(formula(poisson_fit), data = ed_data)
  nb_program <- summary(nb_fit)[["coefficients"]]["program预防项目", ]
  comparison_table <- data.frame(
    模型 = c("Poisson", "quasi-Poisson", "负二项"),
    项目_IRR = c(
      exp(coef(poisson_fit)["program预防项目"]),
      exp(quasi_program["Estimate"]), exp(nb_program["Estimate"])
    ),
    项目标准误 = c(
      summary(poisson_fit)[["coefficients"]]["program预防项目", "Std. Error"],
      quasi_program["Std. Error"], nb_program["Std. Error"]
    ),
    AIC = c(AIC(poisson_fit), NA, AIC(nb_fit))
  )
  knitr::kable(comparison_table, digits = 3,
               caption = "Poisson、quasi-Poisson 与负二项的比较（quasi 无 AIC）")
} else {
  cat("未安装 MASS，跳过负二项拟合；可运行 install.packages('MASS') 后重渲染。\n")
}
Poisson、quasi-Poisson 与负二项的比较（quasi 无 AIC）
模型 项目_IRR 项目标准误 AIC
Poisson 0.776 0.049 5400
quasi-Poisson 0.776 0.057 NA
负二项 0.779 0.061 5234

聚类不是“多加一个分布”就解决 若同一人有多段观察、同一诊所服务多个参与者，独立个体的标准误会偏小。考虑以聚类为单位的稳健方差、GEE、随机效应/混合模型或按设计进行重抽样。负二项的个体层面异质性参数不自动等同于正确处理诊所相关性。

4 零很多，究竟意味着什么？

4.1 先描述零，再猜机制

计数数据出现许多零很常见：低基线率与短随访本身就能产生大量 Poisson 零。对指定协变量，Poisson 的零概率为 P(Y=0)=exp⁡(−μ)P(Y=0)=\exp(-\mu)。因此正确比较是“模型预测了多少零”，而不是仅报告零比例高。

observed_zero <- mean(ed_data$ed_visits == 0)
poisson_zero_pred <- mean(dpois(0, lambda = fitted(poisson_fit)))
zero_table <- data.frame(
  项目 = c("观测零比例", "Poisson 平均预测零比例", "模拟中处于捕获范围外的比例（教学真值）"),
  比例 = c(observed_zero, poisson_zero_pred, mean(ed_data$outside_capture == 1))
)
knitr::kable(zero_table, digits = 3, caption = "零比例的描述：真实研究通常没有最后一行的机制标签")
零比例的描述：真实研究通常没有最后一行的机制标签
项目 比例
观测零比例 0.605
Poisson 平均预测零比例 0.532
模拟中处于捕获范围外的比例（教学真值） 0.318

结构零是“在本研究的事件生成/捕获过程下不可能发生或不可能被记录”的子群；抽样零是“仍在风险或可记录，但这个观察期恰好没发生”。在真实项目中，结构零的定义必须有领域依据，例如确实未覆盖的医疗网络、明确不符合观察资格的人群；不能通过模型事后给每个零贴标签。

4.2 ZIP 模型的组成与似然

ZIP 假定每人先以概率 πi\pi_i 落在结构零组，否则以概率 1−πi1-\pi_i 进入 Poisson 风险组：

P(Yi=0)=πi+(1−πi)e−μi,P(Yi=y>0)=(1−πi)e−μiμiyy!. P(Y_i=0)=\pi_i+(1-\pi_i)e^{-\mu_i},\quad P(Y_i=y>0)=(1-\pi_i)\frac{e^{-\mu_i}\mu_i^y}{y!}.

通常以两个回归式连接协变量：

log⁡(μi)=log⁡(ti)+XiTβ,logit⁡(πi)=ZiTγ. \log(\mu_i)=\log(t_i)+X_i^T\beta,\qquad \operatorname{logit}(\pi_i)=Z_i^T\gamma.

这不是“零的逻辑回归 + 正数的 Poisson”两次独立拟合；零的观测同时可能来自两个部分，必须用联合似然。正数必来自 Poisson 风险组，但一个零的组别不可观测。

可识别性与命名的克制 ZIP 能把数据拟合成两个潜在组，不代表两个组在科学上真实、可稳定区分。若 XX 与 ZZ 完全相同、样本中正数很少或协变量强烈分离，两个部分可能高度不稳定。最好让零部分由明确的捕获/资格机制驱动，并做替代模型与敏感性分析。

5 用 base R 从头估计 ZIP

5.1 设计两个部分和稳定对数似然

为便于教学，计数部分使用项目、吸烟、严重度、农村、年龄和人年；零部分使用就医障碍与农村居住，代表与捕获范围相关的变量。现实研究应在分析计划中说明每个部分的变量、时间顺序和科学依据，而不是根据 AIC 反复筛选。

zip_count_formula <- ~ program + smoking + severity + rural + age10
zip_zero_formula <- ~ access_barrier + rural
X_zip <- model.matrix(zip_count_formula, data = ed_data)
Z_zip <- model.matrix(zip_zero_formula, data = ed_data)
y_zip <- ed_data$ed_visits
offset_zip <- log(ed_data$person_years)

# theta = (beta, gamma)。零值以 log-sum-exp 计算，避免 exp(-mu) 下溢或直接求和失稳。
zip_negloglik <- function(theta, y, X, Z, offset) {
  p_beta <- ncol(X)
  beta <- theta[seq_len(p_beta)]
  gamma <- theta[p_beta + seq_len(ncol(Z))]
  eta_count <- drop(offset + X %*% beta)
  # 只排除超出双精度指数函数数值域的参数，不改变可计算范围内的模型。
  if (any(!is.finite(eta_count)) || any(abs(eta_count) > 700)) return(1e100)
  mu <- exp(eta_count)
  eta_zero <- drop(Z %*% gamma)
  log_pi <- plogis(eta_zero, log.p = TRUE)
  log_one_minus_pi <- plogis(-eta_zero, log.p = TRUE)
  is_zero <- y == 0
  loglik <- numeric(length(y))
  loglik[is_zero] <- log_sum_exp2(
    log_pi[is_zero], log_one_minus_pi[is_zero] - mu[is_zero]
  )
  loglik[!is_zero] <- log_one_minus_pi[!is_zero] +
    dpois(y[!is_zero], lambda = mu[!is_zero], log = TRUE)
  if (any(!is.finite(loglik))) return(1e100)
  -sum(loglik)
}

# Poisson 系数提供计数部分的稳定起点；零部分使用多个起点降低局部最优风险。
poisson_start <- coef(glm(
  ed_visits ~ program + smoking + severity + rural + age10 + offset(log(person_years)),
  family = poisson, data = ed_data
))
poisson_zero_start <- mean(exp(-fitted(poisson_fit)))
pi_start <- (mean(y_zip == 0) - poisson_zero_start) / (1 - poisson_zero_start)
pi_start <- pmin(pmax(pi_start, 0.02), 0.80)
zero_intercept_starts <- unique(c(qlogis(pi_start), qlogis(c(0.05, 0.25, 0.50))))
zip_starts <- lapply(
  zero_intercept_starts,
  function(intercept) c(poisson_start, intercept, rep(0, ncol(Z_zip) - 1L))
)
stopifnot(all(vapply(
  zip_starts,
  length,
  integer(1)
) == ncol(X_zip) + ncol(Z_zip)))

5.2 执行最大似然、Hessian 与标准误

zip_candidates <- lapply(zip_starts, function(start) {
  optim(
    par = start, fn = zip_negloglik, method = "BFGS",
    y = y_zip, X = X_zip, Z = Z_zip, offset = offset_zip,
    hessian = TRUE, control = list(maxit = 2000, reltol = 1e-10)
  )
})
valid_candidate <- vapply(
  zip_candidates,
  function(candidate) candidate$convergence == 0 && is.finite(candidate$value),
  logical(1)
)
if (!any(valid_candidate)) {
  stop("所有 ZIP 优化起点都未收敛；应检查尺度、模型复杂度与数据机制。", call. = FALSE)
}
valid_fits <- zip_candidates[valid_candidate]
zip_optim <- valid_fits[[which.min(vapply(valid_fits, `[[`, numeric(1), "value"))]]

hessian_symmetric <- (zip_optim$hessian + t(zip_optim$hessian)) / 2
hessian_eigen <- eigen(hessian_symmetric, symmetric = TRUE, only.values = TRUE)$values
hessian_condition <- max(hessian_eigen) / min(hessian_eigen)
if (any(!is.finite(hessian_eigen)) || min(hessian_eigen) <= 0 ||
    !is.finite(hessian_condition) || hessian_condition > 1e10) {
  stop("ZIP Hessian 不够正定或条件数过大；当前 Wald 标准误不可靠。", call. = FALSE)
}
zip_vcov <- solve(hessian_symmetric)
zip_se <- sqrt(diag(zip_vcov))

zip_gradient <- vapply(seq_along(zip_optim$par), function(j) {
  step <- 1e-5 * max(1, abs(zip_optim$par[j]))
  upper <- lower <- zip_optim$par
  upper[j] <- upper[j] + step
  lower[j] <- lower[j] - step
  (zip_negloglik(upper, y_zip, X_zip, Z_zip, offset_zip) -
     zip_negloglik(lower, y_zip, X_zip, Z_zip, offset_zip)) / (2 * step)
}, numeric(1))
zip_gradient_max <- max(abs(zip_gradient))
if (!is.finite(zip_gradient_max) || zip_gradient_max > 0.01) {
  stop("ZIP 解的数值梯度仍过大；当前结果不可报告。", call. = FALSE)
}
p_beta <- ncol(X_zip)
zip_beta <- zip_optim$par[seq_len(p_beta)]
zip_gamma <- zip_optim$par[p_beta + seq_len(ncol(Z_zip))]
zip_beta_se <- zip_se[seq_len(p_beta)]
zip_gamma_se <- zip_se[p_beta + seq_len(ncol(Z_zip))]

zip_status <- data.frame(
  指标 = c("收敛代码（0 为成功）", "尝试的起点数", "负对数似然",
           "最大绝对数值梯度", "最小 Hessian 特征值", "Hessian 条件数"),
  数值 = c(zip_optim$convergence, length(zip_starts), zip_optim$value,
           zip_gradient_max, min(hessian_eigen), hessian_condition)
)
knitr::kable(zip_status, digits = 4, caption = "ZIP 优化与 Hessian 检查")
ZIP 优化与 Hessian 检查
指标 数值
收敛代码（0 为成功） 0.00e+00
尝试的起点数 4.00e+00
负对数似然 2.53e+03
最大绝对数值梯度 4.00e-04
最小 Hessian 特征值 1.78e+01
Hessian 条件数 1.71e+02

这里 Hessian 的逆是观测信息矩阵近似，Wald 标准误依赖局部二次近似。对边界参数、很小样本或不稳定混合模型，轮廓似然、bootstrap 或重新设计模型常比盲目报告 Wald 区间更诚实。

zip_count_table <- data.frame(
  计数部分变量 = colnames(X_zip),
  在风险组率比 = exp(zip_beta),
  下限 = exp(zip_beta - 1.96 * zip_beta_se),
  上限 = exp(zip_beta + 1.96 * zip_beta_se),
  P值 = 2 * pnorm(abs(zip_beta / zip_beta_se), lower.tail = FALSE)
)
zip_zero_table <- data.frame(
  零部分变量 = colnames(Z_zip),
  结构零优势比 = exp(zip_gamma),
  下限 = exp(zip_gamma - 1.96 * zip_gamma_se),
  上限 = exp(zip_gamma + 1.96 * zip_gamma_se),
  P值 = 2 * pnorm(abs(zip_gamma / zip_gamma_se), lower.tail = FALSE)
)
zip_count_effects <- zip_count_table[zip_count_table$计数部分变量 != "(Intercept)", ]
zip_zero_effects <- zip_zero_table[zip_zero_table$零部分变量 != "(Intercept)", ]
knitr::kable(zip_count_effects, digits = 3,
             caption = "ZIP 计数部分：条件于非结构零组的率比")
ZIP 计数部分：条件于非结构零组的率比
计数部分变量 在风险组率比 下限 上限 P值
program预防项目 program预防项目 0.77 0.693 0.854 0.000
smoking吸烟 smoking吸烟 1.43 1.273 1.596 0.000
severity severity 1.46 1.380 1.534 0.000
rural农村 rural农村 1.16 1.017 1.322 0.027
age10 age10 1.10 1.053 1.152 0.000 
knitr::kable(zip_zero_effects, digits = 3,
             caption = "ZIP 零部分：属于结构零组的优势比（需有机制依据）")
ZIP 零部分：属于结构零组的优势比（需有机制依据）
零部分变量 结构零优势比 下限 上限 P值
2 access_barrier有障碍 3.78 2.77 5.17 0
3 rural农村 2.12 1.49 3.02 0

表中不展示截距：计数截距指数化后是参考画像每人年的基线率，零部分截距指数化后是参考画像的结构零基线 odds；它们都不是协变量变化的比值。

5.3 两个部分的解释不能混为一谈

  • 计数部分的 IRR：在模型定义的非结构零（风险/可捕获）人群中，协变量变化对应的条件发生率比。
  • 零部分的 OR：协变量变化对应“属于结构零潜在组”的优势比，不是“观测到零”的优势比。
  • 边际均值：E(Yi∣Xi,Zi)=(1−πi)μiE(Y_i\mid X_i,Z_i)=(1-\pi_i)\mu_i，同时受两个部分影响；这通常更接近服务量预测。
  • 零概率：P(Yi=0)=πi+(1−πi)e−μiP(Y_i=0)=\pi_i+(1-\pi_i)e^{-\mu_i}，不是单独的 πi\pi_i。

零部分 OR 常不直观，尤其事件不罕见时；对决策者，优先给出明确人群和暴露时间下的 μ\mu、π\pi、边际均值以及零概率。

6 ZIP 的预测、频数与校准

6.1 为每个人计算四个不同的量

zip_eta_count <- drop(offset_zip + X_zip %*% zip_beta)
zip_mu <- exp(zip_eta_count)
zip_pi <- inv_logit(drop(Z_zip %*% zip_gamma))
zip_mean <- (1 - zip_pi) * zip_mu
zip_p0 <- zip_pi + (1 - zip_pi) * exp(-zip_mu)

pred_summary <- data.frame(
  量 = c("在风险组条件均值 mu", "结构零概率 pi", "边际期望次数", "观测到零的概率"),
  平均值 = c(mean(zip_mu), mean(zip_pi), mean(zip_mean), mean(zip_p0)),
  最小值 = c(min(zip_mu), min(zip_pi), min(zip_mean), min(zip_p0)),
  最大值 = c(max(zip_mu), max(zip_pi), max(zip_mean), max(zip_p0))
)
knitr::kable(pred_summary, digits = 3, caption = "ZIP 的四种预测量不可互换")
ZIP 的四种预测量不可互换
量 平均值 最小值 最大值
在风险组条件均值 mu 1.074 0.103 5.670
结构零概率 pi 0.326 0.187 0.648
边际期望次数 0.716 0.060 3.814
观测到零的概率 0.602 0.203 0.944 
zip_profile_data <- profile[rep(seq_len(nrow(profile)), each = 2), ]
zip_profile_data$access_barrier <- factor(
  rep(c("无明显障碍", "有障碍"), times = 2),
  levels = levels(ed_data$access_barrier)
)
X_profile <- model.matrix(zip_count_formula, data = zip_profile_data)
Z_profile <- model.matrix(zip_zero_formula, data = zip_profile_data)
profile_mu <- exp(log(zip_profile_data$person_years) + drop(X_profile %*% zip_beta))
profile_pi <- inv_logit(drop(Z_profile %*% zip_gamma))
profile_zip <- data.frame(
  项目 = zip_profile_data$项目,
  就医障碍 = zip_profile_data$access_barrier,
  风险组条件均值_mu = profile_mu,
  结构零概率_pi = profile_pi,
  边际期望次数 = (1 - profile_pi) * profile_mu,
  边际期望次数_每100人年 = 100 * (1 - profile_pi) * profile_mu,
  观测零概率 = profile_pi + (1 - profile_pi) * exp(-profile_mu)
)
knitr::kable(profile_zip, digits = 3,
             caption = "指定个体画像在 1 人年内的 ZIP 分量与边际预测；率另换算为每 100 人年")
指定个体画像在 1 人年内的 ZIP 分量与边际预测；率另换算为每 100 人年
项目 就医障碍 风险组条件均值_mu 结构零概率_pi 边际期望次数 边际期望次数_每100人年 观测零概率
1 常规服务 无明显障碍 0.750 0.187 0.610 61.0 0.571
1.1 常规服务 有障碍 0.750 0.465 0.402 40.2 0.717
2 预防项目 无明显障碍 0.577 0.187 0.470 47.0 0.643
2.1 预防项目 有障碍 0.577 0.465 0.309 30.9 0.765 
# 参数模拟保留 beta、gamma 及两个部分之间的完整交叉协方差。
set.seed(20260821)
n_draws <- 1200
coefficient_draws <- matrix(
  rnorm(n_draws * length(zip_optim$par)), nrow = n_draws
) %*% chol(zip_vcov)
coefficient_draws <- sweep(coefficient_draws, 2, zip_optim$par, "+")

beta_draws <- coefficient_draws[, seq_len(p_beta), drop = FALSE]
gamma_draws <- coefficient_draws[, p_beta + seq_len(ncol(Z_zip)), drop = FALSE]
mu_draws <- exp(sweep(
  beta_draws %*% t(X_profile),
  2,
  log(zip_profile_data$person_years),
  "+"
))
pi_draws <- inv_logit(gamma_draws %*% t(Z_profile))
marginal_draws <- (1 - pi_draws) * mu_draws
zero_probability_draws <- pi_draws + (1 - pi_draws) * exp(-mu_draws)

profile_zip$边际均值下限 <- apply(marginal_draws, 2, quantile, 0.025)
profile_zip$边际均值上限 <- apply(marginal_draws, 2, quantile, 0.975)
profile_zip$零概率下限 <- apply(zero_probability_draws, 2, quantile, 0.025)
profile_zip$零概率上限 <- apply(zero_probability_draws, 2, quantile, 0.975)
knitr::kable(profile_zip, digits = 3,
             caption = "ZIP 画像预测及参数不确定性的近似 95% 区间")
ZIP 画像预测及参数不确定性的近似 95% 区间
项目 就医障碍 风险组条件均值_mu 结构零概率_pi 边际期望次数 边际期望次数_每100人年 观测零概率 边际均值下限 边际均值上限 零概率下限 零概率上限
1 常规服务 无明显障碍 0.750 0.187 0.610 61.0 0.571 0.555 0.668 0.543 0.600
1.1 常规服务 有障碍 0.750 0.465 0.402 40.2 0.717 0.347 0.462 0.680 0.755
2 预防项目 无明显障碍 0.577 0.187 0.470 47.0 0.643 0.427 0.515 0.616 0.668
2.1 预防项目 有障碍 0.577 0.465 0.309 30.9 0.765 0.267 0.357 0.733 0.796

这些区间传播了回归系数的不确定性，并保留两个 ZIP 部分之间的交叉协方差；它们不是未来个体真实次数的预测区间，也不包含模型选择、结构假设或外部推广的不确定性。

6.2 观测频数与预测频数

只比较平均数可能掩盖尾部和零的错误。下面把每人的预测概率加总，得到各计数的期望人数；尾部合并为“6 次及以上”以避免稀疏单元。ZIP 的频数概率由混合分布给出。

count_bins <- 0:5
obs_frequency <- c(tabulate(pmin(y_zip, 6) + 1, nbins = 7))
poisson_expected <- sapply(count_bins, function(k) sum(dpois(k, fitted(poisson_fit))))
zip_expected <- sapply(count_bins, function(k) {
  if (k == 0) sum(zip_p0) else sum((1 - zip_pi) * dpois(k, zip_mu))
})
poisson_tail <- nrow(ed_data) - sum(poisson_expected)
zip_tail <- nrow(ed_data) - sum(zip_expected)
frequency_table <- data.frame(
  次数类别 = c(as.character(count_bins), "6 次及以上"),
  观测人数 = obs_frequency,
  Poisson_期望人数 = c(poisson_expected, poisson_tail),
  ZIP_期望人数 = c(zip_expected, zip_tail)
)
knitr::kable(frequency_table, digits = 1, caption = "观测与模型预测的计数频数：零和尾部都应检查")
观测与模型预测的计数频数：零和尾部都应检查
次数类别 观测人数 Poisson_期望人数 ZIP_期望人数
0 1451 1276.1 1445.9
1 498 718.8 506.2
2 257 274.0 257.8
3 120 90.9 113.2
4 47 28.3 46.6
5 15 8.5 18.6
6 次及以上 12 3.4 11.7 
set_count_plot_font()
barplot(
  t(as.matrix(frequency_table[, c("观测人数", "Poisson_期望人数", "ZIP_期望人数")])),
  beside = TRUE, names.arg = frequency_table$次数类别,
  col = c(palette_count["navy"], palette_count["orange"], palette_count["teal"]),
  ylab = "人数", xlab = "ED 就诊次数类别",
  main = "分布校准：不要只比较 AIC 或回归系数"
)
legend("topright", legend = c("观测", "Poisson", "ZIP"),
       fill = c(palette_count["navy"], palette_count["orange"], palette_count["teal"]), bty = "n")
分组柱状图显示零至五次及六次以上 ED 就诊类别的观测人数、Poisson 期望人数和 ZIP 期望人数。

观测频数与 Poisson、ZIP 预测频数的比较

6.3 按预测风险分组校准

校准回答“同样预测为 0.4 次的人群，实际平均是否约为 0.4 次？”它不是证明模型机制正确的检验。这里按 ZIP 边际均值十分位分组，比较观测平均次数、Poisson 平均预测和 ZIP 平均预测；在真实分析中还应预先保留验证样本或做交叉验证。

calibration_group <- cut(
  zip_mean,
  breaks = unique(quantile(zip_mean, probs = seq(0, 1, 0.1))),
  include.lowest = TRUE
)
calibration <- aggregate(
  cbind(观测平均次数 = y_zip, Poisson预测 = fitted(poisson_fit), ZIP预测 = zip_mean) ~ calibration_group,
  FUN = mean
)
calibration$n <- as.integer(table(calibration_group))
knitr::kable(calibration, digits = 3, caption = "按 ZIP 边际预测十分位的均值校准")
按 ZIP 边际预测十分位的均值校准
calibration_group 观测平均次数 Poisson预测 ZIP预测 n
[0.0604,0.243] 0.179 0.224 0.179 240
(0.243,0.319] 0.300 0.333 0.281 240
(0.319,0.407] 0.392 0.405 0.361 240
(0.407,0.49] 0.396 0.495 0.450 240
(0.49,0.584] 0.546 0.580 0.536 240
(0.584,0.696] 0.600 0.645 0.639 240
(0.696,0.845] 0.729 0.768 0.768 240
(0.845,1.04] 0.987 0.905 0.943 240
(1.04,1.37] 1.167 1.120 1.189 240
(1.37,3.81] 1.858 1.680 1.811 240 
set_count_plot_font()
plot(
  seq_len(nrow(calibration)), calibration$观测平均次数, type = "b", pch = 16,
  col = palette_count["navy"], ylim = range(calibration[, c("观测平均次数", "Poisson预测", "ZIP预测")]),
  xlab = "ZIP 边际预测十分位组", ylab = "平均 ED 次数",
  main = "均值校准应结合外部或重抽样验证"
)
lines(seq_len(nrow(calibration)), calibration$Poisson预测, type = "b", pch = 17, col = palette_count["orange"])
lines(seq_len(nrow(calibration)), calibration$ZIP预测, type = "b", pch = 15, col = palette_count["teal"])
legend("topleft", legend = c("观测", "Poisson", "ZIP"), pch = c(16, 17, 15),
       col = c(palette_count["navy"], palette_count["orange"], palette_count["teal"]), bty = "n")
折线图横轴为十个 ZIP 预测分组，纵轴为平均 ED 次数，比较观测值、Poisson 预测和 ZIP 预测。

按预测十分位的观测平均次数与两种模型平均预测

6.4 AIC 有用，但不是裁判

对同一数据、同一响应且有定义良好完整似然的模型，较小 AIC 表示在拟合优度与参数数目之间较好的近似预测折衷。它不检验结构零的科学真实性，不惩罚不良的队列设计，也不能比较 quasi-Poisson。混合模型的非嵌套、边界参数和不同支持集会让常规似然比检验更加棘手。

zip_aic <- 2 * length(zip_optim$par) + 2 * zip_optim$value
aic_table <- data.frame(
  模型 = c("Poisson", if (nb_available) "负二项" else NULL, "ZIP（base R 实现）"),
  AIC = c(AIC(poisson_fit), if (nb_available) AIC(nb_fit) else NULL, zip_aic)
)
knitr::kable(aic_table, digits = 1,
             caption = "AIC 仅作候选模型的一个比较维度；quasi-Poisson 不在表中")
AIC 仅作候选模型的一个比较维度；quasi-Poisson 不在表中
模型 AIC
Poisson 5400
负二项 5234
ZIP（base R 实现） 5086

7 ZIP、hurdle 和其他分析策略

7.1 Hurdle 模型回答的是另一种过程假设

特征 ZIP Hurdle（两部分）
零的来源 结构零或 Poisson 抽样零 所有零都由“是否跨过门槛”过程产生
正数部分 普通 Poisson（可在理论上产生零） 零截断 Poisson/负二项
零与正数的关系 零的潜在来源不可观测 先发生任意一次，再决定正数大小
适合的科学故事 未覆盖系统者 vs 可捕获风险者 接触服务的障碍 vs 接触后使用强度

例如，“是否进入任何参与医院”与“进入后急诊次数”可能支持 hurdle；“部分人始终不在本系统的捕获框内”可能支持 ZIP。两者都应从研究过程出发，而不是从零柱最高这一事实反推。

7.2 设计、聚类与因果推断的边界

若项目按社区实施，个人结果可能在社区内相关；若每人多次年度随访，数据也是纵向的。可考虑 GEE、混合效应负二项/ZIP、带随机效应的两部分模型或以社区为重抽样单位。模型选择须与抽样单位、暴露分配单位和目标推断总体一致。

对于非随机项目评估，还要在建模前明确目标试验：资格、时间零点、治疗策略、随访、结局、混杂控制和估计量。倾向评分、标准化、加权或 g 方法可用于处理已测量混杂；零膨胀模型仍只是结局模型的一部分。参见 因果推断 与 倾向评分匹配 模块。

8 从问题到报告：可审核工作流

8.1 推荐分析路径

  1. 定义事件与人年：明确事件窗口、重复计数规则、进入/退出规则和数据捕获边界。
  2. 描述分布：按暴露组报告人数、人年、总次数、粗率、零比例和高计数尾部。
  3. 预先定义均值结构：基于时间顺序和领域知识选择混杂因素、非线性和交互；不要用结局驱动变量筛选。
  4. 拟合基准 Poisson：使用 offset，并报告 IRR、绝对预测、Pearson 离散度和残差图。
  5. 处理额外变异/相关性：依据目标分别考虑稳健标准误、quasi、负二项、GEE 或混合模型。
  6. 提出并审计零机制：说明什么人为什么可能结构性为零；比较 Poisson 预测零与观测零。
  7. 如使用 ZIP/hurdle：报告两个公式、优化收敛、识别风险、分量预测、频数与校准。
  8. 验证与敏感性分析：比较候选分布、改变均值结构、检查高影响观察、做重抽样或外部验证。
  9. 透明报告边界：说明数据来源、缺失、聚类、未测量混杂、模型依赖性与结果可推广范围。

8.2 报告结果的句式模板

在 N=...N=... 名符合资格的参与者中，累计随访 … 人年，观察到 … 次呼吸系统 ED 就诊（粗率 …/人年；零比例 …）。使用以对数人年为 offset 的 [Poisson/负二项/ZIP] 回归，调整 … 后，预防项目相对常规服务的 [条件 IRR/边际预测差异] 为 …（95% CI …）。[若 ZIP：在风险组计数部分的 IRR 为 …；结构零部分建模的是 …，不应解释为观测零的 OR。] Pearson 离散度为 …；观测—预测频数与分组校准见图 …。结果描述条件关联；由于 …，不能作出/需谨慎作出因果解释。

8.3 常见错误及修正

常见错误 为什么不对 修正
不同随访时长却不设 offset 把暴露机会混入组别差异 核实人年，使用 offset(log(person_years))
把 IRR 说成风险比或绝对减少 估计量与尺度被改变 明确“率比”，另报标准化绝对预测
看到离散度大就自动 ZIP 过度离散有许多来源 先检查均值结构、聚类、NB/quasi 和零机制
把 ZIP 零部分 OR 称作“零的 OR” 观测零来自两种成分 报告 π\pi、P(Y=0)P(Y=0) 与边际均值
用 AIC 宣称发现结构零群体 AIC 不提供机制证据 用领域定义、数据流程和敏感性分析支持
忽略 Hessian/收敛警告 混合模型可处在边界或不可识别 检查参数、起点、模型复杂度和不确定性
用个体独立模型分析诊所聚类数据 标准误和模型结构可能错误 按设计使用 GEE、混合模型或聚类重抽样

9 知识检查与练习

快速自测

  1. 某人被观察 0.5 人年、另一个人被观察 2 人年；为什么简单比较其原始次数可能不公平？
  2. Poisson 的 Pearson 离散度为 2.4，能否仅据此判定数据存在结构零？
  3. ZIP 中 πi=0.30\pi_i=0.30、μi=1.50\mu_i=1.50，该人的边际均值与观测零概率分别是什么？
  4. 若研究问题是“是否首次进入任何急诊”，为什么本页的重复事件率模型未必最合适？
查看答案
  1. 后者有四倍事件观察机会；应以人年为分母，并在率模型中使用 offset。
  2. 不能。它说明残差变异高于 Poisson 工作假设，可能来自遗漏变量、异质性、聚类、非线性或零机制。
  3. 边际均值为 (1−0.30)×1.50=1.05(1-0.30)\times1.50=1.05；零概率为 0.30+0.70e−1.50≈0.4560.30+0.70e^{-1.50}\approx0.456。
  4. 该问题是首次事件/风险或时间到事件问题；应明确删失与时间零点，考虑二项或生存分析，而不是把重复次数当作唯一结局。

练习：修改并审计模型

  1. 将计数部分加入 program:severity 交互，说明交互项的 IRR 在什么条件下解释；比较指定画像的边际预测，而不是仅看 P 值。
  2. 将零部分的 access_barrier 去掉，比较观测—预测零人数、AIC 与实际数据流程的合理性。哪一个证据最重要？
  3. 令所有人年固定为 1，重新拟合。哪些估计量改变？这说明为何应保留原始暴露时间？
  4. 设想数据来自 20 个诊所。写出至少两种会改变推断的相关性处理方案，并说明各自的边际/条件解释。
练习讨论要点
  1. 交互 IRR 修饰的是项目的条件率比；主效应只适用于参考严重度。应画出每个严重度、同一人年和协变量画像的边际均值与区间。
  2. AIC 和零频数都是诊断；最优先的是是否真的存在由就医障碍刻画的捕获范围机制。没有机制依据时，不应把更低 AIC 写成“发现结构零”。
  3. 率模型中的协变量 IRR 在正确模型下可能相近，但均值、截距、拟合、标准误和预测会改变；真实人年承载了可观察事件机会。
  4. 可用人口平均的 GEE（聚类稳健方差）或诊所随机截距的混合模型；前者常解释为边际关联，后者为给定随机效应的条件关联。选择由 estimand 与设计决定。

10 速查表与提交前清单

速查表

目标 关键量/代码 解释提醒
率模型 glm(y ~ x + offset(log(py)), family=poisson) offset 系数固定为 1
IRR exp(coef(fit)) 条件率比，非风险比
Poisson 离散度 sum(residuals(fit,"pearson")^2)/df.residual(fit) 大于 1 需查原因，不是 ZIP 诊断
quasi-Poisson family = quasipoisson 调整方差；无 AIC
负二项 MASS::glm.nb(...) 允许二次均值—方差关系
ZIP 边际均值 (1 - pi) * mu 同时受两个部分影响
ZIP 观测零概率 pi + (1 - pi) * exp(-mu) 不等于 pi
频数检查 逐人预测概率后求和 同时看零、中心与尾部

提交前清单

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60      
##  [5] lattice_0.22-9  cachem_1.1.0    knitr_1.51      htmltools_0.5.9
##  [9] rmarkdown_2.31  stats4_4.6.1    lifecycle_1.0.5 cli_3.6.6      
## [13] grid_4.6.1      sass_0.4.10     jquerylib_0.1.4 compiler_4.6.1 
## [17] tools_4.6.1     nlme_3.1-169    evaluate_1.0.5  bslib_0.12.0   
## [21] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0  MASS_7.3-65