新型コロナウイルス(COVID-19)感染リスク予測分析の基礎

    目次

1.概要

10月以降の本格的なCOVID-19第2波危機を前に,日本では早くも感染が広がってしまった. 現在の感染者数はピークを超えた兆しをみせているものの,社会的距離戦略などの適切な対応策がないままでは,どのような条件で,またいつ感染が拡大をしてしまうのか? その感染リスクを抱えながら,日常生活を送らなければならない. これまでに公開されている多くのサイトは,国内外の累積統計を可視化するにとどまり, 「どのような生活を送れば良いか?」 「どのような具体的な戦略で経営を回復させるか?」 などの意思決定は,私たち国民の自己判断に委ねられているのが現状である. (株)データフォーシーズ AI Labは,弊社の経営理念・ミッション・ビジョンに従い,Data Science × AIのエバンジェリストとして, 国民一人一人の感染リスク・感染予防に関する意識向上と経営改善に関する適切な施策の提案を目指し, 累積統計の可視化に止まらず,統計・AI・機械学習の基礎から透明性の高い科学的根拠を持って,国民の人々に分かりやすく感染リスクの定量化を説明する.

2.導入

中国・武漢を発生起源と考えられる新型コロナウィルス感染症(COVID-19)の世界規模の大流行は,日本を含む世界の約7割でCOVID-19感染が再拡大し[1], 感染拡大が終息する兆しさえ見せていない(図1). Johns Hopkins University & Medicine, Coronavirus Resource Centerが公開するデータから[2][3], 現在までの世界と日本での感染者数・死亡者数・回復者数の累計は下の表1で示される. COVID-19流行は,気象条件などの季節性にも影響されることが示唆されている[4]. 特に,日本は,中国や韓国と隣接していたにも関わらず,1次感染拡大は,アジアの他の地域に比べ非常にゆっくりしていたが(図2.各国での感染拡大の時期の違いを見やすくするため,最大感染確認者数で規格化もので示した), 梅雨という独特の気象条件の影響の可能性も重なり,4月中の自粛時期以上に感染が再拡大した. 8月上旬から中旬にかけ感染者数ピークを迎え,現在は新規感染者数・感染者数は減少傾向を示している. しかし,社会的距離戦略などを踏まえた適切な拡大予防対策がなされないままでは,10月以降の本格的な第2波で8月以上に感染拡大してしまうことが懸念される. 本分析レクチャーサイトでは,国民の感染危機や感染予防の意識向上や,経営改善に関する適切な施策の提案を目指す. 累積統計の可視化に止まらず,統計・AI・機械学習の基礎から透明性の高い科学的根拠を持って,国民の人々に分かりやすく感染リスクの定量化を説明する.

表1:最新COVID-19感染状況
感染確認者数死亡者数回復者数
世界
日本



3.AI・機械学習アルゴリズム

3.1データセット

本サイトでは,海外のCOVID-19感染状況データは,Johns Hopkins University & Medicine, Coronavirus Resource Centerが公開するサイト[2][3]から取得される. 国内のCOVID-19感染状況データは,ジャッグジャパン(株)[5]が公開するデータを利用する.

3.2ベイズ統計学

ベイズ統計学は,トーマス・ベイズ[6]によって提案されたベイズの定理が基盤となり,1700年代に確立された理論である. ベイズの定理により,人々が持つ先験的な確信(原因:例えば,専門家の意見に基づいたりして獲得される事前の情報)を観測データ(結果)と結び付けることで, 推定される確率(逆確率)を求めることができるので[7], しばしば,「確率の更新(アップデート)」として指すこともあり, その応用範囲は,診断検査[8][9]から時系列予測[10]に至るまで非常に広い. この「確率の更新」から病気の罹患率から考えてみる. 例えば,肺炎の罹患率が\(P_h(検査前)=0.2\%\)(テーブル内の前提:肺炎の罹患率が500人に1人)だったとする.

