前回投稿した「新型コロナは本当に脅威か？ Stan で検証した (かった)」は、次の二点で消化不良でした:

結局 COVID-19 の致死率がどんなもんなのか分からない
最終的に Stan がうまく使えていない

という訳で、今回はそのリベンジとして、 Stan を使って「真の致死率」を求めるという試みをやってみます。

前回の記事で説明したことに関しては繰り返し説明していないので、そちらもご参照ください

免責

筆者はベイズモデリングや Stan の専門家ではないため、この記事の内容には誤りや、よりよい手法が存在する可能性があります。

また、間違っている箇所を見つけたら、ご教示いただけると幸いです。

仮定

前回の投稿では、ダイヤモンド・プリンセス号に於ける感染調査 ¹ が、世界で最も正確であるという前提のものに議論を展開しました。今回もその方針を踏襲します。

しかしながら、ダイヤモンド・プリンセス号だけでは症例数が少なくて、年齢ごとの致死率を知ることができません。そこで、中国疾病預防控制中心のデータ ² を利用します。

この2つのデータは母集団がぜんぜん違うので、えいやで仮定を設定してなんとか統合します。ここでは、中国のデータの感染数は、真の感染者全体から、年齢層や致死率に関わらず一定の確率 $c$ で補足され得られたものであると仮定を用いました。

年齢層により症状の傾向が変わり、結果、捕捉された感染者の割合は異なると考えられるため、この仮定は明らかに不正確です。したがって、この記事の結論には、この仮定に起因する誤差が発生します。しかし、国立感染症研究所のデータ ¹ によれば、症状がある場合 (Symptomatic confirmed cases) に対する症状がない場合 (Asymptomatic confirmed cases) の割合はむしろ若年層の方が低くなっています。また、若年層の重症率が意外と高いという報告 ³ も存在します。このことから、直感よりは荒唐無稽な仮定ではないのではないかと思います。

ダイヤモンド・プリンセス号に於ける致死率は、ここで仮定した捕捉率 $c$ を用いて計算した年齢層ごとの致死率に従うものとすることで、この2つのデータを結合します。

一方、死亡者の捕捉については、中国の統計に於いてももれなく行われていると仮定します。

次に、ダイヤモンド・プリンセス号に於ける死亡者は 6 人としました。3月7日に 7 人目の死亡者が確認されていますが、国立感染症研究所のデータ ¹ から2週間以上経過していることから、その時点で感染者でなかった可能性も高いことと、 COVID-19 に関連する死亡かどうかは公表されていないこと、3月7日の時点を基準にするならば感染者の方もデータの時点より増えていたことが根拠です。

最後に、人種による差や、感染者分布の違い (本当は、中国の感染者分布と世界の感染者分布は違うはず)、医療環境の違い (特に、急激に患者が増えると医療崩壊が発生し致死率が高まるはず) を考慮しないことにしました。ただし、中国とダイヤモンド・プリンセス号に於ける感染者数分布の違いは例外的に考慮しています。

仮定内容をまとめると、以下のようになります。

中国に於ける感染者の捕捉率 $c$ は、年齢層に関わらず一定である
中国に於ける死亡者の捕捉率は 100% とする
ダイヤモンド・プリンセス号の統計は正確である (= 捕捉率 100%)
ダイヤモンド・プリンセス号での致死率は、年齢層ごとに、中国での致死率に捕捉率 $c$ を乗じて計算したものに一致する
ダイヤモンド・プリンセス号に於ける死亡者を 6 人とする
致死率の計算に際して、人種による差・感染者分布による差・医療環境の差を考慮しなくてよい

モデル

以下、 $N$ を感染者数、 $D$ を死亡者数、 $p$ を致死率とします。特に断りがない場合、上付き文字で年齢層、下付き文字で場所 ( ${\rm dp}$ = ダイヤモンド・プリンセス号、 ${\rm ch}$ = 中国 ) を示します。上付き文字がないことで全年齢層を表現します。

