Pythonで因果推論したい(CausalImpact)

タイトルの通り、CausalImpactをPythonで試してみます。 下記本を元にし、RからPythonで書き直し、同様の効果が推定できるか試してみます。

導入

pyが最初につくのに注意します。

pip install pycausalimpact

データの作成

statsmodelsのcigarデータを使います。

CausalImpactはデータフレームの一番左をyとし、それ以降が予測に使う変数となるため、その形となるようにpivotと並び替えを行います。

import statsmodels.api as sm  
df = sm.datasets.get_rdataset('Cigar', 'Ecdat').data

#上記本の初期加工を関数化しておく
def create_data(df):
  tmp = df[(df['year']>=70) & (~df['state'].isin([3,9,10,22,21,23,31,33,48])) ]
  tmp['area'] = tmp.state.apply(lambda x:'CA' if x== 5 else 'Rest of US')
  return tmp

use_df = create_data(df)

#ciに使える形に変換
pv = use_df.pivot(index='year',columns='state',values='sales')
cols = [c for c in pv.columns if c != 5]
cols.insert(0,5)
pv2=pv[cols]
pv2.columns = [str(c) for c in cols]

効果検証

from causalimpact import CausalImpact

#介入前、後の区分を用意
#88年~が介入効果がある、という前提
pre = [70,87]
post = [88,92]

ci = CausalImpact(pv2,pre,post)

#可視化
ci.plot()

#結果の数値表示
print(ci.summary())

#各変量の効果の表示
ci.trained_model.summary()

上から、

  • yの実数値と予測値
  • 予測に対しての各時点の差
  • 上記差の累積

がグラフになります。

f:id:esu-ko:20200813234940p:plain

また、summaryをすると、予測値,実測値,95%区間の平均/累積の値および、効果の実数値/相対比率が表示されます。

これは上記本と比較すると値、グラフともに異なっており、このライブラリが使っているstatsmodelsのカルマンフィルターは確定的であり、RのCausalImpactで用いる確率的なカルマンフィルターと違うため、らしいです。