はじめに
新型コロナウイルス感染症(COVID-19)の収束に向けて、緊急事態宣言がなされ、諸外国のロックダウンほどでないにしろ、様々な業種の休業、イベントの中止、リモートワーク、各種自粛といったように、多くの人々の生活に影響が出ています。コロナのような感染症が収束に向かう条件は、実効再生産数（１人の感染者が生み出す２次感染者の総数）が持続的に１未満となることです。以前の記事では、実効再生産数を感染者数のデータから時系列で計算する方法の一つとして、潜伏期間と感染期間を分けて計算する簡便な方法を示しました。
今回の記事では、感染力、すなわち感染してから他者に染す強度は時間によって変化する効果を取り入れた場合にどうなるか、特に感染力のピークが感染後何日目になるか、について推定してみたいと思います。特に、報道などで潜伏期間中にも感染させるのではないか？という疑惑が報告されているので、その真偽の検証にもなるかと思います。
前提知識
本記事の検討の前提となる３つの研究を紹介します。
１つ目は、厚生労働省のクラスター対策班による調査です。クラスター対策班によれば、基本再生産数（基本とつく場合は環境因子を考慮しない、実効とつく場合は政策なども加味するみたいです）は次のグラフにようにサンプル集計されています。
大変ありがたいことに、平均値が1.49と明記されています（ここ重要）。すなわち、殆どの人が感染させない（0人）一方で、11人や12人というように一度に多くの人々に感染させてしまう場合があるのです。なので、R0は裾の広い分布と考えるべきでしょう。
２つ目は、イギリスのImperial College Londonの研究レポートNo.13です。この大学はCOVID-19に関するレポートを早くから公開しており、イギリスの政策にも強い影響を与えているようです。このレポートで注目すべきは、下図のように感染力に関してガンマ分布を使って表現している点です。
平均値は6.5[days]とありますが、ピークは4～5日目くらいにあるように見えます。
３つ目は、アメリカのJohns Hopkins 大学のCOVID-19の潜伏期間に関する報告です。これによると、潜伏期間の中央値は5.1日とのことです。このことから、２つ目の研究のグラフと合わせて考えると、発症日よりも前に感染力のピークがある可能性が示唆されます。
以上の３つの研究をベースに、特に日本の感染者数のデータを使って感染力のピークを検討していきます。
計算モデル
感染力を考慮した基本再生産数を計算するのに、以前の記事をちょっと改良した、SEIRモデルをベースとした、モデルを考えたいと思います。次の図を見てください。
以前の記事では、感染期間(I)にある感染者から等しい感染強度で２次感染が起こるとしていましたが、今回は下記の前提を採用します。
- 潜伏期間(E)中でも感染は起こりうる。
- ２次感染は感染者数と感染強度を積和した感染者数から生じる。
- 感染強度は感染日からの日数に依存した関数で与えられる。
- 潜伏期間(lp)と感染期間(ip)の和は13とする。
また、３の感染強度の関数は２項分布を採用することにします。ガンマ分布に比べパラメータが１つで済むため扱いやすいからです。さて、２項分布に基づく基本再生産数の式を次で定めます。
R_0(y) = r_0 \cdot {}_N C_y \cdot \theta^y (1-\theta)^{N-y}
$r_0$は定数、$\theta$は２項分布のパラメータ、$y$は感染日からの日数、$N$は潜伏期間と感染期間の和(=13)です。感染者から感染者への２次感染は次式のように畳み込みで定義します。$P(t)$は時刻$t$での確定陽性者数です。
P(t+N) = \sum_{y=0}^N R_0(N-y) P(t+y)
これらを変形すると、$r_0$を計算する式が得られます。
r_0 = \frac{P(t+N)}{\sum_{y=0}^N {}_N C_{N-y} \cdot \theta^{N-y} (1-\theta)^{y} P(t+y)}
分母にも分子にも$P(t+N)$がいるのがやや不自然ですが、感染したその日にすぐ２次感染者を生じさせる効果を入れた計算となっています。
さて、この$R_0(y)$に従って、感染強度のパラメータ$\theta$、つまり感染力が最大となる時間$N \theta$を決めたいのですが、以下の作戦で行きます。
- $\theta$を適当に決める。
- $R_0(y)$の式と実データから$r_0$を求める。
- $r_0$の分布の平均値を求める。
- 1から3を繰り返し、$E[r_0; \theta]$が1.49日に最も近い$\theta$を真値と推定する。
ステップ4は、クラスター班の分布が数値データとしてあれば、確率分布間の距離を計算して最適化するなどの方法も考えられたのですが、いかんせん画像でしかデータがないので諦めました…。csvで公開して欲しいですねぇ。（ﾁﾗｯ|дﾟ)
分析対象データ
都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社提供）において、公開されているcsvデータを使わせて頂きました。民間企業がボランティアでこのようなデータ集積と公開にご尽力されているのは大変素晴らしいことだと思います。
Pythonで計算してみる＆結果
今回は、計算しながら結果も同時に見ていきたいと思います。
1. 前準備
ライブラリとフォントの設定をします。
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import locale
#
from scipy.special import comb
font = {'family' : 'IPAexGothic'}
plt.rc('font', **font)
ジャッグジャパンさんのデータから都道府県別の感染者数をカウントする関数です。
def readCsvOfJapanPref(pref : None):
# 下記URLよりダウンロード
# https://jag-japan.com/covid19map-readme/
fcsv = u'COVID-19.csv'
df = pd.read_csv(fcsv, header=0, encoding='utf8', parse_dates=[u'確定日'])
# 確定日, 受診都道府県のみ抽出
df1 = df.loc[:,[u'確定日',u'受診都道府県']]
df1.columns = ['date','pref']
# 受診都道府県で抽出
if pref is None or pref == u'日本国内':
df1 = df1
else:
df1 = df1[df1.pref == pref]
df1 = df1.loc[:,'date']
# 確定日でカウント
df2 = pd.DataFrame( df1.value_counts() ).reset_index()
df2.columns = ['date','P']
df2 = df2.sort_values('date').reset_index(drop = True)
#
return df2
感染強度を考慮して基本再生産数を計算する関数です。ゼロ割対策を入れています。
def binomi(y, N, th):
c = comb(N, y)
p = c * pow(th, y) * pow(1. - th, N - y)
return p
def calcR0(df, keys):
# ２項分布
N = keys['N']
th = keys['th']
bcoef = [binomi(y, N, th) for y in range(N+1) ]
#
nrow = len(df)
getP = lambda s: df.loc[s, 'P' ] if (s >= 0 and s < nrow) else np.NaN
for t in range(1, nrow):
df.loc[t, 'Ppre' ] = sum([ getP(t + y) * bcoef[N - y] for y in range(0, N + 1)])
df.loc[t, 'Pat' ] = getP(t + N )
if df.loc[t, 'Ppre'] > 0.1: # ゼロ割対策
df.loc[t, 'R0' ] = df.loc[t, 'Pat'] / df.loc[t, 'Ppre']
else:
df.loc[t, 'R0' ] = np.NaN
return df
都道府県別の基本再生産数を計算する関数です。
def makeCalcFrame(days):
t_1 = pd.Timestamp(2020,1,24) # 計算開始日
td = pd.Timedelta('1 days')
#
npd = [[t_1 + td * i, 0, 0, 0 ] for i in range(0,days)]
df1 = pd.DataFrame(npd)
df1.columns = ['date', 'Ppre','Pat', 'R0']
#
return df1
def mergeCalcFrame(df1, df2):
return pd.merge(df1, df2, on='date', how='left').fillna(0)
def R0inJapanPref(pref, keys):
#
df1 = makeCalcFrame(100)
df2 = readCsvOfJapanPref(pref)
df = mergeCalcFrame(df1, df2)
#
df = calcR0(df, keys)
#
return df
2. ２項分布のパラメータθと基本再生産数R0平均値との関係
では、感染強度$R_0(y)$に含まれるパラメータ$\theta$ごとに、基本再生産数を都道府県別に計算して、その分布を見てみましょう。
計算のための関数です。パラメータ$\theta$の候補をグリッドでゴリゴリ計算していきます。
def calcDfMean(dflist):
mlist = []
for df, pref in dflist:
R0list = df[(df.R0 >= 0)]['R0']
m = np.mean(R0list)
mlist.append(m)
return np.mean(mlist)
def calcThToDfMean():
preflist = [u'日本国内', u'東京都', u'大阪府', u'愛知県', u'福岡県', u'北海道']
keys = {'N':13}
dc = {}
for y in np.linspace(0, 13, 14):
keys['th'] = y / keys['N']
dflist = [[R0inJapanPref(pref, keys), pref] for pref in preflist]
dm = calcDfMean(dflist)
dc[y] = dm
#
return dc
計算して表示する部分です。数分くらい時間がかかります。
dc = calcThToDfMean()
fig = plt.figure(figsize=(5,5))
ax1 = fig.add_subplot(1,1,1)
ax1.plot(list(dc.keys()), list(dc.values()))
ax1.set_xlabel('Nθ=E[y]')
ax1.set_ylabel('E[R0(y)]')
ax1.set_xticks(np.linspace(0,13,14))
ax1.set_xlim(0,)
ax1.set_ylim(0,)
ax1.grid()
ax1.set_title('The peak day vs the expected R0')
plt.show()
では、計算結果を見てみましょう。
横軸は$N \theta$、すなわち２項分布の平均値を、縦軸は基本再生産数$r_0$の平均値を表しています。さきほどのクラスター対策班のグラフでは、基本再生産数の平均値は1.49でした。
このグラフから、縦軸が1.49の交叉位置を見ると、$N \theta=5$であることが分かります。よって、パラメータ$\theta$を推定値を$5/N = 0.385$としましょう。
3. 推定パラメータθを使ったときの感染強度と基本再生産数
次に、推定パラメータ$\theta=0.385$を使った時の、基本再生産数の分布を見てみましょう。
分布を計算するための関数です。
def makeR0frame(dflist):
dc = {}
for df, pref in dflist:
R0list = df[(df.R0 >= 0)]['R0']
dc[pref] = R0list
return pd.DataFrame(dc)
def showR0bar(ax, df):
color = dict(boxes="Gray",whiskers="DarkOrange",medians="Black",caps="Gray")
df.plot.box(color=color, ax = ax)
def showBCoef(ax, keys):
N = keys['N']
th = keys['th']
ylist = [y for y in range(N+1) ]
bcoef = [binomi(y, N, th) for y in ylist ]
ax.plot(ylist, bcoef)
ax.set_xticks(np.linspace(0,N,N+1))
ax.grid()
ax.set_ylim(0,)
ax.set_xlim(0,)
分布を計算して表示する部分です。
preflist = [u'日本国内', u'東京都', u'大阪府', u'愛知県', u'福岡県', u'北海道']
keys = {'N':13, 'th':5./13. }
dflist = [[R0inJapanPref(pref, keys), pref] for pref in preflist]
fig = plt.figure(figsize=(10,5))
plt.rcParams["font.size"] = 10
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
#
df1 = makeR0frame(dflist)
showR0bar(ax1, df1)
showBCoef(ax2, keys)
ax1.legend()
ax1.set_title("R0")
ax2.set_xlabel('Nθ=E[y]')
ax2.set_ylabel('P(y;θ)')
ax2.set_title('binomi(y;θ)')
plt.show()
では、計算結果を見てみましょう。
左が都道府県別の基本再生産数の分布の箱ひげ図、右が感染強度の２項分布部分、です。
左図をみると、結構裾の長い分布になっていることが分かります。また、外れ値は北海道が最大で17.5付近にあります。パラメータを変えてみれば分かりますが、このパラメータでは都道府県別の違いが少ない印象です。
右図をみると、Imperial College Londonのレポートにあるガンマ分布と似ているような気もします。ピーク位置は５日目なのでややこちら結果の方が後ろでしょうか。
4. 基本再生産数のヒストグラムで答え合わせ
最後に、基本再生産数のヒストグラムを作り、クラスター対策班のヒストグラムと答え合わせをしてみます。
計算するコードです。
def showR0hist(ax, df, pref):
R0list = df[df.R0 >= 0]['R0']
R0list = R0list / R0list.sum() * 100
R0list.hist( bins=30, label=pref, ax=ax)
fig = plt.figure(figsize=(5,5))
ax1 = fig.add_subplot(1,1,1)
for df, pref in dflist:
showR0hist(ax1, df, pref)
ax1.legend()
#
plt.show()
計算結果を見てみましょう。
左が感染者数のデータから推定したR0、右がクラスター対策班の調査で集計したR0です。
形状はそこそこ合っていますね。指数分布で近似できそうな形です。平均値は一致しているので、ヨシ！としましょう。
考察
以上から、日本の各都道府県の感染者データに基づく、感染力のピークの推定に関して、次の傾向が計算結果から導けます。
- 感染力のピークは５日目、すなわち発症当日である（潜伏期間5.1日として）。
- ２項分布の形状から発症２日前から発症２日後の間で感染力が強いと思われる。
- 感染力のピークを５日目として都道府県別データから計算した基本再生産数のヒストグラムは、クラスター対策班の疫学調査のヒストグラムと平均値で一致し、形状もほぼ近い。
さらに言えば・・・
- 発症当日に感染力が強いという結果は、Imperial College Londonのレポートとも一致します。日本のデータとヨーロッパのデータで同じ傾向となったということは再現性がそこそこあるのではないでしょうか。
- クラスター追跡をするうえで、発症２日前から発症２日後を狙って追跡すると効率が良さそうな気がします。
- ２項分布で本当にいいのかは検討が必要でしょう。例えば、発症者のウイルス量を調査して、日々のウイルス量がどのような分布で推移するのかといったデータとも突き合わせてみればいいかもしれません。
- 計算結果で、パラメータ$\theta$が大きいほど基本再生産数$r_0$の平均値と分散が大きくなる傾向がありました。つまり、$\theta$と$r_0$は相関しているので適切に見極めることが大事です。その点、クラスター対策班の集計分布を正解として使えた点は重要です。
参考リンク
下記のページを参考にさせて頂きました。
- 新型コロナウイルスの都道府県別の基本再生産数の推移を計算してみる
- 新型コロナクラスター対策専門家（２）
- 新型コロナクラスター対策専門家が伝える「新型コロナクラスター対策用語定義」
- Report 13 - Estimating the number of infections and the impact of non-pharmaceutical interventions on COVID-19 in 11 European countries
- New Study on COVID-19 Estimates 5.1 Days for Incubation Period
- 【Python】組み合わせ(nCr) 計算の高速化
- 都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社提供）