はじめに

　新型コロナウイルス感染症(COVID-19)の感染の中心地はもはや中国本土から欧米に移行しつつあり、EU域内でも自由な国境の移動を制限したり、外出禁止令が発令されるなどの緊急事態になっています。日本において、新型コロナウイルスの脅威が認識されはじめたのは、2020年2月3日のダイヤモンド・プリンセス号の横浜港入港、2月25日の政府の基本方針表明、2月28日の北海道における緊急事態宣言、などが契機ですが、具体的な対策として2月26日から大規模なイベントの自粛要請、3月2日からの全国の小中高校の臨時休校などがなされて来ました。
　それから約半月が経過しますが、その間、日経平均株価が22000円台から17000円台まで急落するなど、先行きに暗雲が立ち込めています。新型コロナウイルス感染症の影響は健康のみならず、社会経済に深刻な打撃を与えており、いつ収束するのかが極めて重大な関心を集めています。
　そこで、本記事では、１月末から３月頭にかけて、日本における新型コロナウイルスの感染力の指標である、基本再生産数がどのように推移してきたか、計算で求める一つの方法を提示したいと思います。
　計算式やコードは少し長いので、計算結果を先に見て頂いた方が良いかもしれません。

基本再生産数とは？

　基本再生産数は、実は非常に奥が深いようで、京大数理解析研究所講究録や日本数理生物学会レター等を見ると、世代間の変遷を考慮するなど単純ではないようですが、敢えて、一言でいうと、

人口集団の典型的な１人の感染者が、その全感染期間において再生産する２次感染者数の期待値

と定義されるようです。この数値をR0とすると、感染が拡大する条件は$R_0>1$、感染が収束する条件は$R_0 < 1$となります。

　イギリスのボリス・ジョンソン首相が、3月16日の会見で、イギリス国民の約60%が集団免疫を獲得するまでピークをコントロールする趣旨の説明をしたようですが（ニューズウィークの記事）、これは、$R_0 = 2.6$と仮定して、ワクチンや自己免疫などでの集団免疫率$H$が、

(1-H)R_0 < 1

を満たす条件から$H > 1-1/R_0=0.615$と想定されたようです。つまり、集団免疫によって、$R_0$の条件が緩和されるのですね。しかし、ワクチンの開発には臨床試験など１年以上かかる見込みですし、敢えて過剰に感染を抑えないことで集団免疫の獲得をしようという大胆な戦略です。

基本再生産数の計算モデル

　ここで、基本再生産数を計算するために、上記の定義とSEIRモデルをベースとした、簡易モデルを考えたいと思います。次の図を見てください。

上段は日付、下段の棒はある日の感染者集団を表しています。各記号の意味は、

E:感染症が潜伏期間中の者(Exposed)
I:発症者(Infectious)
P:検査陽性者(Positive)
lp:潜伏期間
ip:感染期間

を表しています。つまり、ある日tにおけるIからEの生産は、t以前に感染し発症した状態にあるIの集団からの感染であり、その日tの再生産数は$R_0(t)/ip$である、というモデルです。
　これを式で表すと、

\frac{1}{ip} R_0(t) \times \sum_{s=t+1}^{t+ip} P(s) = P(t+lp+ip), \\
\therefore R_0(t) = \frac{ ip \times P(t+lp+ip)} {\sum_{s=t+1}^{t+ip} P(s)}

となります。以下、この式を使って、公開されている検査陽性者のデータから、基本再生産数$R_0(t)$の時間推移を計算してみましょう。

元になるデータ

　本記事では、都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社提供）において、公開されているcsvデータを使わせて頂きました。都道府県のデータが１つに集約されており、大変使いやすい形になっています。

Pythonで計算してみる

では、Pythonを使って、都道府県別の基本再生産数R0(t)を計算してみましょう。

前提条件

前提条件として、以下を用います。

