感染症のデータ分析
このページでは、確率モデルを使って、感染間隔の分布や、感染症の感染が確認された人数の推移から新規発症者数を推定する方法などを考える。
発症間隔(serial interval)の分布の推定
† Nishiura et al., Int. J. Infec. Dis. 93, pp. 284-286 (2020).
ある人が感染した日を起点に、その人が新たな感染を生じさせるまでの間隔
3, 2, 4, 2, 2, 1, 4, 2, 10, 4, 3, 4, 4, 3, 3, 9, 9, 3, 4, 3, 9, 4, 5, 9, 7, 7, 3, 7
という28の時間間隔の例を用いた解析が行われている。 ここでは、簡単のため、それぞれのデータの重みは等しいと考え、これらの間隔がどのような分布形に従っているのかを推定してみよう。
上記の論文を参考に、その分布がWeibull分布
# coding: utf-8 import math import random import numpy as np import pymc3 as pm import matplotlib.pyplot as plt intervals=[3,2,4,2,2,1,4,2,10,4,3,4,4,3,3,9,9,3,4,3,9,4,5,9,7,7,3,7] with pm.Model() as model_weibull: alpha = pm.Lognormal('alpha') beta = pm.Lognormal('beta') tau = pm.Weibull('tau', alpha=alpha, beta=beta, observed=intervals) trace_w = pm.sample(chains=1, tune=2000, draws=1000) a_m = np.mean(trace_w['alpha']) a_s = np.std(trace_w['alpha']) b_m = np.mean(trace_w['beta']) b_s = np.std(trace_w['beta']) print('alpha=',a_m,'+/-',a_s) print('beta=',b_m,'+/-',b_s) def weibull(x,alpha,beta): return (alpha/beta) * (x/beta)**(alpha-1) * np.exp(-(x/beta)**alpha) x = np.linspace(0,20,100) y = weibull(x,a_m,b_m) plt.plot(x, y, color=(0.0,0.0,1.0), linewidth=1.0) plt.xlabel('TAU (DAYS)') plt.ylabel('F') plt.grid(True) plt.show()
練習: モデルの比較
Weibull分布以外の連続分布(対数正規分布やガンマ分布等)について、同様の推定を行ってみなさい。 さらに、PyMC3のマニュアルを参考に、モデルの「良さ」を比較してみなさい。
ヒント
PyMC3で使われているWAICは、こちらのページで説明されているAICのある種の拡張版である。
発症から確認までの遅延
次に、新型コロナウイルスの感染拡大の様子を想定しながら、報道発表などで知られる「その日の新規感染者数」と、 実際に発症した日とのずれについて考察してみよう。
実際のデータは離散間隔(日毎)で集計される場合が多いので、
ここで、沢山のサンプルから、遅延の分布型はあらかじめわかっているとしよう。
特にそうしたエビデンスがあるわけではないが、以下ではポアソン分布を仮定することにする:
感染確認日からの感染日の確率的な推定
次に、
ナイーブには、連立方程式
上記のlinalg.pinv()
が用意されているので、以下のように計算すればよい:
X = (np.linalg.pinv(A)).dot(Y.T)
確率モデルによって考える
そこで、新規感染者の時系列
内部モデル
内部モデルとして、
観測モデル
この
2020年1月25日から5月9日にかけての、東京都のCOVID-19の感染者数の時系列を例に、
PyMC3でこのモデルを記述してみた例を以下に示す。
オンセットから確認までの平均日数は14日とした(lambd=14
)。
# coding: utf-8 import math import random import numpy as np import pymc3 as pm import theano.tensor as tt from theano.compile.ops import as_op # 2020/1/25〜5/9の東京都の陽性者数(累計) y_cum =[ 2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,6,14,19,19,22,25,25,28,29,29, 32,32,35,36,36,37,39,39,40,44,52,58,64,64,64,67,73,75,77,87,90,90,102, 111,118,129,136,138,154,171,212,259,299,362,430,443,521,587,684,773,889, 1032,1115,1194,1338,1516,1703,1900,2066,2157,2316,2442,2590,2791,2972,3079, 3181,3304,3436,3570,3731,3834,3906,3945,4057,4104,4150,4315,4473,4564,4651, 4709,4747,4770,4809,4845 ] y_cum = np.array(y_cum) y = np.diff(y_cum) N = len(y) M=63 kernel = np.zeros((M,)) lambd=14 tot = 0 for k in range(M): fk = lambd**k / math.factorial(k) tot += fk kernel[k] = fk kernel = kernel/tot kernel = np.roll(kernel, M//2) @as_op(itypes=[tt.dvector],otypes=[tt.dvector]) def convolve(x): return np.convolve(x,kernel,mode='same') model = pm.Model() with model: z_0 = pm.Normal('z_0') z_steps = pm.Normal('z_steps', sigma=0.1, shape=N-1) zs = tt.concatenate([[z_0], z_steps]) z = zs.cumsum() x = pm.Deterministic('x',pm.math.exp(z)) mus = pm.Deterministic('mus', convolve(x)) cases = pm.Poisson('cases', mu=mus, observed=y) trace = pm.sample(chains=1, tune=2000, draws=1000) pm.traceplot(trace)
上記のコードで畳み込みは、関数 convolve()
で行っているが、以下の点に注意が必要である:
- NumPyの畳み込みは、モードによって動作が異なる。ここでは、
mode="same"
を選択。 その際、積和を取る順番が配列の順とは逆転していること、配列要素の対応関係がずれている点に注意。 - モデルの中では関数にオブジェクトとしてデータが渡されるので、関数宣言の前にデコレータ
@theano.compile.ops.as_op
を付けておく必要がある。
続いて、サンプリングの結果を出力するコードと出力例を示す:
import matplotlib.pyplot as plt xm = np.mean(trace['x'],axis=0) xd = np.std(trace['x'],axis=0) days=range(0,len(y)) plt.plot(days, y, color=(0.0,0.0,1.0), linewidth=1.0) plt.errorbar(days, xm, yerr=xd, marker='.', color=(1.0,0.0,0.0), capthick=1, capsize=2, linewidth=0.5) plt.xlabel('DAYS') plt.ylabel('X') plt.grid(True) plt.show()
練習:死亡者数を用いた推定
データは 朝日新聞社のウェブサイトによる。
以下は、2020年2月24日以降の、日本でのCOVID-19による死亡者の累積データである。
# 2020年2月24日から5月19日までの日本の累積死亡者数 death_tolls=[ 1,1,1,1,1,1,1,1,1,1,1,1,1,3,4,5,5,6,6,6,6,6,6,6,7,9,12,15,19,21,22,24,28,29,31, 33,35,36,41,42,43,45,46,49,52,54,56,57,60,63,69,70,73,80,81,85,88,94,98,102,109, 119,136,148,154,161,171,186,277,287,317,334,348,351,376,389,415,432,458,492,510, 521,543,551,557,600,613,621,643,668,687,710,725,744,749,763 ]
これを元に、感染から死亡までの日数や致死率などに仮定を設けて、感染日ごとの感染者数を推定してみなさい。