例えば、 $N_{{\rm dp}}^{60-69}$ は「ダイヤモンド・プリンセス号で確認された 60-69 歳の感染者数」を示します。(実際には、以後、特定の年齢層について言及せず、 $N_{{\rm dp}}^{a}$ と書いて、「ダイヤモンド・プリンセス号で確認された年齢層 $a$ の感染者数」を表すような使い方をしています。)

推定対象として、仮定が正しい場合の推定致死率を $p$ とします。また、この推定致死率 $p$ を年齢層ごとに分けたものを $p^{60-69}$ のように (下付き文字なしで) 表現します。

中国に於ける年齢層 $a$ の感染者に占める死亡者数は二項分布に従い、

$$ D_{{\rm ch}}^a \sim {\rm binomial}(N_{{\rm ch}}^a, p_{{\rm ch}}^a) $$

とモデル化できます。

また、ダイヤモンド・プリンセス号に於ける、感染者に占める死亡者数も同様に、

$$ D_{{\rm dp}} \sim {\rm binomial}(N_{{\rm dp}}, p_{{\rm dp}}) $$

となります。ダイヤモンド・プリンセス号の死亡者の年齢分布が手元にないため、こちらは年齢層ごとにモデリングしていないことに注意してください。

このふたつを仮定の章で言及した捕捉率 $c$ でつなぎ合わせます。

捕捉率に関してはなんの情報もないので、一様分布を仮定します。

$$ c \sim {\rm uniform}(0, 1) $$

$c$ は年齢層や致死率に関わらず一定と仮定したので、 $p_{{\rm ch}}^a$ に対して独立です。だから、年齢層 $a$ の推定致死率は、

$$ p^a = c \times p_{{\rm ch}}^a $$

全年齢では、

$$ p = c \times {1 \over N_{{\rm ch}}} \sum_a{p_{{\rm ch}}^a \times N_{{\rm ch}}^a} = {1 \over N_{{\rm ch}}} \sum_a{p^a \times N_{{\rm ch}}^a} $$

となります。

また、ダイヤモンド・プリンセス号の致死率は推定致死率に従う ($p_{{\rm dp}}^a = p^a$) と仮定したので、 $p_{{\rm dp}}$ は

$$ p_{{\rm dp}} = {1 \over N_{{\rm dp}}} \sum_a{p^a \times N_{{\rm dp}}^a} $$

となります。

実装

役者が揃ったので、 Stan で実装していきます。

入力データ

変数名は、年齢層を 60-69 歳なら _6x (i.e. $N_{{\rm ch}}^{60-69}$ = N_ch_6x)、 80 歳以上を _80_up (i.e. $N_{{\rm ch}}^{80-}$ = N_ch_80_up) と表現していること以外は概ねモデルの章の記号と同様です。

ダイヤモンド・プリンセス号のデータ ¹ は 90 代に関してのデータがありますが、中国のデータ ² は 80 歳以上でまとめられている点に注意してください。

data {
    // ダイヤモンド・プリンセス号に於ける死者数
    int D_dp;

    // ダイヤモンド・プリンセス号に於ける年齢層ごとの感染者数
    int N_dp_0x;
    int N_dp_1x;
    int N_dp_2x;
    int N_dp_3x;
    int N_dp_4x;
    int N_dp_5x;
    int N_dp_6x;
    int N_dp_7x;
    int N_dp_8x;
    int N_dp_9x;

    // 中国に於ける年齢層ごとの死亡者数
    int D_ch_0x;
    int D_ch_1x;
    int D_ch_2x;
    int D_ch_3x;
    int D_ch_4x;
    int D_ch_5x;
    int D_ch_6x;
    int D_ch_7x;
    int D_ch_80_up;

    // 中国に於ける年齢層ごとの感染者数
    int N_ch_0x;
    int N_ch_1x;
    int N_ch_2x;
    int N_ch_3x;
    int N_ch_4x;
    int N_ch_5x;
    int N_ch_6x;
    int N_ch_7x;
    int N_ch_80_up;
}

