连续时间随机流行病模型拟合部分观察发病数据的高效贝叶斯数据增强算法

《Biometrics》：An efficient Bayesian data augmentation algorithm for continuous-time stochastic epidemic models fit to partially observed incidence data

　　

传染病暴发的统计建模依赖发病数据（incidence data），对于理解流行病动力学和指导公共卫生干预至关重要。随机流行病模型（stochastic epidemic models，SEMs）能够考虑流行病演变过程中固有的随机性，在小型暴发或低患病率环境中尤为有用。在SEMs中，流行病被自然地表示为连续时间马尔可夫跳过程（continuous-time Markov jump process，MJP），而监测数据在离散时间点收集并存在检测不足（underdetection）。将SEMs拟合至部分观察到的流行病轨迹需要对所有可能的转移路径进行积分，而这通常是不可解的。个体水平数据增强（individual-level data augmentation，DA）方法已被发展以应对这些挑战。然而，现有方法要么依赖逐受试者更新，可扩展性有限，要么对状态空间SEMs施加通常不切实际的假设，以简化复杂的维度变化问题。研究人员提出了一种高效的DA-马尔可夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）框架，用于将连续时间状态空间SEMs拟合到发病数据。该方法使用块采样器（block sampler）联合更新个体轨迹，并使用潜在MJP的离散时间近似（discrete-time approximation）来提议转移时间。通过模拟研究，研究人员证明与先前的DA方法相比，该方法显著改善了MCMC混合（MCMC mixing）和可扩展性，并在粒子MCMC（particle MCMC）算法表现不佳的场景中仍保持计算稳健性。该方法被应用于2014年塞拉利昂一个农村村庄的埃博拉（Ebola）疫情和2009年加拿大一个偏远社区的H1N1疫情。

传染病暴发的统计建模对于理解流行病动力学、预测疫情演变和指导公共卫生干预至关重要。 compartmental模型（如易感-感染-康复（susceptible-infectious-recovered，SIR）模型）被广泛用于描述疫情演变。与确定性模型相比，随机流行病模型（stochastic epidemic models，SEMs）能够考虑流行病演变固有的随机变异性，这在疾病患病率低或人群规模小时尤为重要，因为个体病程的随机性会显著塑造流行病轨迹，忽略它可能导致错误结论。

对于许多疫情，统计推断主要依赖监测系统以离散时间（如每日或每周）收集的汇总水平发病数据（incidence data）。此类数据仅提供关于连续时间演变流行病过程的不完整信息，且由于无症状病例漏报、监测能力有限或随机噪声，数据质量常因不完美病例检测而受损。状态空间模型（state-space model）框架被用于联合建模潜在流行病过程与给定潜在流行病状态的数据观测过程。潜在流行病过程可由SEMs的连续时间马尔可夫跳过程（continuous-time Markov jump process，MJP）形式表示，其无穷小转移率依赖于状态。将此类连续时间状态空间SEMs拟合到发病数据需要对离散观测时间之间的所有可能潜在连续时间转移路径进行积分，这在即使小型人群中也不可行。

已有多种统计方法试图规避这一问题。基于模拟的方法（包括伪边缘方法、近似贝叶斯计算（Approximate Bayesian Computation）、序贯蒙特卡洛（sequential Monte Carlo）和粒子边缘Metropolis-Hastings（particle marginal Metropolis-Hastings，PMMH）方法）依赖对潜在MJP轨迹的前向模拟来近似不可解的似然，但其缺陷在于需要仔细调参，并可能遭受粒子退化（particle degeneracy）和收敛问题，尤其在观测噪声低或数据高度 informative 时。离散时间近似和线性噪声近似（linear noise approximation）虽可使离散观测时间之间的转移密度可解并支持高效推断，但其假设通常依赖渐近性质，在低患病率或小人群场景中往往不足以推断流行病动力学。

作为互补方法，个体水平数据增强（individual-level data augmentation，DA）通过个体水平流行病轨迹增广潜变量空间，使完全数据似然可解，可在无额外假设下进行精确推断，能灵活整合个体水平协变量与相关结构，且在基于模拟的方法失效时仍可能有效。然而其主要挑战是高效更新高维增广潜状态：早期reversible jump MCMC算法需要确切的个体水平事件时间且计算密集；后来发展的non-centered提议机制虽将全局提议适配于当前参数值，但不受数据指导，在观测信息丰富时效果较差；Fintzi等设计的针对患病率数据的逐受试者数据引导算法，计算成本高且人群规模增大时混合缓慢；近期出现的块采样器假设感染或康复被精确观测，未考虑层次状态空间结构，简化了潜状态空间给定观测数据的复杂维度变化问题，且不清楚如何扩展到简单SIR结构之外。因此，现有DA方法或可扩展性有限，或施加不切实际的假设，这促使研究人员开发新方法。

为解决上述挑战，研究人员提出了一种个体水平DA-MCMC方案，对状态空间SEMs拟合到部分观察发病数据进行精确推断。为高效探索潜变量空间，研究人员使用潜在MJP的离散时间近似在每个M-H（Metropolis-Hastings）步骤形成潜在事件时间的有效块提议。具体而言，在每个观测区间，先由数据引导提议一个聚合水平辅助变量以控制离开易感状态的个体数量，从而确定块大小，随后联合提议这些个体的完整流行病史。该提议机制可紧密近似真实潜在MJP，实现潜在流行病事件时间的有效联合模拟，从而缓解传统逐事件或逐受试者DA-MCMC采样器混合缓慢的问题，并改善全局提议低M-H接受率的不足。模拟研究表明，与先前DA方法相比该方法显著改善了MCMC混合和可扩展性，并在粒子MCMC算法困难的场景中保持计算稳健，为小型至中型疫情的连续时间随机流行病动力学推断提供了实用工具。

