はじめに

アメリカ大統領選挙の選挙予測なんかで有名なアメリカの統計学者ネイト・シルバーが、FiveThirtyEightと言うサイトに、"Coronavirus Case Counts Are Meaningless" (コロナウィルスの発生件数には意味がない) と言う記事を書いています。

この記事の主張は、 「どのように検査が行われているかを十分知らない限りは、報告されているCOVID-19の感染者数は、役に立つ指標にはならない¹。」 と言うものです。直感的に「そりゃそうだ」と思う一方で、ウィルスの感染は指数関数的な現象だから、検査数が(もしかすると恣意的に)絞られていたとしても、再生産率(reproduction ratio; 後で説明)は観測できるはずだ、とか言う主張も見聞きします。また全く別の観点から、検査を絞らなければ医療崩壊すると言うような主張も見受けられます。ネイトの記事には、読者が試してみることのできるシンプルなシミュレーションがエクセルスプレッドシートとして添付されていて、検査ポリシー(検査対象の選び方や検査の総数)と感染対策(social distancingやlockdown)によって、感染者数の報告と真の感染者数がどれほど異なるものになるのかを体感できるようになっています。真の再生産レートと観測される見かけの再生産レートの違いも観察できます。

そのスプレッドシートのモデル²は簡単にパラメータを弄って、結果の違いをみるには便利なのですが、スプレッドシートの一つ一つのセルが $\text{BH32}$ と言った「カラム=アルファベット」と「行=数字」の組み合わせで表されていて、それらの関係がBASIC likeな式で書かれているため、何をどうやって計算しているのか理解するには、あまり向いていません。

そこで、そのスプレッドシートのモデルを Python プログラムに(手作業で)変換して、何を計算しているのかを理解してみました。さらに、現在 (2020年4月22日) の日本の状況をパラメータで与えて試してみました。そしたら、まだ報告感染者数は1万人を超えたばかりだと言うのに、実際の患者数は実に 7５0万人を超え、全人口の６%くらいが感染していることになっちゃいました。これだと、「外に出たら、隣の人は感染者だ」と思っといた方がいいですね。

プログラム

githubに jupyter notebook があります。

感染現象を連続的に起こる現象として捉えると、微分方程式になって素人には手が出しにくいので、とても大雑把な離散的な現象と捉えます。つまり、感染の連鎖を前の世代³から次の世代への連鎖と考え、さらに単純化して世代間の時間差をある定数(例えば5日間)と固定して計算します。

一番大事なパラメータは、再生産率 (reproduction rate; $R$) で、第$n$世代の感染者1名が何人の第$n+1$世代に病気をうつすかを指定するパラメータです。これは病気の感染力によっても変わりますが、社会の人々の行動によっても変わります。クルーズ船の再生産率はとても高かったでしょうし、ロックダウンして外出が全て禁止されてからの武漢の再生産率は低かったでしょう。$R$が1よりも大きいと感染が広がり、1より小さいといずれ病気は消えていきます。

パラメータ設定

病気の拡散に関するパラメータを Disease_spread_params class で表現します。

R_uncontrolled: 何の対策も取らなかったときの再生産率。
R_intermediate: 社会的距離、手洗いの励行、マスク着用など、人々が注意深くなった中間段階の再生産率。
R_lockdown: ロックダウンが施行されたときの再生産率。
cluster: 人口に占める感染者の割合が多くなったときに調整するパラメータ(0=no adjust, 0.5: slightly, 1: moderately, 2: significantly)
begin_internmediate: 無対策の期間が終わって、中間段階に入る世代番号
begin_lockdown: ロックダウンが始まる世代番号
mild: 感染して中等症になる割合
asymptomatic: 感染しても無症状の割合
(感染して重症になる割合は 1-mild-asymptomatic)

人口に関するパラメータを Population_params class で表現します。

population: この国の総人口
initial_infection: 最初(第0世代)の感染者の数
faux_severe: コロナウィルス以外の原因で、COVID-1９に似た重症患者の人口に対する割合
faux_mild: コロナウィルス以外の原因で、COVID-19に似た中等症患者の人口に対する割合
desire_severe: 重症患者の中で、PCR検査を求める人の割合
desire_mild: 中等症患者の中で、PCR検査を求める人の割合
desire_asymptomatic: 無症状の人の中で、PCR検査を求める人の割合

