🐍

機械学習で因果推論~Meta-LearnerとEconML~

2022/07/28に公開約13,500字2件のコメント

はじめに

Meta-LearnerとEconMLについて、Pythonによる実装を交えてまとめました。内容について誤り等ございましたら、コメントにてご指摘いただけますと幸いです。

機械学習を用いた因果推論

機械学習を用いた因果推論手法は大きく分けて下記の2通りが存在します。

  • Meta-Learner系
  • Causal-Tree系

この中でも今回は1つ目のMeta-Learnerについて紹介します。

Meta-Leanrerとは

Meta-Learnerとは、機械学習と因果推論の考え方を掛け合わせて条件付き平均処置効果(CATE: Conditional Average Treatment Effect)を推定する手法の総称です。

条件付き平均処置効果(CATE)とは、平均処置効果(ATE: Average Treatment Effect)をある条件(X=x)に限定して算出したものであり、処置を実施した場合の結果をY(1)、処置を実施しなかった場合を結果Y(0)とすると

\hat{\tau}(x) = E[Y(1)-Y(0)|X=x] = E[Y(1)|X=x] - E[Y(0)|X=x]
で与えられます。X=xについて条件付けした場合の効果を推定するため、個体の属性における効果の異質性、すなわちxによって効果が異なるような非線形な効果も捉えることができます。

ここで、X=xの条件である処置を実施した場合、Y(1)|X=x \:のデータは取れますが、Y(0)|X=x \:のデータは取れません。逆にX=xの条件である処置を実施しなかった場合、Y(0)|X=x \:のデータは取れますが、Y(1)|X=x \:のデータは取れません。

そのため、機械学習を使って取得できないデータの予測値を算出し、推定に利用しようとするのがMeta-Learnerの考え方です。

Meta-Learnerには

  • T-Learner: Two Learner
  • S-Learner: Single Learner
  • X-Learner
  • DA-Learner: Domain Adaptation Learner
  • DR-Learner: Doubly Robust Learner

など、さまざまな手法があります。今回は1つ1つの手法のアルゴリズムには触れず、これから紹介するEconMLというパッケージを用いて、簡易的に条件付き平均処置効果(CATE)を算出する方法について説明します。

EconMLとは

EconMLとは、観察可能なデータから機械学習を用いて条件付き平均処置効果(CATE)を推定するPythonパッケージです。Microsoft社が開発しているもので、さまざまなMeta-Learner系やCausal-Tree系の因果推論手法が実装されています。

EconMLの公式ドキュメントはこちらです。

https://econml.azurewebsites.net/

Pythonによる実装

ダイエット商品の広告効果を例に、Pythonでデータを作成し、Meta-Learner手法を用いて効果を算出してみます。

設定

とある健康食品会社では、ECにてダイエット商品を販売しています。ダイエット商品の販促のため広告を回した結果から広告効果を推定することになりました。

