带有控制干预的随机流行病SEIR模型的统计推断：埃博拉病例研究

总结

为了从1995年刚果民主共和国埃博拉疫情的每日发病率和死亡率时间序列中估计参数，建立了一个传染病随机离散时间易感性-暴露-感染-再覆盖（SEIR）模型。关联时间序列显示出许多低整数以及零计数，需要采用内在随机建模方法。为了捕捉这种模型中隔间种群之间转移的随机性，我们指定了适当的条件二项式分布。此外，还引入了一个相对简单的时变传输率函数，以考虑控制干预的效果。我们发展了马尔可夫链蒙特卡罗方法进行推理，用于探索参数的后验分布。该算法被进一步扩展，以对未观测到的模型状态变量进行数值积分。这提供了一个现实的随机模型，流行病学家可以使用该模型来研究疾病的动力学和控制干预的效果。

1.简介

数学建模已成为了解传染病传播动力学的重要工具。最常用的理论框架是基于将人类宿主人群划分为包括易感人群、已感染但尚未感染（暴露）的人群、传染性人群和康复人群的类别。这些易感感染-再感染（SEIR）模型通常表示为微分方程系统（参见安德森和梅，1991年)其中，隔室之间的流速由特定于疾病自然史的参数决定。人们还认识到，随机建模很重要（参见，例如。，贝利，1975年;Ball、Mollison和Scalia-Tomba，1997年;安德森，1999)尤其是在疾病发病率很低的情况下，传播过程的离散性和随机性不容忽视。这就是本文研究的埃博拉暴发疫情。

传染病模型的推断因以下事实而变得复杂：第一，一个或多个模型变量可能不被观察到（潜在的），第二，数据通常在离散时间点可用，而潜在的真实过程在时间上是连续的。此外，模型参数可能会随着时间的推移而变化，例如，如果在流行病期间采取干预措施来控制疾病的传播。在连续观测数据的情况下，可以获得完整数据的参数估计(贝克尔，1976年)根据不完全数据的性质，还使用鞅方法开发了估计量(安德森和布里顿，2000年)，精确的前向-后向滤波(Fearnhead和Meligkotsidou，2004年)或马尔可夫链蒙特卡罗（MCMC）方法(Gibson和Renshaw，1998年;奥尼尔和罗伯茨，1999年;奥尼尔，2002年;Neal和Roberts，2004年;斯特里夫塔利斯和吉布森，2004年). 所有这些研究都基于这样的假设，即所有事件时间或其子集都可供研究人员使用。不幸的是，很少记录单个事件发生的时间。更常见的是，观察到的数据集是一段时间间隔（如一天或一周）内发生的事件计数的时间序列。

在本研究中，我们开发了在离散时间流行病数据可用于自然历史属于SEIR流行病模型的传染病的情况下的推断方法。当时间间隔等于疾病的发生期时，基于似然的推断（链式肿瘤模型）已被考虑为贝利（1975）和奥尼尔和罗伯茨（1999）在链式肿瘤模型中，假设发生期，即潜伏期和感染期是固定的。在这里，我们放宽了这个假设，允许潜伏期和感染期都是随机的，并且数据可以在距离可能不同于生成期长度的时间点观察。通过引入状态变量转移的概率密度，我们可以建立一个概率离散时间模型，在足够小的区间长度内，该模型可以很好地逼近随机SEIR模型生成的潜在连续时间过程。然后可以根据跃迁密度近似数据的似然函数。

根据可用的测量方法，处理传染病数据的研究人员面临着重要的观察变量和未观察变量的不同情况。第一步是构建待研究疾病的概率模型，并研究可用数据场景下的参数可识别性。MCMC的使用允许对未观测变量的分布进行数值积分。这对于使用标准误差实现参数估计很重要，标准误差实际反映了变量未被观测时增加的不确定性。我们通过将其应用于1995年刚果民主共和国埃博拉疫情的数据集来说明我们的方法的性能。Chowell等人（2004）仅考虑了此次疫情可用数据的一部分，用于估计一组具有时变传输率的确定性SEIR微分方程的参数，以允许进行控制干预。他们的估计方法是基于模拟确定性SEIR方程的解，并确定使观测和模拟累积案例数之间的平方误差之和最小的参数值。优化过程从10个不同的初始参数值开始，所报告的参数估计值是误差平方和最小的估计值。