PCR検査に関するパラメータを Testing_params class で表現します。

initial_tests: 第１世代に対して提供されるPCR検査のキャパシティ
ramp_period: 検査数を増やし始める世代番号
test_growth_rate: 検査数を増やし始めてからの増加率
test_max: 検査の最大数
rationed_test_ratio: 重症 > 中等症 > 無症状の優先順位を守るために確保される検査の割合
false_negative: PCR検査の偽陰性の割合
false_positive: PCR検査の偽陽性の割合
delay: 検査してから結果が出るまでの遅れ (世代数)

これらのパラメータの値は各種論文等からもっともらしい物を選んで、デフォルト値として与えています。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time, re
import datetime

iround = lambda x: int((x * 2 + 1) // 2)

class Disease_spread_params:
    def __init__(self,
                 R_uncontrolled = 2.7,
                 R_intermediate = 1.4,
                 R_lockdown = 0.7,
                 cluster = 0.5,
                 begin_intermediate = 11,
                 begin_lockdown = 15,
                 mild = 0.6,
                 asymptomatic = 0.3,
                 zero_date = '2020-01-01',
                 serial = 5):
        self.R_uncontrolled = R_uncontrolled
        self.R_intermediate = R_intermediate
        self.R_lockdown = R_lockdown
        self.cluster = cluster
        self.begin_intermediate = begin_intermediate
        self.begin_lockdown = begin_lockdown
        self.severe = 1 - mild - asymptomatic
        self.mild = mild
        self.asymptomatic = asymptomatic
        self.zero_date = zero_date
        self.serial = serial

class Population_params:
    def __init__(self,
                 population = 10_000_000,
                 initial_infection = 1,
                 faux_severe = 0.001,
                 faux_mild = 0.025,
                 desire_severe = 1.0,
                 desire_mild = 0.5,
                 desire_asymptomatic = 0.02):
        self.population = population
        self.initial_infection = initial_infection
        self.faux_severe = faux_severe
        self.faux_mild = faux_mild
        self.desire_severe = desire_severe
        self.desire_mild = desire_mild
        self.desire_asymptomatic = desire_asymptomatic

class Testing_params:
    def __init__(self,
                 initial_tests = 1_000,
                 ramp_period = 3,
                 test_growth_rate = 0.5,
                 test_max = 10_000_000,
                 rationed_test_ratio = 0.75,
                 false_negative = 0.20,
                 false_positive = 0.002,
                 delay = 2):
        self.initial_tests = initial_tests
        self.ramp_period = ramp_period
        self.test_growth_rate = test_growth_rate
        self.test_max = test_max
        self.rationed_test_ratio = rationed_test_ratio
        self.false_negative = false_negative
        self.false_positive = false_positive
        self.delay = delay

ここで、実数を整数に丸める関数として、Python組み込みの round() ではなく自前の iround() と言う関数を使っています。Excel の ROUND 関数との互換性を取るためです⁴。

シミュレータ

参考にしたスプレッドシートのセルごと依存関係を見ながら、どのように計算しているかを理解して、実装し直しました。正直言って（謙遜ではなく）汚い。

Simulator class は、パラメータを受け取って instantiate し、run() method を呼ぶと実行され、結果は pandas.DataFrame として返ります。

class Simulator:
    def __init__(self,
                 disease_spread_params,
                 population_params,
                 testing_params):
        self.disease_spread_params = disease_spread_params
        self.population_params = population_params
        self.testing_params = testing_params
        
        self.columns = [ 'date', 'actual.R', 'doubling time in days' ]
        self.date_regex = re.compile('(\d+)-(\d+)-(\d+)')

    def decode_date(self, date):
            match = self.date_regex.fullmatch(date)
            if match:
                y, m, d = match.group(1, 2, 3)
                timestamp = time.mktime((int(y), int(m), int(d), 0, 0, 0, -1, -1, -1))
                return timestamp
            return None
    
    def encode_date(self, timestamp):
        t = time.localtime(timestamp)
        return '{0:04d}-{1:02d}-{2:02d}'.format(t[0], t[1], t[2])

    def run(self, zero_day = '2020-01-01', generations = 40):
        self.generations = generations
        self.df = pd.DataFrame(index = range(0, self.generations))
        t0 = self.decode_date(zero_day)
        self.df['date'] = [ self.encode_date(t0 + d * self.disease_spread_params.serial * 24 * 60 * 60) for d in range(0, self.generations) ]

        self.set_target_R()
        self.compute_actual_infection()
        self.compute_tests()
        
        return self.df

    def set_target_R(self):
        begin_lockdown = self.disease_spread_params.begin_lockdown
        begin_intermediate = self.disease_spread_params.begin_intermediate
        self.df['target_R'] = np.NaN
        
        for i in range(0, self.generations):
            if begin_lockdown != None and i >= begin_lockdown:
                self.df.at[i, 'target_R'] = self.disease_spread_params.R_lockdown
            elif  begin_intermediate != None and i >= begin_intermediate:
                self.df.at[i, 'target_R'] = self.disease_spread_params.R_intermediate
            else:
                self.df.at[i, 'target_R'] = self.disease_spread_params.R_uncontrolled

    def compute_actual_infection(self):
        population = self.population_params.population
        initial_infection = self.population_params.initial_infection
        cluster = self.disease_spread_params.cluster
        serial = self.disease_spread_params.serial

        df = self.df
        df['susceptible'] = np.NaN
        df.at[0, 'susceptible'] = population - initial_infection
        df['new_infection'] = np.NaN
        df.at[0, 'new_infection'] = initial_infection
        df['cumulative_infection'] = np.NaN
        df.at[0, 'cumulative_infection'] = initial_infection
        df['actual_R'] = np.NaN
        df['actual_doubling_time'] = np.NaN
        
        for i in range(1, self.generations):
            df.at[i, 'new_infection'] = iround(df.at[i-1, 'susceptible']*(1-(1-((df.at[i-1, 'target_R']*(df.at[i-1, 'susceptible']/population)**cluster)/population))**df.at[i-1, 'new_infection']))
            df.at[i, 'cumulative_infection'] = df.at[i-1, 'cumulative_infection'] + df.at[i, 'new_infection']
            df.at[i, 'susceptible'] = population - df.at[i, 'cumulative_infection']
            df.at[i-1, 'actual_R'] = df.at[i, 'new_infection'] / df.at[i-1, 'new_infection']
            if df.at[i-1, 'cumulative_infection'] != 0:
                df.at[i-1, 'actual_doubling_time'] = np.inf if df.at[i, 'new_infection'] == 0 else serial * np.log(2) / np.log(df.at[i, 'cumulative_infection']/df.at[i-1, 'cumulative_infection'])

    def compute_tests(self):
        population = self.population_params.population
        ramp_period = self.testing_params.ramp_period
        tests_max = self.testing_params.test_max
        test_growth_rate = self.testing_params.test_growth_rate
        rationed_test_ratio = self.testing_params.rationed_test_ratio
        mild = self.disease_spread_params.mild
        asymptomatic = self.disease_spread_params.asymptomatic
        faux_severe = self.population_params.faux_severe
        faux_mild = self.population_params.faux_mild
        desire_severe = self.population_params.desire_severe
        desire_mild = self.population_params.desire_mild
        desire_asymptomatic = self.population_params.desire_asymptomatic
        false_negative = self.testing_params.false_negative
        false_positive = self.testing_params.false_positive
        delay = self.testing_params.delay
        serial = self.disease_spread_params.serial        

        cumulative_tests_conducted = 0
        cumulative_detected_cases = 0
        
        df = self.df
        df['tests_available'] = 0
        df['new_detected_cases'] = 0
        df['cumulative_detected_cases'] = 0
        
        for i in range(0, self.generations):
            if i  == 0:
                df.at[i, 'tests_available'] = 0
            elif i == 1:
                df.at[i, 'tests_available'] = self.testing_params.initial_tests
            elif i < ramp_period:
                df.at[i, 'tests_available'] = df.at[i-1, 'tests_available']
            else:
                df.at[i, 'tests_available'] = iround(min(tests_max, df.at[i-1, 'tests_available'] * (1 + test_growth_rate)))
            tests_available = df.at[i, 'tests_available']
            rationed_tests =  iround(tests_available * rationed_test_ratio)
            on_demand_tests = tests_available - rationed_tests

            new_infection_severe = iround(df.at[i, 'new_infection'] * (1 - mild - asymptomatic))
            new_infection_mild = iround(df.at[i, 'new_infection'] * mild)
            new_infection_asymptomatic = df.at[i, 'new_infection'] - new_infection_severe - new_infection_mild
            
            population_severe = iround((population - df.at[i, 'new_infection']) * faux_severe) + new_infection_severe
            population_mild = iround((population - df.at[i, 'new_infection']) * faux_mild) + new_infection_mild
            population_asymptomatic = population - population_severe - population_mild
            
            desiring_tests_severe = iround(population_severe * desire_severe * (1 - cumulative_tests_conducted/population))
            desiring_tests_mild = iround(population_mild * desire_mild * (1 - cumulative_tests_conducted/population))
            desiring_tests_asymptomatic = iround(population_asymptomatic * desire_asymptomatic * (1 - cumulative_tests_conducted/population))
            
            alloc_rationed_tests_severe = min(rationed_tests, desiring_tests_severe)
            alloc_rationed_tests_mild = min(desiring_tests_mild, rationed_tests-alloc_rationed_tests_severe)
            alloc_rationed_tests_asymptomatic = min(desiring_tests_asymptomatic, rationed_tests-alloc_rationed_tests_severe-alloc_rationed_tests_mild)

            unfilled_test_demand_severe = desiring_tests_severe - alloc_rationed_tests_severe
            unfilled_test_demand_mild = desiring_tests_mild - alloc_rationed_tests_mild
            unfilled_test_demand_asymptomatic = desiring_tests_asymptomatic - alloc_rationed_tests_asymptomatic
            unfilled_test_demand = unfilled_test_demand_severe + unfilled_test_demand_mild + unfilled_test_demand_asymptomatic
            
            alloc_on_demand_tests_severe = 0 if unfilled_test_demand == 0 else iround(on_demand_tests * unfilled_test_demand_severe / unfilled_test_demand)
            alloc_on_demand_tests_mild = 0 if unfilled_test_demand == 0 else iround(on_demand_tests * unfilled_test_demand_mild / unfilled_test_demand)
            alloc_on_demand_tests_asymptomatic = 0 if unfilled_test_demand == 0 else iround(on_demand_tests * unfilled_test_demand_asymptomatic / unfilled_test_demand)
            
            tests_conducted_severe = alloc_rationed_tests_severe + alloc_on_demand_tests_severe
            tests_conducted_mild = alloc_rationed_tests_mild + alloc_on_demand_tests_mild
            tests_conducted_asymptomatic = alloc_rationed_tests_asymptomatic + alloc_on_demand_tests_asymptomatic
            df.at[i, 'tests_conducted_severe'] = tests_conducted_severe
            df.at[i, 'tests_conducted_mild'] = tests_conducted_mild
            df.at[i, 'tests_conducted_asymptomatic'] = tests_conducted_asymptomatic
            tests_conducted = tests_conducted_severe + tests_conducted_mild + tests_conducted_asymptomatic
            
            cumulative_tests_conducted += tests_conducted
            
            positive_tests_severe = iround(tests_conducted_severe * new_infection_severe / population_severe * (1 - false_negative)) + \
              iround(tests_conducted_severe * (1 - new_infection_severe / population_severe) * false_positive)
            positive_tests_mild = iround(tests_conducted_mild * new_infection_mild / population_mild * (1 - false_negative)) + \
              iround(tests_conducted_mild * (1 - new_infection_mild / population_mild) * false_positive)
            positive_tests_asymptomatic = iround(tests_conducted_asymptomatic * new_infection_asymptomatic / population_asymptomatic * (1 - false_negative)) + \
              iround(tests_conducted_asymptomatic * (1 - new_infection_asymptomatic / population_asymptomatic) * false_positive)
            if i+delay < self.generations:
                df.at[i+delay, 'new_detected_cases'] = positive_tests_severe + positive_tests_mild + positive_tests_asymptomatic

            cumulative_detected_cases += df.at[i, 'new_detected_cases']
            df.at[i, 'cumulative_detected_cases'] = cumulative_detected_cases

            if i > 0 and df.at[i-1, 'new_detected_cases'] > 0:
                df.at[i-1, 'observed_R'] = df.at[i, 'new_detected_cases'] / df.at[i-1, 'new_detected_cases']
                df.at[i-1, 'observed_doubling_time'] = np.inf if df.at[i, 'new_detected_cases'] == 0 else serial * np.log(2) / np.log(df.at[i, 'cumulative_detected_cases']/df.at[i-1, 'cumulative_detected_cases'])

シミュレーション

Covidia シナリオ1

ネイト・シルバーの記事には仮想的な国 Covidia のシミュレーションが乗っています。その結果を可視化したものが以下のグラフです。Covidia国は、人口1千万人で、2020年1月1日に最初の感染者１名が入国することで、感染拡大が始まります。当初の再生産率は$2.7$、中間段階が $1.4$、ロックダウンしてからが$0.7$と仮定しています。

青い線が真の新規感染者数、黄色い線が検査で報告される新規感染者数の推移です。このシナリオでは当初の検査キャパシティは1世代あたり1,000件と少なめなものの、無制限に検査を増やすことで、時間差はあるものの真の感染者数に追従して報告数が伸びているように見えます。

しかしよく見ると、感染者数が増えているフェーズでは２０倍以上の過小評価になり、減少傾向になると逆に過大評価してしまうことがわかります。

日本のシナリオ?

日本はPCR検査数が低く抑えられていると言う批判を受けています。厚生労働省が毎日報告しています⁵が、いまだに検査数が需要に追いついていないと言われています。そこで、日本の人口1億2千6百万人に対して、PCR検査数を低く設定したらどうなるかを試してみました。ついでに、これまでに現実に報告されている陽性確定者数も合わせてプロットしてみました。上のグラフが新規感染者数だったのに対して、こちらは累積感染者数です⁶。

![fig2.png]
(https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223481/193b7bce-f37f-4b9e-727b-223c95b838f5.png)

ここで、緑の線は実際の日本の感染者数⁶です。

このシミュレーションで何かを予測するつもりはないのですが、まず真の値(青線)とシミュレーション結果及び報告されている値(緑線)の違いの大きさをみていただきたいと思います。以下のようなことが観察できます。

上下方向に今日2020/4/20頃の値を読むと、報告されている感染者数 (1万人くらい) に対して、モデルが計算した真の感染者数は700万人を超えています。実に700倍の違いがあります。
報告数(黄色、緑)も真の感染数もログスケールのグラフで直線的に増えているので指数的に増えていますが、傾きが違います。傾きが違うと言うことは、再生産率 $R$ が違うように見えています。つまり、「PCR検査数が絞られていても、再生産数は観測できる」と言う主張は間違っていると言うことです。

では、もし早期から検査のキャパシティーを増強して、無制限に検査を行なっていたとしたらどうなったでしょう。検査キャパシティーの初期値を10,000とし、徐々にしかし無制限に増やした場合の結果が次のグラフです。

この場合でもやはり真の感染者数を捉えることはできませんが、現在の報告者数の10倍程度は観測することができました。
また、報告者数(黄色)の立ち上がりは少し鈍いものの、その傾きは真の感染者数の傾きに近いので、両者の違いは左右方向の差（時間の遅れ）と捉えることができ、正しい再生産率 $R$ が観測できているように思われます。

まとめ

感染者が分かったところで病院や隔離施設がないので、治療という観点では違いがないという指摘は、その通りだと思います。でも、真の感染者数が報告されている10倍、場合によっては100倍も多いのだと知ったら、人々の日々の行動が少し変わるかもしれないとも思います。

感染学でテレビに出てくる専門家の方は、このモデルのずっと緻密なバージョンを作って未来を予測してるのかな。

知らんけど。

The number of reported COVID-19 cases is not a very useful indicator of anything unless you also know something about how tests are being conducted. ↩
ネイト・シルバーは、スプレッドシートをモデルとは読んでいませんでした。何かを予測することを目的とはしていないからです。 ↩
もちろんこの世代は、団塊の世代とかミレニアル世代とかの世代ではなく、最初に感染した人々を第0世代、そこから感染した人々を第1世代、...と呼びます。 ↩
Pythonのround()は偶数への丸めであるのに対して、EXCELのround()は四捨五入。 ↩
例えば、最新の報告では、前日の報告から7826件新たにPCR検査が行われたことがわかります。前の図と違って、累積感染者数のプロットです。 ↩
実データが累積で与えられていて、そちらの方が見慣れていると思ったので。 ↩ ↩²