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

西浦氏の Scientific Reports 論文について（３）活動度は報告率によって変化すべきではない

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


目次

（１）序論

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

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

本稿で Figure 1 について直接言及することはありませんが、補足的事項で、計算した活動度が Figure 1 に与える影響を見るために、重ね合わせたグラフを示しています。

西浦論文は、幾つかの要素から構成されています。
論文 式(4)で計算される活動度は一つの要素です。 また、reporting coverage 報告率も重要な要素であり、西浦論文はこの値を複数用意し、それぞれについて計算を示しています。 他にも論文を構成する要素はあるのですが、本稿が注目するのは今挙げた「活動度」と「報告率」の2つです。

本稿はまず、西浦論文が、

  • 活動度と報告率を独立したものとして扱うべきなのに、
  • 独立したものとして扱っていない、

ということを示します。この点は、明確に西浦論文に示されています。

本稿は次に、西浦論文において活動度と報告率が独立しているべき理由を示します。

ただし、活動度と報告率が独立していないことが、西浦論文の結果に定量的に大きな影響を与えているかどうかは不明です。本稿が述べるのは、（定量的に大きな問題の指摘というより、）定性的に、モデル構成のあるべき姿として問題ではないか、という話です。

また、報告率が論文 式(4)で計算している活動度 h_t に与えている影響をグラフで示します。

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

（２）活動度と報告率の関係 ─ あるべき構成と西浦論文の構成

西浦論文のモデルは、多数の要素から構成されています。

論文で、活動度は h_t と表記されています。式(4)で計算され、式(3)で用いられています。 h_t は、「google mobility の6指標のうち3指標と、その3指標を調整するパラメータから計算される数値」です。

本来、西浦論文のモデルは図3.2の構成になっているべきです。図3.2のポイントは、活動度と報告率とが独立していること、です。ところが西浦論文のモデルは、図3.3の構成になっています。西浦論文のモデルは、図3.3に示したように、活動度が報告率から独立していません。

図3.2　西浦論文のモデルのあるべき構成
図3.2　西浦論文のモデルのあるべき構成

図3.3　西浦論文のモデル構成
図3.3　西浦論文のモデル構成

（３）西浦論文の活動度が報告率から独立していないことの根拠

西浦論文の活動度が報告率から独立していないことは、論文に明記されています。
補足情報 の Table S1 には、4つの報告率毎の ω_h,\ ω_w が記載されています。（西浦論文の活動度に関する表記には不統一があるので、別項で整理しています。）

表3.1　西浦論文の補足情報の Table S1
表3.1　西浦論文の補足情報の Table S1。（クリックすると大きく表示される。）

表3.1 にあるように、活動度（を決めるパラメータ ω_h,\ ω_w ）は報告率毎に異なります。これは、活動度が報告率から独立していないことを示しています。

ω_h,\ ω_w が報告率毎に異なっているのは、西浦論文がこれらを報告率毎に計算したからです。西浦氏ら論文執筆グループがそのような計算方法を選んだのであって、これは必然ではありません。
報告率によって変化しない、（4つの）報告率に共通な ω_h,\ ω_w を計算するような尤度関数を作成し、計算することも選択できます。

西浦氏ら論文執筆グループは、（意識して行ったかは不明ですが、）図3.2 の構成ではなく、図3.3 の構成を選択しています。
この選択は不適切だと思います。

（４）活動度が報告率から独立しているべき理由

活動度と報告率は、以下のものを表しています。

  • 活動度は、人々の活動がコロナに与える影響をモデル化したものです。
  • 報告率は、行政がコロナに感染した人を捕捉する割合を示したものです。

実社会における「人々の活動がコロナに与える影響（活動度）」は、「行政がコロナに感染した人を捕捉する割合（報告率）」と独立していると考えられます。両者の間に特段の関係があるとは考えられません。ですから、西浦論文のモデルが良いモデルであるなら、実社会におけるコロナの感染態様を適切にモデル化したものであるなら、活動度は報告率から独立しているべきです。

本稿の主たる考察は以上です。次節は、以上に関連する参考事項です。