lp:潜伏期間 = 5[day]
ip:感染期間 = 8[day]
上記の図のPは確定診断日とする。

なお、lpとipには諸説ありますが、一つの目安として設定しています。また、計算式の都合上、データの最新の日付から、lp+ip=13[day]以前の$R_0(t)$しか計算できません。グラフ上では、13日前～最新日の期間の値が0になっていても意味のない値ですので、ご注意ください。

コード

まず、ライブラリをインポートします。

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import datetime
import locale

都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社提供）で提供されているcsvから、各都道府県別のデータを読み込み、データフレームに変換する関数です。

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'確定日YYYYMMDD'])
    # 確定日, 受診都道府県のみ抽出
    df1 = df.loc[:,[u'確定日YYYYMMDD',u'受診都道府県']]
    df1.columns = ['date','pref']
    # 受診都道府県で抽出
    if pref is not None: # Noneであれば日本全体
        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')
    return df2

計算用のデータフレームを定義する関数です。列のPpreは計算式の分母用、Patは計算式の分子用です。2020年1月24日を開始日として、days分だけ計算の枠を作ります。

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)

計算式に従って、$R_0$を計算する関数です。計算を簡単にするため、インデックスの範囲外にアクセスしたらNaNを返すラムダ式を作りました。

def calcR0(df, keys):
    lp = keys['lp']
    ip = keys['ip']
    nrow = len(df)
    getP = lambda s: df.loc[s, 'P'] if s < nrow else np.NaN
    for t in range(nrow):
        df.loc[t, 'Ppre'] = sum([ getP(s) for s in range(t+1, t + ip + 1)])
        df.loc[t, 'Pat' ] = getP(t + lp + ip)
        if df.loc[t, 'Ppre'] > 0:
            df.loc[t, 'R0'  ] = ip * df.loc[t, 'Pat'] / df.loc[t, 'Ppre']
        else:
            df.loc[t, 'R0'  ] = np.NaN
    return df

結果をグラフ表示する関数です。

def showResult(df, title):
    # R0=1 : 収束のためのターゲット
    ptgt = pd.DataFrame([[df.iloc[0,0],1],[df.iloc[len(df)-1,0],1]])
    ptgt.columns = ['date','target']
    # show R0
    plt.rcParams["font.size"] = 12
    ax = df.plot(title=title,x='date',y='R0', figsize=(10,7))
    ptgt.plot(x='date',y='target',style='r--',ax=ax)
    ax.grid(True)
    ax.set_ylim(0,)
    plt.show()

各都道府県に対する計算を行う部分です。

def R0inJapanPref(pref, label):
    keys = {'lp':5, 'ip':8 }
    df1 = makeCalcFrame(60) # 60 days
    df2 = readCsvOfJapanPref(pref)
    df = mergeCalcFrame(df1, df2)
    df = calcR0(df, keys)
    showResult(df, 'COVID-19 R0 ({})'.format(label))
    return df

preflist = [[None, 'Japan'], [u'東京都', 'Tokyo'],\
            [u'大阪府', 'Osaka'],  [u'愛知県', 'Aichi'],\
            [u'北海道', 'Hokkaido']]
dflist = [[R0inJapanPref(pref, label), label] for pref, label in preflist]

上の計算結果をまとめて表示する部分です。

def showResult2(ax, df, label):
    # show R0
    plt.rcParams["font.size"] = 12
    df1 = df.rename(columns={'R0':label})
    df1.plot(x='date',y=label, ax=ax)

# R0=1
dfs = dflist[0][0]
ptgt = pd.DataFrame([[dfs.iloc[0,0],1],[dfs.iloc[len(dfs)-1,0],1]])
ptgt.columns = ['date','target']
ax = ptgt.plot(title='COVID-19 R0', x='date',y='target',style='r--', figsize=(10,8))
#
for df, label in dflist:
    showResult2(ax, df, label)
