跳转至

Kalman Filtering and Smoothing for Improving Precision in Horvitz--Thompson Estimation of Infectious Disease Prevalence

作者: Jeongjin Lee, Grzegorz A. Rempala, Patrick M. Schnell
主题: 流行病学
相关性: 6/10
链接: https://arxiv.org/abs/2609.09325

一、领域脉络与小综述

这个方向是什么

本文所处的子方向是传染病监测中的患病率估计，其根本科学问题是：在纵向重复检测制度下（如大学校园每周例行检测），如何从非随机、非代表性的检测数据中无偏且高效地估计人群患病率。该方向的成熟度处于"方法已建立、精度待提升"的阶段——Horvitz–Thompson（HT）估计量已能纠正检测机制带来的选择偏差，但其逐日估计噪声大、方差随时间波动、且在无检测日完全缺失。本文的贡献在于将状态空间模型与 Kalman 滤波/平滑引入这一估计框架，通过借用时间邻域信息来提升精度。

发展脉络（history）

  • 奠基工作：Horvitz & Thompson (1952) 提出逆概率加权估计量，为处理非随机抽样提供了基础工具。这是所有后续工作的统计根基。
  • 主要进展：Schnell et al. (2024) 将 HT 估计量适配到纵向重复检测场景，解决了检测进度表偏差（schedule bias）问题；Lee et al. (2026) 进一步扩展至包含症状驱动检测和接触者追踪的复杂检测机制，通过反事实框架形式化了检测过程的因果结构。这两篇工作确立了"用 HT 估计量纠正选择偏差"这一主流路线。
  • 当前 frontier：HT 估计量虽然无偏，但逐日方差大、精度随时间变化、无检测日无估计。如何在不引入偏差的前提下提升精度、填补缺失日，是当前的前沿问题。
  • 本文的位置：本文不改变 HT 估计量的构造，而是将其视为带噪观测，用状态空间模型（局部线性趋势）描述潜在患病率过程，用 Kalman 滤波/平滑输出更稳定的估计序列。这是"后处理"路线——在保留 HT 无偏性的同时，通过时间结构换取精度。

子线索聚类

  1. 检测机制建模与偏差校正（Schnell et al. 2024; Lee et al. 2026）：核心是理解"谁被检测、为什么被检测"，通过逆概率加权消除选择偏差。这条线索的边界是：估计量无偏但方差大。
  2. 时间序列平滑与状态空间方法（本文）：将逐日估计视为潜在过程的带噪观测，用 Kalman 滤波/平滑提升精度。这条线索的边界是：需要指定状态演化模型（如局部线性趋势），且模型误设可能引入偏差。
  3. 方差估计与置信区间构造（Kaplan & Liu 2024）：针对有偏估计量的置信区间校准问题，本文直接借用了其 CI2/CI5/CI6 构造框架。

这个方向在追问的核心问题

  1. 如何在不引入偏差的前提下提升估计精度？ 主流方法是 HT 加权（无偏但噪声大），本文的答案是时间平滑（有偏但方差小），两者的权衡是核心张力。
  2. 如何处理检测中断（如周末、节假日）？ HT 估计量在无检测日直接缺失，Kalman 滤波通过预测步骤填补，平滑则用后续观测回溯修正。
  3. 如何构造覆盖概率可信的置信区间？ 模型驱动的区间（基于状态协方差）可能因偏差而欠覆盖，本文尝试了多种替代构造。

⚠️ 作者的 framing

作者将缺口 frame 成："HT 估计量虽然无偏，但逐日噪声大、精度随时间变化、无检测日无估计，因此需要一个系统性的后处理框架来提升精度。" 这使得本文成为"显然的下一步"——既然已有无偏估计量，下一步自然是让它更精确。被淡化的竞争路线包括：直接对检测过程建模（而非后处理 HT 估计）、贝叶斯分层模型（可自然处理缺失和先验信息）、以及非参数平滑方法（如局部多项式回归，作者仅在附录中作为 baseline 提及）。值得研究者去查的问题：为什么作者不直接对原始检测数据建模，而是选择对 HT 估计量做后处理？后处理是否丢失了信息？以及，局部线性趋势假设在传染病传播场景下是否合理——疫情爆发期患病率可能呈指数增长而非线性趋势。

