コロナの感染者数を初学者が予測してみた

#はじめに
最近学んでいる数学を利用してコロナに関して分析したいと思ったためこの記事を書いた。
統計学や機械学習などは初学者なため、稚拙な間違いもあると思われるが、適宜教えていただけるとありがたい。

#目的
##動機
統計、数学、機械学習を利用することによって新型コロナウイルスの感染について分析したいと考えた。

##背景
新型コロナウイルスの感染者が増加してから、自粛の必要性がたびたび議論されている。
感染者数を減らすためには、自粛するべきなのは明らかであるが、経済的な影響を考慮すると自粛の線引きは難しい。
そこで私は外出自粛率がどれくらいコロナの感染者数に影響しているのかを数値的に求めたいと感じた。

##仮説
外出を自粛するほど、新型コロナウイルス感染者数は減少する。

##方法
外出自粛率を含んだデータを用いて新型コロナウイルス感染者数を予測するモデルを作成する。
外出自粛率を変化させて新型コロナ感染者数の変化を調べる。

#手順
1. データの用意
2. モデルの作成
3. 学習
4. 精度計算
5. 予測

##データの用意
###データ入手場所
新型コロナ感染者を予測するために以下のデータを6/21に東京都庁、モバイル空間統計より取得した。
1.東京都 新型コロナウイルス陽性患者発表詳細(https://catalog.data.metro.tokyo.lg.jp/dataset/t000010d0000000068)
2.緊急事態宣言前後における全国主要都市の人口変動分析(https://www.nttdocomo.co.jp/utility/demographic_analytics/)
3.増減率の推移(https://mobaku.jp/covid-19)
それぞれのデータをデータ1,データ2,データ3とする。

###データの選択
データ1は東京都おける1/24~6/19の陽性診断のデータである。
今回は陽性数をコロナ感染者数とする。
データ1を用いて日ごとの感染者数、累計感染者数を求めた。

データ2とデータ3はどちらも地域ごとの人口の増加に関するデータである。
人口の増加率が地域ごとに記載されていて、前日比など様々なデータと比べた人口増加率になっている。
データ3は下記の表-1のようなデータであった

**表-1 データ3　地域ごとの人口増加データ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/9074bc6b-687d-bc93-d6f9-0db4578e7c52.png)

今回は自粛率を感染拡大前比率として分析をする。
感染拡大前比率の定義は取得元のNTTドコモより
「2020/1/18～2020/2/14の午後3時時点の平均値と比較した値である。（平日は平日平均と、休日は休日平均との比較）」
とする。
データ3は5/1~6/12の地域ごとの人口増加率のデータであった。
データ2は4/18~4/30の人口増化率の地域ごとのデータであった。
データ2、データ3の東京都に関する感染拡大前比率を利用する。
データ1、データ2、データ3の利用するデータを2020/1/24~2020/6/12分をまとめるたものをデータセットとする。

###欠損値の処理
感染者数の欠損値はなかった。
自粛率の欠損値について考えていきたい。
銀座、渋谷センター街、新宿駅、立川駅は4/17以前の自粛率のデータを取得できなかった。
東京駅、大手町、丸の内、東京駅南、霞が関、品川駅の4/21以前の自粛率のデータを取得できなかった。
羽田空港第1、第2ターミナルの4/23以前の自粛率のデータを取得できなかった。
さらに6/13以降のすべての地域の自粛率も取得できなかった。
つまり上記に記した範囲の自粛率のデータが欠損値となっている。
感染拡大前比率の定義より、感染拡大前の1/24~2/14の自粛率が平均値とほぼ変わりなかったと考えて自粛率を0.00%とする。
4月中の欠損値は4月中の中で一番早い存在する自粛率とする。
つまり、銀座、渋谷センター街、新宿駅、立川駅などは4/18の自粛率と同じであると仮定する。
そのほかの4月中の自粛率の欠損値も同様の扱いとする。
残りの欠損値は2/15~3/31と6/13以降の自粛率のデータである。
2/15~3/31の自粛率は１次関数的に自粛率が増加していたと仮定した。
つまり2/14の自粛率0.00%から4/1まで1次関数のように自粛率が変化していったと仮定した。
6/13以降のデータは後述する目標値の設定上欠損値の処理は行わない。

###目標値の設定
今回は感染者数を予測したい。
今回はある日から2週間後の感染者数を目標値とする。
そうすると、6/6~6/19のの2週間分の目標値は欠損値となってしまう。
今回はこの2週間分のデータは学習には利用しないものとする。
欠損値処理後、目標値を設定した後のデータセットが下記の表-2である。

**表-2 欠損値処理後、目標値設定後のデータセット**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/c9a0b2d1-5821-1744-9ee1-60adda216812.png)

###データの基本統計量
欠損値処理後、目標値処理後の基本統計量は下記の表-3、表-4になった。

**表-3 num,num_sum,targetの基本統計量**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/3773d329-ed5e-2674-e728-1f1ecff46e2f.png)

