解读：基于马尔科夫链蒙特卡洛方法对武汉市2019-nCoV爆发的潜在国内和国际传播的预测：一项模型研究-优快云博客

本文链接：https://blog.youkuaiyun.com/zx_zhang01/article/details/104390212

论文ID

原名： Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study
翻译名：对源自中国武汉的2019-nCoV爆发的潜在国内和国际传播的预测和预测：一项模型研究
期刊：The Lancet
IF:59.102
发表时间：2020.02
作者：Joseph T Wu*, Kathy Leung*, Gabriel M Leung
数据来源：中国疾病预防控制中心

文章简介

背景

自2020年1月23日开始，新型冠状病毒肺炎（2019-nCoV）开始逐渐引起了人民的关注。大量的科研人员开始了对新型冠状病毒（2019-nCoV）肺炎的研究。nature上的一篇论文指出截止到2月4日至少发表了77篇关于新型冠状病毒（2019-nCoV）肺炎的学术论文。本篇文章描述到：作者根据从武汉出口到中国大陆以外城市的病例数来估算武汉市的流行病规模，并建模预测流行病在国内和全球范围内造成的公共卫生风险（包括社会和非公共卫生风险），并提出来相关的药物预防干预措施。

结果

文章使用了马尔科夫链蒙特卡洛方法进行估算基本再生数 $R_0$ 为2.68.基本再生数在一定意义上表征了2019-nCoV的传播能力， $R_0$ 值越大，流行病的传播能力越强。基本再生数为2.68是什么概念呢？如下图所示。
在这里插入图片描述
通过将易感者——潜伏者——感染者——恢复者（SEIR）纳入metapopulation模型对中国主要城市的新型冠状病毒肺炎（2019-nCoV）进行模拟。经过模拟文章得到了截止2020年1月底估计有75,815 个人（95％CrI37304–130330）已经被感染。其中作为中国的交通重市，重庆，北京，上海，广州和深圳分别进入了 461（227-805），113（57-193）， 98（49–168），111（56–191）和 80（40–139）感染个人。
在这里插入图片描述

方法

文章首先基于国外城市的病例数，推断出了2019-nCoV的基本再生数和武汉市从2019年12月1日到2020年1月25日的八法规模；然后，进而估计出了从武汉市转移到中国大陆其他城市的病例规模。最后预测了2019-nCoV在中国大陆内部和外部的传播。从文章分析步骤来看，文章与流行病学的论文的思路基本一致。
文章中的研究病例定义为了：有症状的个体，或者进行住院治疗或者两者同时加上武汉旅游史。

模型建立

通过使用susceptible——eposed——infectious——recovered（SEIR）model对武汉流行病进行估算。
$\begin{aligned} \frac{\mathrm{d} S(t)}{\mathrm{d} t}=-\frac{S(t)}{N}\left(\frac{R_{0}}{D_{I}} I(t)+z(t)\right) &+L_{I, W}+L_{C, W}(t) \\ &-\left(\frac{L_{W, I}}{N}+\frac{L_{W, C}(t)}{N}\right) S(t) \end{aligned}$
$\begin{aligned} \frac{\mathrm{d} E}{\mathrm{d} t}=\frac{S(t)}{N}\left(\frac{R_{0}}{D_{l}} I(t)+z(t)\right) &-\frac{E(t)}{D_{E}} \\ &-\left(\frac{L_{W_{l}}}{N}+\frac{L_{w}(t)}{N}\right) E(t) \end{aligned}$
$\frac{\mathrm{d} I(t)}{\mathrm{d} t}=\frac{E(t)}{D_{E}}-\frac{I(t)}{D_{I}}-\left(\frac{L_{W I}}{N}+\frac{L_{W, C}(t)}{N}\right) I(t)$
其中 $S (t)$ , $E (t)$ , $I (t)$ , $R (t)$ 分别代表在时间t处易感，潜伏，感染，恢复个体的数量。
$D_E$ 和 $D_I$ 分别为平均潜伏期和传染期时间。 $Z (t)$ 代表在武汉华南贸易市场关闭之前人畜共病的案例数，此后为0。
假设旅行行为不受疾病影响，因此国际案例中是根据：
$\lambda(t)=\frac{L_{W, I}}{N}(E(t)+I(t))$
因此第d天国际案例出口的预期数量为：
$\lambda_{d}=\int_{d-1}^{d} \lambda(u) \mathrm{d} u$
两者之和的似然函数为：
$L\left(R_{0}\right)=\prod_{d=D_{s}}^{D_{e}} \frac{e^{-\lambda_{d} \lambda_{d} x_{d}}}{x_{d} !}$
文章使用马尔科夫链蒙特卡洛方法结合吉布斯采样（Gibbs sampling）（统计学中用于马尔科夫蒙特卡洛（MCMC）的一种算法，用于在难以直接采样时从某一多变量概率分布中近似抽取样本序列。该序列可用于近似联合分布、部分变量的边缘分布或计算积分（如某一变量的期望值）。某些变量可能为已知变量，故对这些变量并不需要采样。）使用后验方法表示点估计，使用95%的置信区间表示统计的不确定性。