【因果推斷】優惠券政策對不同店鋪的影響

這次依然是用之前rossmann店鋪競賽的數據集。
之前的數據集探索處理在這里已經做過了，此處就不再贅述了CSDN鏈接
數據集地址：競賽鏈接
這里探討數據集中Promo2對于每家店鋪銷售額的影響。其中，Promo2是一個基于優惠券的郵寄活動，發送給參與商店的顧客。每封信里都有幾張優惠券，大部分是所有產品的一般折扣，有效期為三個月。所以在這些優惠券到期之前，我們會給客戶發新一輪的郵件。具體的變量解釋可以看這里：Promo2釋義
此處單純是為了作因果推斷的練習。由于筆者是因果推斷的初學者，有些概念與處理可能會有疏漏錯誤，還請發現的讀者不吝賜教。

一、數據可視化

首先，我們從單純地數據可視化的角度來看總體均值以及店鋪自身使用優惠券政策前后的比較。

import warnings
warnings.filterwarnings("ignore")import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
data = pd.read_csv("train.csv")
data = data[data["Open"]==1]
store_info = pd.read_csv("store.csv")#假設一年的第一周是包含1月4日的那周
data = data[data["Sales"]>0]
data["Date"]=pd.to_datetime(data["Date"])
data["Year"]=[i.year for i in data["Date"]]
data["Month"]=[i.month for i in data["Date"]]
def get_first_day_of_week(year, week):if np.isnan(year):return np.nanyear = int(year)week = int(week)first_thursday = datetime(year, 1, 4)  while first_thursday.weekday() != 3:  # weekday() returns 0 (Monday) through 6 (Sunday)  first_thursday += timedelta(days=1)  first_day_of_week = first_thursday - timedelta(days=first_thursday.weekday() - 1)  target_date = first_day_of_week + timedelta(weeks=week - 1)  return datetime.date(target_date)
res = []
for i in range(store_info.shape[0]):res.append(get_first_day_of_week(store_info.loc[i,"Promo2SinceYear"],store_info.loc[i,"Promo2SinceWeek"]))
store_info["Promo2StartTime"] = pd.to_datetime(res)
data = data.merge(store_info,on="Store")
data["Post"] = (data["Date"]>data["Promo2StartTime"]).astype("int8")
data["StateHoliday"] = (data["StateHoliday"]!="0").astype("int8")
## 先簡單比較一下每天·店的均值
#均值比較
Promo2_date=data[data["Post"]==1].groupby(["Year","Month"])["Sales"].mean()
NoPromo2_date=data[data["Post"]==0].groupby(["Year","Month"])["Sales"].mean()
Promo2_date.index = ["-".join([str(j) for j in i]) for i in Promo2_date.index]
NoPromo2_date.index = ["-".join([str(j) for j in i]) for i in NoPromo2_date.index]
fig, ax = plt.subplots()
plt.plot(Promo2_date,label="has Promo2")
plt.plot(NoPromo2_date,label="does not have Promo2")
ax.xaxis.set_visible(False)
plt.legend()
plt.show()

在這里插入圖片描述
圖中藍線代表實行了優惠券政策的店鋪每個月的月均銷量，而橙線則代表控制組的月均銷量，可以看出控制組的銷量明顯高于優惠券政策的店鋪。

## 去除那些沒有“促銷前”信息的店鋪
tmp = data[~pd.isna(data["Promo2StartTime"])].groupby(["Store"])[["Promo2StartTime","Date"]].min()
drop_stores = list(tmp[tmp["Date"]>=tmp["Promo2StartTime"]].index)
for i in drop_stores:data = data[data["Store"]!=i]
data.index = range(data.shape[0])
data_promo2 = data[~pd.isna(data["Promo2StartTime"])].reset_index(drop=True)
isin_thirty_days=[abs(i).days<=30 for i in data_promo2["Date"]-data_promo2["Promo2StartTime"]]
data_promo2 = data_promo2[isin_thirty_days]
promo2_before=data_promo2[data_promo2["Date"]<data_promo2["Promo2StartTime"]].groupby(["Store"])["Sales"].mean()
promo2_after=data_promo2[data_promo2["Date"]>=data_promo2["Promo2StartTime"]].groupby(["Store"])["Sales"].mean()
fig, ax = plt.subplots(figsize=(10,10))
num_bars = 15
idx = range(num_bars)
plt.bar(idx,promo2_before[0:num_bars].values,width=0.35,label="before Promo2")
plt.bar([i+0.35 for i in idx],promo2_after[0:num_bars].values,width=0.35,label="after Promo2")
ax.xaxis.set_visible(False)
plt.legend()
plt.show()

