はじめに

最初は最近話題(？)のCausalImpactについて書こうと思っていたのですが、その基礎となる差分の差分法（difference in difference）についての知識が不足していたため、この記事では差分の差分法について試したことを記述していきます。

差分の差分法とは

書いたことは下記資料の抜粋なので、より確実な情報が得たい方は下記の資料を参照してください。

参考資料：

以下の図を用いて説明します。

f:id:rmizutaa:20190501195352p:plain

出典：https://ja.wikipedia.org/wiki/%E5%B7%AE%E5%88%86%E3%81%AE%E5%B7%AE%E5%88%86%E6%B3%95

Pが処置群、Sが対照群で、それぞれの点は集団を示します。例えばTime1でP群のみに投薬を行い、Time2でその結果を検証するというケースがイメージしやすいかと思います。差分の差分法ではPとSが平行トレンドであることが条件となるため、Time2において仮にP1に処置が行われなかったとしたらその値はQになると想定されます。その仮定を用いることで、Time2におけるP2-Qを処置効果と捉えることができるようになります。

仕組みとしては結構簡単なのですが、差分の差分法を使う条件の処置群と対照群が平行トレンドであることが割とネックになります。この条件が結構厳しく、実験としてデザインされておらず気軽にデータを取り直せない環境下で処置群に対して平行トレンドをもつ対照群をどう選択するの？ということがこの手法を用いる上での課題になると思われます。(実施する可能性があるならデータ取得の部分にも関われということかもしれませんが)

計算するには下記の式で回帰分析を行えば良いみたいです。

$y= \beta _ {0}+ \beta _ {1} T+ \beta _ {2} S+\beta _ {3}(T \cdot S) + \beta _ {4} X +\varepsilon$

ここでTが処置群フラグ、Sが期間フラグ(図のTime1と2の違い)、T・Sは処置群*期間フラグで、この係数 $\beta _ {3}$ が求めたい処置効果になります。Xはその他の共変量で、ここを調整することで平行トレンド仮定を満たすような形に近づけることができるようです。

実験

kaggleのRecruit Restaurant Visitor Forecastingのレストランの来店人数のデータを使います。

ある店舗では立地的にお盆は来店人数が増えるという仮説があり、差分の差分法を用いてそうでない店舗と比較してどれくらいの来店人数の向上がみられるかを確認したいという課題を考えます。

とりあえずこの仮定に合いそうな2店舗を抽出します。期間はざっくり2016/5/1~8/20とします。抽出した2店舗の来客人数推移をプロットしたものを下記に示します。

f:id:rmizutaa:20190501195425p:plain

ここでstore1が処置群、store2が対照群で、store1ではお盆による集客効果があるように見えます。

データを回帰分析にかけられる形にするために前処理を行います。ここで、データに明らかな週の周期性があるためその情報も加えます。

#データ読み込み
df=pd.read_csv("air_visit_data.csv",parse_dates=["visit_date"])

#元のデータは欠損があるため日時のベースを作成
basedate=pd.DataFrame(list(pd.date_range('2016-05-01', '2016-08-20')))
basedate.columns=["visit_date"]

#店舗を指定、日時を限定
store1=df[df["air_store_id"]=="air_e55abd740f93ecc4"]
store1=store1[(store1["visit_date"]>="2016-05") & (store1["visit_date"]<="2016-08-20")]
store1["visitors"]=store1["visitors"].interpolate() #欠損を線形補間
store2=df[df["air_store_id"]=="air_1653a6c513865af3"]
store2=store2[(store2["visit_date"]>="2016-05") & (store2["visit_date"]<="2016-08-20")]
store2["visitors"]=store2["visitors"].interpolate() #欠損を線形補間

df=pd.concat([store1,store2])

#store1が処置群
store1["treatment"]=1
store2["treatment"]=0

#お盆(8/11-8/20)が処置期間
df["postperiod"]=(df["visit_date"]>="2016-08-11").astype(int)

#処置効果ダミー
df["treatmentXpostperiod"]=df["treatment"]*df["postperiod"]

#曜日ダミーの追加
weekdays=pd.get_dummies(df["visit_date"].dt.weekday).rename(columns=lambda x:f"weekday_{x}")
df=pd.concat([df,weekdays],axis=1)

この結果データはこのような形になります。

f:id:rmizutaa:20190501195445p:plain

線形回帰を行います。

import statsmodels.api as sm
X=df[['treatment', 'postperiod','treatmentXpostperiod','weekday_0', 'weekday_1', 'weekday_2',
       'weekday_3', 'weekday_4', 'weekday_5', 'weekday_6']]
X["coef"]=1 #切片
y=df["visitors"]
mod = sm.OLS(y,X)
res = mod.fit()
print(res.summary())

f:id:rmizutaa:20190501195534p:plain

処置効果であるtreatmentXpostperiodの係数は13.5245と比較的それっぽい値が出ました。立地的にこの店舗ではお盆の効果により14人ほど来店する人数が増えるようです。調整済み決定係数が0.355と当てはまりは悪いですが、係数のp値は0.05以下なのでそれなりの確からしさはあるかと思います。

平行トレンド仮定に基づくなら変数に使用した曜日と時間情報が同じであれば2店舗のトレンドは同様になるはずですが、平行になるのは数ではなく率ではないかという疑いもあり、そこの厳密性は少し欠けるかと思います。

また今回は時系列データを取り扱ったのですが、8/10以前と8/11以後という分け方しておらず、細かい前後関係を見られていないのは少し気になるところです。