手元には下記のデータがあります。

  • iさんのダイエットへの意識の高さを表す指標x_i
    • 一様分布(-1, 1)に従う
  • iさんが広告を見たかどうかを表すダミー変数T_i
    • 下記のモデルによって決定されるとする
      T_i = \left\{\begin{array}{ll} 1 & (x_i + noise_i > 0): 広告を見た \\ \\ 0 & (x_i + noise_i \leq 0): 広告を見てない \end{array}\right.
    • noise_iは標準正規分布に従う
  • iさんの購入金額Y_i
    • x_iT_iの影響を受け、真のモデルは下記のように表される
      Y_i = 1000(\: [\: 3 + \tau_i T_i + 3x_i + noise_i \:] \: )
    • ただし[x]は、n \leq x < n+1を満たす整数nを表す
    • また、効果\tau_ix_iによって異なる
      \tau_i = \left\{\begin{array}{ll} 1 & (x_i \leq 0) \\ \\ 2 & (0 < x_i \leq 0.5) \\ \\ 3 & (0.5 < x_i \leq 1) \end{array}\right.
    • noise_iは一様分布(0, 1)に従う

データの作成

Pythonで上記の設定を満たすデータを作成します。

# 必要なライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
# グラフをJupyter上に描画
%matplotlib inline

# データ数
size = 1000
# シードの設定
np.random.seed(0)

# ダイエットへの意識の高さ
x = np.random.uniform(-1, 1, size)
noise = np.random.randn(size)

# 広告ダミー
_T = x + noise
T = np.where(_T>0, 1, 0)

# ダイエットへの意識の高さによって広告効果が異なる
t = np.zeros(size)
for i in range(size):
    if x[i] < 0:
        t[i] = 1
    elif x[i] < 0.5:
        t[i] = 2
    else:
        t[i] = 3
        
# 売上
noise = np.random.uniform(0, 1, size)
Y = np.clip(t*T + 3*x + 3 + noise, 0, 10).astype("int") * 1000

# 売上のヒストグラムを描画
plt.hist(Y)
plt.xlabel("ダイエット商品の売上(単位:円)")
plt.show()

(出力結果)

Meta-Learnerによる効果検証

今回は下記5つのMeta-Learner手法で効果を推定します。

  • T-Learner: Two Learner
  • S-Learner: Single Learner
  • X-Learner
  • DA-Learner: Domain Adaptation Learner
  • DR-Learner: Doubly Robust Learner

T-Learner

T-Learnerは推定に利用する機械学習モデルを選択するだけで実装できます。
今回はランダムフォレスト回帰をモデルに選択します。

# 必要なライブラリのインポート
from sklearn.ensemble import RandomForestRegressor
from econml.metalearners import TLearner

# モデルの構築
models = RandomForestRegressor(max_depth=3, random_state=0)
T_learner = TLearner(models=models)
T_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = T_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

多少のブレはありますが、比較的よく効果を推定できているようです。

S-Learner

S-LearnerもT-Learnerと同様、推定に利用する機械学習モデルを選択するだけで実装できます。
今回はランダムフォレスト回帰をモデルに選択します。

# 必要なライブラリのインポート
from sklearn.ensemble import RandomForestRegressor
from econml.metalearners import SLearner

# モデルの構築
overall_model = RandomForestRegressor(max_depth=3)
S_learner = SLearner(overall_model=overall_model)
S_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = S_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

目感にはなりますが、T-Learnerよりも推定の精度が少し落ちたような気がします。

X-Learner

X-Learnerでは、売上Y(1), Y(0)を予測するモデルに加えて、傾向スコアを算出するためのモデルも選択します。今回は売上を予測するモデルにランダムフォレスト、傾向スコアを算出するモデルにロジスティック回帰を選択します。

# 必要なライブラリをインポート
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LogisticRegression
from econml.metalearners import XLearner

# モデルの構築
models = RandomForestRegressor(max_depth=3, random_state=0)
propensity_model = LogisticRegression()
X_learner = XLearner(models=models, propensity_model=propensity_model)
X_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = X_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

T-Learnerよりもかなり良く推定できているようです。

DA-Learner: Domain Adaptation Learner

DA-Learnerでは下記の3つのモデルを選択する必要があります。

  • データを予測するモデル
  • 傾向スコアを算出するモデル
  • 効果の推定に用いるモデル

今回は

  • データを予測するモデル: ランダムフォレスト
  • 傾向スコアを算出するモデル: ロジスティック回帰
  • 効果の推定に用いるモデル: ランダムフォレスト

を選択します。

# 必要なライブラリをインポート
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LogisticRegression
from econml.metalearners import DomainAdaptationLearner

# モデルの構築
models = RandomForestRegressor(max_depth=3, random_state=0)
final_models = RandomForestRegressor(max_depth=3, random_state=0)
propensity_model = LogisticRegression()
DA_learner = DomainAdaptationLearner(models=models, 
                                     final_models=final_models, 
                                     propensity_model=propensity_model)
DA_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = DA_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

こちらも非常によく推定できているようです。

DR-Learner: Doubly Robust Learner

DR-Learnerでは、DA-Learnerと同様に下記の3つのモデルを選択する必要があります。

  • データを予測するモデル
  • 傾向スコアを算出するモデル
  • 効果の推定に用いるモデル

今回は

  • データを予測するモデル: ランダムフォレスト
  • 傾向スコアを算出するモデル: ロジスティック回帰
  • 効果の推定に用いるモデル: ランダムフォレスト
# 必要なライブラリをインポート
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LogisticRegression
from econml.dr import DRLearner

# モデルの構築
outcome_model = RandomForestRegressor(max_depth=3, random_state=0)
pseudo_treatment_model = RandomForestRegressor(max_depth=3, random_state=0)
propensity_model = LogisticRegression()
DR_learner = DRLearner(model_regression=outcome_model, 
                       model_propensity=propensity_model, 
                       model_final=pseudo_treatment_model)
DR_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = DR_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

こちらも良く推定できているようです。

補足

先ほども少し触れましたが、現実社会においては広告効果を推定したい場合、広告を見た人もそうでない人もその多くが商品を購入しないというシーンがほとんどです。

先ほどの例で売上Yの真のモデルを

Y_i = 1000(\: [\: \tau_i T_i + 3x_i + noise_i \:] \: )
とし、Pythonで実装すると下記のようになります。

# 必要なライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
# グラフをJupyter上に描画
%matplotlib inline

# データ数
size = 1000
# シードの設定
np.random.seed(0)

# ダイエットへの意識の高さ
x = np.random.uniform(-1, 1, size)
noise = np.random.randn(size)

# 広告ダミー
_T = x + noise
T = np.where(_T>0, 1, 0)

# ダイエットへの意識の高さによって広告効果が異なる
t = np.zeros(size)
for i in range(size):
    if x[i] < 0:
        t[i] = 1
    elif x[i] < 0.5:
        t[i] = 2
    else:
        t[i] = 3
        
# 売上
noise = np.random.uniform(0, 1, size)
Y = np.clip(t*T + 3*x + noise, 0, 10).astype("int") * 1000
# Y = np.clip(t*T + 3*x + 3 + noise, 0, 10).astype("int") * 1000

# 売上のヒストグラムを描画
plt.hist(Y)
plt.xlabel("ダイエット商品の売上(単位:円)")
plt.show()

(出力結果)

このデータにおいて、たとえばDR-Learnerを実装してみます。

# 必要なライブラリをインポート
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LogisticRegression
from econml.dr import DRLearner

# モデルの構築
outcome_model = RandomForestRegressor(max_depth=3, random_state=0)
pseudo_treatment_model = RandomForestRegressor(max_depth=3, random_state=0)
propensity_model = LogisticRegression()
DR_learner = DRLearner(model_regression=outcome_model, 
                       model_propensity=propensity_model, 
                       model_final=pseudo_treatment_model)
DR_learner.fit(Y, T, X=x.reshape(-1, 1))

# 効果の推定
tau = DR_learner.effect(x.reshape(-1, 1))

# 推定された効果の可視化
plt.scatter(df[["x"]], tau, alpha=0.3)
plt.hlines(1000, -1, 0, linestyles='--', color="red")
plt.hlines(1000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(2000, 0, 0.5, linestyles='--', color="red")
plt.hlines(2000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.hlines(3000, 0.5, 1.0, linestyles='--', color="red")
plt.hlines(3000, -1, 1, linestyles='--', color="red", alpha=0.3)
plt.xlabel("ダイエットへの意識の高さ")
plt.ylabel("広告効果")
plt.ylim(0, 4000)
plt.show()

(出力結果)

このとき、ダイエットへの意識が低いグループの広告効果が過小評価されていることが分かります。

参考文献

おわりに

最後まで読んでいただきありがとうございました。他にも「Python×データ分析」をメインテーマに記事を執筆しているので、参考にしていただけたら幸いです。内容の誤り等がございましたら、コメントにてご指摘くださいませ。

他にも下記のような記事を書いています。ご一読いただけますと幸いです。

Discussion

こんにちは!
卒論で操作変数のデータ分析をやろうと思って、調べていたら記事を見つけました。
わかりやすくまとまっていて、非常に参考になります。

すごい細かいところなのですが、「おわりに」の見出しの下の「操作変数を用いた効果検証」と「回帰不連続デザインを用いた効果検証」のURLを押すと404に飛ばされます。
うとしんさんのプロフィールから「操作変数を用いた効果検証」と「回帰不連続デザインを用いた効果検証」の記事が閲覧できたので、多分URLがちょっと変わっちゃってるのではないかと思います。

もし問題ないようなら、無視しちゃって大丈夫です。

ご連絡いただきありがとうございます!
修正しておきます!

ログインするとコメントできます