在本文中，我们引入了一个完全概率模型，因此使用所有可用数据进行刚性的类似于hood的推理，同时考虑了未观察变量的不确定性以及报告过程中的错误。更一般地说，我们提供了一个模型，流行病研究人员可以使用该模型真实地模拟埃博拉疫情的随机动力学，以研究控制干预的效果和其他生物感兴趣的问题。

我们将在下一节中描述数据和建模方法第3节，基于模型的可能性进行地址推断。我们首先考虑这样一种情况，即每天都会观察到所有相关变量。事实上，易感个体感染的数量是一个潜在的过程，估计算法必须考虑这一点，以及一些变量报告不足的情况。将埃博拉数据应用于埃博拉病毒，可以得出表征该疾病自然史的参数估计值，并使我们能够调查发生的控制干预措施的有效性。

2.数据和模型

2.1数据

埃博拉出血热，俗称埃博拉病毒通过身体接触感染者的体液、分泌物、组织或精液传播(Chowell等人，2004年). 这种疾病的特征是最初的流感样症状，迅速发展为呕吐、腹泻、皮疹和内外出血。一旦接触，个体会经历大约6.3天的潜伏期，之后会感染，估计会持续3.5到10.7天(Breman等人，1977年). 大多数病例在首次感染后10天内死亡，该疾病的死亡率为50-90%(世界卫生组织，2003年). 有两种已知的埃博拉病毒株，即埃博拉-苏丹和埃博拉-Zaire，以它们最初被发现的国家命名(世界卫生组织，2003年). 在这里，我们研究了1995年刚果民主共和国埃博拉疫情的数据，该数据由Khan等人（1999）数据由两个时间序列组成（参见图1)记录的时间为3月1日至7月16日，即按症状出现日期统计的埃博拉每日病例数，共291例，以及埃博拉死亡日数，共236例。据记载，第一例于1995年1月6日发病，最后一例于7月16日死亡，共发现316例，死亡率为81%。5月9日，在对一些早期病例的标本进行检测时，埃博拉病毒被确认为病原体。立即采取了控制措施。其中包括使用防护服、积极监测和社区教育(Khan等人，1999年). 疫情持续了约200天，疫情开始后约130天采取了控制措施。3月1日之前疫情的确切开始时间和演变尚不清楚。此外，从316例确诊病例的总数中可以推断出，在给定的时间序列中，25例患者的症状出现日期和80例患者从感染类别中移除的日期没有报告。

2.2型号

考虑一个时间间隔(t吨，t吨+小时]，其中小时表示测量时间点之间的长度，此处小时=1天。让B类(t吨)表示被感染的易感个体的数量，C类(t吨)症状出现日期的病例数，以及D类(t吨)在此时间间隔内从感染类别中移除（死亡或恢复）的病例数。此外，让τ*表示流行病灭绝的时间点，即人口中没有暴露或感染个体的第一个时间点。让B类={B类(t吨)}^τ*_t吨=0表示的时间序列B类(t吨)从疫情开始到结束并定义C类和D类类似地。我们对随机连续时间SEIR模型使用离散时间近似（参见Gibson和Renshaw，1998年). 定义S公司(t吨),E类(t吨),我(t吨)，以及R（右）(t吨)作为当时人群中易感、暴露、传染性和迁移个体的数量t吨分别是。给定初始条件S公司（0）=秒₀，E类(0) =e（电子）₀，我(0) =一以及人口规模N个，离散化随机SEIR模型由

哪里

是具有二项式Bin的随机变量(n个，第页)概率分布：

