いろいろなことを考察する

「42万死亡推計」は、2020-03-19 専門家会議の資料と同じ計算


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



西浦氏が 2020-04-15 に行った「42万死亡推計」の計算は、約一ヵ月前に行われた 2020-03-19 専門家会議の資料で西浦氏が示した計算と同じものです。西浦氏がそう述べています。

この点、必ずしも全ての人がご存じではないと思いますので説明します。


(1)序論

西浦氏は、2020-04-15 に非公開の記者会見を開き、「42万死亡推計」を発表しました。
「42万死亡推計」は、発表の約一ヵ月前、2020-03-19 に行われた専門家会議の資料で西浦氏が示した計算と同じものです。

2020-03-19 専門家会議
https://www.kantei.go.jp/jp/singi/novel_coronavirus/taisaku_honbu.html
の資料
https://www.kantei.go.jp/jp/singi/novel_coronavirus/senmonkakaigi/sidai_r020319.pdf
の9ページ、図6左と右の計算です。

西浦氏が発表の一ヵ月半後に「42万死亡推計」を説明した Newsweek 記事 [*1-1] などで、そう説明しています。


[*1-1] 紙媒体(2020-06-02 発売の、Newsweek 2020-06-09 号)と 2020-06-11 公開の電子記事 の「「8割おじさん」の数理モデルとその根拠」。どちらも github 資料 を参照している。

(2)2つで同じグラフを示している

  • 「42万死亡推計」を説明する Newsweek 記事 [*1-1] と、
  • 2020-03-19 専門家会議の資料とは、

同じグラフを示しています。

西浦氏が「42万死亡推計」を説明する Newswek 記事 [*1-1] が示すグラフがこちらです。

図1 Newsweek 2020-06-11 記事に掲載のグラフ
[図2-1] Newsweek 2020-06-11 記事に掲載のグラフ

2020-03-19 専門家会議の資料、9ページの図6の左がこちらです。

図2 2020-03-19 専門家会議の資料、9ページの図6左
[図2-2] 2020-03-19 専門家会議の資料、9ページの図6左

[図2-2] には、

図6.大規模流行時に想定される10万人当たりの新規感染者数(左)と重篤患者数(右)

と説明があります。

2つのグラフは同じものです。したがって、2つの計算は同じものです。

(3)Newsweek 記事の本文でも、2つの計算が同じだと説明している

[*1-1] の、上に [図2-1] として示した図の付近に、以下の記述があります。

SIRモデルを活用したシミュレーションで42万人死亡と試算したものが、上のグラフだ(3月19日の専門家会議の提言書内に示したもの)。

「42万死亡推計」のグラフは、2020-03-19 専門家会議の資料のグラフと同じだ、と説明しています。

(4)結論:2つの計算は同じもの

西浦氏は、(2)に示したように [*1-1] のグラフにおいて、また(3)に示したように [*1-1] の本文において、「42万死亡推計」と「2020-03-19 専門家会議の資料で示した計算」が、同じ計算だと述べています。

更新履歴

  • 2024-03-22
    公開。


「42万死亡推計」の計算方法では、「死亡者数」は「感染者総数」にほぼ比例する


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



目次


(1)序論

「42万死亡推計」の計算方法では、「死亡者数」は「感染者総数(感染したことのある人の数)」に「ほぼ比例」します。

これは意外なことではないので、同じ様にお考えであれば、本エントリを読んで頂く必要はありません。

「42万死亡推計」の問題点を述べる中でこの「ほぼ比例する」という特徴を使うのですが、「ほぼ比例する」が自明とまでは言えないかと思ったので、こうしてご説明しています。

本エントリでは SIR モデル内での計算には触れません。触れる必要はないからです。
(本エントリは、twitter のスレ (まだ書いてません。書いたら貼ります。) と同じ論旨ですが、本エントリの方が詳しいはずです。)

図だと、こうなります。

図1 「42万死亡推計」計算の概要
[図1] 「42万死亡推計」計算の概要

(2)「42万死亡推計」の前提

計算に際して西浦氏は、以下などを前提としました [*2-1]。

  • (a) 一度感染したら再感染はしません。
  • (b) 全国民を年齢で3つの世代(子供 (0~14)、成人 (15~64)、高齢 (65~))に分けました。
  • (c)「流行終了まで、何も感染症対策はしない」との現実にはあり得ないシナリオを採用。
  • (d) 基本再生産数(皆に免疫がない状態で、1人の感染者が感染させる人数の平均値)の値は不明だが、2.5 と仮定。

この他の主な前提としては、

  • (e) 世代毎に感染しやすさが違うと想定(感受性の異質性の導入)。
  • (f) 世代毎の感染させやすさの違いは、データ不足により組み込めなかった(感染性の異質性は導入しない)。
  • (g) 計算開始時には、成人世代に10人の感染者がいると仮定。

があります。


[*2-1] 紙媒体(2020-06-02 発売の、Newsweek 2020-06-09 号)と 2020-06-11 公開の電子記事 の「「8割おじさん」の数理モデルとその根拠」。どちらも github 資料 を参照している。

(3)各世代の感染者総数の計算

これら前提を組み込んだ SIR モデルを使って、流行終了までの感染の推移を計算します。 すると、子供、成人、高齢それぞれの、感染者総数(=感染したことのある人の数)を得ることができます。

こうして、
(高齢の感染者総数)
(成人の感染者総数)
(子供の感染者総数)
が計算できました。

この3つを合計すると、(全体の感染者総数) になります。

図1 「42万死亡推計」計算の概要(「全体の感染者総数」まで)
[図2] 「42万死亡推計」計算の概要(「全体の感染者総数」まで)

(4)全体の死亡者数の計算

西浦氏は、子供の死亡者数をゼロと考えた [*2-1] ので、以下では子供を省きます。

西浦氏は、(高齢の感染者死亡率) と (成人の感染者死亡率) という定数を用意していたので、計算で得た感染者総数とこれら定数があると、(高齢の死亡者数) と (成人の死亡者数) が計算できます。

(高齢の死亡者数)  = (高齢の感染者総数)  × (高齢の感染者死亡率)
(成人の死亡者数)  = (成人の感染者総数)  × (成人の感染者死亡率) ...... (式1)

(全体の死亡者数) は、(高齢の死亡者数) と (成人の死亡者数) の合計です。

図1 「42万死亡推計」計算の概要(「全体の死亡者数」まで)
[図3] 「42万死亡推計」計算の概要(「全体の死亡者数」まで)

(5)結論:全体の死亡者数は、全体の感染者総数にほぼ比例する

(式1) で計算しているので、
(高齢の死亡者数)は、(高齢の感染者総数)に正確に比例します。
(成人の死亡者数)は、(成人の感染者総数)に正確に比例します。

したがって、
「全体の死亡者数」は、「全体の感染者総数」にほぼ比例します。

世代毎ならば、死亡者数と感染者総数は正確に比例します。
全体の死亡者数は、全体の感染者総数に正確には比例せず、ほぼ比例となります。

以上、

「42万死亡推計」の計算方法では、「死亡者数」は「感染者総数(感染したことのある人の数)」に「ほぼ比例」します。

をご説明しました。

更新履歴

  • 2024-03-19
    公開。


西浦氏が示す「42万死亡推計」のパラメータの不自然さ


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