在這里插入圖片描述
選取部分店鋪在發行優惠券前后的銷售額均值對比，發現大多數情況下，銷售額會隨著降低，不過也有部分店鋪在發行優惠券后銷售額反而提升了。

## 作圖平行趨勢
promo2_bydate = data.groupby(["Promo2StartTime","Year","Month"])["Sales"].mean().reset_index()
nopromo2 = data[pd.isna(data["Promo2StartTime"])].groupby(["Year","Month"])["Sales"].mean().reset_index()
nopromo2.index = [str(i)+"-"+str(j) for i,j in zip(nopromo2["Year"],nopromo2["Month"])]
fig, ax = plt.subplots()
sns.lineplot(nopromo2.sort_values(["Year","Month"])["Sales"],label="no_promo")
for i in promo2_bydate["Promo2StartTime"].unique()[[0,10,20,15]]:tmp = promo2_bydate[promo2_bydate["Promo2StartTime"]==i].sort_values(["Year","Month"]).reset_index()tmp.index = [str(i)+"-"+str(j) for i,j in zip(tmp["Year"],tmp["Month"])]sns.lineplot(tmp["Sales"],label=str(i)[0:10])ax.vlines(x=str(i)[0:7].replace("-0","-"),ymin=50,ymax=9000,ls="dashed")ax.xaxis.set_visible(False)

在這里插入圖片描述
選取幾個不同時間開始發行優惠券的店鋪的歷史銷售作圖，可以發現大多數發行優惠券店鋪的整體歷史銷售和不發行優惠券的那店鋪趨勢大致相同。

二、分組雙重差分

從計量經濟學的角度來看，本身我們獲得的數據是一份面板數據（Panel Data），所以可以使用雙重差分（DID）的方法來判斷優惠券政策是否顯著影響了店鋪的銷售額。但現在有一個問題，我們每一家店鋪的優惠券政策（干預）實際上并非同時發生的。這并不符合傳統DID的假設，因為對于不同的店鋪而言，優惠券政策的影響效果是不同的。因此我們需要使用一個更為靈活的模型考慮。我們期望對于不同的時間而言，每個店鋪有不同的效應。為了保證不出現梯度爆炸，我將時間分組為每年每周，而由于這一數據集中，干預或許并非是隨機發生的，因此需要將其他的混淆因子也加入到模型中，能夠使得模型解釋性更好。

import statsmodels.formula.api as smf
## 按照群組分組進行OLS以作DIDdata["Promo2StartTime"] = data["Promo2StartTime"].fillna(pd.to_datetime("2100-01-01"))
data["Post"] = (data["Date"]>data["Promo2StartTime"]).astype("int8")
data["StateHoliday"] = (data["StateHoliday"]!="0").astype("int8")
data["WeekOfYear"] = [i.weekofyear for i in data["Date"]]
data_for_ols = data.groupby(["Store","Year","WeekOfYear"]).agg({"Sales":"sum","Promo":"sum","StateHoliday":"sum","SchoolHoliday":"sum","Promo2":"sum","Post":"sum"}).reset_index()
data_for_ols = data_for_ols.merge(store_info[["Store","StoreType","Assortment","CompetitionDistance","Promo2StartTime"]],on="Store")
boolize = ["Promo","StateHoliday","SchoolHoliday","Promo2","Post"]
for col in boolize:data_for_ols[col] = pd.Series([1 if i>0 else 0 for i in data_for_ols[col]]).astype("int8")data_for_ols["Sales"] = np.log(data_for_ols["Sales"])
formular = "Sales~Promo2:Post:C(Promo2StartTime):C(Year):C(WeekOfYear)+Promo+StateHoliday+SchoolHoliday+StoreType+Assortment+CompetitionDistance"
twfe_model = smf.ols(formular,data=data_for_ols).fit()
df_predict = data_for_ols[data_for_ols["Post"]*data_for_ols["Promo2"]==1].reset_index()
df_predict["Promo2"] = 0
df_predict["predict"] = twfe_model.predict(df_predict)
print("參數個數：",len(twfe_model.params))
print("估計的ATT：",(df_predict["Sales"]-df_predict["predict"]).mean())

在這里插入圖片描述
最終我們使用OLS最小二乘建模，獲得了一個有3754個參數的回歸模型。對于銷售額而言，平均干預效果為負數。

三、干預異質性檢驗