参数β(t吨)、1/ϱ和1/γ分别是与时间相关的传播率、平均潜伏期和平均感染期。Mode and Sleeman（2000）对于具有三个疾病阶段或分区的SIR模型，导出了（5，6）中规定的二项式密度。其基本思想是，个体从疾病的前一阶段到下一阶段的转变被视为相应的人群分区之间的随机运动。在每一段时间里，一个人要么呆在隔间，要么走到下一个隔间。假设个人在隔间内花费的时间长度与特定隔间比率λ呈指数分布(t吨)，然后将逗留时间延长一段时间的概率小时是exp（-λ(t吨)小时)因此，离开的概率为1−exp（−λ(t吨)小时). 二项分布(5)假设每个Bernoulli试验对一个隔间的所有成员都是独立的和相同的，则这些试验的总和得出的结果。此外，需要注意的是，特定于分区的指数率是⁠以及γ（对于易感、暴露和感染的隔间）导致停留在隔间内的概率，如(6)（请参见Mode and Sleeman，2000年). 由此可见，潜伏期和感染期的指数分布近似于具有平均值1的相应几何分布/第页_C类和1/第页_R（右）分别是。对所有最新信息有条件t吨，二项式随机变量B类(t吨),C类(t吨)，以及D类(t吨)都是独立的。该模型进一步假设人口规模N个保持不变，个体混合均匀。

为了解释控制干预，我们假设传输参数β(t吨)在引入控制措施的时间点之前保持不变，之后呈指数衰减。这可以表述为

哪里t吨_*是引入控制措施的时间点，β是初始传输速率，以及q个>0是β的速率(t吨)衰变为t吨>t吨_*.Chowell等人（2004年）考虑传输速率的指数衰减，但引入第三个参数来进一步表征衰减。我们发现，对于实际样本大小，似然函数的几何结构不允许识别该参数，因为其估计值与指数相关q个。我们的传输速率模型(7)因此，是一种更节省的参数化。请注意，除非该疾病可以治愈，否则干预不会影响γ，而埃博拉并非如此。基本复制编号R（右）₀定义为当引入大量易感人群时，主要病例在其感染期内产生的次要病例的平均数量(Diekmann、Heesterbeek和Metz，1990年). 常量R（右）₀因此可以测量疫情的初始增长率，对于上面的模型，可以表明R（右）₀=β/γ (Hethcote，2000年). 此外，Chowell等人（2004）定义与时间相关的有效的生殖数每次感染病例的继发病例数t吨.因为S公司(t吨) ≈N个，因此是与时变传输速率成比例的函数(7)。时间点R（右）₀(t吨)假设值小于1表示控制措施何时在控制疫情方面有效。

中指定的流行病模型(1)–(6)以及接触率模型(7)参数向量θ={β，q个，ϱ，γ}，我们希望通过对初始条件、种群大小的了解以及对{B类，C类，D类}或其子集。有效值的时间演变R（右）₀(t吨)然后根据估计参数导出。

3.推断

如果初始条件、种群大小和向量{B类，C类，D类}以一定的时间间隔观察小时，那么这里我们说数据是完成。请注意{S公司(t吨),E类(t吨),我(t吨)}在这种情况下，完全由应用方程确定(1)–(3)根据给定的初始条件。因为B类(t吨),C类(t吨)，以及D类(t吨)是条件独立的，数据的可能性可以近似为

哪里克₁，克₂、和克_三代表中规定的二项式跃迁密度(5)和(6)以θ和所有最新信息为条件t吨.θ的最大似然（ML）估计，随后R（右）₀和R（右）₀(t吨)，可以通过最大化获得(8).

与实际情况一样，我们进一步假设B类未观察到。在本文中，我们在贝叶斯MCMC方法的框架内进行数据插补，这使我们能够在未观测过程的概率分布上进行数值积分。因为疫情一直持续到最后E类(τ*) = 0,我（τ*）=0，易感个体的最终数量由⁠疫情的最终规模，定义为最终感染该疾病的总人数，在本例中由米=S公司(τ*) −秒₀.系列{我(t吨)}也被称为可以通过以下方式从观测变量重建方程式（3）.系列{S公司(t吨),E类(t吨)}然而，这取决于未观察到的B类在马尔可夫链的每一次扫描中，我们现在对随机过程进行插补B类并重建{S公司(t吨)}和{E类(t吨)}使用(1)和(2).增加的可能性B类，C类、和D类是我(B类，C类，D类 | θ）如(8)将似然乘以先验π（θ），得到比例常数的后验分布

我们希望从中取样。从条件分布π中依次采样的MCMC算法(B类 | C类，D类，θ）和π（θ | C类，D类，B类)从所需的π（θ，B类 | C类，D类). 该算法的总体结构如下：

