こんにちは。レッジのデータサイエンティストの松本です。今回が2回目の投稿になります。
レッジでは、クライアント先に常駐してデータ・ドリブンな課題解決に取り組んだり、ダイナミックプライシングアルゴリズムを受託開発したり、クライアント先へのBI導入の推進など、幅広くデータ利活用に関わる業務に取り組んでいます。
さて、近頃の話題といえば新型コロナウィルス（COVID-19）で持ちきりですが、忘れていませんか？毎年冬に流行しているヤツの存在を。
そう、インフルエンザです。新型コロナウィルス関連の報道に隠れてしまって、今年はインフルエンザの報道を見かける機会がほとんどありませんが、ヤツは一体どうなったのか気になりませんか？
確かに新型コロナウィルスへの感染対策として、ソーシャルディスタンス（社会的距離）やマスク・手洗い等のエチケットの徹底、さらには人の移動の激減によって、インフルエンザを含む感染症の発症数が例年よりも抑えられていることは容易に想像できます。 しかし、例年よりも具体的にどのくらいの人たちがインフルエンザに感染しなかったのかはどのニュースメディアでも取り上げられていないように思います。
そこで、今回の記事では新型コロナウィルスの影響によって、どれくらいインフルエンザの感染者数が減ったのかをPythonを使って検証してみます。
効果検証のアイデア
日本国民が新型コロナウィルスを意識するようになって以降の期間を介入1があった期間とみなします。
効果検証は、介入があったデータ群と介入がなかったデータ群との比較を基本としますが、今回に関して言えば、タイムマシンとか別の世界線に行かない限り介入がなかったデータ群を得ることはできません2。
したがって、2019年以前のインフルエンザ感染者数のデータをもとにした予測結果、つまり、介入がなかったとしたらこういう結果になったであろうという反実仮想的な結果を非介入のデータ群とみなすことにします。
このように、 介入が行われていないデータ群をもとに、現実には存在しない非介入のデータ群を補完する手法はCausalImpactと呼ばれます。 効果検証の進め方をまとめると以下になります。
- 2019年以前のインフルエンザ感染者数のデータをもとに時系列予測モデルを構築
- 時系列予測モデルを使って2020年1月から3月のインフルエンザ感染者数を予測
- 予測結果と実際に観測されたインフルエンザ感染者数との差分を算出
検証
以下の環境で検証を行いました。
- Python: 3.7
- pandas: 1.0.1
- fbprophet: 0.6
データの読み込み
WHO（世界保健機関）の公式ホームページで、各国のインフルエンザ発症数のデータをダウンロードできます。 https://apps.who.int/flumart/Default?ReportNo=12
今回は以下の条件に合致するデータをダウンロードし、使用します。
| 検索条件 | 選択項目 |
|---|---|
| Select by | Country, area or territory |
| Filter by | Japan |
| Year from | 2010 |
| Week from | 1 |
| Year to | 2020 |
| Week to | 52 |
また、効果検証対象のカラムは、A型インフルエンザ合計感染者数を意味するINF_Aとします。
前処理
今回の検証で使用するライブラリを読み込みます。
import pandas as pd from fbprophet import Prophet import matplotlib.pyplot as plt import seaborn as sns from datetime import date sns.set()
次に、ダウンロードしたデータをpandas.DataFrameで読み込み、時系列予測ライブラリProphet3でfitできるようにデータ型を変換し、カラム名を変更します。
# https://apps.who.int/flumart/Default?ReportNo=12からダウンロードしたcsv flu = pd.read_csv("FluNetInteractiveReport.csv") # 週開始日を文字列から日付型に変換 flu.loc[:, 'SDATE'] = pd.to_datetime(flu['SDATE']) # Prophetで利用できるようにカラム名変更 flu_all = flu.rename(columns={'SDATE': 'ds', 'INF_A': 'y'})
2020年1月1日より前を学習データ、2020年1月1日以降を効果検証用のデータとします。
# 学習データと効果検証用データに分割 train = flu_all.query("ds < '2020-01-01'") actual = flu_all.query("ds >= '2020-01-01'")
予測
学習データをもとにProphetで時系列予測モデルを構築します。 なお、今回はCausalImpactという効果検証の手法に焦点を当てています。 したがって、時系列予測モデルの精度向上のためのチューニングはせず、デフォルトで実行しています。
# Prophetで学習 model = Prophet() model.fit(train) # 週単位で2020-01-01から20週間先までの予測対象日を生成 future = model.make_future_dataframe(periods=20, freq='w') # 予測 pred = model.predict(future) # 予測結果を可視化 def plot_prophet( model, predicted, xlabel='ds', ylabel='y', figsize=(10, 6) ): """ Prophet学習済みモデルをもとに予測結果を可視化するメソッド。 Parameters ------- model : Prophet Prophet学習済みモデル。 predicted : pd.DataFrame Prophetによる予測結果。 xlabel : str xlabel。 ylabel : str ylabel。 figsize : tuple figsize。 Returns ------- figure : figure object Note ------- https://github.com/facebook/prophet/blob/master/python/fbprophet/plot.py """ fig = plt.figure(figsize=figsize) ax = fig.add_subplot(111) model_t = model.history['ds'].dt.to_pydatetime() predicted_t = predicted['ds'].dt.to_pydatetime() ax.plot(model_t, model.history['y'], 'k.', label='actual') ax.plot(predicted_t, predicted['yhat'], ls='-', c='#0072B2', label='predicted') ax.fill_between(predicted_t, predicted['yhat_lower'], predicted['yhat_upper'], color='#0072B2', alpha=0.2) # Specify formatting to workaround matplotlib issue #12925 locator = AutoDateLocator(interval_multiples=False) formatter = AutoDateFormatter(locator) ax.xaxis.set_major_locator(locator) ax.xaxis.set_major_formatter(formatter) ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2) ax.axvspan(model_t.max(), predicted_t.max(), color="gray", alpha=0.3) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) ax.legend(loc='upper right') ax.set_title('actual and predicted by Prophet') fig.tight_layout() plot_prophet(model, pred)
効果検証期間（2020年1月から3月）を灰色で色付けしています。 2020年1月1日から20週間先までの予測結果が可視化されていることを確認できました。
効果検証
「効果検証のアイデア」で説明した通り、CausalImpactを用いて効果検証を行います。 まず、Prophetでpredictした結果から、効果検証の期間のデータを抽出します。
# 予測結果のうち2020-01-01以降を抽出 pred_future = pred.query("ds >= '2020-01-01'") # Prophetの仕様で年を跨いだweeklyの加算で不具合があるため、観測データに揃えるように補正 pred_future.loc[:, 'ds'] = actual.ds
時系列予測モデルProphetによる予測結果と、実際に観測されたインフルエンザ感染者数を可視化します。
def plot_effect( actual, predicted, xlabel='date', ylabel='Number of infected', figsize=(10, 6) ): """ 効果検証期間の予測結果と観測結果を可視化するメソッド。 Parameters ------- actual : pd.DataFrame 観測結果。 predicted : pd.DataFrame 予測結果。 xlabel : str xlabel。 ylabel : str ylabel。 figsize : tuple figsize。 Returns ------- figure : figure object """ fig = plt.figure(figsize=figsize) ax = fig.add_subplot(111) ax.plot(actual.ds, actual.y, label='actual') ax.plot(predicted.ds, predicted.yhat, label='predicted') ax.axvspan(actual.ds.min(), date(2020,4,1), color="gray", alpha=0.3) locator = AutoDateLocator(interval_multiples=False) formatter = AutoDateFormatter(locator) ax.xaxis.set_major_locator(locator) ax.xaxis.set_major_formatter(formatter) ax.legend() ax.set_title('difference between predicted and actual') fig.tight_layout() plot_effect(actual, pred_future)
効果検証期間（2020年1月から3月）を灰色で色付けしています。
予測結果では1月末にインフルエンザ感染者数がピークを迎えていますが、実際に観測された感染者数は1月初日から減少傾向にあります。
グラフからも予測結果と実際の観測結果との間でかなり数値に乖離があることがわかります。
次に、2020年1月から3月で合計何人のインフルエンザ感染者数が減ったのかを確認します。
pred_sum = pred_future.query("ds < '2020-04-01'").yhat.sum() actual_sum = actual.query("ds < '2020-04-01'").ALL_INF.sum() print(f'predicted total:{pred_sum}, actual total: {actual_sum}') print(f'diffrence between predicted and result:{pred_sum - actual_sum}')
predicted total:4915.167388577098, result total: 2091.0 diffrence between predicted and result:2824.167388577098
上記より、1月から3月のインフルエンザ感染者数の予測結果と観測結果の差分は2,824人となりました。
つまり、新型コロナウィルスの影響によって、インフルエンザ感染者数が2,824人低減したということになります。
考察
今回の検証から、新型コロナウィルスの影響によって、インフルエンザ感染者数が低減していることがわかりました。
ただし、その原因は単純にソーシャルディスタンス（社会的距離）やマスク・手洗い等のエチケットに対して国民の意識が向上したり、緊急事態宣言によって人の移動が激減したことだけではないでしょう。
おそらく、新型コロナウィルス感染者による病床数の逼迫や、通院・入院による新型コロナウィルスへの感染の恐れといった国民の心理的要因によって、本来はインフルエンザに感染しているものの自宅療養している方も一定数いると思われます。
今回はインフルエンザ感染者数の時系列データのみの単純な効果検証でしたが、他の要因を考慮した多変量解析によって、より詳細な効果検証が可能になります。
