差分の差分法（DID）とCasual Impact

# 分析用にEcdatからCigarデータを取得
install.packages('Ecdat')
library('Ecdat')

# 確認用
data(Cigar)

# データの出力
out_path = '../02_work/01_dataset/01_cigar.csv'
write.csv(Cigar, out_path, row.names = FALSE)

from google.colab import drive
drive.mount('/content/drive')

# 必要なパッケージのインストール
!pip install tfcausalimpact pandas==1.5.0 numpy
# matplotlibでの日本誤文字化け対策用
!pip install japanize-matplotlib

import pandas as pd
from causalimpact import CausalImpact
import matplotlib.pyplot as plt
import japanize_matplotlib

# データの読み込み
file_path = '../02_work/01_dataset/01_cigar.csv'
df_raw = pd.read_csv(file_path)

# 1988年以降にたばこ税が50セント以上上昇した州を除外
skip_state = [3,9,10,22,21,23,31,33,48]
df = df_raw[(~df_raw['state'].isin(skip_state)) & (df_raw['year'] >= 70)]
df = df.reset_index(drop=True)

# 目的変数
# カリフォルニア州の売上のみを抜き出す
y = df[df.state == 5][['year', 'sales']]
y = y.set_index('year')
y.columns = ['y']

# 説明変数
# カリフォルニア州以外の州の売上データを抜き出す
df[df.state != 5][['year', 'state', 'sales']]

# 説明変数について、売上が州ごとに１カラムとなるように横持ち変換
df_other_usecol = df[df['state'] != 5][['year', 'state', 'sales']]
X = pd.pivot_table(df_other_usecol, values='sales', index='year', columns='state')
# カラム名を州番号だけではなくPrefixをつける
X = X.add_prefix('STATE_')

df_ci = pd.concat([y, X], axis=1)

# #########################################
# 州別売上推移の可視化
# #########################################
fig, ax = plt.subplots(figsize=(10, 6))
columns = df_ci.columns
for column in columns:
    ax.plot(df_ci.index, df_ci[column], label=column)
# 88年に破線のラインを追加
ax.axvline(x=88, color='gray', linestyle='--', label='intervention')
title = '州別たばこ売上推移'
x_lab = '年'
y_lab = '売上'
ax.set_title(title)
ax.set_xlabel(x_lab)
ax.set_ylabel(y_lab)
ax.legend(loc='upper left', bbox_to_anchor=(1, 1), ncol=2)

# 画像の出力
save_path = '../02_work/02_CI実行_CA禁煙キャンペーン効果検証/02_州別売上推移.png'
plt.savefig(save_path, dpi=150, bbox_inches='tight')

plt.show()

# #########################################
# Casual Impact
# #########################################
# indexの年を削除して0からの通し番号にする
# （pre_period, post_periodの行番号が施策実施前/後期間のデータに合致するように）
df_ci = df_ci.reset_index().drop(columns=['year'])

# 介入前期間と介入後期間をindex番号で指定
pre_period = [0, 17]
post_period = [18, 22]

# 施策実施前/後の期間をindex番号で指定してCasual Impact実施
ci = CausalImpact(df_ci, pre_period, post_period, 
                  model_args={'fit_method': 'hmc'})
ci.plot()

ci_summary = ci.summary()
print(ci_summary)
# Posterior Inference {Causal Impact}
#                           Average            Cumulative
# Actual                    77.3               386.5
# Prediction (s.d.)         86.91 (7.45)       434.55 (37.24)
# 95% CI                    [69.79, 98.99]     [348.96, 494.94]

# Absolute effect (s.d.)    -9.61 (7.45)       -48.05 (37.24)
# 95% CI                    [-21.69, 7.51]     [-108.44, 37.54]

# Relative effect (s.d.)    -11.06% (8.57%)    -11.06% (8.57%)
# 95% CI                    [-24.95%, 8.64%]   [-24.95%, 8.64%]

# Posterior tail-area probability p: 0.09
# Posterior prob. of a causal effect: 90.61%

# For more details run the command: print(impact.summary('report'))

ci_report = ci.summary('report')
print(ci_report)

#  Analysis report {CausalImpact}

# During the post-intervention period, the response variable had
# an average value of approx. 77.3. In the absence of an
# intervention, we would have expected an average response of 86.91.
# The 95% interval of this counterfactual prediction is [69.79, 98.99].
# Subtracting this prediction from the observed response yields
# an estimate of the causal effect the intervention had on the
# response variable. This effect is -9.61 with a 95% interval of
# [-21.69, 7.51]. For a discussion of the significance of this effect,
# see below.


# Summing up the individual data points during the post-intervention
# period (which can only sometimes be meaningfully interpreted), the
# response variable had an overall value of 386.5.
# Had the intervention not taken place, we would have expected
# a sum of 434.55. The 95% interval of this prediction is [348.96, 494.94].


# The above results are given in terms of absolute numbers. In relative
# terms, the response variable showed a decrease of -11.06%. The 95%
# interval of this percentage is [-24.95%, 8.64%].
# This means that, although it may look as though the intervention has
# exerted a negative effect on the response variable when considering
# the intervention period as a whole, this effect is not statistically
# significant and so cannot be meaningfully interpreted.


# The apparent effect could be the result of random fluctuations that
# are unrelated to the intervention. This is often the case when the
# intervention period is very long and includes much of the time when
# the effect has already worn off. It can also be the case when the
# intervention period is too short to distinguish the signal from the
# noise. Finally, failing to find a significant effect can happen when
# there are not enough control variables or when these variables do not
# correlate well with the response variable during the learning period.


# The probability of obtaining this effect by chance is p = 9.39%.
# This means the effect may be spurious and would generally not be
# considered statistically significant.