【因果推断】优惠券政策对不同店铺的影响

这次依然是用之前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.nan
    year = 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以作DID

data["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来衡量，笔者有空也会对此进行学习。