Go to Qiita Advent Calendar 2024 Top
LoginSignup
45
20

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@akeyhero(Akihiro Katsura)

新型コロナは本当に脅威か？ Stan で検証した (かった)

Last updated at Posted at 2020-03-15

【告知】 続編で COVID-19 の致死率を求めています。

目的

  • COVID-19 の致死率が本当に高いのか検証する
  • 確率的プログラミング言語 Stan を触ってみる

確率的プログラミング言語とは

統計学的な関係性 (モデル) を記述したらいい感じに解いてくれるプログラミング言語です。実装としては、マルコフ連鎖モンテカルロ法によるサンプリングを行うようです。

手続き型プログラミング言語とは役割も手法も全く異なっているため、それを置き換える存在ではありません。

Stan は確率的プログラミング言語の一つです。

準備

データ

ダイヤモンド・プリンセス号の感染状況を使います。
なぜなら、世界でこれよりも詳しく感染状況の調査が行われた空間は存在せず、最も正確な検証ができると期待されるからです。

項目 データ 日時 ソース
感染者数 696人 3/5 時事ドットコム 1
感染者数 (有症状) 301人 2/20 国立感染症研究所 2
感染者数 (無症状) 318人 2/20 国立感染症研究所 2
死亡者数 7人 3/7 読売新聞 3

年齢別でのインフルエンザとの比較も行いたいので、 CDC の統計 4 を加工して国立感染症研究所 2 のデータと比較できるようにします。階級の変換は単純にその幅の比率によって行いました。
CDC の統計は有症状のインフルエンザ感染数についての数値しか存在しないので、 COVID-19 に関しても有症状の感染数データを使います。

階級 CODIV-19 有症状感染数 インフルエンザ致死率 (推定)
0-9 0 0.0050%
10-19 2 0.0063%
20-29 25 0.0206%
30-39 27 0.0206%
40-49 19 0.0206%
50-59 28 0.0614%
60-69 76 0.4465%
70-79 95 0.8315%
80-89 27 0.8315%
90-99 2 0.8315%
Total 301 0.0962%

実行環境

Stan をインストールします。
Stan はコマンドライン単独でも動かせますが、ラッパーを使う方がいろいろと楽です。

今回は pip で簡単に導入できる、 Python インターフェイスの PyStan を使います。
numpy と cython 、さらにグラフの表示に scipy と matplotlib が必要なので、それも一緒に入れておきます。

Anaconda を使う場合、

$ conda create -n dp-corona-stan python=3.7 numpy cython scipy matplotlib pystan
$ conda activate dp-corona-stan

とりあえずジャブから

感染者数 (696人) と死亡者数 (7人) から致死率を Stan で推定します。

import pystan

model = pystan.StanModel(model_code='''
data {
    int N; // 感染者数
    int D; // 死亡者数
}
parameters {
    real<lower=0, upper=1> p;
}
model {
    D ~ binomial(N, p);
}
''')

data = {
    'N': 696,
    'D':   7
}

fit = model.sampling(data=data, chains=4)
print(fit)

StanModel に渡している文字列が Stan のコードです。

この記事では Stan の書き方に関して深入りしませんが、簡単に説明すると、 data に指定されたデータを入力し、 parameters に指定した変数を model を満たすように最適化するという意味のコードになっています。

感染者のうち何人死亡したか、という事象は、統計的には二項分布 (binomial distribution) に従うことになります。すなわち、 model の記述、

    D ~ binomial(N, p);

は、 感染者 N がそれぞれ p の確率で死亡するとき、 D 人死亡した、という状況をモデル化したものです。

このコードを実行すると、全感染者に対する致死率 p を推定できます。

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
p      0.01  9.8e-5 4.1e-3 5.0e-3 8.4e-3   0.01   0.01   0.02   1724    1.0
lp__ -44.22    0.02   0.73  -46.3  -44.4 -43.95 -43.76  -43.7   1531    1.0

corona_death_estimation_stan_simple.png

p は 1% くらいと推定されます。
世間では 2% と言われたりもしていますが、このデータから、 2% というのは 95% 信頼区間 ( 5.0e-3 〜 0.02 ) の当落線上であり、 その可能性は棄却はできないがそんなに高くないだろうと推定できます。

ただし、感染者の半分は無症状なので、有症状の感染者に対する致死率と基準を置き直せば確かに 2% くらいになります。 (ちなみに、インフルエンザは 1/3 が無症状と言われているらしいです 5。)

年齢を考慮する

しかしながら、ダイヤモンド・プリンセス号の乗員・乗客の半分以上は 60 歳以上でした。これは世間の人口分布とは違うので、実態に即した致死率を得るためには、その点を考慮する必要があります。