张力

未见明显对立引用。但存在一个隐含张力：HT 估计量的设计目标是无偏，而 Kalman 滤波/平滑的目标是最小均方误差——后者必然引入偏差。作者在模拟中展示了 RMSE 的大幅下降，但代价是估计量不再无偏，这一权衡在文中未被系统讨论。

二、最核心、最简单的例子 / 数学问题

第一步：符号、模型、可观测数据

符号：

  • \(p_t\)：第 \(t\) 天的潜在患病率（latent prevalence），是要估计的 estimand，不可直接观测。
  • \(v_t\)：潜在斜率（latent slope），表示患病率从 \(t\) 到 \(t+1\) 天的预期变化，是辅助状态变量。
  • \(\alpha_t = (p_t, v_t)^\top \in \mathbb{R}^2\)：潜在状态向量。
  • \(\hat{p}_t\)：第 \(t\) 天的 HT 估计量，是可观测数据（由检测数据计算得出）。
  • \(R_t\)：\(\hat{p}_t\) 的日别方差，由 delete-a-group jackknife 估计，视为已知。
  • \(Q = \text{diag}(Q_{\text{level}}, Q_{\text{slope}})\)：状态演化噪声的协方差矩阵，是要估计的超参数。
  • \(F = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}\)：状态转移矩阵。
  • \(H = \begin{pmatrix} 1 & 0 \end{pmatrix}\)：观测矩阵。
  • \(w_t \sim N(0, Q)\)：状态演化噪声。
  • \(e_t \sim N(0, R_t)\)：观测噪声。
  • \(\tilde{\alpha}_{t|s}\)：基于截至 \(s\) 时刻观测对 \(\alpha_t\) 的条件期望。
  • \(P_{t|s}\)：对应的条件协方差矩阵。

模型（状态空间模型）：

  • 状态方程：\(\alpha_t = F \alpha_{t-1} + w_t\)，即 \(p_t = p_{t-1} + v_{t-1} + \eta_t\)，\(v_t = v_{t-1} + \zeta_t\)。这表示患病率的变化由斜率驱动，斜率本身随机游走。
  • 观测方程：\(\hat{p}_t = H \alpha_t + e_t = p_t + e_t\)，即 HT 估计量是潜在患病率的无偏但带噪观测。

可观测数据：研究者实际能观测到的是每日的 HT 估计量 \(\hat{p}_t\) 及其方差估计 \(R_t\)。潜在状态 \((p_t, v_t)\) 不可观测，只能通过 Kalman 滤波/平滑推断。

第二步：最小内核

最小内核：假设只有 \(T=3\) 天，每天都有 HT 估计 \(\hat{p}_1, \hat{p}_2, \hat{p}_3\) 及其方差 \(R_1, R_2, R_3\)。忽略斜率（即设 \(v_t \equiv 0\)），模型退化为：

\[p_t = p_{t-1} + \eta_t, \quad \eta_t \sim N(0, Q_{\text{level}})\]
\[\hat{p}_t = p_t + e_t, \quad e_t \sim N(0, R_t)\]

要解决的数学问题：给定 \(\hat{p}_1, \hat{p}_2, \hat{p}_3\)，如何估计 \(p_2\)（第 2 天的患病率）？

核心思想：HT 估计 \(\hat{p}_2\) 是无偏的，但方差 \(R_2\) 可能很大。Kalman 滤波的做法是：

  1. 预测：用第 1 天的滤波估计预测第 2 天：\(\tilde{p}_{2|1} = \tilde{p}_{1|1}\)（因为 \(v \equiv 0\)），预测方差 \(P_{2|1} = P_{1|1} + Q_{\text{level}}\)。
  2. 更新：将预测与观测 \(\hat{p}_2\) 加权平均：
    \[\tilde{p}_{2|2} = \frac{R_2}{P_{2|1} + R_2} \cdot \tilde{p}_{2|1} + \frac{P_{2|1}}{P_{2|1} + R_2} \cdot \hat{p}_2\]