ベイズ
表2:検査マトリックス
前提肺炎の罹患率\(P_h(検査前)\):500人に1人
検査陽性陰性
罹患\(82\%\)\(18\%\)
非罹患\(7\%\)\(93\%\)


この時,図3で示される確率を持つ検査を受けて陽性と診断されてしまった時,実際の肺炎を罹患している確率\(P_h(検査後)\)は次のように求めることができる:
  1. 検査前の前提だった肺炎の罹患率が事前確率とする.
  2. 肺炎を罹患が「原因」となって,陽性と判定される「結果」に関する条件付確率は,テーブルの\(82\%\)である.
  3. 検査での陽性判定が「原因」となって肺炎を罹患した「結果」となる同時確率は,\[P(罹患\cap陽性)=P(陽性|罹患)P(罹患)\]より,\[0.82 \times 0.002 = 0.00164\]となる.
  4. 肺炎にならなかった場合の非罹患率は,\[事前確率 = 1-0.002 = 0.998\]である.
  5. 肺炎を罹患しないが「原因」となって,陽性と判定される「結果」の条件付確率は,テーブルの\(7\%\)である.
  6. 検査での陽性判定が「原因」となって,肺炎にならなかった「結果」となる同時確率は,\[0.07 \times 0.998 = 0.06986\]となる.
  7. 検査で陽性判定を受けて,実際に肺炎になってしまった罹患率は,\[ \cfrac{ 0.82 \times 0.002 }{ 0.082 \times 0.002 + 0.07 \times 0.998}=0.023\]で,\(2.3\%\)となる.
これは,検査前の肺炎の罹患率が,上の表の検査(肺炎の罹患が「原因」で陽性判定を受けた「結果」の条件付確率)を通して,その罹患率が約10倍に更新されたことになる. ベイズの定理は,1つの疾患だけでなく,複数の疾患を想定する場合や,複数の疾患を複数の診断法で診断する場合にも適用することができる[8] 実際,ベイズの定理を用い,中国における感染拡大時でのCOVID-19感染検査の信頼性を評価することができる[11].前述のように,過去の疾病の罹患率などの事前に獲得された情報をCOVID-19感染状況の観測データと結び付けることで,感染リスクやをCOVID-19との合併による重症化のリスクの定量化も可能であると考えられる.

\(k\)-平均法

世界各国のCOVID-19時系列を,下記の最適化問題を解く,非階層型クラスタリングアルゴリズムの1つ\(k\)-平均法[12]を用い, COVID-19時系列に潜むデータ構造を見出す:\[ \argmin_{V_1,\ldots,V_k} \sum_{i=1}^n \min_j \| x_i − V_j \|^2 \] 上記の各クラスタ平均\(V_j\)(\( j = 1, \ldots, k\))からの距離が最小になるようにクラスタ数\(k\)に分類することができるようになる. 図4は,2020年1月22日から2020年6月1日までのクラスタ数が\(k=5\)の時の東アジアから南アジアにかけての第1次感染拡大の推移を可視化した結果である. この結果が示すように,各クラスタが任意の拡大期間で分類されるように,その平均\(V_j\)(\( j = 1, \ldots, k\))からの距離を最小にしているがわかる. また,日本は,中国と隣接しているにも関わらず,台湾や韓国,東南アジア諸国に比べ,第1次感染拡大が異様に遅いことは分析上大変興味深い.

時系列フィッティング

