外出自粛効果を考慮してコロナウィルスの感染者数を予測してみた

初めに

　私は感染症の専門家ではありませんので、理解した上で読んでください。

　2019年12月から中国湖北省武漢で発生した新型コロナ肺炎(Covid-19)が、日本にも広がり感染者数が増加しています。2020年4月7日には、日本でも緊急事態宣言が発令され、外出自粛が呼びかけられています。在宅勤務などが広がっていますが、まだ外出自粛が十分でないという報道も見かけられます。
　そのためどのくらい外出自粛すればよいのか、外出自粛により感染者数の増加がどう変化するのか、予測したいと思います。

予測アルゴリズム

　今回は、機械学習の回帰により感染者数を予測します。
　一口に回帰といっても、さまざまな方式のアルゴリズムがありますが、今回はscikit-learnのサポートベクトル回帰を使用します。

データ

感染者数データ

　KaggleのCOVID-19 Complete Dataset (Updated every 24hrs)のcovid_19_clean_complete.csvを使用させてもらいました。このデータは全世界のデータが入っていますが、日本の感染者数のみを使用します。
　日本については、2020/1/22から2020/4/9までのデータがあります。
　感染者数と、日々の新規感染者数は以下の図の通りです。

外出自粛効果

　外出自粛により日本全国でどのくらい人々の外出自粛しているのか、良いデータありません。しかし東京都の新型コロナウイルス感染症対策サイトに、都営地下鉄の利用者数の推移のデータがあります。東京都の都営地下鉄で、週毎という限られたデータではありますが、このデータを使用します。
　しかしながら、このデータはPDFで配布されています。そのため、CSV形式に手で入力しました(泣)。新型コロナウイルス感染症対策サイトで配布されている感染者データは、CSV形式で配布されているので、都営地下鉄の利用者数もCSV形式で配布してもらいたかったです。
　以下に都営地下鉄利用者数の増減率のグラフを示します。

日	新宿7時	新宿8時	新宿9時	渋谷7時	渋谷8時	渋谷9時	東京7時	東京8時	東京9時
2020/1/31	1.88%	-2.96%	0.39%	0.57%	-4.58%	-1.86%	-1.49%	-1.93%	0.44%
2020/2/7	0.18%	-1.03%	2.06%	-0.56%	-4.05%	1.65%	1.15%	0.84%	1.97%
...	...	...	...	...	...	...	...	...	...

評価方法

シンプルな予測

Kaggleの感染者数データから、以下のようなデータを抽出します。

Index	日	感染者数
1	2020/1/22	2
2	2020/1/23	2
3	2020/1/24	2
4	2020/1/25	2
5	2020/1/26	4
6	2020/1/27	4
7	2020/1/28	7
8	2020/1/29	7
9	2020/1/30	11
10	2020/1/31	5
...	...	...
79	2020/4/9	4667

説明変数Xは、[0, 1, ,2, 3, ...., 79]です。
予測対象Yは、感染者数[2, 2, 2, 2, 4, 4, 7, ...]です。

2020/1/22から2020/4/9のうち、前方の9割を学習データ、後方の1割をテストデータとしました。
学習データで回帰モデルを学習し、テストデータで評価しました。学習データ部分はほぼぴったりあっていますが、テストデータ部分は少し乖離があります。
単純な回帰だけでは、正しく予測できません。
さらに将来10日分も予測しました。

外出自粛を考慮した予測

シンプルな予測ではKaggleのデータのみ使用しましたが、さらに都営地下鉄の利用者数の減少率を追加して予測を行います。自粛の効果が現れるのは、約2週間後と言われています。そのため、感染者数と、2週間前の都営地下鉄利用者減少率を結合して、予測を行います。

2020/1/22から2020/4/9のうち、前方の9割を学習データ、後方の1割をテストデータとしました。
学習データで回帰モデルを学習し、テストデータで評価しました。学習データ部分とテストデータ部分も、ぴったりあっています。
さらに将来10日分も予測し、赤い太線で示します。
今後、感染者数は急激に増加していきます。

