新型コロナ新規陽性者数と都営地下鉄の利用者数の関連分析
１．目的
SIRモデル等では接触率を変化させることで、将来の状況をシミュレーションするわけですが、これは「接触率と新規感染者数には因果関係がある」という前提の話です。
定性的に考えれば、この前提は正しいとは思いますが、定量的に示している記事などを見たことがありません。
ということで、本稿の目的は当該前提をデータを使って定量的に考察してみる、です。
２．方法
時系列解析で外生変数の影響を分析する、という目的に合致しているモデルを探していたら、どうやら「状態空間モデル」がよさげなので、今回はこのモデルを用いて考察してみます。
状態空間モデルについては、以下のブログを参照しました。
logics-of-blue.com
また、各データの取得元は以下の通りです。
・新規感染者数
toyokeizai.net
・外生変数（都営地下鉄の利用者数の推移）
メインの出勤時間帯を7時30分～9時30分と考え、当該時間帯のデータを利用。
stopcovid19.metro.tokyo.lg.jp
３．結果と考察
言語はPythonで行ないますが、その前にエクセルで「都営地下鉄の利用者数の推移」を加工します。
当該データは週単位であるため、日単位にするために「線形補完→Grevilleの3次13項式」によりなめらかな曲線にします。
また、感染から発症まで約２週間のタイムラグが存在しているそうなので、２週間前数値を現在値として使用することにします。
計算式などは割愛しますが、結果は以下の通り。徐々に利用者数が増加してますね（当該数値は、ある時点を1とした時のrelativity）。
準備ができたので、Pythonで分析していきます。
# 基本のライブラリを読み込む import numpy as np import pandas as pd from scipy import stats # グラフ描画 from matplotlib import pylab as plt import seaborn as sns %matplotlib inline import warnings warnings.simplefilter('ignore') # グラフを横長にする from matplotlib.pylab import rcParams rcParams['figure.figsize'] = 15, 6 # 統計モデル import statsmodels.api as sm #data import df = pd.read_csv(prefectures.csv') #加工（東京都に限定する等） tokyo = df[df['prefectureNameJ']=='東京都'] tokyo['ymd'] = '' for i in range(len(tokyo)): tokyo.iloc[i,9] = str(tokyo.iloc[i,0]) + '/' + str(tokyo.iloc[i,1]) + '/' + str(tokyo.iloc[i,2]) tokyo.ymd = pd.to_datetime(tokyo.ymd) tokyo = tokyo.set_index("ymd") tokyo['testedPositive'] = tokyo['testedPositive'].astype('float64') tokyo['NewPositive'] = 0.0 for i in range(len(tokyo)-1): tokyo.iloc[i+1,9] = tokyo.iloc[i+1,5] - tokyo.iloc[i,5] tokyo['ymd_num'] = 0.0 for i in range(len(tokyo)): tokyo.iloc[i,10] = tokyo.iloc[i,0]*10000 + tokyo.iloc[i,1]*100 + tokyo.iloc[i,2] tokyo = tokyo[tokyo['ymd_num']>=20200317] # 日付形式にする ts = tokyo['NewPositive'] #Subway Passengerをインポート df2 = pd.read_csv(SubwayPassengers_input.csv') df2.ymd = pd.to_datetime(df2.ymd) df2 = df2.set_index("ymd")
ライブラリやデータの準備ができたので、早速モデリングをしてみます。
モデルとしては「季節変動あり+外生因子ありのローカル線形トレンドモデル」を採用しています。
# 季節変動あり+外生因子ありのローカル線形トレンドモデル mod_season_trend_x = sm.tsa.UnobservedComponents( ts, 'local linear trend', seasonal=7, exog=x.values ) # まずはNelder-Meadでパラメタを推定し、その結果を初期値としてまた最適化する。2回目はBFGSを使用。 res_season_trend_x = mod_season_trend_x.fit( method='bfgs', maxiter=500, start_params=mod_season_trend_x.fit(method='nm', maxiter=500).params, ) # 推定されたパラメタ一覧 print(res_season_trend_x.summary()) # 推定された状態・トレンド・季節の影響の描画 rcParams['figure.figsize'] = 15, 20 fig = res_season_trend_x.plot_components()
結果は以下の通り。
外生変数の係数は232.8916ですが、95%CIは0を含んでいます。
ということで、統計的には「何とも言えない」という結論なわけですが、係数はプラスなので、定性的な感覚とは合致する結果となりました。
もちろん「接触率」を「都営地下鉄の利用者数」だけで完全に表現することはできないわけで、今回のモデルでは接触率を表す要素は他の変数で（無理矢理）表現しようとしています（どの変数が頑張っているかはわかりませんが）。
Unobserved Components Results
====================================================================================
Dep. Variable: NewPositive No. Observations: 115
Model: local linear trend Log Likelihood -512.541
+ stochastic seasonal(7) AIC 1035.082
Date: Sat, 11 Jul 2020 BIC 1048.447
Time: 11:42:22 HQIC 1040.500
Sample: 03-17-2020
- 07-09-2020
Covariance Type: opg
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
sigma2.irregular 299.7816 70.591 4.247 0.000 161.426 438.137
sigma2.level 143.6784 77.113 1.863 0.062 -7.461 294.817
sigma2.trend 0.8645 1.179 0.733 0.464 -1.447 3.176
sigma2.seasonal 12.0896 9.912 1.220 0.223 -7.337 31.516
beta.x1 232.8916 185.786 1.254 0.210 -131.243 597.026
===================================================================================
Ljung-Box (Q): 46.12 Jarque-Bera (JB): 11.90
Prob(Q): 0.23 Prob(JB): 0.00
Heteroskedasticity (H): 0.67 Skew: 0.21
Prob(H) (two-sided): 0.23 Kurtosis: 4.58
===================================================================================
PCR検査の実施基準が変わっているっぽいので、将来予測はあまり意味がないのですが、状態空間モデルでは外生変数を用いて将来予測ができるようなので、念のために計算してみました。
#利用比率の値（7/31まで） df3 = pd.read_csv(SubwayPassengers_pred.csv') pred_x = df3.iloc[:,1].as_matrix().reshape(22,1) # 予測 pred = res_season_trend_x.predict('2020-03-17', '2020-07-31', exog=pred_x) pred.index = pred.index -1 # 実データと予測結果の図示 rcParams['figure.figsize'] = 15, 6 plt.plot(ts) plt.plot(pred, "r")
案の定、嘘くさい予測になりました＿|￣|○