Stan では便利のために加工したデータを定義しておけるので、全年齢層の値を計算しておきます。

transformed data {
    // ダイヤモンド・プリンセス号に於ける感染者数
    int N_dp;

    // 中国に於ける感染者数
    int N_ch;

    N_dp = N_dp_0x + N_dp_1x + N_dp_2x + N_dp_3x + N_dp_4x + N_dp_5x + N_dp_6x + N_dp_7x + N_dp_8x + N_dp_9x;
    N_ch = N_ch_0x + N_ch_1x + N_ch_2x + N_ch_3x + N_ch_4x + N_ch_5x + N_ch_6x + N_ch_7x + N_ch_80_up;
}

推定対象パラメータ

推定対象のパラメータを指定します。

parameters {
    // 捕捉率 c
    real<lower=0, upper=1> c;

    // 中国に於ける致死率
    real<lower=0, upper=1> p_ch_0x;
    real<lower=0, upper=1> p_ch_1x;
    real<lower=0, upper=1> p_ch_2x;
    real<lower=0, upper=1> p_ch_3x;
    real<lower=0, upper=1> p_ch_4x;
    real<lower=0, upper=1> p_ch_5x;
    real<lower=0, upper=1> p_ch_6x;
    real<lower=0, upper=1> p_ch_7x;
    real<lower=0, upper=1> p_ch_80_up;
}

改めてこう書いてみると、中国に於ける致死率は普通に $D_{{\rm ch}} \over N_{{\rm ch}}$ で計算できるから推定対象ではないかと思われる方もいらっしゃると思います。

実は、中国に於ける致死率はまだこのモデルでは確定していません。例えばサイコロを 2 回振って 2 回とも 6 が出たとして、そのサイコロは 6 が 100% 出るサイコロだとは言えないように、感染者数と死亡者数が与えられても直ちに致死率を確定することはできないのです。

こういう不確定のものを不確定のまま簡単に扱えるのがベイズモデリングの良さであり、 Stan の良さですね。

さて、データと同じように、 parameters も値の変換を行えるので、やっておきます。

transformed parameters {
    // 年齢ごとの推定致死率
    real<lower=0, upper=1> p_0x;
    real<lower=0, upper=1> p_1x;
    real<lower=0, upper=1> p_2x;
    real<lower=0, upper=1> p_3x;
    real<lower=0, upper=1> p_4x;
    real<lower=0, upper=1> p_5x;
    real<lower=0, upper=1> p_6x;
    real<lower=0, upper=1> p_7x;
    real<lower=0, upper=1> p_80_up;

    // 推定致死率
    real<lower=0, upper=1> p;

    // ダイヤモンド・プリンセス号に於ける致死率
    real<lower=0, upper=1> p_dp;

    p_0x    = c * p_ch_0x;
    p_1x    = c * p_ch_1x;
    p_2x    = c * p_ch_2x;
    p_3x    = c * p_ch_3x;
    p_4x    = c * p_ch_4x;
    p_5x    = c * p_ch_5x;
    p_6x    = c * p_ch_6x;
    p_7x    = c * p_ch_7x;
    p_80_up = c * p_ch_80_up;

    p = (p_0x    * N_ch_0x +
         p_1x    * N_ch_1x +
         p_2x    * N_ch_2x +
         p_3x    * N_ch_3x +
         p_4x    * N_ch_4x +
         p_5x    * N_ch_5x +
         p_6x    * N_ch_6x +
         p_7x    * N_ch_7x +
         p_80_up * N_ch_80_up
         ) / N_ch;

    p_dp = (p_0x    * N_dp_0x +
            p_1x    * N_dp_1x +
            p_2x    * N_dp_2x +
            p_3x    * N_dp_3x +
            p_4x    * N_dp_4x +
            p_5x    * N_dp_5x +
            p_6x    * N_dp_6x +
            p_7x    * N_dp_7x +
            p_80_up * (N_dp_8x + N_dp_9x)
            ) / N_dp;
}

