新型コロナ新規陽性者数と都営地下鉄の利用者数の関連分析
１．目的
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")
案の定、嘘くさい予測になりました＿|￣|○
新型コロナ観察記2020/7/11
１．最新データに基づくアップデート
最新のデータ（7/10更新）に基づいて、時系列分析データをアップデートしました。
２．方法
引き続きstatsmodelsを使用します。
東京都の新型コロナウィルス時系列データを時系列分析してみる
１．背景
時系列分析の勉強をしている最中に新型コロナウィルスが発生したので、ホットな話題ということもあり、時系列分析をぶつけてみた次第です。
２．方法
・使用データは東洋経済オンラインのサイトから取得しました。
toyokeizai.net
・言語はPython。時系列分析はstatsmodelsを使用します。tsa.seasonal_decomposeによるいわゆる乗法モデルベースの基本成分分解ですね。
３．結果と考察
まずはパッケージとデータのインポートから。
import statsmodels.api as sm import pandas as pd import numpy as np import requests import io from matplotlib import pylab as plt %matplotlib inline # グラフを横長にする from matplotlib.pylab import rcParams rcParams['figure.figsize'] = 15, 6 import warnings warnings.simplefilter('ignore') df = pd.read_csv('prefectures.csv')
データの中身。とりあえず、目的である東京だけ見てみる。
tokyo = df[df['prefectureNameJ']=='東京都'] tokyo.tail() year month date prefectureNameJ prefectureNameE testedPositive peopleTested discharged deaths 4524 2020 6 15 東京都 Tokyo 5592 16351.0 4887.0 314.0 4571 2020 6 16 東京都 Tokyo 5619 16351.0 4936.0 316.0 4618 2020 6 17 東京都 Tokyo 5633 61112.0 4978.0 317.0 4665 2020 6 18 東京都 Tokyo 5674 63083.0 5017.0 317.0 4712 2020 6 19 東京都 Tokyo 5709 64844.0 5071.0 320.0
新規感染者数を分析したいのですが、累積の感染者数しかありません。
なので、手作業で新規感染者数を作りにいきます。
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] plt.figure(figsize=(15, 8)) plt.plot(tokyo.NewPositive)
よくみるグラフが出てきました。ぱっと見、上昇トレンドかなと。
次に自己相関係数と偏自己相関係数を見てみましょう。
fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(tokyo.NewPositive, lags=30, ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(tokyo.NewPositive, lags=30, ax=ax2)
15日前との偏自己相関が大きい。しかもマイナスの偏自己相関になりました。
ということは約2週間前の結果が小さくなると、当日の結果が大きくなる関係にあります。
人々はある日の新規感染者数を見て警戒度を変えるので上記のような関係になるということでしょうか。
すなわち、2週間前の新規感染者数が小さい→油断する→当日の新規感染者数が増える、ということなのかもしれないです。
ただし、トレンドがある場合に偏自己相関係数をそのまま解釈してよいのか、勉強不足故、よくわかりません。
トレンドがある時系列なので、周期性を観察するために１階差分をとって自己相関係数をみてみます。
diff = tokyo.NewPositive.diff() diff = diff.dropna() fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(diff, lags=30, ax=ax1)
まあ、7日ですかね。
次に、成分分解をしてみます。
fig, (ax1,ax2,ax3) = plt.subplots(3,1, figsize=(15,8)) seasonal_decompose_res = sm.tsa.seasonal_decompose(tokyo.NewPositive, freq=7) seasonal_decompose_res.trend.plot(ax=ax1) seasonal_decompose_res.seasonal.plot(ax=ax2) seasonal_decompose_res.resid.plot(ax=ax3)
上昇トレンドにみえますね。
特段何もしなければ、今後も上昇トレンドは続いていくのか、注視が必要でしょう。