そこで、 COVID-19 の致死率がインフルエンザ相当であれば、ダイヤモンド・プリンセス号規模の感染があった場合、どの程度の死者が発生するのかを推定し、実際の死者数 (7 人) と比較します。

準備: データを Python で表現

データ で用意したデータを py ファイルに落としておきます。

S_xx が年代別の有症状の感染者数、 p_flu_xx が年代別の有症状インフルエンザ感染者数に対する致死率の推定値です。

# data.py
data = {
    'S_0x':  0,
    'S_1x':  2,
    'S_2x': 25,
    'S_3x': 27,
    'S_4x': 19,
    'S_5x': 28,
    'S_6x': 76,
    'S_7x': 95,
    'S_8x': 27,
    'S_9x':  2,
    'p_flu_0x': 0.000050,
    'p_flu_1x': 0.000063,
    'p_flu_2x': 0.000206,
    'p_flu_3x': 0.000206,
    'p_flu_4x': 0.000206,
    'p_flu_5x': 0.000614,
    'p_flu_6x': 0.004465,
    'p_flu_6x': 0.008315,
    'p_flu_7x': 0.008315,
    'p_flu_8x': 0.008315,
    'p_flu_9x': 0.000962
}

まず NumPy で検証する

NumPy には二項分布に従う乱数を発生させる機能があります。その機能を使えば簡単にシミュレーションが行なえます。

import numpy as np
from data import data

N = 10000
MIN_DEATH = 7

sample = (np.random.binomial(data['S_0x'], data['p_flu_0x'], N) +
          np.random.binomial(data['S_1x'], data['p_flu_1x'], N) +
          np.random.binomial(data['S_2x'], data['p_flu_2x'], N) +
          np.random.binomial(data['S_3x'], data['p_flu_3x'], N) +
          np.random.binomial(data['S_4x'], data['p_flu_4x'], N) +
          np.random.binomial(data['S_5x'], data['p_flu_5x'], N) +
          np.random.binomial(data['S_6x'], data['p_flu_6x'], N) +
          np.random.binomial(data['S_7x'], data['p_flu_7x'], N) +
          np.random.binomial(data['S_8x'], data['p_flu_8x'], N) +
          np.random.binomial(data['S_9x'], data['p_flu_9x'], N))
probability = float(sum(sample >= MIN_DEATH)) / N

print('Average # of deaths: %.2f' % (sample.mean()))
print('# of deaths >= %d: %.2f%%' % (MIN_DEATH, probability * 100))

こちらも深入りはしませんが、簡単に説明します。

np.random.binomial(n, p, N) は、確率 p で発生する事象を n 回試行するして発生回数を記録する作業を N 回繰り返します。例えば np.random.binomial(2, 0.5, 3) だったら array([2, 0, 1]) とかが返ってきます。

NumPy の配列は普通に足し算すると要素ごとの足し算になるので、これにより全階級に亘っての発生回数が得られます。

また、 NumPy の配列に対して比較演算子を使うと、要素ごとに真偽判定してできた配列を返します。 例えばnp.array([2, 0, 1]) > 1 をすると array([ True, False, False]) が返ってきます。

Python の sum 関数は True を 1、 False を 0 と見なすので、結局、 float(sum(sample >= MIN_DEATH)) / N により、試行 N 回中、死亡者数が MIN_DEATH を超えた割合を計算できます。

結果は

Average # of deaths: 1.67
# of deaths >= 7: 0.21%

のようになりました。

COVID-19 の死亡率がインフルエンザ並だったとしたら、 301 人の有症状感染者に対して 7 人も死亡するというのは相当にありえないと言えそうです。

※ 有症状感染者のデータは 2/20 時点のもので、 7 人目の死亡者が確認された時点 (3/7) では更に感染数が増えていましたが、それを考慮しても COVID-19 並に死亡する確率が 1% を超えることにはならなさそうです。

Stan で検証してみる

二項分布で解いてみる

最適化の対象が確率から死亡者数に変わるのと、年齢の階級を考慮する点が違うだけなので、さっきの Stan コードと同じようなものを書けばいけるはずです。

年齢別の推定死亡数を d_xx とすると、