为什么成立：这是高斯线性状态空间模型下的精确条件期望。权重由方差决定——若 \(R_2\) 大（HT 估计噪声大），则更多权重放在预测上；若 \(P_{2|1}\) 大（历史信息少），则更多权重放在观测上。关键点：滤波估计 \(\tilde{p}_{2|2}\) 的方差 \(P_{2|2} = \frac{P_{2|1} R_2}{P_{2|1} + R_2} \leq R_2\)，即严格小于 HT 估计的方差——这就是精度提升的来源。

平滑：滤波只用 \(t \leq 2\) 的信息，平滑则用全部 \(T=3\) 天的信息。第 2 天的平滑估计 \(\tilde{p}_{2|3}\) 会利用第 3 天的观测来回溯修正，方差进一步降低。

一般情形的推广：加入斜率 \(v_t\) 后，状态变为二维，但滤波/平滑的递归结构不变。缺失日（如周末）的处理方式是：不执行更新步骤，仅执行预测步骤，让状态继续演化。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：如何将每日 Horvitz–Thompson 患病率估计序列转化为更稳定、更精确的患病率轨迹，并处理无检测日的缺失问题。
  2. 核心工具/方法：联合局部线性趋势状态空间模型 + Kalman 滤波（实时估计）与 Kalman 平滑（回顾性估计），观测方差用 delete-a-group jackknife 估计，过程方差用创新似然（innovation likelihood）最大化估计。
  3. 主要结论：在模拟和俄亥俄州立大学 2020 年秋季 SARS-CoV-2 监测数据中，Kalman 滤波和平滑相对原始 HT 估计大幅降低 RMSE、保留主要时间趋势，并在无检测日提供 HT 无法给出的估计和区间。

关键设定与假设

  • 状态空间模型：联合局部线性趋势模型（式 1-5），状态 \(\alpha_t = (p_t, v_t)^\top\)，状态转移矩阵 \(F = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}\)。相比标准随机游走（仅 \(p_t\)），增加了斜率项 \(v_t\)，允许患病率趋势持续上升或下降，更灵活。
  • 观测方程：\(\hat{p}_t = p_t + e_t\)，\(e_t \sim N(0, R_t)\)。关键假设：HT 估计量是潜在患病率的无偏观测，且观测方差 \(R_t\) 已知（由 jackknife 估计）。这比标准 Kalman 滤波假设观测方差已知更强——实际中 \(R_t\) 是估计的，但文中未讨论这一不确定性对推断的影响。
  • 过程方差估计：\(Q_{\text{level}}, Q_{\text{slope}}\) 通过最大化 Gaussian innovation likelihood 估计（式 21-22）。关键细节：滤波使用截至 \(t\) 时刻的数据估计 \(Q\)，保证实时性；平滑使用全序列估计 \(Q\)，是回顾性的。
  • 缺失观测：当 \(\hat{p}_t\) 或 \(R_t\) 缺失时，滤波跳过更新步骤，仅做预测；平滑则通过后续观测回溯重构缺失时段。
  • 置信区间：默认区间基于状态协方差矩阵（式 24-25），但作者承认可能欠覆盖，因此考虑了 Kaplan & Liu (2024) 的 CI2/CI5/CI6 构造。CI2 用 HT 标准误作为尺度，CI5 用校准临界值，CI6 用 HT 与滤波/平滑估计的凸组合。

主要结果

  • 模拟研究：基于 Lee et al. (2026) 的检测机制模拟 21 天疫情，比较 HT、KF、KS 的 RMSE 和覆盖率。核心结果：KF 和 KS 相对 HT 大幅降低 RMSE（具体数值需查原文图 1-3），KS 通常优于 KF；在缺失日（第 10-11 天），KF 和 KS 仍能给出估计，而 HT 完全缺失。默认区间欠覆盖，CI2/CI5 改善覆盖率但区间更宽，CI6 在覆盖率和区间长度间取得较好平衡。
  • 真实数据：俄亥俄州立大学 2020 年秋季 11,335 名住校本科生监测数据。KF 和 KS 在无检测日（周末、劳动节、退伍军人节）给出 HT 无法提供的估计；在观测日，KF/KS 区间比 HT 更窄。KF 和 KS 轨迹在大多数观测日几乎重合，但在缺失时段后 KS 会利用后续观测修正，与 KF 出现差异。
  • 与简单平滑器的比较（附录 B）：在完整数据下，3 天和 5 天中心移动平均与 KS 的 RMSE 相近，但移动平均需要预先选择窗口宽度，且无法自然处理缺失或日别方差；KS 的优势在于数据自适应地确定平滑强度，并显式利用日别方差信息。

