ビジネスで利用する因果推論実践 ~アンチパターンを添えて~

こんにちは。次世代システム研究室のY. O.です。

前回の記事では、筆者が共に働いている、スマホアプリを運営するグループ会社のチームで取り入れている因果推論手法について紹介しました。
今回は、ABテストのようなコントロールされた環境下で施策を実施することができないケースにおいて利用される種々の因果推論手法について、実際にデータを用いながら結果を考察していきたいと思います。また、それぞれの手法で陥ってはいけないアンチパターンの再現も行ってみました。

1. 本記事で行うデータ解析の概要

ビジネスケース

• 本施策の目的
  
  • ユーザーの愛着度を上げるキャンペーン施策を行い、アプリを開く日数を多くする
• 本施策の実施方法
  
  • 愛着度が低いユーザーほどキャンペーンが適用される確率が高くなるように、全ユーザーの20%にキャンペーンを適用する
• 本施策の測定指標
  
  • ユーザーあたりのアプリを開く日数

データ解析のオーバービュー

1. データ生成
   
   • 今回のケースに合致するデータをkaggleで数時間探したが、見つからなかった（マーケティング系のデータって全然ない…）
   • ユーザーがアプリを開くかどうかを定式化し、その式をもとにデータを生成した
2. 複数手法による施策評価
   
   • 差の差分法・IPW推定・傾向スコアマッチングを用いて、生成データ内の真のキャンペーン効果を正しく測定できるかを評価
   • それぞれの評価において、各種手法の弱点をつくような考察を実施

2. 実験データ作成

import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler

import matplotlib.pyplot as plt
import seaborn as sns

def generate_engagement_dataset(n_users=10000):
    '''
    ユーザーの特徴量とengagementを生成する
    '''
    # データ生成時に利用する、特徴量間の相関係数を設定する
    cor_12 = 0.4
    cor_13 = -0.5
    cor_14 = 0.3

    # engagementを定義する真の回帰係数を設定する
    b_0 = 0
    b_1 = -2
    b_2 = -1.7
    b_3 = 1.4
    b_4 = 1

    # engagementにおける変動係数の定義
    # （重畳させるガウシアンノイズの標準偏差が、engagementの平均値の何倍になるかを定義する）
    eps = 0.3


    user_feat = pd.DataFrame()
    user_id = [("0"*len(str(n_users))+str(i+1))[-len(str(n_users)):] for i in range(n_users)]
    user_feat["user_id"] = user_id

    seed_1 = np.random.beta(2, 2, n_users)
    seed_2 = np.random.beta(1, 5, n_users)
    seed_3 = np.random.beta(3, 2, n_users)
    seed_4 = np.random.beta(5, 3, n_users)
    x_1 = seed_1
    x_2 = seed_2 + np.sqrt(1 - cor_12**2) * x_1
    x_3 = seed_3 + np.sqrt(1 - cor_13**2) * x_1
    x_4 = seed_4 + np.sqrt(1 - cor_14**2) * x_1

    scaler = MinMaxScaler()
    user_feat["x1"] = scaler.fit_transform(x_1.reshape(-1, 1))
    user_feat["x2"] = scaler.fit_transform(x_2.reshape(-1, 1))
    user_feat["x3"] = scaler.fit_transform(x_3.reshape(-1, 1))
    user_feat["x4"] = scaler.fit_transform(x_4.reshape(-1, 1))

    user_feat["eng_baseline"] = b_0 + b_1 * user_feat["x1"] + b_2 * user_feat["x2"] + b_3 * user_feat["x3"] + b_4 * user_feat["x4"]

    eng_mean = user_feat["eng_baseline"].mean()
    user_feat["eng_baseline"] += np.random.randn(n_users) * eps*eng_mean
    user_feat["eng_baseline"] = 1/(1+np.exp(-1 * user_feat["eng_baseline"]))

    return user_feat

user_feat = generate_engagement_dataset(n_users=10000)

print(user_feat.shape)
display(user_feat.head())

sns.pairplot(user_feat[["x1","x2","x3","x4","eng_baseline"]])
plt.show();

def apply_treatment(e, mean_engagement):
    apply_ratio = 0.2
    assert mean_engagement>=apply_ratio, "error"

    apply_prob = (1-e) * apply_ratio / (1 - mean_engagement)
    apply_flg = np.random.binomial(1, apply_prob)

    return apply_flg

effect = 0.05
mean_engagement = user_feat["eng_baseline"].mean()
user_feat["treated"] = user_feat["eng_baseline"].apply(lambda x : apply_treatment(x, mean_engagement))
user_feat["eng_treatment"] = np.clip(user_feat["eng_baseline"] + effect * user_feat["treated"], 0, 1)