モデル

モデルの章で書いたままなので特に言うことがないですね。この表現力が Stan の良さです。

model {
    c ~ uniform(0, 1);

    D_ch_0x    ~ binomial(N_ch_0x,    p_ch_0x);
    D_ch_1x    ~ binomial(N_ch_1x,    p_ch_1x);
    D_ch_2x    ~ binomial(N_ch_2x,    p_ch_2x);
    D_ch_3x    ~ binomial(N_ch_3x,    p_ch_3x);
    D_ch_4x    ~ binomial(N_ch_4x,    p_ch_4x);
    D_ch_5x    ~ binomial(N_ch_5x,    p_ch_5x);
    D_ch_6x    ~ binomial(N_ch_6x,    p_ch_6x);
    D_ch_7x    ~ binomial(N_ch_7x,    p_ch_7x);
    D_ch_80_up ~ binomial(N_ch_80_up, p_ch_80_up);

    D_dp ~ binomial(N_dp, p_dp);
}

実行

データ

以下のデータを入力しました。
これらは仮定の章で記載したデータに基づきます。

data = {
    # ダイヤモンド・プリンセス号に於ける死者数
    'D_dp': 6,

    # ダイヤモンド・プリンセス号に於ける年齢層ごとの感染者数
    'N_dp_0x':   1,
    'N_dp_1x':   5,
    'N_dp_2x':  28,
    'N_dp_3x':  34,
    'N_dp_4x':  27,
    'N_dp_5x':  59,
    'N_dp_6x': 177,
    'N_dp_7x': 234,
    'N_dp_8x':  52,
    'N_dp_9x':   2,

    # 中国に於ける死者数
    'D_ch_0x':      0,
    'D_ch_1x':      1,
    'D_ch_2x':      7,
    'D_ch_3x':     18,
    'D_ch_4x':     38,
    'D_ch_5x':    130,
    'D_ch_6x':    309,
    'D_ch_7x':    312,
    'D_ch_80_up': 208,

    # 中国に於ける年齢層ごとの感染者数
    'N_ch_0x':     416,
    'N_ch_1x':     549,
    'N_ch_2x':    3619,
    'N_ch_3x':    7600,
    'N_ch_4x':    8571,
    'N_ch_5x':   10008,
    'N_ch_6x':    8583,
    'N_ch_7x':    3918,
    'N_ch_80_up': 1408,
}

実行環境

前回と同様、 PyStan を使って実行しました。

全体のソースコードはこちら

結果

parameters をたくさん設定したので、実行するとなにやら長大な結果が得られます。

             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
