SIR模型中繁殖率的不确定性

，30，并计算t处的基本繁殖数和有效繁殖数。图2中报告了这些路径，复制数等于S=100。我们观察到，在SIR的离散化版本中，即使是基本生产数量，即根据风险人口规模调整的数量，在流行病期间也不是恒定不变的。图2：基本生育率和有效生育率的演变0。81.01.21.40 100 200 300 400 500次（天）R0R0类型：R*0，t（基本）R0，t（有效）我们还观察到，有效繁殖数的最终水平与其起始值相等。事实上，对于大t，风险人口的规模与最终规模一致，然后是R(∞) = 空调也是。图2给出了未来路径N（t）在100天内的演化情况。实际上，总和可以被截断，这种截断可能会对R的评估产生影响。图3提供了分别在30天、60天和100天内计算的繁殖率的演变。图3：截断条件下繁殖率的演变0。751.001.250 100 200 300 400 500次（天）R0运行*0，t（30）R0，t（30）R*0，t（60）R0，t（60）R*0，t（100）R0，t（100）3估计。1挑战SIR模型的估计，以及更普遍的任何流行病学模型的估计，由于三个主要原因有点具有挑战性。i） SIR模型是一种具有混沌特性的非线性动力学模型[参见Harko，Lobo，Mack（2014）]。这意味着，从中长期来看，参数值a、c的微小变化，尤其是估计误差，可能会对过程的演变产生重大影响。众所周知[见Allen（1994）]确定的离散时间版本的SIR模型满足群体免疫力。

然而，在我们的随机框架中，群体免疫水平以及达到免疫水平的时间对a、c值和初始条件非常敏感。ii）疾病的演变是非平稳的，如图1所示。如果R0,0>1，受感染个体的比例增加到一个峰值，然后下降到一个渐近平稳状态。这种非平稳性使得很难将估计量的性质分析为观测日期T的函数。此外，在疫情开始时，T通常很小，在20-60天之间。iii）与前一点相比，横截面尺寸n非常大，当n趋于一致，T固定时，我们期望一个渐近理论。然而，命题2显示了二元分布B[N（t）的关键作用- 1） ，p（t）]和B[N（t）]- 1） ，c]，t=1，T对于无症状分析，重要的不是n，而是边缘计数sn（t）-1） ，N（t）-1). 虽然易感人群通常非常多，至少在疾病开始时是如此，但感染者的数量却少得多。然而，对于大N（t- 1） ，N（t）- 1） ，我们可以应用二项分布的标准渐近结果。也就是说，通过泊松分布或高斯分布来逼近它的可能性。因此，B（N，p）的近似值是p（Np），如果N→ ∞, P→ 0，这样的Np→ λ>0，或N[Np，Np（1- p） ]如果N→ ∞, p固定。在我们的框架中，p（t）=c和p（t）都很小。两种方法之间的选择取决于N（t）的大小-1） p（t），N（t）-1） p（t），t=1，T也就是说，新感染者和新康复者的数量分别为。粗略地说，如果它们小于45-50，比如说，可以使用泊松近似，否则使用高斯近似。

但在疫情开始和结束时，N（t）和N（t）都相当小，而在疫情高峰期则更大。因此，近似值将取决于观测日期。他们还依赖于感兴趣人口的规模。例如，如果我们想考虑多伦多的一个亚群体，比如说，比75.3.2机械模型年长的男性，那么这个规模就更小。文献的主要部分是基于一个确定性动态模型，该模型隐含地假设了通过频率对应物近似理论转换概率的可能性。也就是说，使用高斯近似。更准确地说，在假设A.1下，我们有：-1^p（t）=p[^p（t）-1） ]p（t-1). （3.1）因此如果^p（t）~ p（t），我们得到p（t）s的以下确定性动力学模型：p（t）=p[p（t-1） ]p（t-1). （3.2）这通常被称为机械模型[参见Breto等人（2009）和附录1，因为它与连续时间SIR模型的联系]。3.3（近似）最大似然估计在我们的框架中，对数似然函数L（a，c）可以分解为和L（a，c）=L（a）+L（c）。这使我们能够通过分别关注（观察到的）转换矩阵的第一行和第二行，分别估计A和c（见附录2）。不同的对数似然函数也可参见Zhang et al.（2020）的分析，该分析仅限于钻石公主号游轮上的流行病分析。例如，基于二项分布的真实分布，或基于泊松或高斯近似的近似分布。3.3.1二项式对数概率我们有：L（a）=TXt=1{N（t）log[1- a^p（t）- 1） ]+N（t）log[a^p（t）- 1） ]}，（3.3）L（c）=TXt=1{N（t）log（1）- c） +N（t）log c}。