那么對于不同的店鋪而言，優惠券政策是否有不同影響？為了研究這一問題，筆者使用了元機器學習（Meta-Learner）的技巧，構建了一個X-Learner：
在這里插入圖片描述
X-Learner分為兩個階段：第一個階段訓練1個傾向得分模型（也就是圖中最左側，根據協變量X來預測T的模型），同時訓練2個模型回應(respond)模型（也就是圖中左側2個并列的模型），分別將受到干預和沒收到干預的部分分開使用協變量預測因變量的模型。到了第二階段，我們使用第一階段的回應模型預測出兩個反事實結果：對于沒有受到干預的店鋪，如果當初受到干預的話會有多少銷售額（圖中的 $Y^*|T=0$ ）以及對于受到了干預的店鋪如果當初沒有受到干預的話會有多少銷售額（圖中的 $Y^*|T=1$ ）。將 $Y 1$ 和 $Y 0$ 相見，我們就獲得了條件平均處理效果(CATE)。我們將CATE作為應變量，再使用協變量X訓練2個模型（也就是圖中帶有綠色CATE的部分）。使用這2個模型按照傾向得分倒數的權重預測最終的CATE。（因為傾向得分越高表示越有可能受到干預，倒數反而代表對應樣本更稀有，對于非隨機實驗而言參考價值更大）

## X-learner## PS Model
Y_PS = "treatment"
train = pd.concat([train_T1,train_T0],axis=0)
train_ps = pd.get_dummies(train)
X_PS = [i for i in train_ps.columns if i != "Sales" and i != "treatment"]
ps_model = LogisticRegression(penalty='none')
ps_model.fit(train_ps[X_PS],train_ps[Y_PS])## Predict_Model
X_lgbm = [i for i in data_t1.columns if i != "Sales" and i != "treatment"]
Y = "Sales"
model_t1 = LGBMRegressor()
model_t0 = LGBMRegressor()model_t0.fit(train_T0[X_lgbm],train_T0[Y],sample_weight=1/ps_model.predict_proba(pd.get_dummies(train_T0)[X_PS])[:, 0])
model_t1.fit(train_T1[X_lgbm],train_T1[Y],sample_weight=1/ps_model.predict_proba(pd.get_dummies(train_T1)[X_PS])[:, 1])## second_stage
tau_hat_0 = model_t1.predict(train_T0[X_lgbm])-train_T0[Y]
tau_hat_1 = train_T1[Y]-model_t0.predict(train_T1[X_lgbm])model_tau_t1 = LGBMRegressor()
model_tau_t0 = LGBMRegressor()model_tau_t0.fit(train_T0[X_lgbm],tau_hat_0)
model_tau_t1.fit(train_T1[X_lgbm],tau_hat_1)
# estimate the CATE
valid = pd.concat([valid_T1,valid_T0],axis=0)ps_valid = ps_model.predict_proba(pd.get_dummies(valid)[X_PS])[:, 1]cate_valid = valid.assign(cate=(ps_valid*model_tau_t0.predict(valid[X_lgbm]) +(1-ps_valid)*model_tau_t1.predict(valid[X_lgbm]))
)
sns.displot(cate_valid["cate"])

在這里插入圖片描述
可以看出，使用X-Learner預測的結果中，在驗證集大部分的店鋪在受到了干預之后，銷售量明顯是負數；只有小部分的店鋪銷售額會受到正面影響。
接下來將驗證集中各個CATE排序，求出累計的干預效果，以此做出QINI曲線。

sorted_cate = cate_valid["cate"].values[np.argsort(-cate_valid["cate"].values)]
cumulative_cate = np.cumsum(sorted_cate)/np.arange(1,len(sorted_cate)+1)
cumulative_fractions = np.arange(1, len(cumulative_cate) + 1) / len(cumulative_cate)
plt.plot(cumulative_fractions,cumulative_cate,"r",label="QINI Curve")
plt.plot([0,1],[cumulative_cate[0],cumulative_cate[-1]],"k--",label="Random Assignment")
plt.title("QINI Curve")
plt.legend()
plt.show()

在這里插入圖片描述
可以看到，在一開始，的確有些店鋪的銷售額是受到了優惠券政策的正向影響；但是到了之后，優惠券政策依然是負向影響。

四、總結

本文使用了競賽用的銷售數據進行了銷售額與優惠券政策的因果推斷。實際上更多的是因果推斷方法論的學習。對于優惠券政策而言，銷售額很有可能并非是真正的干預目標，而且店鋪是否要發現優惠券，也需要考量除了數據集以外的其他因素。而對于銷售額的干預影響，很多公司都會做Uplift Model來衡量，筆者有空也會對此進行學習。