c            0.21  8.7e-4   0.08   0.08   0.15    0.2   0.25   0.38   7792    1.0
p_ch_0x    2.4e-3  3.0e-5 2.3e-3 7.1e-5 7.1e-4 1.7e-3 3.2e-3 8.3e-3   5850    1.0
p_ch_1x    3.6e-3  3.3e-5 2.6e-3 4.3e-4 1.7e-3 3.0e-3 4.9e-3   0.01   6364    1.0
p_ch_2x    2.2e-3  9.7e-6 8.0e-4 9.5e-4 1.6e-3 2.1e-3 2.7e-3 4.0e-3   6894    1.0
p_ch_3x    2.5e-3  6.3e-6 5.7e-4 1.5e-3 2.1e-3 2.5e-3 2.9e-3 3.7e-3   8234    1.0
p_ch_4x    4.6e-3  8.7e-6 7.4e-4 3.2e-3 4.1e-3 4.5e-3 5.0e-3 6.1e-3   7259    1.0
p_ch_5x      0.01  1.3e-5 1.1e-3   0.01   0.01   0.01   0.01   0.02   7562    1.0
p_ch_6x      0.04  2.4e-5 2.0e-3   0.03   0.03   0.04   0.04   0.04   7096    1.0
p_ch_7x      0.08  5.2e-5 4.3e-3   0.07   0.08   0.08   0.08   0.09   6650    1.0
p_ch_80_up   0.15  1.1e-4 9.2e-3   0.13   0.14   0.15   0.15   0.17   7552    1.0
p_0x       4.8e-4  7.1e-6 5.3e-4 1.2e-5 1.3e-4 3.1e-4 6.5e-4 1.9e-3   5530    1.0
p_1x       7.4e-4  8.5e-6 6.3e-4 7.0e-5 3.0e-4 5.8e-4 9.8e-4 2.4e-3   5462    1.0
p_2x       4.5e-4  3.2e-6 2.4e-4 1.3e-4 2.8e-4 4.0e-4 5.7e-4 1.1e-3   5637    1.0
p_3x       5.2e-4  2.7e-6 2.3e-4 1.8e-4 3.4e-4 4.8e-4 6.4e-4 1.1e-3   7336    1.0
p_4x       9.4e-4  4.6e-6 3.9e-4 3.6e-4 6.5e-4 8.8e-4 1.2e-3 1.9e-3   7248    1.0
p_5x       2.7e-3  1.2e-5 1.0e-3 1.1e-3 1.9e-3 2.5e-3 3.3e-3 5.0e-3   7484    1.0
p_6x       7.4e-3  3.1e-5 2.8e-3 3.0e-3 5.4e-3 7.0e-3 9.0e-3   0.01   7796    1.0
p_7x         0.02  6.9e-5 6.1e-3 6.6e-3   0.01   0.02   0.02   0.03   7770    1.0
p_80_up      0.03  1.3e-4   0.01   0.01   0.02   0.03   0.04   0.06   7937    1.0
p          4.7e-3  2.0e-5 1.8e-3 1.9e-3 3.5e-3 4.5e-3 5.8e-3 8.6e-3   7853    1.0
p_dp         0.01  4.7e-5 4.2e-3 4.5e-3 8.3e-3   0.01   0.01   0.02   7865    1.0
lp__        -4215    0.06   2.35  -4220  -4216  -4214  -4213  -4211   1641    1.0

重要なところを抜き出すと

             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
c            0.21  8.7e-4   0.08   0.08   0.15    0.2   0.25   0.38   7792    1.0
p          4.7e-3  2.0e-5 1.8e-3 1.9e-3 3.5e-3 4.5e-3 5.8e-3 8.6e-3   7853    1.0

となります。

すなわち、仮定が正しいならば、中国での感染者の捕捉率は 20% くらいで、それを考慮した致死率は 0.5% くらいであり、 95% 信頼区間においては 0.19% 〜 0.86% と推定されました。

前述の通り、実際にはこの仮定は正確ではないため、 95% 信頼区間なんて何も信頼できないのですが、乗客の過半が 60 歳以上だったダイヤモンド・プリンセス号に於ける致死率が実測 1% くらいなので、この推定は大きく外してはいないんじゃないかと期待しています。

まとめ

今回は前回と違って Stan を使って新型コロナウイルスの致死率を推定することができました。

前よりは Stan の良さも表現できたんじゃないかと思います。 (前ができてなさすぎた)

結果として、強めの仮定を置いたので眉に唾を塗っておく必要があるものの、中国の感染状況に対応する推定致死率として 0.19% 〜 0.86% を得ました。

余談

前回の投稿で、二項分布をベータ分布に変換して推定を行いましたが、今回も、中国に於ける年齢層 $a$ の感染者に占める死亡者数のモデリングを

$$ D_{{\rm ch}}^a \sim {\rm binomial}(N_{{\rm ch}}^a, p_{{\rm ch}}^a) $$

ではなく、

$$ p_{{\rm ch}}^a \sim {\rm beta}(D_{{\rm ch}}^a + 1, N_{{\rm ch}}^a - D_{{\rm ch}}^a + 1) $$

に置き換えて、同じように推論できます。興味がある方は試してみてください。