🔎 结论是否比证明窄

  • Lemma 1 和 Lemma 2 的证明是严格的：Lemma 1 证明滤波后验方差 \(\leq R_t\)（式 46），Lemma 2 证明平滑方差 \(\leq\) 滤波方差（式 53）。这两个引理在固定 \(Q\) 和 \(R_t\) 的条件下成立，证明是干净的。
  • 但"精度提升"的结论比证明宽：Lemma 1-2 证明的是条件方差的序关系，而非 RMSE 的序关系。RMSE 包含偏差项，而 KF/KS 是有偏的（因收缩和模型误设）。模拟中 RMSE 的降低是经验结果，没有理论保证。作者在讨论中承认了这一点，但未给出偏差的界。
  • 过程方差估计的不确定性未被纳入：\(Q\) 被估计后当作已知，区间未反映这一不确定性。这在短序列（如 21 天模拟）中可能造成欠覆盖，作者通过 CI2/CI5/CI6 部分缓解，但未给出理论保证。
  • "实时可用"的结论有条件：滤波的实时性依赖于 \(Q\) 的估计只用截至 \(t\) 的数据，但文中未讨论 \(Q\) 估计的收敛速度——在疫情早期（数据少）时，\(Q\) 估计可能不稳定，滤波性能可能下降。

四、开放问题

  1. 偏差-方差权衡的理论刻画：KF/KS 引入的偏差有多大？能否给出偏差的界，或在什么条件下 RMSE 严格优于 HT？这扎根于 Lemma 1-2 只证明了方差序关系、未涉及偏差的事实。要确认这是否为真 gap：去读近期关于 shrinkage estimator 的偏差-方差权衡文献（如 Donoho 1995 后的 minimax shrinkage 理论），看是否已有可迁移的框架。

  2. 过程方差估计的不确定性传播：\(Q\) 的估计误差如何影响滤波/平滑区间？能否构造考虑 \(Q\) 不确定性的区间（如通过 parametric bootstrap 或 Bayesian 方法）？这扎根于文中将 \(Q\) 视为已知的简化处理。要确认：去读状态空间模型的 Bootstrap 滤波文献（如 Stoffer & Wall 1991），看是否有现成方法。

  3. 模型误设的稳健性：局部线性趋势假设在疫情指数增长期是否合理？如果真实患病率呈指数增长，KF/KS 是否会因模型误设而产生系统性偏差？这扎根于作者仅在讨论中提及"模型选择"而未系统研究。要确认：去读关于状态空间模型误设稳健性的文献（如 Qu & Chen 2012 的 robust filtering），看是否有可借鉴的框架。

  4. HT 估计量的相关性结构：本文假设观测误差 \(e_t\) 独立，但同一批个体可能被重复检测，HT 估计量在不同日期间可能相关。忽略这一相关性是否会导致区间过窄？这扎根于文中未讨论 \(e_t\) 的时间相关性。要确认：去读纵向抽样下 HT 方差估计的文献（如 Berger 2004 的 variance estimation for longitudinal surveys），看相关性对推断的影响是否已被研究。

  5. 非线性/非高斯扩展：患病率有界（0 到 1），但 Kalman 滤波假设高斯状态和观测，可能产生越界估计。能否用粒子滤波或 logistic 状态空间模型改进？这扎根于文中未处理有界性问题。要确认：去读关于有界状态空间模型的文献（如 Smith 2015 的 non-Gaussian state space models），看是否有现成的扩展。

提醒：要确认上述问题是否是真 gap，建议去读同子领域（传染病监测中的状态空间方法）近期约 5 篇论文的引言——如果多篇都指向同一个问题，说明是共识性 gap；如果互相打架，则可能是机会所在。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论