**表-4 地域ごとの自粛率の基本統計量**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/4a455ea5-00bd-a311-d2f6-34a5f6123dcf.png)


##モデルの作成
今回は感染者数を予測したいので、数値予測のモデルを作成する。
数理モデルはニューラルネットワークを利用した。
活性化関数はrelu関数、optimizerはadam、lossはMSE、metricsはMeanAbsoluteErrorを使用した。
中間層は今回は2層とした。
中間層の1層、２層のノード数は64とした。
出力層のノード数は1である。
構造の図式は下記の図-1である。

**図-1 ニューラルネットワーク構造**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/2328616e-12bb-eaf9-ba17-5c8e48c0f304.png)

##学習データ、テストデータ分割
作成したデータのうちの75%を無作為抽出して学習データとした。
残りの25%をテストデータとする。

##学習
精度を上げるために、学習データを平均0、分散を1にするための標準化を行った。
学習データの数が少ないため、K-分割交差検証を使用して学習をした。
k-分割交差検証の分割数は20とした。
学習するときのエポック数は30、バッチサイズは32とした。
横軸をk-分割回数、縦軸をval_lossとしたグラフが下記のグラフ-1である。

**グラフ-1 ニューラルネットワークのval_loosの変化のグラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/b780ac99-8582-b9d3-5c1f-2b3ffcf71b54.png)

val_lossが大きくなる時の検証データの目標値の分布を知るりたいと考えた。
上記のグラフ-1においてval_lossの値が1000以上の値をとった8回目と15回目のval_y(検証利用した目標値)のデータをが表-5、表-6である。

**表-5 k-分割交差検証8回目の検証データの目標値**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/114b1281-fdca-5c69-d0cc-313f95f57009.png)

**表-6 k-分割交差検証15回目の検証データの目標値**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/9fb663ed-1ffe-e1b0-38ac-9817733962b7.png)


##精度計算
モデル1のテストデータの目標値を横軸、縦軸を予測値とした散布図が下記のグラフ-2である。
（赤い線はy=xのグラフ）

**グラフ-2 目標値と予測値の散布図**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/97d18683-8453-0c4b-af97-eb0a6ee42fff.png)

1つのテストデータのみで評価したメトリクスを精度として扱うのはリスクがあると考えた。
そのリスクとはテストデータのデータ数が少ないため、テストデータの分布がメトリクスに大きく影響を与えてしまう可能性である。

このリスクを回避するために、学習データ、テストデータの分割を10回行った。
そして各分割に対して別々のモデルを用意し、学習させメトリクスを計算した。
この10個のメトリクスの値の平均値を、平均メトリクスと定義する。
1回目の分割に対応するモデルの名前をモデル1として、そのほかも同様とした。

この各モデルに対するメトリクスの値の注意点は使用するテストデータも変わっていると同時に、学習データも変わっていることである。
その場合、通常は学習データが違うモデルは全く異なるモデルになる。
しかし今回学習データとテストデータの分割は3:1なので、学習データの内容は大きくは変わらない。
よって、モデル1~モデル10の中身も大きく変わらないと仮定した。
仮定より、平均メトリクスの値は作成したモデルの精度評価として利用してよいと考えた。