• 1

  初始化B类。这可以通过设置B类(0) =米，而向量中的所有其他位置都用零填充。对于任何初始B类系列{S公司(t吨),E类(t吨)}使用重建(1)和(2)分别是。

• 2

  初始化参数向量θ。

• 三

  更新B类从B类 | C类，D类，θ并计算新序列{S公司(t吨)}和{E类(t吨)}来自(1)和(2).

• 4

  从θ更新θ | C类，D类，B类.

• 5

  重复步骤(3)和(4)直到链收敛后获得所需的样本。

3.1从B中取样 | C、 D，θ

一种自然的求婚方式B类将对每个样品进行取样B类(t吨)根据其在每个时间点的条件二项分布。此外，检查提案是否与观察到的疫情最终规模和长度一致⁠、和我(τ*) +E类(τ*) = 0. 不幸的是，此类提案很少出现这种情况，因此该方案效率不高。相反，我们通过使用以下更新方案明确地将我们的建议限制在观测到的消光时间上。更新当前实现B类，选择一个时间点t吨′（满足B(t吨′）>0）从{0均匀随机，小时，第2页小时， … , τ*−小时}并设置B类(t吨′) =B类(t吨′) − 1. 然后，选择均匀随机设置⁠.更新系列{S公司(t吨),E类(t吨)}使用(1)和(2)，并检查是否满足上述所有条件。如果是，那么新配置B类'被视为B类我们接受B类'有概率在马尔可夫链的每次迭代中B类可以通过将过程迭代多次（大约为最终大小的10%）来更新。这可以改善MCMC算法的混合和收敛性(Neal和Roberts，2004年).

3.2从θ取样 | C、 D、B

我们使用随机游走方案更新θ的每个元素，其中高斯扰动的方差被调整为总接受率在20%到40%之间(罗伯茨和罗森塔尔，2001年). 初步分析表明B类在链之间引入负相关性q个我们发现，如果我们使用以二元条件分布模式为中心的二元正态建议，对它们的先前值进行联合独立采样，收敛性会得到改善（参见，例如。，Chib和Greenberg，1994年). 独立的伽马先验值被分配给θ中的每个参数，即π（ζ）～Γ（ν_ζ, λ_ζ)，式中ζ=β，q个、γ和ϱ，其中γ(一，b条)指具有参数的伽马分布一和b条，平均值一/b条、和方差一/b条².

3.3模拟数据的应用

首先，我们通过将算法应用于模型中的模拟疫情时间序列来测试和演示算法的性能(1)–(6)具有参数值秒₀= 5,364,500,e（电子）₀= 1,一= 0, β= 0.2,q个= 0.2, ϱ= 1/5 = 0.2, γ= 1/7 ≈ 0.143,小时=1天，干预时间t吨_*= 130. 这些参数值尽可能根据数据进行调整，因为它们是受早期埃博拉研究和已知埃博拉自然历史的影响(Breman等人，1977年;Khan等人，1999年;Chowell等人，2004年). 模拟疫情时间序列的最终规模为米=297，终止于τ*=172。