（3.4）a的最大似然估计是一阶条件的解：-TXt=1N（t）^p（t）-1)1 - ^a^p（t）-1)+^aTXt=1N（t）=0，（3.5），没有封闭形式的表达式。c的最大似然估计量为：^c=TXt=1N（t）/TXt=1N（t）-1） =TXt=1N（t）-1） TXt=1N（t-1） ^p（t）. （3.6）这是注明日期的过渡频率的加权组合。3.3.2泊松近似对数概率：LP（a）∝TXt=1{N（t）log[aN（t- 1） ^p（t）- 1)] - 安（t）- 1） ^p（t）- 1） （3.7）LP（c）∝TXt=1{N（t）log[N（t）- 1） c]- N（t）- 1） c]}。（3.8）我们得到了泊松近似最大似然（AML）估计量，其封闭形式为：^aP=nTXt=1N（t）/TXt=1[N（t）-1） N（t）-1） ，（3.9）^cP=TXt=1N（t）/TXt=1N（t）-1） =^c.（3.10）第一个公式表明，^aPis是日期估计的过渡系数的加权平均值：^at=N（t）/[N（t）]-1） ^p（t）-1） ]中，权重与N（t）成比例-1） ^p（t）-1).我们推导了初始生殖数的相应估计量的解析公式：^R0，P=nTXt=1N（t）TXt=1N（t）-1） TXt=1[N（t-1） N（t）-1） ]TXt=1N（t）。（3.11）如果xt=1N（t）非零，也就是说，如果观察到回收率，则可以使用该公式。3.3.3高斯近似对数概率：LG（a）∝ -TXt=1log（a^p（t-1)[1 - a^p（t）-1)]) -TXt=1N（t-1） [p（t）-a^p（t）-1） [a^p（t）-1)[1 - a^p（t）-1） [，（3.12）LG（c）∝ -Tlog[c（1- c） ]-TXt=1N（t-1） [p（t）-c] c（1）-c） 。（3.13）3.3.4不可行高斯近似对数似然通过替换方差a^p（t）获得近似对数似然-1)[1 -a^p（t）-1） ]通过估计^p（t）[1- ^p（t）]。我们得到：吊耳（a）=-TXt=1N（t）-1） （p（t）-a^p（t）-1） [p（t）（1）- ^p（t）]. （3.14）我们得到了^aUG的闭式表达式，它对应于a:^aUG=TXt=1（N（t）的不可行广义最小二乘（GLS）估计量-1） ^p（t）-1)/[1-^p（t））/TXt=1N（t）-1） ^p（t）-1） ^p（t）[1- ^p（t）].

（3.15）3.3.5泊松/高斯近似对数概率当n大、p小和np大时，泊松分布p（np）可近似为高斯分布n（np，np）。因此，与3.3.3中的近似值相比，方差中的项被忽略。我们有：LP G（a）∝ -TXt=1log[a^p（t-1)] -TXt=1（N（t-1） [p（t）-a^p（t）-1） [a^p（t）-1） ）（3.16）LP G（c）∝ -T log c-TXt=1N（t）-1） [p（t）-c] c（3.17）当n趋于一致时，这可能不一致。然后AML估计是2次多项式方程的正解，即：TTXt=1{N（t-1） ^p（t）-1） }a+a-TTXt=1{N（t-1） ^p（t）}=0，ttxt=1N（t）-1） nc+c-TTXt=1{N（t-1） ^p（t）}=0。总之，我们得到的a，c和初始复制数R0，0=a/c的AML估计数与（近似）对数似然数一样多。这也解释了R0,0的不同近似值，即使适用于同一系列的总计数。3.4反洗钱估计器的性质反洗钱估计器的性质可通过第4节中所述的蒙特卡罗方法得出。它们的渐近性质取决于泊松或高斯渐近（取决于哪一种最合适）以及所选的估计量。例如，当高斯渐近条件满足时，我们可能选择了PoissonAML估计。在这种情况下，B（N，p）很好地由N[Np，Np（1）逼近-p） ，它被p（Np）取代，p（Np）与N（Np，Np）相近。

因此，我们没有使用正确的高斯近似，忽略了inp项。为了说明，我们考虑以下两种情况：i）当泊松渐近有效时，泊松AML估计量^aP的行为。ii）当高斯渐近有效时，二项式ML估计量^a的行为。3.4.1 Poisson AML和Poisson渐近性让我们考虑T=1的情况，即对聚集的两个观察。以下主要结果适用于任何特定情况。然后我们有：^aP=nN（1）/N（0）N（0），^cP=N（1）/N（0），^R0，P=^aP/^cP=N（1）N（1）nN（0）。（3.18）在[N（0），N（0）]的条件下，估计量^apa和^cpa是独立的，因此N（0）N（0）N^aP~ P[aN（0）N（0）N]，N（0）^cP~ P[cN（0）]。我们推断：E^aP=a，E^cP=c，这表明泊松估计量对于T=1是无偏的。它们的方差为：V^aP=anN（0）N（0），V^cP=cN（0）。（3.19）在实践中，N（0）\'N和N（0）对于泊松渐近有效性来说是相当小的（<30或40）。因此，V（^aP）和V（^cP）都不是小的，即使对于大的n，我们也不能期望在泊松渐近下，对于nlarge，V（^aP）和V（^cP）的一致性。此外，在疫情刚开始时，受感染的个体尚未恢复，这意味着N（1）=0。我们推断：^R0P=^aP/^cP=^aP/0=∞. 这说明在暴发的初始阶段，基本生产比率缺乏准确性。备注1：无偏性特性适用于T=1的情况。如果t=2，我们有：^aP=n[n（1）+n（2）]n（0）n（0）+n（1）n（1）我们推断出它在日期1时的期望：E（^aP）=nN（1）+aN（1）n（1）n（0）n（0）+n（1）n（1），通过迭代期望，E（^aP）=nEN（1）+aN（1）N（1）N（0）N（0）+N（1）N（1）,这是计数N（1），N（1），N（1）的复杂非线性函数的期望。3.4.2二项式ML和高斯渐近这是适用大数定律和中心极限定理的标准渐近理论。

样本频率倾向于它们的理论对应：^pjk（t）→ pjk（t），^pj（t）→ pj（t），j，k=1，2，3，当n趋于一致时。ML估计器倾向于真实参数值：^a→ a、 ^c→ c、 ^R=^a/^c→ a/c，速度1/√n、 ^a，^c是渐近独立的，渐近正态的，它们的方差由以下公式一致估计：^V（^a）=（TXt=1N（t）^p（t）-1)[1 - ^a^p（t）-1)]+^aTXt=1N（t））-1，（3.20）^V（^c）=^c（1）- ^c）TXt=1N（t-1). （3.21）备注2：高斯渐近也可以应用于其他AML估计，如泊松AML。在这种情况下，isstill的Poisson AML估计是一致的，渐近正态的。然而，由于近似对数似然是错误的，其渐近方差是通过一个三明治公式获得的，该公式涉及信息矩阵的两个表达式[seeHuber（1967）]。4蒙特卡罗研究即使可以使用高斯渐近，我们也不知道它们在确定不同参数a、c、R的置信区间时是否准确，我们对第3节中介绍的一些估计器进行蒙特卡罗分析。我们将设计计算如下：N（0）=3000000，N（0）=1001000，T=20，c=0.07，R=2这对应于在[0，T]期间计算的估值。注意，边际计数的过程是马尔可夫过程。因此，它也适用于在（t，t+t）上计算的RollingSimator，例如，t处的边际计数是N（0），N（0）的固定计数。这解释了为什么我们在设计中允许N（0）的大值。