各学習済みモデルのメトリクスは下記のようになった。
モデル1  :13.404623985290527
モデル2  :12.988094329833984
モデル3  :16.170719146728516
モデル4  :10.167283058166504
モデル5  :16.49835205078125
モデル6  :11.155875205993652
モデル7  :16.170719146728516
モデル8  :7.813872337341309
モデル9  :13.984491348266602
モデル10 :12.988094329833984

これらの平均値は約13.1342となった。

モデル1（1回目の分割のデータセットを利用して学習したモデル)をこれ以降では使用する。
モデル1の学習データの基本統計量、モデル1におけるテストデータの基本統計量が下記の表-7、表-8である。

**表-7　モデル1の学習データの基本統計量**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/b32c0be3-f0cc-9a20-ddb7-dd82de6946e6.png)

**表-8 モデル1におけるテストデータの基本統計量**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/ba52ce0f-9373-1a29-df46-92685ef4c2b2.png)


##予測
今回の目的は自粛率の効果を検証することである。
つまり自粛率を変化させて、感染者数の増減を調べたい。
予測に必要データは、感染者数、累計感染者数、自粛率である。

予測用のデータとして感染者数、累計感染者数数は6/28~7/2のデータを利用する。
6/28の感染者数は60人、累計感染者数は6115人である。(https://www.asahi.com/articles/ASN6X6T24N6XUTIL001.html朝日新聞より)
各地域での自粛率は0%～-99%まで-1%刻みで変化させたデータを作成して予測する。
つまり下記の表-3のようなデータを利用して2週間後の感染者数者数を予測する。

**表-3 6/28における自粛率の影響を予測するためのデータ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/ff2116dd-b92c-2d13-e546-8999e77f3571.png)
モデルを利用する前に学習データの平均値、分散値を利用して、このデータの標準化を行った。

これらのデータを使用して2週間後のコロナ感染者数を予測した。
予測に利用したモデルはモデル1である。
横軸を自粛率(%)、縦軸を予測値(6/28から2週間後のコロナ感染者数)としたグラフが下記のグラフ-3である。

**グラフ-3 6/28のデータを用いた自粛率による感染者の変動グラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/39bfd779-c7f3-4cce-b76c-163e01dbb1a5.png)

同様の作業を6/29~7/2の範囲で行った。
それが下記のグラフ4~グラフ7である。
各日の感染者数は順番に、58、54、57、107人であった。

**グラフ-4 6/29のデータを用いた自粛率による感染者の変動グラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/b188ae02-932c-96df-4977-79a861546159.png)

**グラフ-5 6/30のデータを用いた自粛率による感染者の変動グラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/9998e72b-8c6a-bac9-33e5-26f307c4cd16.png)

**グラフ-6 7/1のデータを用いた自粛率による感染者の変動グラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/8622af05-2726-b1b5-00fd-7679281ac4fe.png)

**グラフ-7 7/2のデータを用いた自粛率による感染者の変動グラフ**
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/631346/1de8ee7b-dfab-ac99-ab6a-537e07edef77.png)

#考察
##モデル評価
今回自粛率、コロナ感染者数、累計感染者数を用いて2週間後のコロナ感染者数を予測するモデルを作成した。
メトリクスの値が約13のモデルを作成することができた。
今回メトリクスはMAEを利用した。
MAEの値は作成したモデルの予測値と目標値の平均誤差である。

目標値の値の範囲によって予測の精度が変わっていることがグラフ-2よりわかる。
目標値が小さい範囲ではy=xの直線に近い位置に点が分布しているため、精度は高いと考える。
しかし目標値が大きくなってしまうと誤差が大きいことがグラフ-2よりわかる。
表-4、表-5からも2週間後の目標値が大きいとval_lossの値が大きくなってしまうことがわかる。
まとめると、予測値が小さい時は精度が高く、予測値が大きい時は精度が低くなっている。

今回のモデルの誤差の原因はいくつか考えられる。
1番誤差に影響があったのは検査数の変化であると考える。
今回作成したモデルは検査数の増加を考慮せずに感染者数の予測をしている。
検査数が増加すれば、自然と感染者数も増加するがそれが今回考慮されいないため、誤差が生まれてしまったと考えられる。