关键技术方法方面，方法以连续时间随机SIR模型为核心，以MJP描述易感者S(t)、感染者I(t)、康复者R(t)的转移，个体轨迹Z_j(t)由感染时间τ_j^I与康复时间τ_j^R刻画；观测模型为二项分布Y_k|X_k,φ～Binomial(X_k,φ)。块采样器基于离散时间SIR马尔可夫链的Poisson转移概率在每个观测区间提议真实感染增量X_k，再联合提议感染时间（截断指数分布）与康复时间。参数更新采用部分折叠Gibbs采样器（partially collapsed Gibbs sampler），对φ积分，并以r₀=βP/γ的对数变换重参数化。算法用C++与R包Rcpp实现；模拟数据从SIR与SEIR的MJP生成每日发病数据，与逐受试者数据增强（sBDA）算法及PMMH算法（分别以τ-leaping算法与Gillespie算法模拟粒子）比较；实例数据源自2014年塞拉利昂Kailahun农村村庄埃博拉疫情（P=1161）和2009年加拿大Manitoba偏远社区H1N1疫情（P=3300）。

研究结果如下。模拟研究在多种人群规模P、再生成数r₀与病例检测率φ场景下，每个场景生成1000个数据集。三种算法在大多数场景中推断一致，均能无偏恢复真实参数值；计算性能评估显示，DA-MCMC所需CPU时间在所有测试场景中显著少于PMMH，尤其在小型暴发中优势明显；当φ接近1、观测噪声降低时，DA-MCMC计算效率改善而PMMH效率下降。P=500时DA-MCMC在所有场景的multi-ESS/CPU时间更高；P=1000时在φ=0.7场景更高。与sBDA相比，DA-MCMC在所有场景计算效率均显著提高。误设稳健性研究表明，当数据由r₀分四个阶段变化的SEIR流行病生成时，PMMH在约五分之一的重复中不收敛，而DA-MCMC保持计算稳健并给出合理参数估计。

在2014年塞拉利昂农村村庄埃博拉疫情应用中，以SEIR模型拟合周度确诊与疑似病例数，DA-MCMC与PMMH估计高度一致。平均潜伏期估计为7.83天（90% CrI：4.92–12.12天），平均传染期估计为8.19天（90% CrI：5.15–12.89天）；病例检测率φ估计为0.88，提示回顾性病例追踪近乎彻底；r₀估计为1.10（90% CrI：0.73–1.63），接近1表明村庄疫情规模有限，干预与社区行为调整有效控制了传播。在2009年加拿大Manitoba偏远社区H1N1疫情应用中，以SEIR模型拟合每日记录病例，结果再次显示两算法高度一致。平均潜伏期估计为2.33天（90% CrI：1.45–3.31天），平均传染期估计为2.37天（90% CrI：1.50–3.52天），r₀估计为1.55（90% CrI：1.26–1.85），约57%（90% CrI：34%–70%）人群在调查期结束时已感染或暴露，病例检测率仅0.11。针对第二波流行（2009年5月9日之后）的敏感性分析给出更高的r₀估计2.17（90% CrI：1.73–2.69）。

讨论部分指出，所提DA-MCMC算法在所有测试场景中计算时间少于PMMH算法；在涉及数百例感染的小型至中型暴发中，ESS/CPU时间与PMMH相当，但在更大暴发中效率较低。与Fintzi等（2017）的sBDA相比，通过联合采样个体流行病路径块，MCMC混合与计算效率得到显著提高。然而，该方法继承了个体水平DA方法的基本局限：增广参数空间维度随暴发规模扩大与模型复杂化而增加，拟合非常大规模暴发时MCMC混合可能显著变慢，这限制了逐受试者DA方法在超大规模人群中的实用价值；不过连续时间SEM精确推断的优势在小型至中型暴发中最为明显，此类场景中常微分方程、随机微分方程和离散时间SEM等近似模型可能导致有偏估计与错误的不确定性刻画。其他局限包括：观测模型采用简单的二项分布，现实监测可能需要负二项分布或移动平均模型等更复杂结构；仅依赖发病数据推断时模型参数可能出现弱识别性或后验多模态，需先验信息或其他数据来源；研究仅考察封闭且均匀混合人群、参数恒定的情形。在结论部分，论文指出：本文提出了高效的贝叶斯数据增强算法，为将连续时间SEMs拟合至存在不完美病例检测的发病数据这一具有挑战性的推断任务提供了解决方案；模拟研究表明所提出的DA-MCMC算法能够在多样的流行病场景中产生有效推断；该方法在所有测试场景中均比PMMH算法耗费更少的计算时间；在拟合涉及数百例感染的小型至中型暴发时，其计算效率（ESS/CPU时间）与PMMH算法相当，尽管在更大暴发中效率较低；模拟研究表明其计算稳健性优于PMMH算法；该方法的实际应用通过对2014年塞拉利昂农村村庄的埃博拉疫情和2009年加拿大偏远社区的H1N1疫情的分析得到说明；尽管本文仅展示了SIR和SEIR模型拟合发病数据的方法，该框架还可推广至其他SEMs及患病率数据等数据类型的拟合。