基于嵌套粒子滤波的随机SIR模型统计推断方法及其在新冠数据中的应用

《Biostatistics》：A filtering approach for statistical inference in a stochastic SIR model with an application to Covid-19 data

本文针对传统流行病学模型难以捕捉感染率随机波动和未检测感染病例的问题，提出了一种离散时间随机SIR模型，采用贝叶斯方法和嵌套粒子滤波技术，在部分可观测环境下实现了对系统状态和参数的联合估计。通过奥地利新冠疫情数据的实证分析，验证了该方法在有效再生数估计和疫情预测方面的优越性，为公共卫生决策提供了可靠的技术支持。

　　

在传染病建模领域，传统确定性模型往往难以充分捕捉疾病传播过程中的关键特征。当人口规模较小时，疾病传播过程存在明显的随机性；感染率可能受到病毒传染性变化、环境条件、季节性和政策措施等多种随机因素的影响；此外，由于许多感染病例未被检测到，分析师无法完全观测到流行病学系统的所有组成部分。这些现实挑战催生了对更复杂建模方法的需求。

在这篇发表于《Biostatistics》的研究中，Katia Colaneri等人开发了一个随机且部分可观测的SIR模型，其中传播率和真实感染人数都是潜在变量。该模型通过隐马尔可夫模型（HMM）框架进行构建，包含两个组成部分：潜在马尔可夫过程（状态过程）和可观测过程（测量过程）。状态过程由SIR模型的舱室和传播率组成，而观测过程则由新确诊病例数给出。

研究人员采用贝叶斯方法进行统计推断，目标是近似状态变量和模型参数的后验分布。他们使用了随机滤波领域的先进技术，特别是将Crisan和Miguez（2018）提出的嵌套粒子滤波（NPF）适配到当前研究设定中。该方法通过两个嵌套的粒子滤波层：一个"外层"滤波，用于近似给定观测值的θ后验分布；一组"内层"标准bootstrap滤波，每个对应外层生成的一个样本，产生在给定观测值和θ样本条件下状态后验测度的近似。

在模型设定方面，研究考虑了两种版本：无隔离的基本HMM版本和包含隔离措施的扩展版本。在扩展版本中，检测呈阳性的人立即从感染者池中移除，这反映了隔离措施或自我隔离预防措施。未被检测到的感染者在感染康复后转移到移除状态。这种设定使得状态变量的演化依赖于观测值，因此该模型不再属于标准HMM类别。

关键技术方法包括：离散时间随机SIR模型构建、隐马尔可夫模型框架、贝叶斯推断方法、嵌套粒子滤波算法、后验预测分布预测技术。研究使用奥地利2020年5月1日至2022年6月15日的COVID-19感染数据进行分析，数据来源于奥地利健康与食品安全局（AGES）的公开数据。

模型设定与理论基础

研究首先建立了离散时间随机SIR模型的基本框架。关键变量包括：S_n（t_n时刻易感者数量）、I_n（t_n时刻可产生新感染的感染者数量）、I_n⁺（[t_n, t_n+1)区间内新感染人数）、I_n^-（从感染者中移除的人数）、R_n（移除者数量）、P_n（t_n时刻新报告感染数）和β_n（感染率）。

对数传播率Ψ_n = log(β_n)被建模为一阶自回归过程：Ψ_n = Ψ_n-1 + κ(μ - Ψ_n-1) + σZ_n，其中Z_n为独立标准正态随机变量。新感染I_n⁺被建模为泊松随机变量，其期望值与感染者数量和总人口中易感者比例成正比。

隔离扩展模型

在隔离扩展模型中，感染者的减少包括两部分：检测阳性被隔离者（P_n）和自然康复者（γI_n），即I_n^- = P_n + γI_n。观测过程保持不变，仍为P_n ～ Binomial(?I_n?, q)，但隔离引入的状态-观测依赖关系使模型超出了标准HMM框架。

有效再生数计算

研究推导了有效再生数R_n的计算公式：R_n = β_nS_n/[N(γ + q - γq)]。这一公式考虑了感染者的平均传播能力和移除前的预期时间，其中移除时间是康复时间τ^rec和隔离时间τ^quar的最小值，两者分别服从参数为γ和q的几何分布。

实证分析结果

将NPF方法应用于奥地利COVID-19数据后，研究获得了感染率β_n和有效再生数R_n的滤波估计。与奥地利健康机构AGES公布的官方估计相比，NPF方法得到的R_n估计显示出更高的变异性和更明显的峰值，这可能更接近真实情况。

参数θ = (κ, σ, μ)的后验分布随时间推移越来越集中。与模拟研究结果一致，参数μ的后验均值估计波动最小。从数据第二年开始，μ的后验均值估计似乎有所增加，这可能暗示参数实际上是时变的，特别是新病毒变体的出现可能开启了新的疫情阶段。

预测与模型检验

研究讨论了基于后验预测分布的预测和模型检验方法。通过NPF算法生成未来阳性检验的预测分布，可以计算各种分位数预测。使用Rosenblatt（1952）的经典结果进行正式的统计检验：如果预测分布设定正确，则随机变量?_j := F?_nj(P_nj+Δ)是独立同分布的标准均匀随机变量。

实证预测结果显示，预测分布具有高度偏态特征。以2021年12月20日（Omicron变体引发的感染波开始前夕）的预测为例，模型很好地预测了Delta变体病毒引发的感染波的下降，但由Omicron变体病毒引发的高感染数在预测期结束时达到了预测分布的90%分位数。到2022年1月17日，模型已经学习了不同的疫情阶段，实际病例数介于F_n,Δ的中位数和75%分位数之间。

研究结论与意义

本研究提出的离散时间随机SIR模型能够比传统确定性模型更贴近现实世界流行病动态，特别是在观测有限或延迟的环境中。NPF方法的递归特性使其能够有效纳入新信息来更新后验分布，这一特点在处理预测时特别有用，是NPF相对于迭代滤波等竞争算法的主要优势。

NPF的另一个优势是能够检测模型参数的变化点，例如由新病毒变体出现触发的变更，前提是这些变化相对较少发生。然而，对于随时间快速或连续变化的参数，NPF方法无法准确估计，因为该算法不是为处理高频参数变化而设计的。

研究也指出了模型的若干简化假设，如舱室间的转移可能导致非整数值，这在大型人群中不会影响结果，但在小群体应用中可通过将转移建模为泊松过程来改进。虽然研究重点在于随机传播率和未检测感染的影响，但模型可扩展纳入疫苗接种或死亡等额外舱室。

这项研究为传染病建模领域提供了重要的方法论贡献，特别是在部分可观测环境下进行状态和参数估计方面。通过结合随机SIR模型和先进的滤波技术，研究为公共卫生决策提供了更可靠的分析工具，有助于在未来的流行病应对中制定更有效的防控策略。