（５）活動度は報告率からどのような影響を受けているのか

論文 式(4) で計算される活動度 h_t が、どのように報告率の影響を受けるのかを実際に計算します。

西浦論文が用いたのは、google mobility の、6個ある指標のうち、"retail_and_recreation", "workplaces", "residential" の3個です。（ただし西浦論文は、それぞれ "community", "workplace", "household" と表記しています。この付近の表記には他にも問題があるので、別に整理しています。）
西浦論文は、式(4) の計算のために google mobility の指標を加工し、Figure S4 に示しています。

図3.4　西浦論文の補足情報 Figure S4
図3.4　西浦論文の補足情報 Figure S4

私は、図3.4 を再現しました。

図3.5　図3.4の再現
図3.5　google mobility データを加工して、図3.4を私が再現したもの
西浦論文 式(4) の α_t^{community},\ α_t^{house},\ α_t^{work}

図3.5は、google mobility データを加工して、図3.4を私が再現したものです。よく再現できていると思います。加工の詳細は（１０）に。 図3.5を描くのに用いたデータはこちら。

再現した活動度データと、西浦論文の補足情報 Table S1 を使うと、論文 式(4)を計算できます。

式(4)は、こちら。

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

再現した図3.5 が 式(4)の α_t^{community},\ α_t^{house},\ α_t^{work} です。 ω^{community} は定義から 1 であると本文に書いてあり、 ω^{house},\ ω^{work} は Table S1 に（ ω_h,\ ω_w の表記で）記載されています。

reporting coverage ω_h ω_w
1 −0.490 −0.083
0.5 −0.483 −0.088
0.25 −0.467 −0.096
0.125 −0.427 −0.118

表2 西浦論文の補足情報 Table S1 のうち、式(4) に関係する部分

これらを用いると h_t が計算できます。

図3.6　論文 式(4)のグラフ
図3.6　論文 式(4)を計算してグラフにしたもの
西浦論文 式(3)、式(4) の h_t

図3.6 は論文 式(4) の h_t です。描くのに用いたデータはこちら。4本の曲線が描かれているのは、表2 に示されている4つの reporting coverage 報告率毎の ω_h,\ ω_w を使って、論文 式(4) を計算したからです。

これが、報告率 reporting coverage が、活動度 h_t に与えている影響です。

ただし、論文 式(3)の左辺 R_{ab,t} に対して、図3.6 が示している報告率によるパラメータの変化が与える影響を考える際は、図3.6 以外の要因にも注意する必要があります。 R_{ab,t} が報告率から受ける影響は、（報告率によって変化する）活動度 h_t 経由のもの（＝図3.6）だけではなく、他の（報告率によって変化する）パラメータ経由のものもあるからです。

ややこしく書きましたが要するに、図3.6 からすぐに何かを言えるということはないと思います。

次節以降は、補足的な事項です。

（６）（参考グラフ）図3.7

図3.6 に示した論文 式(4) の活動度 が、西浦論文 Figure 1 にどのように影響したのかは気になるところです。論文 Figure 1 と 図3.6 を重ね合わせた図を作成しましたので、ご参考までに示します。こちらです。

図3.7　西浦論文 Figure 1 と図3.6 の一部を重ね合わせたグラフ
図3.7　西浦論文 Figure 1 と図3.6 の一部を重ね合わせたグラフ

図3.7 の重ね合わせに使ったのは、図3.6 の一部、報告率 1 の曲線だけです。西浦論文 Figure 1 が報告率 1 における曲線なので、この曲線だけを使いました。

ご説明するまでもなく、Figure 1 の形状に影響を与えるのは、図3.6 の活動度だけではありません。この他、感染やワクチンによる免疫獲得や、市中におけるデルタ株の比率なども影響を与えています。
ただしこれらは、活動度ほど「細かく」上下には動きません。（例えばワクチン接種による免疫保持者の割合を示した、補足情報 Figure S2 を参照。上下しますが、ゆったりとした動きです。）したがって、Figure 1 曲線の細かな上下動に大きく寄与しているのは、活動度ではないかと思われます。