2つ目の誤差の要因は、時系列を考慮していないことが考えられる。
今回作成したモデルは過去の感染者数をnum_sum(累計感染者数)として利用した。
つまり前日との感染者数の変化などを考慮してない。
感染者数の増加スピードなどの時系列データを予測するときに利用していれば精度が上がる可能性がある。

3つ目ににデータ数の少なさが誤差を生んでしまったと考えられる。
学習に利用したデータ数は100個、説明変数の数は14個である。
データ数、説明変数の数がともに少なかったため誤差が生まれてしまったと考えられる。
データ数を増やすためには、東京都だけでなく全国や外国のデータを学習させる方法が考えられる。

誤差を抑えるためには、ハイパーパラメータの設定などを工夫する方法も考えられる。

##仮説検証
今回自分が立てた仮説は、「外出を自粛するほど、新型コロナウイルス感染者数は減少する。」であった。
仮説を検証するために5日分のデータを用いて自粛率を変化させた結果の感染者数のグラフを作成した。
それがグラフ-3～7である。
これらはどれも似た概形を示した。

予測値の値はすべて30以下であったので予測値が小さいとして、モデルの精度はこの状況においては良いと仮定する。
グラフ-3~7を見てみると-70%~-40%の範囲では自粛するほど感染者が大きく減少することがわかる。
-35%~0%の間でも同様に自粛するほど感染者は少し減少する。
しかし、-35%~-40%辺りでは自粛率が減少すると感染者数も増加してしまっている。
-70%~-99%の間も同様に自粛するほど感染者数が大きく増加してしまっている。
これは一般的なイメージで考えるとおかしいと感じる。
直感的に考えると、自粛率が減少すると人同士の接触が少なるなるので感染者数も減少するはずである。

しかし、大きな間違いにここで気づいた。
この自粛率を変化させた場合の予測は、現実社会のシュミレーションになっていないことである。
学習データの自粛率というのは世間のコロナへの不安感などの情報で決まるものである。
よって自粛率が低いということはコロナ感染の可能性が高いと考えられる。
その世間の不安感が自粛率に考慮したうえでのモデルが今回作成したモデルであると考えられる。

よってそもそもの自粛率というものを自分で仮定していいものではないため仮説の検証は不可能である。

#結論/課題
・テストデータと比べたとき、今回作成したモデルは予測値が小さい時は精度が高い
・精度を高めるために、検査数の考慮、時系列分析、データ数の増加が必要
・自粛率は世間のコロナへの不安などが絡んだうえのデータのため、自粛率を自分で変えてシュミレーションを行い、感染者数を予測することは今回のモデルではできない
・今回のモデルでは仮説の検証は不可能

#感想
AI、ニューラルネットワーク、データサイエンスの難しさを自分の手を動かして初めて感じました。
そして、なんでもニューラルネットワークにぶち込めばいいものでもないことも学びました。
ニューラルネットワークの中身がブラックボックスなため、原因の究明というものが難しかったです。
自粛率というものを勝手に自分で仮定してはいけないという点に気づいたときは驚きました。
変数として自分で変えていい場合と変えてはいけない場合を意識していきたいと思います。
分析の中で難しかった点は下の2点でした。

・コロナの検査数を考慮していないこと
・k-分割交差検証で学習させたモデルを利用して予測するときは予測に使うデータは何を使って標準化すればいいのかがわからなった

特に検査数の変化は専門家の方にとってもなかなか分析をするうえでめんどくさそうな要因でありそうだなと初学者ながら感じました。
途中で検査数が変わるためそれを考慮したほうがいいことには気づいたのですが、方法がわかりませんでした。
時系列分析をすればより良いモデルが作れそうだなと感じたので、時系列分析の勉強も進めていきたいと思います。

コロナが収まるのがいつかわからない状況ですが、頑張っていきましょう!

ここまで読んでいただきたありがとうございました。
何か疑問点、問題点がありましたら指摘していただけると嬉しいです。

--追記--
江崎貴裕　(2020),データ分析のための数理モデル入門　本質をとらえた分析のために　ソシム株式会社
を読んで、この記事において理解志向型モデリングと応用志向型モデリングの区別ができていなかったことに気づきました。
詳しく知りたい方は上記の本を読んでください、すごくおもしろかったです。