上記のベイズ統計学の頻度主義的な確率の概念を用い,COVID-19感染者数や回復者数,死亡者数等の観測された変数にノイズがのっていることを仮定することで,COVID-19感染時系列の多項式曲線フィッティングを考える[13]. これは,実際の感染状況時系列データから再構築される時系列モデルが,どの程度,以下で記述される感染症数理モデルに当てはめることができるのかをみるためである. 本サイトでは,スパースコーディングを用いた感染状況の時系列データのフィッティングを説明する[14][15]. スパースコーディングの目的とは,Johns Hopkins University & Medicine, Coronavirus Resource Centerが公開するサイト[2][3]から取得したCOVID-19感染者数データ\(\vec{x_i}=(x_i(t_0), \ldots, x_i(t_l))\)(ここで,\(t_0\))を2020年1月22日とし,\(t_l\)を最新日とする)を,「辞書」と呼ばれるm個の基底関数の集合\(\Phi=[\phi_0,\phi_1, \ldots, \phi_{m-1}]\)のできるだけ少ない線形結合で時系列データを復元することである. ここでは,次の1次元の離散コサイン変換を基底関数の集合を「辞書」として用いる: \[ \phi_j(t)=F(j)\cos\left[ \cfrac{\pi}{m} \left( t+\cfrac{1}{2} \right)j \right] \]ここで, \[ F(j) = \begin{cases} \cfrac{1}{\sqrt{m}}, & \mbox{for } j = 0\\ \sqrt{\cfrac{2}{m} }, & \mbox{for } j \not= 0 \end{cases} \] (参考文献[15]を参照する). そして,実感染者数の時系列データ\(\vec{x_i}\)と辞書\(\Phi\)を用い,次の最適化問題 \[w^* =\argmin_{w}\|x_i − \Phi w \|^2_2 + \lambda \| w \|_1 \]を解くことで,下記に定式化されるような,スパース(ほとんどの要素が0)な\(w^*\)を求めることである: \[\vec{x_i} \simeq \Phi w^* \] 図5で、求めたパラメータによる感染者数の時系列(図5の赤線)のフィルタ結果(図5の緑点)を図4に示す. また,このような手法で得られる時系列データにフィッティングした曲線を平均とした確率分布から,将来の感染確率の推移を推定することができるようになるかもしれない.

感染症数理モデル

COVID-19がどのように感染拡大するか,その時の感染者数,或いは,流行の期間などを予測したり,基本再生産数(疫学的には,感染症への感受性をもつ全個体の集団を仮定し,1感染個体から直接生み出される感染個体数の平均と考えられ,集団内の感染症の蔓延する速度を推定するものではないことに注意する)などの様々な疫学的パラメータを推定したりするためには,また,異なる公衆衛生上の介入が流行の結果にどのような影響を与えるかを示すためには,数理モデルを用いることが非常に有効であると考えられている. 代表的な感染症数理モデルの1つとして,1927年にW.O. KermackとA.G. McKendrickによって提案されたSIRモデルが挙げられる[16]. このSIRモデルから推定される基本再生産数\(R_0\)は,\(R_0 > 1\) が感染症の集団内で蔓延を表し,その値が大きいほど流行を抑えることが難しく,\(R0 < 1\) は,感染症が蔓延しないことを表す.ここでは,Johns Hopkins University & Medicine, Coronavirus Resource Centerが公開するサイト[2][3]から取得されるCOVID-19時系列データをSIRモデル[16]に当てはめられると仮定し,[17][18]を参考に,相図分析による基本再生産数\(R_0\)の推定方法を説明する. SIRモデルは,次のように記述される: \begin{aligned} {\frac {dS}{dt}}(t) & = -\beta S(t)I(t)\\ {\frac {dI}{dt}}(t) & = \beta S(t)I(t)-\gamma I(t)\\ {\frac {dR}{dt}}(t) & = \gamma I(t) \end{aligned} このとき,\(S(t)\),\(I(t)\),\(R(t)\)は,それぞれ,感受性保持者数,感染者数,回復者数を表す. 上式での,\(\beta\),\(\gamma\)はそれぞれ感染率,回復率であり,2つのパラメータを用い,基本再生産数は,\[R_0=\cfrac{\beta}{\gamma}\]として表される. COVID-19感染状況では,\(S\)が未知であるが,感染が一旦終息した時期の感染確認者の累計\( {\bar C}\)(ほぼ一定値をとる)を総人口と仮定し, \begin{aligned} S(t) & = {\bar C}-R(t)-I(t) \end{aligned} とする. 図6の上段の\((S,I)\)相図の軌道分析をする.初期値を\(S(0)=1\),\(I(0)=0\)とすると,\(S\)の減少とともに\(I\)が増加し,ピーク後急速に減少していく軌道が描かれる. 2次感染拡大がなければ,本来は,\(S(\infty)=I(\infty)=0\)の終息軌道を描く. そのような軌道上のピーク,即ち,\(\cfrac{dI}{dS}=0\)(即ち,\(\cfrac{dI}{dt}=0\))となる変曲点から基本再生産数\(R_m\)を求められる: \begin{aligned} 0 & = \gamma I_m ( {\frac{\beta}{\gamma} } S_m - 1 )\\ R_m &= {\frac{1}{S_m}} \end{aligned} 実際の総人口は,ほぼ一定値をとる累積感染確認者数よりも多いので,基本再生産数は,\(R_0 \lesssim R_m \)となると考えられる. しかし,今年のCOVID-19大流行は,図2-図5で示されるように2次感染拡大が起きてしまった.図6の中段の\((I,R)\)相図や下段の\((R,S)\)相図を含む3次元\((S,I,R)\)相図から, 日本におけるCOVID-19流行は,将来,\((S(\infty), I(\infty), R(\infty))=(0,0,1)\)となるような終息軌道を描くことが予想される.\((S,I,R)\)の3変数は,感染終息の特徴量となりうる.終息状態からの距離から終息確率を算出できるかもしれない.

