如何计算基本再生数R0?

今天不写计量经济学了,介绍一种计算传染病基本再生数R_{0}的方法。

R_{0}的计算方法有很多,我把它们大致分为两类:一类是通过数学推导直接给出R_{0}的计算公式,基于现有数据直接算出R_{0}的值,比如Shi Zhao等人发表在《International Journal of Infectious Diseases》上的论文;还有一类则建立模型,将R_{0}作为模型的一个参数,而后基于数据资料将参数估计出来,近期有代表性的就是Joseph T. Wu等人发表在《The Lancet》杂志上的那篇论文。

本文主要介绍一种基于数学推导的R_{0}计算方法,详细推导过程可见Wallinga and Lipsitch (2007) 这篇文章。这种计算方法首先假定了每天新增的被感染人数仍然在指数增长,所以更适合被用来计算传染病出现早期的R_{0}值。

基本再生数R_{0}(basic reproductive number)的值表示“一个(典型的)病例在一个全都是易感者的人群中直接造成的新感染者的平均人数”。

与其相对的是另一个再生数,被称为“有效再生数(effective reproductive number)”,一般记为R,它与R_{0}含义类似,只不过“全都是易感者的人群”被替换成了“有部分免疫者存在的人群”,显然它的实际应用更加广泛。本文要介绍的算法对两种再生数都兼容,因为推导过程不含有关于人群是不是完全由易感者组成的假定,因此我直接用R来表示。


首先给出计算公式:

R=(1+r\bar{T_L})(1+{r\bar{T_I}}) (1)

其中r表示指数增长的增长率,如果用 b(t) 表示第t天的新增被感染者人数,则:

b(t)=b(t-\Delta t)e^{r\Delta t} (2)

显然,如果我们放心大胆地假设每天新增的被感染人数是指数增长的,那么你将有100种方法来增长率r,或者从前几天推,或者从第一天推,或者...反正自己算吧。

\bar{T_{L}}\bar{T_{I}} 都来自对现实情况的观察, \bar{T_{L}} 表示被感染者的潜伏期的平均长度, \bar{T_{I}} = \bar{SI} - \bar{T_{L}} 。下标L表示Latent,I表示Infectious。SI表示Serial Interval,也叫Generation Time,它的定义是“一个感染者被感染的时间和他/她感染的下一个人被感染的时间的间隔”。

这个定义可能有点绕口。Serial Interval这个名词在人口学等学科中也在被使用,它也可以表示一位女性的出生和她生孩子这两件事之间的时间间隔,对照一下可能就更好理解SI在我们的场景下的含义了,它在整个计算中非常重要,后文还会提到。

顺便提一下,《The Lancet》上那篇文章使用的 \bar{SI}\bar{T_{L}} 分别为8.4天和6天,参考了当年的非典以及中东呼吸综合征爆发时的数据;CDC发在NEJM上的那篇论文估计的 \bar{SI}\bar{T_{L}} 分别为7.5和5.2天。


那么(1)式是怎么来的呢,使用的时候有什么要注意的呢?事实上,R的计算公式的得出有一个重要基础,那就是对SI的分布的假定。如果对SI的分布的假定不同,我们会得出不同的计算R的公式。

感兴趣的朋友可以认真读一下Wallinga and Lipsitch (2007) 这篇文章以了解整个推导过程。这里我就不再从头说起了,只说(1)式的直接来源,也就是下式:

R=\frac{1}{M(-r)} (3)

M(·)这样的函数叫做“moment generating function”,即“动差生成函数”或者“矩量母函数”。如果我们记SI的概率密度函数为g(t),其中t表示从被感染开始经过的天数,那么M(-r)可以被展开为:

M(-r)=\int^{\infty}_{t=0}e^{-rt}g(t)dt (4)

(4)式看似复杂,但如果我们已知了(或假设了)g(t)的函数形式,也就是SI的分布,那么它马上可以被化简为我们喜闻乐见的简单形式(可以参考这篇文档blog.csdn.net/STcyclone)。

现在我们考虑以下三种SI的分布:

1)假设SI服从均值为 \bar{SI} 的指数分布,那么它的函数表达式即为:

f(x)=\lambda e^{-\lambda x},\space x>0 (5)

将其带入(4)式并把 \lambda 替换成 \frac{1}{\bar{SI}} ,我们可以得到:

R=1+r\bar{SI} (6)

2)进一步地,如果考虑存在潜伏期 T_{L} ,这时候我们假设SI的两个部分 T_{L}T_{I} 各自服从着均值分别为 \bar{T_{L}}\bar{T_{I}} 的指数分布,那么SI就成为了两个两个指数分布的“叠加”,它服从一个我也叫不出名字的分布,但带入整合之后的结果是非常简洁的:

R=(1+r\bar{T_{L}})(1+r\bar{T_{I}}) (1)

这就是开头给出的计算公式。至于为什么我们要假设 T_{L}T_{I} 各自服从指数分布,或者为什么是指数分布,这些问题会比较复杂。简单说,人们常用的传染病模型如SIR、SEIR模型也都有着一些假设,我们这样做正好与这些模型的假设实现了一致。当然,你也可以假设SI服从正态分布或者别的分布,也都有人这样干。

3)如果我们得到了 \bar{SI} 并假设SI服从指数分布,且我们承认潜伏期的存在,但没有得到 \bar{T_{L}}\bar{T_{I}} ,该怎么办呢?可以假设 T_{L}T_{I} 是独立同分布的,即二者总体来说在SI中各占一半,这样一来SI其实就服从了gamma分布(指数分布是gamma分布的一种特例)。Shi Zhao等人发表在《International Journal of Infectious Diseases》上的那篇论文就是假设了SI服从gamma分布。


以上这个方法就介绍完了,总结一下,使用这个方法有几个要点:

1)一开始要假设每天新增的被感染人数是指数增长的

2)新增被感染者数量的增长率r要能被算出来,当然,在假设1)的背景下r非常好算

3)要有对SI的分布的假设,当然,在假设了相应分布后,有关分布的参数要能被获得

之前提到了还有一种通过建立模型,将R_{0}设置为模型参数,并用数据资料估计出参数的方法,这类方法虽然同样受制于相应模型的一些基本假设,但总的来说会更加灵活。比如《The Lancet》上的那篇文章在建模时将武汉市当时的人口流入、流出同时考虑在内。不过,对传染病进行建模是一项危险的工作,它对数据资料的完整性和质量要求太高,100个人可以建立出100种模型,但真正有参考价值的寥寥无几。


参考文献:

Li, Q., Guan, X., Wu, P., Wang, X., Zhou, L., Tong, Y., ... & Xing, X. (2020). Early transmission dynamics in Wuhan, China, of novel coronavirus–infected pneumonia. New England Journal of Medicine.

Marc, Lipsitch, Ted, Cohen, Ben, & Cooper, et al. (2003). Transmission dynamics and control of severe acute respiratory syndrome. Science, 300 (5627), 1966-1970.

Wallinga, J., & Lipsitch, M. (2007). How generation intervals shape the relationship between growth rates and reproductive numbers. Proceedings of the Royal Society B: Biological Sciences, 274(1609), 599-604.

Wu, J. T., Leung, K., & Leung, G. M. (2020). Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. The Lancet. 2020; S0140-6736(20)30260-9.

Zhao, S., Lin, Q., Ran, J., Musa, S. S., Yang, G., Wang, W., ... & Wang, M. H. (2020). Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak. International Journal of Infectious Diseases.

编辑于 02-15