data {
    int S_0x;
    int S_1x;
    int S_2x;
    int S_3x;
    int S_4x;
    int S_5x;
    int S_6x;
    int S_7x;
    int S_8x;
    int S_9x;
    real p_flu_0x;
    real p_flu_1x;
    real p_flu_2x;
    real p_flu_3x;
    real p_flu_4x;
    real p_flu_5x;
    real p_flu_6x;
    real p_flu_7x;
    real p_flu_8x;
    real p_flu_9x;
}
parameters {
    // upper = S_xx + 1 so that S_xx can be 0
    int<lower=0, upper=S_0x+1> d_0x;
    int<lower=0, upper=S_1x+1> d_1x;
    int<lower=0, upper=S_2x+1> d_2x;
    int<lower=0, upper=S_3x+1> d_3x;
    int<lower=0, upper=S_4x+1> d_4x;
    int<lower=0, upper=S_5x+1> d_5x;
    int<lower=0, upper=S_6x+1> d_6x;
    int<lower=0, upper=S_7x+1> d_7x;
    int<lower=0, upper=S_8x+1> d_8x;
    int<lower=0, upper=S_9x+1> d_9x;
}
transformed parameters {
    int d;
    d = d_0x + d_1x + d_2x + d_3x + d_4x + d_5x + d_6x + d_7x + d_8x + d_9x;
}
model {
    d_0x ~ binomial(S_0x, p_flu_0x);
    d_1x ~ binomial(S_1x, p_flu_1x);
    d_2x ~ binomial(S_2x, p_flu_2x);
    d_3x ~ binomial(S_3x, p_flu_3x);
    d_4x ~ binomial(S_4x, p_flu_4x);
    d_5x ~ binomial(S_5x, p_flu_5x);
    d_6x ~ binomial(S_6x, p_flu_6x);
    d_7x ~ binomial(S_7x, p_flu_7x);
    d_8x ~ binomial(S_8x, p_flu_8x);
    d_9x ~ binomial(S_9x, p_flu_9x);
}

※ Stan は配列を扱えるのでもっとスッキリ書けますが、階級を配列の添字に落とし込むのが何となく嫌だったので愚直に列挙してます。

実行すると

ValueError: Failed to parse Stan model 'anon_model_fecb1e77228fe372ef4eb9bc4bcc8086'. Error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Parameters or transformed parameters cannot be integer or integer array;  found int variable declaration, name=d_0x
 error in 'unknown file name' at line 26, column 35
  -------------------------------------------------
    24: parameters {
    25:     // upper = S_xx + 1 so that S_xx can be 0
    26:     int<lower=0, upper=S_0x+1> d_0x;
                                          ^
    27:     int<lower=0, upper=S_1x+1> d_1x;
  -------------------------------------------------

integer ではダメって怒られました。。。

StanとRでベイズ統計モデリング読書会Ch.9 6 の 20 ページ目には、 Stan の弱点として parameters に int が使えないことが挙げられていました。
どうも本当に使えないみたいです。

Ω＼ζ゜) チーン♪

binomial は左辺に int を求めるため、 d_xx たちを int ではなく real で宣言しておく裏技も使えません。

Ω＼ζ゜) チーン♪

ベータ分布で無理くり解いてみる

二項分布がダメだったので、ベータ分布を試しに使ってみました。

マサカリが飛んでくる覚悟で言えば、ベータ分布は二項分布の逆みたいなやつで、二項分布に従う事象に於ける発生確率 p の確率分布です。くわしくはこちら。

    d_xx ~ binomial(S_xx, p_flu_xx);

と書いていたものを

    p_flu_xx ~ beta(d_xx + 1, S_xx - d_xx + 1);

に書き直し、 d_xx の型を real にしてみます。

import pystan
from data import data

model = pystan.StanModel(model_code="""
data {
    int S_0x;
    int S_1x;
    int S_2x;
    int S_3x;
    int S_4x;
    int S_5x;
    int S_6x;
    int S_7x;
    int S_8x;
    int S_9x;
    real p_flu_0x;
    real p_flu_1x;
    real p_flu_2x;
    real p_flu_3x;
    real p_flu_4x;
    real p_flu_5x;
    real p_flu_6x;
    real p_flu_7x;
    real p_flu_8x;
    real p_flu_9x;
}
parameters {
    // upper = S_xx + 1 so that S_xx can be 0
    real<lower=0, upper=S_0x+1> d_0x;
    real<lower=0, upper=S_1x+1> d_1x;
    real<lower=0, upper=S_2x+1> d_2x;
    real<lower=0, upper=S_3x+1> d_3x;
    real<lower=0, upper=S_4x+1> d_4x;
    real<lower=0, upper=S_5x+1> d_5x;
    real<lower=0, upper=S_6x+1> d_6x;
    real<lower=0, upper=S_7x+1> d_7x;
    real<lower=0, upper=S_8x+1> d_8x;
    real<lower=0, upper=S_9x+1> d_9x;
}
transformed parameters {
    real d;
    d = d_0x + d_1x + d_2x + d_3x + d_4x + d_5x + d_6x + d_7x + d_8x + d_9x;
}
model {
    p_flu_0x ~ beta(d_0x + 1, S_0x - d_0x + 1);
    p_flu_1x ~ beta(d_1x + 1, S_1x - d_1x + 1);
    p_flu_2x ~ beta(d_2x + 1, S_2x - d_2x + 1);
    p_flu_3x ~ beta(d_3x + 1, S_3x - d_3x + 1);
    p_flu_4x ~ beta(d_4x + 1, S_4x - d_4x + 1);
    p_flu_5x ~ beta(d_5x + 1, S_5x - d_5x + 1);
    p_flu_6x ~ beta(d_6x + 1, S_6x - d_6x + 1);
    p_flu_7x ~ beta(d_7x + 1, S_7x - d_7x + 1);
    p_flu_8x ~ beta(d_8x + 1, S_8x - d_8x + 1);
    p_flu_9x ~ beta(d_9x + 1, S_9x - d_9x + 1);
}
""")