print(user_feat.shape)
display(user_feat.head(10))
print("施策対象ユーザー数：\t{:,}".format((user_feat["eng_baseline"] != user_feat["eng_treatment"]).sum()))
print("施策対象外ユーザー 平均engagement：\t{:.4f}".format((user_feat.loc[user_feat["eng_baseline"] == user_feat["eng_treatment"], "eng_baseline"]).mean()))
print("施策対象ユーザー 平均engagement：\t{:.4f}".format((user_feat.loc[user_feat["eng_baseline"] != user_feat["eng_treatment"], "eng_baseline"]).mean()))
print("施策により上昇した平均engagement：\t{:.4f}".format((user_feat.loc[user_feat["treated"] == 1, "eng_treatment"] - user_feat.loc[user_feat["treated"] == 1, "eng_baseline"]).mean()))ね

L = 60
a = 0.2
e_1 = 0.0

df_days = pd.DataFrame(np.arange(L), columns=["date"])
df_days["join_key"] = 1
df_days["seasonal"] = ((np.cos(2 * np.pi / L * df_days["date"]) - 1) / 2 + np.random.randn(len(df_days)) * e_1) * a + 1

user_history = user_feat[["user_id", "eng_baseline", "eng_treatment"]].copy()
user_history["join_key"] = 1

user_history = pd.merge(df_days, user_history, on='join_key', how='outer')
user_history.drop(columns=["join_key"], inplace=True)

user_history["before_after"] = 0
user_history.loc[user_history["date"]>=L/2, "before_after"] = 1

user_history["eng"] = user_history["eng_baseline"]
user_history.loc[user_history["before_after"]>=1, "eng"] = user_history.loc[user_history["before_after"]>=1, "eng_treatment"]

user_history["prob"] = np.clip(user_history["eng"] * user_history["seasonal"], 0, 1)
user_history["observe"] = np.random.binomial(1, user_history["prob"])

user_history = pd.merge(user_history, user_feat[["user_id", "treated"]], on="user_id", how="left")

x = user_history.groupby("date")[["observe"]].mean()
x.plot()
plt.show();

3. 差の差分法（DID）

基本形

data_0 = user_history[user_history["treated"]==0].groupby("date")[["observe"]].mean()
data_0.reset_index(drop=False, inplace=True)
data_1 = user_history[user_history["treated"]==1].groupby("date")[["observe"]].mean()
data_1.reset_index(drop=False, inplace=True)

plt.plot(data_0["date"], data_0["observe"], label="control")
plt.plot(data_1["date"], data_1["observe"], label="treated")
plt.xlabel("date")
plt.legend()
plt.show();

did_result = user_history.groupby(["treated", "before_after"])[["observe"]].agg(["count","mean","var"])
display(did_result)

アンチパターン

user_history["prob"] = np.clip(user_history["eng"] * user_history["seasonal"] - user_history["date"] * user_history["eng_baseline"]*0.2/L, 0, 1)

4. IPW推定

基本形

import statsmodels.api as sm

# ロジスティック回帰モデルの作成
df = user_feat[["x1","x2","x3","x4","treated"]].copy()
df = sm.add_constant(df)

logit_mod = sm.Logit(df[["treated"]], df[["const","x1","x2","x3","x4"]], family=sm.families.Gamma())
logit_res = logit_mod.fit(method="newton")
print(logit_res.summary())


# モデルの出力結果を傾向スコアとしてテーブル化
res = pd.DataFrame()
res["user_id"] = user_feat["user_id"]
res["propensity"] = logit_mod.endog - logit_res.resid_response

ipw_data = user_history[["user_id", "before_after", "treated", "observe"]].copy()
ipw_data = ipw_data[ipw_data["before_after"]==1]
ipw_data = pd.merge(ipw_data, res, on="user_id", how="inner")

# 定義に従ってIPW推定量を計算
cond_e1 = ipw_data["treated"]==1
E_1 = ((1 / ipw_data.loc[cond_e1, "propensity"])*ipw_data.loc[cond_e1, "observe"]).sum() / (1 / ipw_data.loc[cond_e1, "propensity"]).sum()

cond_e0 = ipw_data["treated"]==0
E_0 = ((1 / (1 - ipw_data.loc[cond_e0, "propensity"]))*ipw_data.loc[cond_e0, "observe"]).sum() / (1 / (1 - ipw_data.loc[cond_e0, "propensity"])).sum()

ipw_data["residual"] = 0
ipw_data.loc[ipw_data["observe"]==1, "square_error"] = ipw_data["propensity"] * np.power(ipw_data["observe"] - E_1, 2)
ipw_data.loc[ipw_data["observe"]==0, "square_error"] = ipw_data["propensity"] * np.power(ipw_data["observe"] - E_0, 2)

se = np.sqrt(
    ipw_data.loc[ipw_data["observe"]==1, "square_error"].sum()/np.power((ipw_data.loc[ipw_data["observe"]==1, "propensity"]).sum(), 2) + 
    ipw_data.loc[ipw_data["observe"]==0, "square_error"].sum()/np.power((ipw_data.loc[ipw_data["observe"]==0, "propensity"]).sum(), 2)
    )

print(E_1 - E_0)
print(se * 1.96)

アンチパターン

logit_mod = sm.Logit(df[["treated"]], df[["const","x1","x2","x3"]], family=sm.families.Gamma())

5. 最後に