#
ax.grid(True)
ax.set_ylim(0,10)
plt.show()

計算結果

　それでは、いよいよ計算結果を見てみましょう。

日本全体の基本再生産数の推移

2月10日までは最大10付近までの大きな流入があるようです。これは、恐らく、春節（1月24日～1月30日）に来日した中国旅行者の影響ではないでしょうか。
その後は、一旦3付近まで上昇するものの、1を辛うじて超える付近（平均して1.5付近）で推移しているようです。

東京都の基本再生産数の推移

2月6日ごろまでは最大16の高い値でしたが、2月7日以降はおおむね落ち着いているようです。
ただし、2月18日～2月23日や3月1日など、時折3を超える山があり、クラスターの発生や、海外からの流入があった可能性を示唆しています。

大阪府の基本再生産数の推移

2月17日～2月24日に最大26の大きなピークが見られます。これは恐らくクラスターの発生でしょう。
実際、高槻市のページによると、大阪市内の４つのライブハウスで2月15日から2月24日にかけて行われたライブの参加者から、感染クラスターの発生が確認されています。

愛知県の基本再生産数の推移

2月12日までの、恐らく海外渡航者による感染が見られます。
その後、2月18日から2月26日にかけておおむね3付近で推移する山が見られます。これは、報道されているように、ジムや福祉施設で起きたクラスターと思われます。
クラスターの山と山が隣接してしまっているので、クラスター連鎖が心配されます。

北海道の基本再生産数の推移

2月8日から2月10日にかけて最大80もの鋭いピークが見られます。これは恐らく、さっぽろ雪祭り（2月4日～2月11日）で発生したクラスターではないかと思われます。
2月16日以降は、概ね1を下回る水準で推移しており、収束に向かっているものと思われます。

全体を比較

2月18日までに主に海外からの流入を主とした感染、3月1日までに国内でのクラスターの散発的な発生があったものと思われます。
3月1日以降は、全般的に落ち着いてきており、まだ1を超える水準であるものの、ギリギリ収束に向かっていると言えるのではないでしょうか。

考察

以上から、日本の各都道府県の感染者データに基づく、基本再生産数R0の推定に関して、次の傾向がシミュレーションから導けます。

2月中旬までの海外からの流入が最もR0が高く、その後国内でのクラスターと思われる影響で散発的に跳ね上がる傾向が見られる。
それ以外の期間では、特に北海道では1を下回る期間をキープ出来ており、概ね良好な状態。
計算結果のクラスター疑惑日と、報道されているクラスター発生日（大阪のライブハウス、さっぽろ雪まつり）がよく一致しており、計算モデルからクラスター発生を当てられる可能性が高い。（もっとも、積極的疫学調査による発見の効果もあると思われます。）

さらに言えば・・・

シミュレーション結果からは、基本再生産数$R_0$は一定ではなく、かなり変動しています。
基本再生産数$R_0$を上昇させる主要因は、海外からの流入と、クラスター発生であることが数値的にも裏付けられます。これらを抑制する対策は、やはり海外渡航制限と、クラスターの要因になりうるイベントや施設の自粛であると思われます。
全体的な傾向では、3月1日以降は収束傾向にあり、日本全体でR0が1に近づいているため、現状の対策が功を奏しているようにも見えます。
$R_0$は陽性者数と陽性者数の比率なので、検査されていない隠れ感染者が検査陽性者の一定倍いたとしても影響は軽微でしょう。

参考リンク

下記のページを参考にさせて頂きました。
都道府県別新型コロナウイルス感染者数マップ（ジャッグジャパン株式会社提供）
京大数理解析研究所講究録
 日本数理生物学会レター
 ニューズウィークの記事
 SEIRモデル
 高槻市のページ

直接引用してませんが、大変有用なリンクを追加しておきます。
COVID-19 reports
全国クラスターマップ

新型コロナウイルスの都道府県別の基本再生産数の推移を計算してみる