自粛の効果が予測にどう表れるのか、予測モデルで評価しました。
都営地下鉄利用者増減率が影響するのは2週間後になるので、もし2週間前(3月27日)から都営地下鉄利用者増減率が実際よりも40%低下していたとしたら、どうなるかシミュレーションしてみました。4月3日の地下鉄利用者は、時間帯により異なりますが約40%減少していました。更に40%削減なので、約80%削減した場合のシミュレーションです。

　2週間前の3月27日時点で、地下鉄利用者数を約80%削減できていたら、感染者数の増加ペースをだいぶ減らせる予測となりました。こういう予測を通して外出自粛の効果が目に見えると、外出を自粛しようという気持ちになりますね。ただし、効果が見えてくるのは2週間以上後なので、じっと我慢する必要がありますね。

参考文献

Kaggle : https://www.kaggle.com/imdevskp/corona-virus-report
東京都の新型コロナウイルス感染症対策サイト:https://stopcovid19.metro.tokyo.lg.jp/

プログラム

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import print_function
import numpy as np
import pandas as pd
from datetime import datetime
import copy
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.model_selection import GridSearchCV
from sklearn.svm import SVR

def read_confirmed():
    data = pd.read_csv('covid_19_clean_complete.csv')
    data_japan = data[data.loc[:, 'Country/Region']=='Japan']
    data_japan = data_japan.reset_index(drop=True)
    data_japan = data_japan.loc[:, ['Date', 'Confirmed', 'Deaths', 'Recovered']]
    data_japan['New_confirmed'] = (data_japan['Confirmed'] - data_japan['Confirmed'].shift(1)).fillna(0)
    return data_japan

def read_subway(data_japan):
    subway = pd.read_csv('20200409_subway.csv')
    data_japan1 = pd.merge(data_japan, subway, left_on='Date', right_on='Date3', how="outer")
    data_japan1 = data_japan1.drop(["Date1", "Date2", "Date3"], axis=1)
    data_japan1 = data_japan1.loc[:78, :]
    data_japan1 = data_japan1.fillna(0)
    data_japan1["x"] = data_japan1.index
    return data_japan1

def plot1(data_japan):
    # x=data_japan.loc[:, 'Date']
    x = np.arange(len(data_japan))
    plt.plot(x, data_japan.loc[:, 'Confirmed'], label='confirmed')
    plt.plot(x, data_japan.loc[:, 'New_confirmed'], label='New_confirmed')
    plt.title('Japan')
    plt.legend()
    plt.savefig('japan_confirmed.png')
    plt.cla()

def predict_svr1(data_japan):
    y = data_japan['Confirmed']
    x = np.arange(len(data_japan)).reshape((-1,1))
    Xtrain, Xtest, Ytrain, Ytest = train_test_split(x, y, test_size=0.10, shuffle=False)
    svm_confirmed = SVR(shrinking=True, kernel='poly',gamma=0.01, epsilon=1,degree=5, C=0.1)
    svm_confirmed.fit(Xtrain, Ytrain)
    Ytrain_pred = svm_confirmed.predict(Xtrain)
    Ytest_pred = svm_confirmed.predict(Xtest)

    # 将来の予測
    Xtest2 = np.arange(Xtest[-1]+1, Xtest[-1]+11).reshape((-1, 1))
    Ytest_pred2 = svm_confirmed.predict(Xtest2)

    # プロット
    plt.plot(np.arange(len(data_japan)), data_japan.loc[:, 'Confirmed'], label="confirmed", color='blue')
    plt.plot(Xtrain, Ytrain_pred, '--', label="train_pred", color='red')
    plt.plot(Xtest, Ytest_pred, label="test_pred", color='red', linewidth=1)
    plt.plot(Xtest2, Ytest_pred2, label="pred2", color='red', linewidth=3)
    plt.legend()
    plt.title('Japan')
    plt.savefig('japan_redict1.png')
    plt.cla()