fit = model.sampling(data=data, chains=4)
print(fit)

実行して得られた結果がこちら:

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
d_0x    0.1  1.4e-3   0.09 2.3e-3   0.03   0.07   0.14   0.35   4361    1.0
d_1x   0.12  1.7e-3   0.12 2.5e-3   0.03   0.08   0.16   0.43   4474    1.0
d_2x   0.19  2.6e-3   0.18 5.0e-3   0.06   0.14   0.27   0.67   4976    1.0
d_3x    0.2  3.1e-3   0.19 4.4e-3   0.06   0.14   0.27   0.71   3641    1.0
d_4x   0.19  2.5e-3   0.19 3.5e-3   0.05   0.13   0.26   0.69   5315    1.0
d_5x   0.25  3.2e-3   0.23 8.0e-3   0.08   0.18   0.34   0.85   5100    1.0
d_6x   0.93    0.01   0.73   0.04   0.36   0.77   1.33    2.8   4387    1.0
d_7x   1.06    0.01    0.8   0.04   0.43    0.9   1.51   3.01   4028    1.0
d_8x   0.54  7.0e-3   0.48   0.01   0.18   0.42   0.76   1.82   4729    1.0
d_9x   0.17  2.4e-3   0.16 4.3e-3   0.05   0.12   0.24   0.58   4580    1.0
d      3.73    0.02   1.25   1.66   2.84   3.59   4.51   6.57   4367    1.0
lp__  -1.64    0.08   2.58  -7.78  -3.13  -1.28   0.24   2.37   1048    1.0

d が死亡者数ですが、 50% が 3.59 になっており、 NumPy により得られた 1.67 より明らかに多いです。

Ω＼ζ゜) チーン♪

なぜでしょう。

これはおそらく d が real として定義されているために起こっている問題です。本来は整数値しか取れない d_xx が小数を取れることで、現実の問題から乖離してしまっているようです。

実際、 d_xx を floor で整数化して計算してみると

       mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
...
d      1.82    0.03   1.29   0.01   1.01   1.69    2.5   4.93   1710    1.0

のようにそれっぽい値になりました。

※ 単に floor を取るだけだと勾配がなくなってうまく収束しないので、 inv_logit とかを使って、うまくスムージングしてやる必要があります。

しかしながら、ここまでやるなら NumPy を使えばいいですね。確率モデリングとしてもはや正しいとは言い難いですし。

この問題を Stan で解くなにかうまい方法はないものでしょうか。

まとめ

部分的ながら、 Stan を使って新型コロナウイルスの危険性を推定することができました。

ダイヤモンド・プリンセスの状況から、 COVID-19 はどうやら従来のインフルエンザと比べても、致死率の点で凶悪なものであることは相当に間違いなさそうです。加えて、インフルエンザは既に抗体を持っている人がたくさんいる一方で、 COVID-19 の抗体を持っている人はほとんどいない点も脅威です。個人のためというより社会全体のために、十分な警戒をしたほうがいいでしょう。

しかしながら、インフルエンザでも今回の規模で感染すると 1.67 人程度の死者が出るという推定から、 COVID-19 をここまで恐れるのであれば、普段のインフルエンザももっと警戒してはいいのではないかということを、個人的見解ながら、申し上げておきたいと思います。

続編では COVID-19 の致死率を求める課題にもチャレンジしています。ぜひご覧ください。

付録

ソースコード: https://github.com/akeyhero/dp-corona-stan

  1. https://www.jiji.com/jc/article?k=2020030501355&g=soc ↩

  2. https://www.niid.go.jp/niid/ja/diseases/ka/corona-virus/2019-ncov/2484-idsc/9422-covid-dp-2.html ↩ ↩2 ↩3

  3. https://www.yomiuri.co.jp/national/20200307-OYT1T50228/ ↩

  4. https://www.cdc.gov/flu/about/burden/2018-2019.html ↩

  5. https://en.wikipedia.org/wiki/Influenza#Signs_and_symptoms ↩

  6. https://www.slideshare.net/KojiKosugi/stanrch9 ↩

45
20
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
45
20

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?