（７）（参考グラフ）図3.8

本稿の内容を把握して頂いたなら、図3.6 に示した論文 式(4) の（報告率毎の）活動度 h_t に、論文 式(3) の（報告率毎の） p×r を掛算したグラフにご興味があるかも知れません。このグラフから読み取れる特記すべき事項はないと考えますが、ご参考として示しておきますと、こちらです。

図3.8　図3.6（論文 式(4)の h_t）に、報告率毎の p×r を掛けたもの
図3.8　（参考グラフ）図3.6（論文 式(4)の h_t）に、報告率毎の p×r を掛けたもの

図3.8 のグラフの形状は、図3.6 と似ていますが、違いも少しあり、報告率 0.125 の曲線と、その他の曲線との相対的位置がやや近づいたように見えます。グラフ縦軸が変化していますが、これは本質的問題ではありません。

（８）論文の活動度に関する表記の不統一

西浦論文の活動度 h_t に関する表記には不統一があるので、整理します。
h_t は論文 式(4)で計算され、式(3)で用いられています。式(4) の計算では、「google mobility の6指標のうち3指標（"retail_and_recreation", "workplaces", "residential"）と、それらを調整するパラメータ」を使います。

google mobility の3指標を指し示すはずの論文の表記は不統一で、以下の表のような対応関係があります。

google mobility での表記 retail_and_recreation residential workplaces
論文本文での表記 community household workplace
式(4)での活動度の表記 α^{community} α^{house} α^{work}
式(4)でのパラメータの表記 ω^{community} ω^{house} ω^{work}
補足情報 Table S1 での
パラメータの表記		 (表記なし *1) ω_{h} ω_{w}

表3 論文の活動度に関する表記の対応
(*1 Table S1 に community に関する項目はありません。この値は常に 1 と論文本文で説明されているので、記載を省略したと思われます。）

西浦論文の表記は不統一ですが、本稿では論文の表記を踏襲しています。

（９）論文の補足情報の所在

論文の補足情報の所在の詳細は、論文
https://www.nature.com/articles/s41598-023-44942-6
をご確認頂いた方がいいですが、以下に概要をご説明します。

論文の下部には、補足情報 Supplementary Information として幾つかの url が示されていて、そのうちの1つに本稿で何度か言及した Table S1 が入っています。

Table S1 が入っているファイルの url は
https://static-content.springer.com/esm/art:10.1038/s41598-023-44942-6/MediaObjects/41598_2023_44942_MOESM1_ESM.docx
になりますが、私の環境でこの url をクリックすると、Microsoft Word のファイルである 41598_2023_44942_MOESM1_ESM.docx が、ブラウザのデフォルトのダウンロードフォルダにダウンロードされます。そのファイルを開くと、ほぼ最下部に Table S1 があります。

（１０）google mobility の加工してパラメータαを得る方法の詳細

西浦論文は、google mobility を加工してパラメータ α^{community},\ α^{house},\ α^{work} を得ていますが、加工の細部は示されていません。私はこれを以下のように推測しました。図3.5が図3.4を非常によく再現していることから、私の推測は正しいと思われます。

α = (google\ mobility\ 指標の7日間移動平均) / 100 + 1

ここで、移動平均を計算する際には、python で、pandas.DataFrame.rolling(7, center=True).mean() を使いました。デフォルトではない center=True を指定している点を特に注意します。

（９９）更新履歴

  • 2023-12-30
    公開。
  • 2024-01-04
    別稿「西浦氏の Scientific Reports 論文について（２）」に、「更新履歴 2024-01-04」の大幅な修正をしました。これに合わせて本稿（７）の関連する箇所を修正しました。
  • 2024-01-10
    図3.7 を鮮明なものに差し替えました。グラフの内容は同じです。
  • 2024-01-11
    図3.5 の凡例の順序が 図3.4 の凡例と不一致だったのを修正しました。
  • 2024-01-23
    図3.5 と 図3.6 のデータを github に置き、リンクを貼りました。
  • 2024-01-25
    （１０）を加筆し、ここへの言及を追記しました。