参考文献

[1]日本経済新聞「感染再拡大、世界の7割 〜経済両立可能な対策探る」(2020年8月2日付)
[2]Johns Hopkins University & Medicine, Coronavirus Resource Center
[3]https://github.com/CSSEGISandData/COVID-19
[4]Mohammad M Sajadi, et al. “Temperature, Humidity, and Latitude Analysis to Estimate Potential Spread and Seasonality of Coronavirus Disease 2019 (COVID-19).” JAMA Netw Open. 2020; 3:e2011834.
[5]都道府県別新型コロナウイルス感染者数マップ(ジャッグジャパン株式会社提供)- Dashboard&Map of COVID-19 Japan Cas
[6]https://en.wikipedia.org/wiki/Thomas_Bayes
[7] J. P. Klein, M. L. Moeschberger(著)打破守(訳)「Survival Analysis-2nd Edition 生存時間解析」(丸善出版,2012年)
[8] 森實敏夫(著)「新版 入門 医療統計学-Evidenceを見出すために」(東京図書,2016年)
[9] 折笠秀樹「ベイズ統計の応用領域 としての診断検査」 Jpn. Pharmacol. Ther. (薬理 と治療) vol.47 no.8, (2019)
[10] S.Kozawa, T.Akanuma, T.Sato, Y.D.Sato, K.Ikeda, T.N.Sato ``Real-time prediction of cell division timing in developing zebrafish embryo” Scientific Reports Vol. 6, 32962, (2016)
[11] A. Chen “Bayes' Rule, Unreliable Diagnostic Testing, And Containing COVID-19” towards data science Sharing concepts, ideas and codes (February 16th 2020)
[12]H. Steinhaus “Sur la division des corps matériels en parties” (French). Bull. Acad. Polon. Sci. 4 (12): 801–804 (1957).
[13] C. M. ビショップ(著)元田浩等(訳)「パターン認識と機械学習 上 - ベイズ理論による統計的予測」(丸善,2007)
[14]Sparse coding with a precomputed dictionary
[15]M. Mai, M. D. Shattuck, C. S. O'Hern “Reconstruction of Ordinary Differential Equations From Time Series Data" arXiv:1605.05420 (2016)
[16]W. O. Kermack and A. G. McKendrick “A Contribution to the Mathematical Theory of Epidemics”. Proc. Roy. Soc. of London. Series A 115 (772): 700-721 (1927). doi:10.1098/rspa.1927.0118. JFM 53.0517.01
[17]Malaria Math “Phase Plane Analysis of Kermack-McKendrick Model." (2017)
[18]M. Martcheva “An Introduction to Mathematical Epidemiology.” (Springer, 2010).