我们将上一节中描述的MCMC算法应用于模拟数据，假设B类并使用三个参数集对先验分布进行推断。第一个{（ν_ζ, λ_ζ); ζ=β, ϱ, γ,q个}={（2，10），（2，10.），（2,14），（2,10）}，第二个{（20，100。在前两种情况下，分布的平均值是真实的参数值，在后一种情况下先验值信息更丰富。在第三种情况下，每个分布的平均值为0.5，这与真实值相差很大。我们将把这个案例称为未输入的先前案例。优先分配q个在前两种情况下，由于很难预先知道干预的效果，因此选择的信息量较小。在每种情况下，在放弃10000个老化期后，使用每100次链迭代取1000个样本来获得后验分布。在马尔可夫链的每次迭代中B类已更新30次。通过对不同起始值的一系列运行，并检查自相关函数，评估马尔可夫链的收敛性。在所有三种情况下，马尔可夫链似乎在老化期后收敛。参数的后验均值和标准偏差报告于表1后验平均值与从完整数据中获得的最大似然估计值非常一致。可以预期，β、ϱ和q个更大，表明该算法包含了增加的不确定性水平，因为有效地，大约30%的完整数据是不可用的。参数γ的不确定性不受潜在数据的影响，因为涉及该参数的可能性分量仅取决于{我(t吨)}，可从观测数据中获得。表1表明估计方法对先验信息的选择不是很敏感，因为不同的选择导致参数估计与ML估计一致。同样的估计算法也被应用于输入大约10%的暴露个体感染的未知日期和大约25%的感染病例恢复的未知日期。此处未报告的结果表明，将这些未知日期与B类得出的后验均值与表1.

4.从埃博拉观察数据推断

4.1估算和结果

上述算法应用于实际埃博拉病毒数据集，我们在其中插补了未观察到的过程B类以及未知的日期，其中25例感染，80例被切除。我们假设最初有一个接触者，第一个病例在被诊断之前接触了6天（近似平均潜伏期）。参数约束由0<β<1给出；q个>0和以下Chowell等人（2004），我们从1<1/ϱ<21和3.5<1/γ<10.7中选择初始条件，以及它们的估计种群大小N个=5363500，如上述模拟所示。先验分布和抽样方法的选择与上一节完全相同。图2显示了参数的估计后验分布。后验均值和标准偏差报告于表2如模拟研究所示，对于不同的先验分布选择，结果没有显著差异。我们估计平均感染期的后验密度平均约为7天，后验标准差约为半天。我们还估计平均潜伏期的后验密度平均为10天，后验标准差约为1天。这些估计值均高于通过以下方法获得的关于这两个参数的5.5天左右的最小二乘估计值Chowell等人（2004）标准误差非常低，约为5小时。我们对传输参数β的估计为0.21，小于通过Chowell等人（2004）而我们对这个参数的后验标准偏差大约是他们报告的标准误差的三倍。我们对R（右）₀约为1.4，这与通过以下方法获得的值1.8大致相符Chowell等人（2004年）.我们估计的标准偏差R（右）₀然而，比标准误差高出两倍。总之，我们的参数后验均值与Chowell等人（2004）尽管我们使用的数据是原来的两倍，但我们的后验标准偏差远大于相应的标准误差。对于非统计学家来说，这可能被认为是一个缺点，但重要的是要注意到，这里的分析是完全概率的，首先是关于疾病传播的随机性和离散性，其次是关于所有未观察到的事件和变量。因此，存在额外的不确定性，这反映在参数的后验标准偏差中。上一节的模拟研究已经非常清楚地表明，潜在变量的存在会增加参数估计的后验标准差。

预计有效R（右）₀(t吨)控制干预5天后，曲线下降到小于1的值。疫情结束时，有效生殖数低至0.00268。这清楚地表明，通过采取干预措施，疫情得到了有效控制。使用基于估计参数后验均值的模拟，我们比较了有控制干预和无控制干预的情况。我们模拟了系统在某个时间点进行干预的情况下的500次疫情t吨_*=130和另外500个流行病，我们通过设置β(t吨) =β. 结果表明，通过干预，疫情持续时间从约950天减少到约200天，最终规模从约350万例（约占总人口的三分之二）减少到约300多例。

4.2模型诊断

例如，可以使用基于以下建议的差异度量的后验预测模型检查技术来评估模型拟合度Gelman、Meng和Stern（1996）这里我们选择最终尺寸作为差异度量。在马尔可夫链的每次迭代中，并给定其当前状态，我们计算流行病最终规模的结果。这产生了最终尺寸的后验预测分布，如图所示图3a如果模型符合数据，则观察到的差异度量不应位于该分布的尾部，在我们的情况下第页值0.441并不表示任何缺乏拟合的证据。模型的拟合也可以从以下方面看出图3b，其中观察到的累积病例数落在模型预测的95%后验可信区间内。

5.讨论

我们考虑了一种基于SEIR系统随机离散时间近似的概率方法，目的是对埃博拉疫情进行建模，并引入了一种用于参数估计的MCMC估计算法。由于模型还包含了在确定病原体时引入的控制干预，因此产生了其他参数。另一个不确定性来源是随机变量之一，即暴露在外的易感个体数量未被观察到。对于具有潜在变量的模型，适当的程序是对其概率分布进行积分，而MCMC提供了一种可行的数值算法。增加的不确定性反映在受潜在过程影响的任何参数的后验分布的更大范围内。虽然这是一种完全不同的方法，但我们已经将我们的结果与通过以下方法获得的参数估计值进行了比较Chowell等人（2004）他们应用最小二乘估计，目的是通过症状发作将SEIR微分方程拟合到观察到的日常病例中。与…对比Chowell等人（2004），我们还使用了可用的每日死亡率数据。这些与去除率相对应，但80个未观察到的去除率除外，我们的分析也可以解释这一点。我们发现，虽然模型参数的后验平均值仍具有可比性，但差异显著。当我们的后验标准偏差超过Chowell等人（2004）尽管我们使用的数据是原来的两倍，但还是有很多人在使用。然而，这并不奇怪，因为如上所述，对传输过程的概率性质所做的假设完全不同。我们对R（右）₀埃博拉疫情数据为1.36Chowell等人（2004）我们得出结论，这些干预措施在控制该病方面是成功的。我们估计，该流行病的无阻碍传播将持续约五倍的时间，并影响到三分之二的人口。然而，由于该病的死亡率如此之高，所以不能忽视的是，即使更早采取控制干预措施，也可以挽救更多的生命。

在开发上述估算方法时，我们做出了一些假设，如果该方法在其他数据场景中证明有用，则可以放宽这些假设。例如，假设从开始到结束都观察到疫情，这意味着最终的规模米这使得MCMC方法的实现更加简单。有可能从开始到某个时间点都能观察到疫情T型< τ*. 在这种情况下，最终尺寸未知，MCMC方法可以修改如下。让米_T型是（0，T型]. 此参数未被观察到，但我们知道因此，我们可以在MCMC的每次迭代中，对每个进行采样B类(t吨)根据相应的二项式分布，检查满足上述条件以及以下所有其他条件B类与观察到的疫情一致。提议的B类'然后可以像以前一样接受或拒绝。