本エントリは、西浦氏が示す「42万死亡推計」に関するパラメータの値に不自然な点があることを示します。私は、西浦氏が示す「42万死亡推計」のパラメータに2%ほどの誤りがある(https://sarkov28.hatenablog.com/entry/2024/03/16/092400)と指摘し、正しいと考えるパラメータを示していますが、私が正しいと考えるパラメータには、本エントリが指摘する不自然さはありません。本エントリは、「2%ほどの誤りがある」との私の指摘が正しいことの傍証になっています。


(1)序論

西浦氏は、2020-04-15 の非公開の会合で、専門家個人として「42万死亡推計」を発表しました。
私は「42万死亡推計」に関して幾つかの問題を指摘していて、それはこちら(https://sarkov28.hatenablog.com/entry/2022/03/29/160915)の (2-1) にまとめてあります。

本エントリが扱うのは、「西浦氏が示しているパラメータに不自然な点がある」という件です。

(2)では、本エントリが扱う資料や数値を示します。(3)では、西浦氏が示している計算に隠れている不自然な値のパラメータについて説明します。(4)では結論を述べます。(5)では、不自然な値のパラメータがある理由について考察します。

(2)本エントリで扱う資料や数値など

2020-04-15 に発表された「42万死亡推計」についての説明としては、新聞や雑誌の記事があります。[*2-1] に挙げたのは、その中でも数値を詳しく報じているものです。

また西浦氏は、「42万死亡推計」について、2020-06-02 発売の Newsweek の記事 [*2-2] で説明しました。ここで西浦氏は github 資料 も示しています。

言うまでもないことですが、[*2-1] と [*2-2] で説明されている計算は、同じものであるはずです。

西浦氏は国民を3つの世代に分けて計算しました。3つの世代は、年齢で (0-14), (15-64), (65-) です。西浦氏は github 資料 でこの3世代をそれぞれ、c, a, e と書いています。(それぞれ、child, adult, elder の略だと思われます。)

3つの世代のうち、世代 c の扱いは他とは違います。github 資料 に I_c とあるように「世代 c にも感染者がいる」のですが、[*2-2] によれば「世代 c の重症者数はほぼゼロ、死亡者数はゼロ」というのが西浦氏の考えです。本エントリでは重症者数や死亡者数において世代 c のことを考慮しません。死亡者数を検討する上で支障がないからです。

[*2-1] [*2-2] の資料には以下のような数値があります。死亡者数や重症者数は流行終了までの値です。

名称 出所
(2a) 基本再生産数  R_0 2.5 注 [*2-1]
(2b) 世代 a の重症者数 20万1301人 注 [*2-1]
(2c) 世代 e の重症者数 65万2066人 注 [*2-1]
(2d) 全世代の重症者数 85万3367人 注 [*2-1]
(2e) 世代共通の重症者死亡率 約49% 注 [*2-1]
(2f) 死亡者数 約41.8万人 注 [*2-1]
(2g) 世代 a の感染者死亡率 0.0015 注 [*2-2]
(2h) 世代 e の感染者死亡率 0.0100 注 [*2-2]

[表2-1] [*2-1] [*2-2] の資料の「42万死亡推計」に関する数値

なお、(2g)(2h) は [*2-1] にも記載がありますが、[*2-2] を出所としたのはこちらは計算コード中の記述なので、示されている数値をそのまま計算に使うことが確実な、概数では「ない」値のはずだからです。


[*2-1] 2021-05 発刊の別冊 Newton。2020-04-15 の毎日記事にも同様の記述がある。2020-04-15 毎日「新型コロナ 85万人が重篤に 「対策なしなら」厚労省班試算」念の為のアーカイブはこちら
[*2-2] 紙媒体(2020-06-02 発売の、Newsweek 2020-06-09 号)と 2020-06-11 公開の電子記事 の「「8割おじさん」の数理モデルとその根拠」。どちらも github 資料 を参照している。

(3)西浦氏が示している計算の不自然な値のパラメータ

西浦氏は、「42万死亡推計」において「重症者数が出てくる計算」と「重症者数が出てこない計算」という2つの計算を示しています。この2つは同じ計算であるはずですが、同じ計算だとすると、西浦氏が不自然な値のパラメータを使った事になります。

(3-1) 重症者数が出てくる「42万死亡推計」計算

「42万死亡推計」を発表した非公開の会合において西浦氏は、メディアに

  • 重症者数
  • 重症者死亡率

を説明した、との記述が共著書や他の報道などにあります。

実は42万人という数は、僕の口からは言っていなんです。(略)「じゃあ、85万人重症でその約半分が死亡」ならいいんですかと聞きますと、「それならよい」ということでした。(「理論疫学者・西浦博の挑戦」[*3-1-1] 第3章9)

ここでは「85万人重症」「約半分が死亡」と表記されていますが、実際に説明されたのはもっと詳しい数値であり、[表2-1] の一部の (2b)~(2e) だと思われるので、再掲します。

名称 出所
(2b) 世代 a の重症者数 20万1301人 注 [*2-1]
(2c) 世代 e の重症者数 65万2066人 注 [*2-1]
(2d) 全世代の重症者数 85万3367人 注 [*2-1]
(2e) 世代共通の重症者死亡率 約49% 注 [*2-1]

[表3-1] 「42万死亡推計」を報じた記事による西浦氏の説明内容


[*3-1-1] 西浦博、川端裕人「理論疫学者・西浦博の挑戦」2020-12、中央公論新社

(3-2) 重症者数が出てこない「42万死亡推計」計算

西浦氏が示している github 資料 が示す計算には、重症者数が出てきません。

github 資料には、[表2-1] の一部 (2g)(2h) があるので、再掲します。

名称 出所
(2g) 世代 a の感染者死亡率 0.0015 注 [*2-2]
(2h) 世代 e の感染者死亡率 0.0100 注 [*2-2]

[表3-2] github 資料 の「42万死亡推計」に関する数値

(3-3) 「42万死亡推計」の概要

(3-1)と(3-2)はどちらも「42万死亡推計」ですから、同じ計算でなければいけません。西浦氏は(3-1)で重症者数を計算しています。(3-2)に示された「計算コード」から考えて、西浦氏は

  • (3-3a) github 資料 のSIR方程式を使って、感染者総数(計算期間中に感染したことのある人数)を計算し、
  • (3-3b) 感染者総数に感染者重症化率を掛算して重症者数を計算し、
  • (3-3c) 重症者数に重症者死亡率(約49%)を掛算して死亡者数を計算した、

のだと思われます。重症者数の計算以後は、世代 a, e それぞれで行っています。

(3-4) 世代 a の感染者重症化率は不自然な値になる

西浦氏は感染者重症化率を明示していませんが、上に挙げた公開情報から計算することができます。これが不自然な値になります。まず世代 a について論じます。明確にするために数値は ( ) に入れて表記します。

表2-1 に 世代 a の (感染者死亡率) が示されています。0.0015 です。
表2-1 に(世代共通の)(重症者死亡率) が示されています。約49%です。

この2つを分数の形で書くと以下になります。

 (死亡者数)/(感染者総数) = 0.0015,\ \ \ \ (死亡者数)/(重症者数) = (約49%) ...... (式3-1)

感染者総数というのは「計算期間中に感染したことのある人の数」です。「感染者数」と書くと、「その時点で感染している人」と紛らわしいのでこの表記にします。

(式3-1)から、(感染者重症化率) を計算することができます。

 \displaystyle (死亡者数) = 0.0015 × (感染者総数) = (約49%) × (重症者数) ...... (式3-2)

 (感染者重症化率) = (重症者数)/(感染者総数) = 0.0015/(約49%) ...... (式3-3)

(式3-3)第3項の  (重症者死亡率) = (約49%) を単に  49% として計算すると、世代 a においては

 (感染者重症化率) = 0.00306122 ...... (式3-4)

となります。この (感染者重症化率) は中途半端で不自然な値です。(重症者死亡率) が「49%」ではなく「約49%」であることを考慮し、

 0.485 ≦ (重症者死亡率)\ <\ 0.495 ...... (式3-5)

としても、

 0.00303030\ <\ (感染者重症化率) ≦ 0.00309278 ...... (式3-6)

となるだけで不自然さは解消されません。

なお、ここで有効数字を6桁にしているのは、この値を使って西浦氏が計算したはずの重症者数が有効数字6桁だから([表2-1] の (2b))です。

(3-5) 世代 e の感染者重症化率も不自然な値になる

世代 e に関しても同様の結果になります。

 \displaystyle (死亡者数) = 0.0100 × (感染者総数) = 約49% × (重症者数) ...... (式3-6)

から、途中を省略して、世代 e では

 (感染者重症化率) = 0.0204082 ...... (式3-7)

となり、やはり中途半端で不自然な値です。世代 a と同様に範囲を検討しても、

 0.0202020\ <\ (感染者重症化率) ≦ 0.0206186 ...... (式3-8)

となるだけで不自然なままです。

(4)結論

西浦氏は「42万死亡推計」の計算として2つの説明を示しています。この2つが同じ計算であることから、西浦氏は感染者重症化率として、世代 a で 0.00306122、世代 e で 0.0204082、という不自然な値のパラメータを使った事を導くことができます。

この不自然なパラメータは、西浦氏が示す「42万死亡推計」の説明に誤りがあることを示唆します。次項では私が別エントリで行った検討を元にして、この誤りについて考察します。

(5)考察:世代 a や世代 e の感染者重症化率が不自然な値になる理由

西浦氏は、「42万死亡推計」の計算として、重症者数と死亡者数を示しているのですから、これらを計算するためのパラメータを用意していたと考えるしかありません。

西浦氏による github 資料 と本エントリでの検討によれば、西浦氏が重症者数と死亡者数を計算するために用意したパラメータは以下になります。

名称 出所
(5a) 世代 a の感染者重症化率 0.00306122 本エントリ (式3-4)
(5b) 世代 e の感染者重症化率 0.0204082 本エントリ (式3-7)
(2g) 世代 a の感染者死亡率 0.0015 注 [*2-2]
(2h) 世代 e の感染者死亡率 0.0100 注 [*2-2]

[表5-1] github 資料 による、感染者数から重症者数と死亡者数を計算するためのパラメータ

一方、私の 別エントリ での検討によると、西浦氏が重症者数と死亡者数を計算するのに使ったパラメータは以下です。

名称 出所
(5c) 世代 a の感染者重症化率 0.003 別エントリ
(5d) 世代 e の感染者重症化率 0.02 別エントリ
(2e) 世代共通の重症者死亡率 約49% 注 [*2-1]

[表5-2] 私の 別エントリ での検討による、感染者数から重症者数と死亡者数を計算するためのパラメータ

2つの表を比較すると、西浦氏が用意していたパラメータとして合理的なのは [表5-2] の方です。[表5-1] は合理性を欠いているだけでなく、このパラメータで計算した死亡者数は、「42万」や「41.8万」と2%ほどの相違を生じます。[表5-2] のパラメータで計算した死者数は、「42万」や「41.8万」になります。

本エントリが問題視している不自然さが発生した理由は、西浦氏による github での情報開示が誤っているからだと思われます。

上で示してきたことと実質同等ですが、参考のため、私の検討による感染者重症化率や感染者死亡率を、西浦氏による値と比較しておきます。

名称 西浦氏の
github 資料
による
私の
別エントリでの
検討による
世代 a の感染者重症化率 0.00306122 0.00300000
世代 e の感染者重症化率 0.0204082 0.0200000
世代 a の感染者死亡率 0.0015 0.00147
世代 e の感染者死亡率 0.0100 0.0098

[表5-3] 感染者死亡率の比較(西浦氏によるものと、私の検討によるもの)

左右には、2%ほどの相違があります。この相違が、[github 資料] で死者数を計算すると2%ほどの相違が生じる原因です。
私の検討による感染者死亡率の計算式は、以下です。
 (感染者死亡率)=(感染者重症化率)×(重症者死亡率)=(感染者重症化率)×(49%)

更新履歴

  • 2024-03-18
    公開。本エントリに関係する説明を twitter に書くかも知れません。書いたら冒頭付近にリンクを書いておきます。


西浦氏が github で示している「42万死亡推計」のパラメータの誤り(改訂版)


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



本エントリは、以前に書いた ブログエントリ の改訂版です。パラメータを修正する際の方針を一部で変更した結果、議論で出てくる数値の一部が約0.036%変化しました。また全体を改稿しました。今回の変更後も、「西浦氏が公開しているパラメータの一部に約2%の誤りがある」との論旨は変わりません。変更した内容や詳細は(11)に書きました。


(1)序論

西浦氏は、2020-04-15 の非公開の会合で、専門家個人として「42万死亡推計」を発表しました。
私が「42万死亡推計」に関して指摘している問題は、こちら(https://sarkov28.hatenablog.com/entry/2022/03/29/160915)の (2-1) にまとめてあります。

本エントリが扱うのは上記とは別の問題であり、「西浦氏が開示している「42万死亡推計」のパラメータの一部に約2%の誤りがある」という件です。

「42万死亡推計」が国民に大きな影響を与えることを西浦氏は発表前に理解していました [*1-1]。国民に大きな影響を与える計算であれば、それが専門家個人としての発表だったとしても、西浦氏は計算の詳細情報を「42万死亡推計」の発表と同時に開示すべきでした。計算根拠の説明として必要ですし、国民による検証を仰ぐためでもあります。

西浦氏は、「42万死亡推計」発表から一ヵ月半も経過してから、やっと計算の詳細を開示しました [*1-3]。「42万死亡推計」は、2020-03 月からの第1波での行動制限との関わりで示された [*1-4][*1-6] のですから、緊急事態宣言が全国で解除された 2020-05-25 の後になってようやく開示しても、時宜を逸していて不適切です。

西浦氏による詳細開示は時宜を逸していましたが、それでも事後的な検証の役には立つはずでした。ところがここでも問題が起きました。開示された情報(パラメータ)を使っても「42万死亡推計」が再現できないのです。 計算結果は約2%異なった値になります。

この2%の差が生じる原因は、西浦氏が開示したパラメータに誤りがあるからです。以下に示すのは、西浦氏のパラメータが間違っていると合理的に考えられる理由の説明です。この誤りが意図的なものならば極めて重大な問題です。私はこの誤りが過失によるものだとは思えません。

少なくとも西浦氏は、開示したパラメータで「42万死亡推計」を計算したことは一度もありません。一度でも計算すれば、約2%の差が生じることはすぐに分かります。西浦氏が計算したのなら、約2%結果が異なるのを知った上で開示したことになります。いずれにしても極めて不適切な情報開示の姿勢です。

私がこの件を最初に指摘したのは、2021-07-08 でした(ブログtwitter)。これ以降、この件について的を射た異論はありません。これは私の指摘に問題がないことの傍証にはなると思います。


[*1-1] 「もちろん、最悪の被害想定を語ることで、ざわざわすることは覚悟していました。」([*1-2] 第3章9)
[*1-2] 西浦博、川端裕人「理論疫学者・西浦博の挑戦」2020-12、中央公論新社
[*1-3] 紙媒体(2020-06-02 発売の、Newsweek 2020-06-09 号)と 2020-06-11 公開の電子記事 の「「8割おじさん」の数理モデルとその根拠」。どちらも github 資料 を参照している。
[*1-4] 「外出自粛に代表される行動制限によって、感染被害を軽減できることを市民に理解してもらうのが狙いという。」(2020-04-15 毎日 [*1-5])
[*1-5] 2020-04-15 毎日「新型コロナ 85万人が重篤に 「対策なしなら」厚労省班試算」念の為のアーカイブはこちら
[*1-6] 「何も対策をしなかったら、流行してこれくらいの規模の死亡者が出るという危機が目前に迫っていると、政権に腹を括ってもらうのが狙いでした。」(週刊文春 2020-04-30号)

(2)本エントリが示すこと

「西浦氏のパラメータが誤っていること」と「私の計算が正しいこと」の直接で詳細な根拠として、私はプログラム
https://colab.research.google.com/drive/1ZvY-HmOvh3AS7OnjJdFh0lBw_OfTX5cU?usp=sharing
を公開しています。本来はこのプログラムを見て頂くことで、ご説明は済んでいます。(このプログラムについては(9)の結論に少し説明を書きます。)しかしプログラミング言語である python をお使いにならない方へのご説明として、これはどうなのかと思います。

また、プログラムを見なくても、公開されている情報を検討したり、数値計算の途中で出てくる値を公開情報と照合すると、西浦氏のパラメータがおかしいと判断できる傍証があるので、それを以下でご説明します。

繰り返すようですが、本エントリで字数を使ってご説明しているのは傍証に過ぎません。説明の中核にあるのはプログラムですので、もし「傍証では足りない」とお考えであれば、プログラムの方をご参照下さい。

以下は傍証をどのようにご説明するのかです。

(3)では、本エントリの検討に用いる事項や数値を示します。(4)では、西浦氏の示しているパラメータに不自然な点があることを示します。これはパラメータのどこかに誤りがあることの傍証です。

(5)では、西浦氏のパラメータを使った数値計算の結果が、「42万」にならないことを示します。(6)では、私の数値計算と「42万死亡推計」を報じる記事の数値とを使って新しいパラメータを計算し、それが自然な値になることを示します。(4)で示した(西浦氏が示している)パラメータは不自然な値なのに、(6)で得る新しいパラメータは自然な値になります。これは新しいパラメータが正しいことの傍証です。

(7)では、西浦氏が示しているパラメータと新しいパラメータとを比較します。(8)では本エントリで得た新しいパラメータでの計算が「42万」になることを確認します。(9)では本エントリの結論を述べます。(10)で示すのは、「42万死亡推計」に関する、西浦氏への質問事項です。(11)では以前に示していた同主旨のエントリから本エントリへの変更点を説明します。

(3)本エントリで用いる数値など

「42万死亡推計」は、国民を3つの世代に分けて計算されました [*1-3]。3つの世代は、年齢で (0-14), (15-64), (65-) です。西浦氏は github 資料 でこの3世代をそれぞれ、c, a, e と書いています。(それぞれ、child, adult, elder の略だと思われます。)

3つの世代のうち、世代 c の扱いは他とは違います。github 資料 に I_c とあるように「世代 c にも感染者がいる」のですが、[*1-3] によれば「世代 c の重症者数はほぼゼロ、死亡者数はゼロ」というのが西浦氏の考えです。本エントリでは重症者数や死亡者数において世代 c のことを考慮しません。死亡者数を検討する上で支障がないからです。

本エントリでは、「42万死亡推計」に関して最も詳細だと思われる以下の数値を使います [*3-1] [*1-3]。死亡者数や重症者数は流行終了までの値です。

名称 出所
(3a) 基本再生産数  R_0 2.5 (略)
(3b) 世代 a の重症者数 20万1301人 注 [*3-1]
(3c) 世代 e の重症者数 65万2066人 注 [*3-1]
(3d) 全世代の重症者数 85万3367人 注 [*3-1]
(3e) 重症者死亡率 約49% 注 [*3-1]
(3f) 死亡者数 約41.8万人 注 [*3-1]
(3g) 世代 a の感染者死亡率 0.0015 注 [*1-3]
(3h) 世代 e の感染者死亡率 0.0100 注 [*1-3]

[表3-1] 本エントリで用いる「42万死亡推計」に関する数値

なお、西浦氏の示す数値の有効数字の取扱いには不適切なところがあります。表3-1 の基本再生産数  R_0 が 2.5 で有効数字2桁であるのに対し、(3b) では有効数字6桁であるなどです。本エントリは、正しい有効数字の扱いではなく、西浦氏の計算の検証が目的なので、可能な範囲で有効桁数の多い数値で議論します。


[*3-1] 2021-05 発刊の別冊 Newton。同様の記述は [*1-5] の 毎日記事 にもある。

(4)西浦氏が示している計算には不自然なパラメータが隠れている

西浦氏は、「42万死亡推計」において、2つの計算を示しています。「重症者数が出てくる計算」と「重症者数が出てこない計算」の2つです。この2つは同じ計算であるはずですが、同じ計算だとすると西浦氏は感染者重症化率で不自然な値のパラメータを使ったことになります。この件については、別エントリ で検討しました。

(5)西浦氏の github のパラメータを使った数値計算は「42万」にならない

西浦氏が github に示している SIR 方程式とパラメータを使って、私が作成したプログラム数値計算すると、以下の要領で全世代の死者数を計算することができます。有効数字は6桁としています。感染者総数は、「計算期間中に感染したことのある人の数」です。

計算方法 計算結果
(5a) github の SIR 方程式とパラメータ
による私の数値計算
世代 a の感染者総数
=67100300人
(5b) github の SIR 方程式とパラメータ
による私の数値計算
世代 e の感染者総数
=32603300人
(5c) 世代 a の
感染者総数(67100300人)
と、感染者死亡率(0.0015)
を掛算
世代 a の死亡者数
=100650人
(5d) 世代 a の
感染者総数(32603300人)
と、感染者死亡率(0.0100)
を掛算
世代 e の死亡者数
=326033人
(5e) 世代 a の死亡者数(100650人)
世代 e の死亡者数(326033人)
を合計
全体の死亡者数
=42.6683万人

[表5-1] 西浦氏の github のパラメータを使った「42万死亡推計」の計算

(5e)で得た全体の死亡者数「42.6683 万」は、どこで四捨五入しても「42万」や「41.8 万」になりません。ここには、約2%の差があります。( (42.6683 / 41.8) - 1 = 0.0208 = 2%

私の数値計算か、西浦氏のパラメータか、どちらかに誤りがあります。なお(4)に示したように、西浦氏の示している説明に誤りがないとすると不自然であることが分かっています。(4)で明らかになっている不自然さは、西浦氏のパラメータに誤りがあることの傍証になっています。

(6)私の数値計算と別冊 Newton 記事から得られる新しいパラメータは自然な値になる

(5)で数値計算した、世代 a, e の感染者総数と、別冊 Newton や毎日記事 [*1-5] の重症者数の値を使うと、(4)で不自然な値になった、世代 a, e の感染者重症化率を計算することができます。有効数字は6桁としています。感染者総数は、「計算期間中に感染したことのある人の数」です。

計算方法 計算結果
(6a) github のSIR 方程式とパラメータ
による私の数値計算
世代 a の
感染者総数
=67100300人
(6b) github のSIR 方程式とパラメータ
による私の数値計算
世代 e の
感染者総数
=32603300人
(6c) 表3-1
世代 a の重症者数(201301人) /
世代 a の感染者総数(67100300人)
世代 a の
感染者重症化率
=0.00300000
(6d) 表3-1
世代 e の重症者数(652066人) /
世代 e の感染者総数(32603300人)
世代 e の
感染者重症化率
=0.0200000

[表6-1] 私の数値計算と別冊 Newton 記事の重症者数から得られる感染者重症化率

私の数値計算を用いて新しく得たパラメータ、

  • 世代 a の感染者重症化率=0.00300000
  • 世代 e の感染者重症化率=0.0200000

は、有効数字6桁でキリのいい値になっています。偶然にゼロが連なっているのではなく、西浦氏が選んだ

  • 世代 a の感染者重症化率=0.003
  • 世代 e の感染者重症化率=0.02

という値を、私が本項で採用した方針で探り当てたのだと思われます(本項の方針=「私の数値計算による感染者総数」と「別冊 Newton の重症者数」から「感染者重症化率を計算する」という方針)。

(7)「西浦氏のパラメータ」と「新たなパラメータ」の比較

(4)の 別エントリ で得たパラメータは、西浦氏が示している情報が全て正しいとした場合のものでした。(6)で得たパラメータは、私の数値計算の感染者総数と記事の重症者数が正しいとした場合のものでした。この2つを比較します。有効数字6桁です。

別エントリ で計算
(重症者死亡率49%と
 github の感染者死亡率
 による)
(6)で計算
(私の数値計算の感染者総数と
 記事の重症者数
 による)
世代 a の
感染者重症化率
0.00306122 0.00300000
世代 e の
感染者重症化率
0.0204082 0.0200000

[表7-1] 西浦氏の github のパラメータによる感染者重症化率と、新しい感染者重症化率

  • 私の数値計算を使った値が、西浦氏の github のパラメータを使った値と2%ほど異なっていること、
  • 西浦氏の github のパラメータを使ったパラメータの値が不自然であること、
  • 私の数値計算の結果を使った新しいパラメータの値が自然であること、

が分かります。 新しいパラメータこそが、西浦氏が「42万死亡推計」を手元で計算した際に使ったパラメータです。

西浦氏が github に示していたのは感染者重症化率ではなく、感染者死亡率でした。これを新しいパラメータを使って計算すると、以下になります。

西浦氏の
github の値
(6)での
計算による値
世代 a の
感染者死亡率
0.0015 0.00147
世代 e の
感染者死亡率
0.0100 0.0098

[表7-2] 西浦氏の github の感染者死亡率と、新しい感染者死亡率

両者には、約2%の差があります。
 (0.0015 / 0.00147) - 1 = 0.0204 = 2%
 (0.0100 / 0.0098)  - 1 = 0.0204 = 2%

(5)では github のパラメータを使った「42万死亡推計」の計算が、2%ほどの差を生じていました。(5)における2%ほどの差は、ここに示した感染者死亡率の2%の差に起因していたことになります。(感染者死亡率が約2%変化すれば、死亡者数は約2%変化します。)

(8)新しいパラメータを使った数値計算は「42万」になる

新しいパラメータは、「私の数値計算の感染者総数」と「別冊 Newton の重症者数」を使って得たものです。この値を使えば死亡者数が「42万」や「41.8万」になることは分かっているので、本項は論理的には不要な蛇足であり、確認のための計算です。

計算方法 計算結果
(8a) github の SIR 方程式とパラメータ
による私の数値計算
世代 a の感染者総数
=67100300人
(8b) github の SIR 方程式とパラメータ
による私の数値計算
世代 e の感染者総数
=32603300人
(8c) 世代 a の感染者総数(67100300人)と、
世代 a の感染者重症化率(0.00300000)の掛算
世代 a の重症者数
=201301人
(8d) 世代 e の感染者総数(32603300人)と、
世代 e の感染者重症化率(0.0200000)の掛算
世代 e の重症者数
=652066人
(8e) 世代 a の重症者数(201301人)と、
世代 e の重症者数(652066人)
の足し算
全体の重症者数
=853367人
(8f) 全体の重症者数(853367人)と、
重症者死亡率(49%)の掛算
全体の死亡者数
=418150人

(8f)で得られた 418150人 を四捨五入すると、「42万」や「41.8万」になります。

(9)結論

西浦氏が githubNewsweek 2020-06-09 号記事([*1-3])で示したパラメータの一部に約2%の誤りがあり、計算結果の死亡者数も「42万」「41.8万」から約2%ほど異なった値になります。

「西浦氏が github など [*1-3] で示しているパラメータで計算すると、死亡者数が「42万」「41.8万」にならないこと」
「本エントリ(6)で新しく計算したパラメータで計算すると、死亡者数が「42万」「41.8万」になること」
を示すプログラム
https://colab.research.google.com/drive/1ZvY-HmOvh3AS7OnjJdFh0lBw_OfTX5cU?usp=sharing
を作成しました。(「新しいパラメータ」と書きましたが、これは西浦氏が「42万死亡推計」を計算した際に使ったはずのパラメータです。)

プログラムを検討しなくても、新しく計算したパラメータの方が正しいと判断できる傍証があるので、(6)に示しました。また (4)で示した別エントリ は、そもそも(私のプログラムによる数値計算と関係なく)西浦氏が示しているパラメータに不自然な点があることを示していて、これは西浦氏が示しているパラメータに誤りがあることの傍証になります。

(10)「42万死亡推計」に関する西浦氏への質問

「42万死亡推計」に関して、西浦氏に質問したい事項を列挙します。この中には、私が既に指摘してきた事項や、示してきた数値が含まれています。

  • (10a) github 資料 に示した「計算コード」での検算を行ったのか。
    (行っていれば、本エントリが指摘している2%のパラメータの誤りに西浦氏が気付かないはずはありません。)
  • (10b) 「計算コード」での計算は、何日目で計算を停止したのか。
    (西浦氏は、横軸100日目までのグラフを示しています。しかし「計算を停止した日付が100日」なのか「計算はもう少し長い期間行ったが、右の方はグラフとして示す意味に乏しいので100日までのグラフとしたのか」は明示されていません。私は念の為として計算停止を200日としています。)
  • (10c) 世代 a、世代 e の、感染者総数(計算期間中に感染したことのある人の数)は何人か。
    (この基本的な数値を西浦氏は明示していません。この関連で西浦氏が示しているのは、全世代での感染者総数の有効桁数2桁の概数だけです。重症者数を [*1-5] で有効数字6桁で示しているのですから、こちらも同程度の精度で、世代毎に示すべきです。この正確な値が分かると、私が 別エントリ で検討した不自然なパラメータの問題が明らかになります。)
  • (10d) 世代 a、世代 e の感染者重症化率は幾つか。世代 a、世代 e の感染者死亡率は github 資料によれば、0.0015、0.1000 だが、これは正しいのか。
    ((6)に示したように、それぞれ、0.003、0.02、0.00147、0.0098 であるはずです。)
  • (10e) 2020-03-19 に、新型コロナの基本再生産数  R_0 を、値が分からないとしながらも唯一の値として 2.5 を示した、その時点における理由は何か。
    (おそらく西浦氏に説明するに足る理由はありません。あるなら既に説明しているはずです。西浦氏は幾つか  R_0=2.5 の理由のような話をしていますが、それらは事後的な説明であり、2020-03-19 時点における理由の説明はありません。)
  • (10f) 2020年春の日本の新型コロナの  R_0 を値は現在も分からないのではないか。この点の確認。
  • (10g) 現在でも 2020年春の日本の新型コロナの  R_0 の値が分からないのなら、今後の感染症においても同様の問題が発生すると思われるが、この点の確認。
    (今後の感染症においても、基本再生産数の値が分からないまま感染症対策を行う可能性があります。分からないなら「 R_0 は分からない」と分かりやすく明示した上で計算などを示すべきです。)

質問したい事項はまだあるのですが、別のエントリで議論する事項に関係するのでそちらを書いてからにします。

(11) 以前に示していた同主旨のエントリからの変更点

以前に示していた同主旨のエントリ から本エントリへの変更点を説明します。

まず共通点を述べます。 本エントリの(6)[表6-1] は、以前のエントリの(4)の表1 と実質的に同等です。世代 a、世代 e の感染者重症化率が共通していることになります。
本エントリと以前のエントリで、別冊 Newton や毎日記事 [*3-1] における重症者数が正しいと考えていることも共通しています。

次に変化した点を述べます。 変化したのは、重症者死亡率(約49%)と、死亡者数(約41.8万人)についてです。

以前と今回の違いを表にします。

記事などの記述 以前のエントリ 今回のエントリ
重症者死亡率 約49% 48.9824% 49%
死亡者数 約41.8万人 418000人 418150人

[表11-1] 以前と今回のエントリにおける違い

この付近には、次の式

 (重症者数)×(重症者死亡率)=(死亡者数) ...... (式11-1)

が成立していて、この式に、
(重症者数)=853367人
(重症者死亡率)=約49%
(死亡者数)=約41.8万
の数値を入れると、

 853367×(49%)=(41.8万) ...... (式11-2)

となります。(式11-2) の左辺は 418150 なので、(式11-2) は厳密には成立しません。

以前の私は、この式の左辺が丁度 41.8万 になるように「約49%」の (重症者死亡率) の値を決定し、(重症者死亡率)=48.9824% としていました。しかしこの値は、(式11-2) を成立させるとしても、不自然です。本エントリで、西浦氏の感染者重症化率が不自然だと批判しているのと同様に不自然です。したがって不適切です。

よって今回私は、(式11-2) の左辺は  853367×(49%) のままとし、式の値は 418150 なので、右辺の41.8万にぴったり一致はしないが、右辺は本来は「約41.8万」なのだからぴったり一致しなくても問題はない、と考え方を変えました。

本エントリ冒頭に書いた「議論で出てくる数値の一部が約0.036%変化しました」は、以下を指しています。

重症者死亡率:
48.9824%から49%に変化。変化率は  (49/48.9824 - 1) × 100 = 0.036%
死亡者数:
418000人から、418150人に変化。変化率は  (418150 / 418000 - 1) × 100 = 0.036%

このように議論に出てくる数値の一部は約0.036%変化しましたが、以前と今回で論旨に変化はありません。

更新履歴

  • 2024-03-18
    公開。本エントリに関係する説明を twitter に書くかも知れません。書いたら冒頭付近にリンクを書いておきます。


西浦氏の Scientific Reports 論文について(6)論文の実効再生産数と活動度の一部期間での一致


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



目次


(1)序論

西浦氏らのグループは、2023-10 に Scientific Reports にコロナワクチンに関する論文を発表しました。
後に「コロナワクチンで死者9割以上減 京都大チームが推計」と共同が報じた論文です。(以下、西浦論文、と書きます。)
Evaluating the COVID-19 vaccination program in Japan, 2021 using the counterfactual reproduction number

本稿は、西浦論文 Figure 2B と、西浦論文 式(4)  h_t の一部とが、ある期間にはほとんど一致すること、その後の両者の上下関係、これらについての大筋での検討を行います。

西浦論文については、複数のエントリで書いています。

(1-1) 図6.1と図6.2

図6.1 (a) 西浦論文の実効再生産数 Figure 2B。(b) 西浦論文の活動度 式(4) h_t。(c) (a)と(b)を重ね合わせたグラフ。
図6.1 (a) は、西浦論文の実効再生産数 Figure 2B。(b) は、西浦論文の活動度 式(4)  h_t の一部(reporting coverage=0.25)。(c) は (a) と (b) との重ね合わせ。(c) の左方で (a) と (b) はほとんど一致している。この期間の実効再生産数の変動が、活動度の変動でほとんど説明できることを示している。

図6.1(c) の左部分で、実効再生産数が活動度とほとんど一致しているのが目を引きます。この図6.1(c) に説明を加筆したのが図6.2です。

図6.2 西浦論文の実効再生産数 Figure 2B と、西浦論文の活動度 式(4) h_t との重ね合わせとその説明
図6.2 西浦論文の実効再生産数 Figure 2B と、西浦論文の活動度 式(4)  h_t の一部 との重ね合わせとその説明

図6.2の期間(1)では、実効再生産数  R_t が活動度  h_t とほとんど一致しています。この期間の実効再生産数の変動が、活動度の変動でほとんど説明できることを示しています。期間(2)では、実効再生産数が活動度より大きくなります。活動度以外の感染拡大要因があったことになりますが、これはデルタ株の比率増加によると思われます。期間(3)では、実効再生産数が活動度より小さくなります。デルタ株による感染拡大効果を上回る感染抑制要因があったことになりますが、これは主にワクチンによる免役だと思われます。ただし、論文からワクチンの効果を定量的に見積もるのは困難です。

(2) では図6.1(c)を作図した意図を説明します。(3) では図6.1(c)の作図方法を説明します。(4) では、図6.2の期間(1)~期間(3)などについて分かる範囲で説明します。(5) では本稿の結論を述べます。(6) は、補足事項です。

(1-2) 本稿の考察の限界

期間(1)で実効再生産数  R_t と活動度  h_t がほとんど一致するのは、(4) に示すように少し考えれば意外なことではなく、むしろ当たり前です。しかし私はこの一致を興味深いと思ったので本エントリを書いてご説明しています。この一致をご説明できるのは、西浦論文 式(3)の重要パラメータの値が期間(1)ではかなり特定できる、という特殊な事情があるからです。

(4) では、図6.2の期間(2)や期間(3)についても検討しますが、これらの期間について言えることはかなり限定されます。というのも、これら期間での西浦論文 式(3)の重要パラメータの値は、大雑把にしか把握できないからです。

西浦論文の重要なパラメータである実効再生産数は、論文 式(3)によって他のパラメータから計算されます。本稿は、これらパラメータ同士がどういう関係にあるのかを検討します。検討は入手しやすい情報のみに基づいていて、前提とする西浦論文のパラメータの妥当性は検討しない上、一部のパラメータは論文での詳細開示がないことから、検討は限定的なものになります。

また論文が対象とする期間の一部(図6.2の期間(1))では定量的な議論ができていますが、他の期間で行っているのはできているとしても定性的な議論です。例えば、図6.2で期間(2)と期間(3)の境界の、縦線(C)の位置を他の資料から見積もり、論文 Figure 2B における位置と照合することはできていません。

こうした事情があるので、本稿が西浦論文の計算の一般的な妥当性を支えることはありません。

(1-3) 比較対象となる実効再生産数は「描かれていない」

なお、図6.1(c)図6.2 を私が作図した主たる意図は、青線で示した活動度と実効再生産数の比較です。ところが本来比較すべき実効再生産数は "Early schedule" と "Late schedule" の中間にある「描かれていない曲線」です。この事情は (6-1) に書きました。

(2)図6.1(c)を作図した意図

私は 別稿 で、西浦論文 式(4)の  h_tグラフ を示し、さらにその一部を 西浦論文 Figure 1重ね合わせたグラフ を示しました。この重ね合わせグラフを示したのは、 h_t の細かな上下が Figure 1 の細かな上下に反映されているのを示すためでした。

この重ね合わせも興味深いものでしたが、 h_t を重ね合わせるグラフは、西浦論文の Figure 1 ではなく、Figure 2B にすべきでした。なぜなら、 h_t は活動度なので、その影響が直接及ぶのは新規感染者数(やその予測値)を描いた Figure 1 ではなく、実効再生産数を描いた Figure 2B だからです。

(3)図6.1(c)の作図方法

図6.1(c) は、以下の要領で作図しました。

  • (3a) まず横方向に、2枚のグラフの日付を合わせました。
  • (3b) 次に縦方向に、Figure 2B の左方のできるだけ広い範囲に論文 式(4)  h_t の一部が重なるように、 h_t を上下に移動し、上下に伸縮しました。

(3b) でパラメータ  h_t を上下に移動し、上下に伸縮させたパラメータを、 h_t^{'} とします。すると両者の間には以下の関係があります。

 h_t^{'} = w_1 h_t + w_2 ...... 式(1)

 w_1,\ w_2 は定数です。( h_t を上下に移動することは、 w_2 を加算することに対応しています。 h_t を上下に伸縮することは、 w_1 を掛算することに対応しています。)

これまで、図6.1(c) や図6.2 で実効再生産数に重ね合わせた活動度を、 h_t と書いてきましたが、このように作図したのですから、正確に言うと重ね合わせたのは( h_t ではなく) h_t^{'} です。

なお重ねる範囲として Figure 2B の左方を選んだのは、この時期の実効再生産数が活動度  h_t^{'} にほとんど一致することが、計算モデル上、予想されるからです。詳細は (4-1-1) に。

(4)図6.2の期間(1)~期間(3)などについての検討

図6.2 を再掲します。

図6.2 [西浦論文の実効再生産数 Figure 2B](https://github.com/sarkov28/storage/blob/master/2023-10_Nishiura_Scientific_Reports/Figure_2B.png?raw=true) と、西浦論文の活動度 式(4) h_t との重ね合わせとその説明
図6.2 西浦論文の実効再生産数 Figure 2B と、西浦論文の活動度 式(4)  h_t の一部 との重ね合わせとその説明

図6.2で、縦線(A)から(B)の期間が期間(1)から期間(2)への移行期で、縦線(C)が期間(2)から期間(3)への移行期です。

以下で実効再生産数  R_t と活動度  h_t の関係を検討する際、西浦論文 式(3)

 \displaystyle R_{ab, t} = \left(1 - l_{a, t} - \frac{\sum_{k=1}^{t-1}{i_{a, k}}}{n_a}\right) K_{ab}\ p\ h_t\ d_t\ c_t ...... 論文 式(3)

を以下のように読み替えた式を使います。

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) K\ p\ h_t\ r\ u_t\ c_t ...... 式(2)

論文 式(3)から年齢層の要素を除いたのが式(2)です。ここで  R_t は全体の実効再生産数、 l_t は全体のワクチン免疫比率、 i_k は全体の新規感染者数、 n は全人口、(したがって、 \sum_{k=1}^{t-1}{i_k}/n は全体の既感染者割合)、 K は次世代を定める定数です。式(2)は年齢層の要素を除いているので厳密には成立しませんが、大筋では成立しています。

(4-1) 実効再生産数の動きの説明がある程度可能な前半

(4-1-1) 実効再生産数が活動度にほとんど一致する、期間(1)

図6.2 の期間(1)は、2021-03-04~2021-05-10 ごろの期間です。期間(1)では、実効再生産数  R_t が青線の活動度  h_t にほとんど一致しています。この期間の式(2)

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) K\ p\ h_t\ r\ u_t\ c_t ...... 式(2)

では、以下のように見なすことができます。

これらから式(2)は、

 R_t = K\ p\ r\ h_t ...... 式(3)

となります。(式(3)は、期間(1)でのみ成立します。)

期間(1)では、 R_t h_t の間に式(3)のシンプルな関係があることが分かりました。(3) の方法で  h_t を変形した  h_t^{'} R_t にほとんど一致する 図6.1図6.2 が得られたのは、 R_t h_t にこうした関係があったからです。

なおこの時期、Figure 2B の4つのシナリオの実効再生産数もほとんど一致しています。この時期にはワクチン接種者がほとんどいないので、異なる接種状況を想定する各シナリオが、実効再生産数に影響しないのです。

(4-1-2) 実効再生産数と活動度の、よりシンプルな関係式

式(3)( R_t = K\ p\ r\ h_t)であり、式(1)( h_t^{'} = w_1 h_t + w_2)であり、期間(1) ではほとんど  R_t = h_t^{'} ですから、概ね、

 h_t^{'} = K\ p\ r\ h_t = w_1 h_t + w_2 ...... 式(4)

と書くことができます。式(4)が複数の  t で成立することから、

 (w_1,\ w_2) = (K\ p\ r,\ 0) ...... 式(5)

を導けます。これは期間(1) の状況から得た値ですが、式(1)( h_t^{'} = w_1 h_t + w_2)のパラメータ  w_1,\ w_2 は研究期間全体で一定(=そのように作図している)なので、式(5)は研究期間全体で成立します。

これを使うと式(1)は、

 h_t^{'} = K\ p\ r\ h_t ...... 式(6)

と書くことができ、式(6)を用いると、式(2)

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) K\ p\ h_t\ r\ u_t\ c_t ...... 式(2)

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) h_t^{'} u_t\ c_t ...... 式(7)

となります。式(7)は、式(2)よりも見通しが良くなっています。(式(7)は研究期間全体で大筋で成立します。)

本項での検討を言葉で書くと、

私が図6.1(c)で行った作図が満たす式(6)を、図6.2 期間(1)でグラフがほとんど一致していることなどから導き、式(6)を式(2)に代入し、式(7)を得た、

となります。

(4-1-3) 期間(1)から期間(2)への移行期の、縦線(A)から縦線(B)

図6.2 の期間(1)(2021-03-04~2021-05-10 ごろ)の間、実効再生産数はほとんど活動度のみによって上下します。しかし、4月末の縦線(A)のころから、図6.2 の実効再生産数のピンク線が、活動度の青線の僅かに上で推移するようになり、ピンク線の幅程度の差が生じています。

これは4月末からのゴールデンウィークが連休として処理されていることによると思われます。ただし、西浦論文 補足情報 Table S1 のパラメータ  e にある通り、連休による感染拡大の効果は1日あたり2%以下なので、それほど大きくはありません。

(4-1-4) 実効再生産数が活動度を上回る期間(2)の前半

図6.2 の 2021-05-20 ごろの縦線(B) を過ぎたところから、実効再生産数のピンク線は、活動度の青線をはっきり上回るようになります。このころ、上で求めた式(7)

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) h_t^{'} u_t\ c_t ...... 式(7)

の各パラメータは、

という値なので、およそ、

 \displaystyle R_t = h_t^{'} u_t ...... 式(8)

となり、増加中のデルタ株比率 を反映した  u_t の分だけ、 R_t h_t^{'} より大きくなります。(デルタ株比率は、0~1ですが、パラメータ  u_t は 1~1.5 だと思われます。別稿(2-11) に関連があります。)

この間、(感染を抑制する)ワクチン免疫比率も上昇しますが、(感染を拡大する)デルタ株比率の上昇が先行するため、実効再生産数  R_t は活動度  h_t^{'} より上を推移することになります。

実効再生産数の動きをある程度説明できるのは、この7月末ごろまでです。

(4-2) 実効再生産数の動きの説明が困難な後半

デルタ株比率のパラメータ  u_t は、8月初めごろから上昇が鈍くなり始め、8月末には明らかに鈍くなります。したがってこの時期の実効再生産数と活動度の関係への  u_t の寄与は小さくなり、Figure S2 の年齢層別のワクチン免疫比率 の寄与が大きくなってきます。

ワクチン免疫比率は、8月に入ると年齢層により上昇するものと下降するものが混ざるようになります。ただしこの時期でも、下降する年齢層は70歳以上に限られること、下降の速度は上昇する年齢層の速度を上回ることから、全年齢層でのワクチン免疫比率  l_t は上昇しているでしょう。ただし全年齢層の免疫比率がどの程度の水準になっているかは見積もりにくいです。このため、 l_t が重要な役割果たす式(7)

 \displaystyle R_t = \left(1 - l_t - \frac{\sum_{k=1}^{t-1}{i_k}}{n}\right) h_t^{'} u_t\ c_t ...... 式(7)

の様子も見積もりにくくなります。したがって、実効再生産数  R_t と活動度  h_t^{'} の関係をこれまでの時期のように説明することは難しくなります。

図6.2 R_t h_t^{'} との関係からは、

  • (4-2a) 図6.2 期間(2)の後半(7月末ごろ~8月半ばごろ)は、 R_t h_t^{'} が近づいています。したがって、 u_t による感染拡大効果を  l_t による感染抑制効果(少し細かく書くと、 (1 - l_t) による効果)が打ち消しつつあるはずです。
  • (4-2b) 図6.2 縦線(C)(8月半ばごろ)付近では  R_t h_t^{'} が交差しています。この時期は、 u_t の効果と  h_t^{'} の効果が「釣り合う」状態になっているはずです。
  • (4-2c) 図6.2 期間(3)(8月半ばごろ~)は、 R_t h_t^{'} を下回り、さらにその差が徐々に拡大しています。したがって、 u_t の効果を、 l_t の効果が上回り、その差を拡大しているはずです。

図6.2 から得られる(4-2a)~(4-2c)の「・・・はずです」との推測を、 u_t を示している Figure S5 や、 l_t の年齢層別の値である Figure S2 の動きと照合することもある程度は行えるかも知れませんが、本稿での検討はここまでにしておきます。

(5)結論

西浦論文の実効再生産数のグラフ Figure 2B に、論文 式(4)の活動度  h_t の グラフを重ね合わせた 図6.1 と、説明を加筆した 図6.2 を作成しました。

図6.2 期間(1) では、活動度は実効再生産数にほとんど一致し、実効再生産数の動きを活動度だけでほとんど説明できることが分かりました。

また、図6.2 期間(2) の途中までは両者の関係をある程度説明できました。その後は、特にワクチン接種による感染予防効果がどの程度なのかが見積もりにくいため、説明が困難になりました。図6.2 期間(2) はデルタ株比率の上昇による感染拡大効果が強い時期だと思われ、期間(3) はデルタ株の効果をワクチンの効果が上回る時期だと思われるのですが、このことを論文のパラメータの動きと照合して説明することはほとんどできませんでした。

(6)補足事項

(6-1) 活動度と比較すべき実効再生産数について

西浦論文の Figure 2B は、4つのワクチン接種シナリオ(早いスケジュール "Early schedule"、遅いスケジュール "Late schedule"、高い接種率 "Elevated coverage"、ワクチンなし "Without vaccination")別の、実効再生産数です。

図6.3 西浦論文の Figure 2B
図6.3 西浦論文の Figure 2B

本来、論文 式(4)  h_t と比較すべきなのは、これら4つのシナリオの実効再生産数ではなく現実と同じスケジュールで接種した場合の実効再生産数ですが、それは論文中に示されていません。

比較すべき「現実と同じスケジュールで接種した場合の実効再生産数」の線は、Figure 2B の「早いスケジュール "Early schedule"」(緑線)と「遅いスケジュール "Late schedule"」(灰色線)との中間付近にあるはずです。この線は Figure 2B に描かれていないので、分かりにくく、比較しにくいです。

(6-2) 重ね合わせに使った西浦論文 式(4) h_t について

私が 別稿 で示した、西浦論文 式(4)の  h_t はこちらです。

図4 西浦論文 式(4)の h_t
図4 西浦論文 式(4)の  h_t

別稿では、西浦論文 Figure 1 に  h_t を重ね合わせていました。このため、使ったのは図4のうち報告率 reporting coverage=1 のグラフであり、横軸範囲は 2021-02-17 から 2021-11-30 でした。

本稿では、西浦論文 Figure 2B に重ね合わせるので、図4のうち 報告率 0.25 のグラフ を使い、横軸範囲は 2021-03-01 から 2021-11-30 にしました。グラフの曲線の色は、Figure 2B の他の線と区別しやすい色にしました。

グラフを作成したデータは、こちらの csv 形式のデータ の、報告率 reporting coverage=0.25 の、期間 2021-03-01~2021-11-30 です。重ね合わせの作図方法は、(3) で説明しました。

(99)更新履歴

  • 2024-01-23
    公開
  • 2024-01-30
    (1-2) から (1-3) を分離しました。
    (1)の記述の一部を (1-2) に移動すると共に、(1-2) に加筆しました。


西浦氏の Scientific Reports 論文について(5)論文 Figure 1 を構成しているパラメータの検討


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



目次


(1)序論

西浦氏らのグループは、2023-10 に Scientific Reports にコロナワクチンに関する論文を発表しました。
後に「コロナワクチンで死者9割以上減 京都大チームが推計」と共同が報じた論文です。
本稿では Scientific Reports 論文を西浦論文と書きます。
Evaluating the COVID-19 vaccination program in Japan, 2021 using the counterfactual reproduction number

図5.1 西浦論文の Figure 1
図5.1 西浦論文の Figure 1

本稿は、西浦論文 Figure 1 を構成している主なパラメータについて検討します。

Figure 1 では、予測値が実測値に概ね一致して見えます。これに対し、「一致しているのは、論文執筆者が裁量でパラメータの値を調整して「作った」からなのでは」との疑問があります。このため、以下では、「論文執筆者の裁量が入る余地がどの程度あるか」に注意しながら各パラメータを検討します。

西浦論文については、複数のエントリで書いています。

(2) Figure 1 を構成しているパラメータに論文執筆者の裁量余地はあるか

まず (2-01) に Figure 1 を構成しているパラメータを挙げ、その後の各項で各パラメータについて検討していきます。 この際、各パラメータに、西浦論文 Figure 1 をフィットさせるための論文執筆者の裁量が入る余地がどの程度あるかを見ていきます。

(2-01) パラメータのリスト

  •  l_{a,t}: 年齢層  a の時刻  t におけるワクチンによる感染予防効果。
  •  v_{a,t}: 年齢層  a の時刻  t におけるワクチン接種者数。
  •  h_t(論文 式(1)): 接種から時間  t が経過した時のワクチンの感染予防効果。
  •  i_{a,t}: 年齢層  a の時刻  t における新規感染者数。
  •  R_{ab,t}: 時刻  t において、一人の年齢層  b から生じる年齢層  a の新規感染者の平均数。実効再生産数。
  •  g_t: 感染から時間  t が経過した時の世代時間の確率密度関数
  •  K_{ab}: 年齢層  b から年齢層  a への次世代行列(の要素)。
  •  p: 調整パラメータ。
  •  h_t(論文 式(3)): 時刻  t における活動度。
  •  d_t: 時刻  t におけるデルタ株による感染性の増強効果。
  •  c_t: 時刻  t における休日による感染性の増強効果。

このうち、 h_t に関しては、不適切なことに論文 式(1)と論文 式 (3) で同じ文字が別の意味に使われています。

(以下の各項の見出しで、本来なら「 l_{a,t}」などと書いてあるべきところが「l_{a,t}」のようになっているのは、見出しに数式表現を使うと目次で期待通りに表示されなかったためです。もし対策をご存じでしたら、https://twitter.com/sarkov28 までお知らせ頂けると助かります。)

(2-02) l_{a,t}: 年齢層 a の時刻 t におけるワクチンによる感染予防効果

西浦論文は、 l_{a,t} で、年齢層  a の時刻  t におけるワクチンによる感染予防効果を組み入れています。このパラメータは、

  • (2-03)  v_{a,t} 年齢層  a の時刻  t におけるワクチン接種者数
  • (2-04)  h_t(論文 式(1)) 接種から時間  t が経過した時のワクチンの感染予防効果

から計算されています。(2-03)(2-04)の検討を踏まえると、 l_{a,t} に論文執筆者の裁量が入る余地は、あってもさほど大きくないのではと思います。

(2-03) v_{a,t}: 年齢層 a の時刻 t におけるワクチン接種者数

 v_{a,t} は、各年齢層毎の実際のワクチン接種数に依っていると説明されているので、ここに論文執筆者の裁量が入る余地はほとんどないと思われます。

(2-04) h_t(論文 式(1)): 接種から時間 t が経過した時のワクチンの感染予防効果

 h_t は接種者本人への感染予防効果を示すパラメータであり、西浦論文はこれを参考文献 53, 54, 19 から得たと説明しています。しかし西浦論文が  h_t をどのようにモデル化したのかの詳細は示されていません。数式もグラフも示されていないので、これはやや不可解です。

上述のように  v_{a,t} と本項  h_t から  l_{a,t} を計算しているのですが、このうち、

ところが、(2-04) の  h_t についてはなにもありません。(西浦論文 参考文献 19 の補足情報 Figure S6 が  h_t のグラフなのかも知れません。)
この説明欠落はやや気になるのですが、一般的に考えれば、 h_t に Figure 1 をフィットさせるための論文執筆者の裁量が入る余地は、あってもさほど大きくないと思われます。なぜなら、このパラメータ  h_t t が、 R_{ab,t} t とは異なる時間軸だからです。前者は接種からの時間であり、後者はカレンダー時刻です。また、 h_t には少なくとも、接種後しばらくは単調増加し、その後は単調減少するとの制約があるからです。

(2-05) i_{a,t}: 年齢層 a の時刻 t における新規感染者数

 i_{a,t} は、実測された新規感染者数に関する情報を感染時点に換算することで得ているようです。 この換算に(確率分布のパラメータのある程度の選択など)多少の裁量はあり得ますが、Figure 1 の形状を調整するような細かな裁量が入ることはないと思われます。

(2-06) R_{ab,t}: 時刻 t において、一人の年齢層の感染者 b から生じる年齢層 a の新規感染者の平均数。実効再生産数

 R_{ab,t} は、これまで述べてきたパラメータや、(2-07)~(2-12) のパラメータから、論文 式(3)で計算されています。

 \displaystyle R_{ab, t} = \left(1 - l_{a, t} - \frac{\sum_{k=1}^{t-1}{i_{a, k}}}{n_a}\right) K_{ab}\ p\ h_t\ d_t\ c_t ...... 論文 式(3)

各パラメータを検討したところ、私は、 R_{ab,t} に論文執筆者の裁量が入る余地は、あってもさほど大きくないと思います。

(2-07) g_t: 感染から時間 t が経過した時の世代時間の確率密度関数

 g_t は、論文 式(2)に出てくる変数です。

 \displaystyle i_{a,t} = \sum_{b=1}^{10}\ \sum_{τ=1}^{t-1}\ R_{ab,t}\ i_{b,t-τ}\ g_{τ} ...... 論文 式(2)

このパラメータで Figure 1 の形状を調整できるとは私には思えません。パラメータ  g_t t が、 R_{ab,t} t とは異なる時間軸だからです。前者は感染からの時間であり、後者はカレンダー時刻です。

(2-08) K_{ab}: 年齢層 b から年齢層 a への次世代行列(の要素)

 K_{ab} について、私はあまり調べていません。 K_{ab} はカレンダー時刻で変化するパラメータではないため、このパラメータで Figure 1 の形状を調整をすることは困難ではないかと思います。

(2-09) p: 調整パラメータ

パラメータ  p は、正に調整のためにあるのですが、これはカレンダー時刻では変化しない、一定の値を取るパラメータです。(一定であることは、西浦論文 補足情報 Table S1 からも分かります。)したがってこのパラメータでカレンダー時刻での変化を調整することは困難だと思われます。
 p については、別稿(2)(3)(4) でも触れています。

(2-10) h_t(論文 式(3)): 時刻 t における活動度

パラメータ  h_t については、別稿(3) で検討しています。他に、別稿(2)(4) でも触れています。

論文 式(3)の  h_t は、論文 式(4)

 h_t = ω^{community}\ α_t^{community} + ω^{house}\ α_t^{house} + ω^{work}\ α_t^{work} ...... 論文 式(4)

で計算されるのですが、このうち3つの  αgoogle mobility データを加工したものです。(別稿で再現しました。図4 と図5。)この加工には、(6つある google mobility の指標のうちどれを使うかなどに、)ある程度の裁量があるかも知れませんが、これは Figure 1 の形状を調整できるようなものではありません。また3つの  ω のうち1つの値は定義により決まっていて、残り2つの値は尤度関数(論文 式(8))を用いて最尤推定されます。したがってこの  h_t に裁量が入る余地は一般的にはないと考えられます。

(2-11) d_t: 時刻 t におけるデルタ株による感染性の増強効果

パラメータ  d_t は、 d_t = r\ u_t だと論文に説明されています。
このうち  r は、論文 式(8)を用いて最尤推定されるパラメータなので、ここに裁量が入る余地は一般的にはないと考えられます。
 u_t は市中におけるデルタ株の比率を示す指標です。西浦論文における  u_t の説明は細部が不明ですが、西浦論文 補足情報 Figure S5 は、実測の「 \circ」に沿うようにロジスティック曲線をフィットさせた説明されています。
図5.2 西浦論文 補足情報 Figure S5
図5.2 西浦論文 補足情報 Figure S5(クリックすると拡大します)

このフィットには図にあるように多少ずれてますが、実測に縛られているので、ここでの裁量はあっても小さいと思われます。
論文によると、Figure S5 を 1~1.5 になるよう調整して  u_t とし、さらにパラメータ  r を掛算して、 d_t を得ているようです。この処理ならば、パラメータ  d_t に Figure 1 の形状を調整する裁量が入る余地があるようには見えません。

なお、パラメータ  r 付近に無駄なパラメータがあるという件を 別稿(2) で、関連を別稿 (3)(4) で検討しています。

(2-12) c_t: 時刻 t における休日による感染性の増強効果

西浦論文のパラメータ  c_t の説明は分かりにくいのですが、このパラメータは以下の値を取ります。

 c_t = e\ \ (if\ \  t = 対象日)
 \ \ \ = 1\ \  (else)

このうちパラメータ  e は論文 式(8)を用いて最尤推定されます。したがってパラメータ  c_t に裁量が入る余地はないと思われます。

(3) 結論

西浦論文 Figure 1 を構成している主なパラメータについて、特に論文執筆者が Figure 1 の形状を調整するような裁量があるかを検討しました。各パラメータの裁量の余地はないか、あっても小さいと思われます。

検討は項目ごとに行ったので、複数のパラメータを組み合わせによる、Figure 1 の形状の調整については検討していません。しかし、一般論になりますが、大きな裁量のない複数のパラメータを組み合わせて、Figure 1 の形状の調整を行うことは、相当に困難だと思われます。

(99)更新履歴

  • 2024-01-23
    公開。


西浦氏の Scientific Reports 論文について(4)Table S1 のパラメータは再計算されるべき


私がネット上でしていることの まとめ は、こちらに。
https://sarkov28.hatenablog.com/entry/2022/03/29/160915



目次


(1)序論

西浦氏らのグループは、2023-10 に Scientific Reports にコロナワクチンに関する論文を発表しました。
後に「コロナワクチンで死者9割以上減 京都大チームが推計」と共同が報じた論文です。(以下、西浦論文、と書きます。)
Evaluating the COVID-19 vaccination program in Japan, 2021 using the counterfactual reproduction number

図4.1 西浦論文の Figure 1
図4.1 西浦論文の Figure 1

本稿で Figure 1 について直接言及することはありません。
西浦論文については複数のエントリで書いています。

本稿では、このうち「論文について(2)」と「論文について(3)」の検討を踏まえ、西浦論文の補足情報 Table S1 に記載されているパラメータは再計算されるべきであることを示します。

(2)西浦論文の補足情報 Table S1 のパラメータは再計算されるべき

西浦論文の補足情報 Table S1 に示されている( e 以外の)パラメータは、再計算されるべきです。

パラメータのうち、 p,\ r が再計算されるべきなのは、この2つのパラメータが、実質的に1つのパラメータ  p×r だからです。2つのパラメータとしてではなく、1つのパラメータとして再計算されるべきです。
この件については、別稿である 論文について(2)論文 式(3)、式(8)の無駄なパラメータ で検討しました。ここではこの件を、西浦論文に示されている式の定義から導出しました。また、西浦論文が計算したパラメータが、別稿で導入した1つのパラメータ  p×r と比べて大きくバラついていることを示しました。

Table S1 のパラメータ、 ω^{house},\ ω^{work} が再計算されるべきなのは、この2つが報告率によらない(=全ての報告率に共通の)値であるべきだからです。現状の  ω^{house},\ ω^{work} は報告率で変化しています。この2つが報告率によらないことを前提とした尤度関数を構成し、再計算すべきです。
この件については、別稿である 論文について(3)活動度は報告率によって変化すべきではない で検討しました。ここでは、 ω^{house},\ ω^{work} が報告率によって変化すべきでない理由も示しました。また、西浦論文 Figure 1 と、計算した活動度  h_t とを重ね合わせたグラフを参考までに示しました。

なお、 ω^{house},\ ω^{work} を全ての報告率に共通の値として計算すると、一般には、西浦論文 Figure 1 に示されているフィッティングの状況は、現状より悪化すると思われます。しかしそれはむしろ、西浦論文のモデルの本来の姿だと思われます。現状からの悪化が観察されるなら、その差はむしろいわゆるオーバーフィッティングではないでしょうか。

(99)更新履歴

  • 2023-12-30
    公開。
  • 2024-01-04
    別稿「西浦氏の Scientific Reports 論文について(2)」に、「更新履歴 2024-01-04」の大幅な修正をしました。これに合わせて本稿のタイトルを修正し、本文の内容も変更しました。修正前の本稿のタイトルは「西浦氏の Scientific Reports 論文について(4)Table S1 のパラメータは、全て再計算されるべき」でした。