def predict_svr2(data_japan):
    y = data_japan["Confirmed"]
    x = data_japan[['Shinjyuku_7', 'Shinjyuku_8', 'Shinjyuku_9', 'Shibuya_7', 'Shibuya_8', 'Shibuya_9', 'Tokyo_7', 'Tokyo_8', 'Tokyo_9', 'x']]
    Xtrain, Xtest, Ytrain, Ytest = train_test_split(x, y, test_size=0.10, shuffle=False)
    Xtest = Xtest.reset_index(drop=True)
    svm_confirmed = SVR(shrinking=True, kernel='poly',gamma=0.01, epsilon=1,degree=5, C=0.1)
    svm_confirmed.fit(Xtrain, Ytrain)
    Ytrain_pred = svm_confirmed.predict(Xtrain)
    Ytest_pred = svm_confirmed.predict(Xtest)

    # 将来の予測1
    Xtest2 = copy.deepcopy(Xtest)
    last = len(Xtest) - 1
    Xtest2.loc[:, 'x'] = np.arange(Xtest.loc[last, 'x'] + 1, Xtest.loc[last, 'x'] + 1 + len(Xtest))
    last = len(Xtest) - 1
    Xtest2['Shinjyuku_7'] = Xtest.loc[last, 'Shinjyuku_7']
    Xtest2['Shinjyuku_8'] = Xtest.loc[last, 'Shinjyuku_8']
    Xtest2['Shinjyuku_9'] = Xtest.loc[last, 'Shinjyuku_9']
    Xtest2['Shibuya_7'] = Xtest.loc[last, 'Shibuya_7']
    Xtest2['Shibuya_8'] = Xtest.loc[last, 'Shibuya_8']
    Xtest2['Shibuya_9'] = Xtest.loc[last, 'Shibuya_9']
    Xtest2['Tokyo_7'] = Xtest.loc[last, 'Tokyo_7']
    Xtest2['Tokyo_8'] = Xtest.loc[last, 'Tokyo_8']
    Xtest2['Tokyo_9'] = Xtest.loc[last, 'Tokyo_9']
    Ytest_pred2 = svm_confirmed.predict(Xtest2)

    # 将来の予測2
    Xtest3 = copy.deepcopy(Xtest)
    Xtest3.loc[:, 'x'] = np.arange(Xtest.loc[last, 'x'] + 1, Xtest.loc[last, 'x'] + 1 + len(Xtest))
    reduce = -0.1
    num = len(Xtest3)
    diff = reduce / (num - 1)
    Xtest3.loc[:, 'Shinjyuku_7'] = [Xtest.loc[last, 'Shinjyuku_7'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Shinjyuku_8'] = [Xtest.loc[last, 'Shinjyuku_8'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Shinjyuku_9'] = [Xtest.loc[last, 'Shinjyuku_9'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Shibuya_7'] = [Xtest.loc[last, 'Shibuya_7'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Shibuya_8'] = [Xtest.loc[last, 'Shibuya_8'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Shibuya_9'] = [Xtest.loc[last, 'Shibuya_9'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Tokyo_7'] = [Xtest.loc[last, 'Tokyo_7'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Tokyo_8'] = [Xtest.loc[last, 'Tokyo_8'] + diff * i for i in range(num)]
    Xtest3.loc[:, 'Tokyo_9'] = [Xtest.loc[last, 'Tokyo_9'] + diff * i for i in range(num)]
    Ytest_pred3 = svm_confirmed.predict(Xtest3)

    # プロット
    plt.plot(np.arange(len(data_japan)), data_japan.loc[:, 'Confirmed'], label="confirmed", color='blue')
    plt.plot(Xtrain['x'], Ytrain_pred, '--', label="train_pred", color='red')
    plt.plot(Xtest['x'], Ytest_pred, label="test_pred", color='red', linewidth=1)
    plt.plot(Xtest2['x'], Ytest_pred2, label="pred2", color='red', linewidth=3)
    plt.plot(Xtest2['x'], Ytest_pred3, label="pred3", color='green', linewidth=3)
    plt.legend()
    plt.title('Japan')
    plt.savefig('japan_redict2.png')

def main():
    # 感染者データリード
    data_japan = read_confirmed()
    data_japan.to_csv('data_japan.csv', index=False)
    plot1(data_japan)
    # 予測1
    predict_svr1(data_japan)

    # 感染者データリード
    data_japan = read_subway(data_japan)
    # 予測2
    predict_svr2(data_japan)

if __name__ == '__main__':
    main()