这里开发的估算框架还可以扩展到其他更复杂的场景，例如，异质混合人口、转移密度的其他分布假设（参见莱利等人，2003年用于描述包含gamma-分布等待时间的模型）和其他类型的不完整数据。该模型的概率设置假设潜伏期和感染期是几何分布的，因为埃博拉文献中没有证据表明有任何替代的分布形式。尽管我们在这项研究中没有严重遇到参数可识别性问题，很明显，潜在变量可能会影响这一点。如果未观察到一个或多个变量，则可对完整数据进行可接受精度估计的参数可能无法识别。另一个重要的数据缺失来源是，随着时间的推移，对过程的测量过于稀疏。例如，如果我们假装小时如果是1周，我们发现参数估计是不精确的，因为它们有相当大的后验标准差。在这种情况下，将未观察到的过程归入（有规律或不规则间隔的）离散时间点之间的方法(Elerian、Chib和Shephard，2001年;达勒姆和加兰特，2002年)可以大大改进参数估计。上述模型的一个重要而现实的特点是，过渡密度的参数，以及条件均值和方差，都是随时间变化的。然而，随着小时在任何时间间隔内小时假设保持不变。因此，这种恒定速率近似的质量取决于系统的未知参数，因此需要一些专家判断来断言疾病的动力学在一个时间单位的长度内不应发生显著变化小时我们有充分的理由相信小时=1天足够小。目前，这似乎是人们希望报告此类疾病的最小时间单位。此外，如果我们有时间单位小于一天的数据，那么引入的SEIR模型必须包含更多细节，例如，感染性接触的强度不会在夜间和白天保持不变。

本文所述方法的实现是使用现代标准计算机，使用OX语言完成的Doornik（1999）根据作者的要求，感兴趣的读者可以获得代码。每次马尔可夫链运行的计算时间都小于45分钟。

致谢

作者感谢华威大学（华威研究生奖学金）和国家海外研究计划（ORS）资助莱科内的研究，并感谢两位匿名审稿人，他们的评论大大改进了手稿。

参考文献