はじめに
本稿の目的は、以下リンク(以下チュートリアル)先チュートリアルを学習し、Pythonのライブラリであるcausalmlを用いたアップリフトモデリングに関する理解を深めることです.
私個人としては、ゆくゆくは以下のリンクのコンテンツに挑戦したく、本稿の位置付けはそのための基礎トレーニングです.
環境など
- 私用Macパソコン(M2チップ)
- conda仮想環境
- VisualStudioCode
ライブラリのインストール
公式サイト(以下リンク)を参考にターミナル上で以下のコマンドを実行しました.
$ conda create --name envforcausalml python=3.8
$ conda install -c conda-forge causalml
参考
消去法的に上記手法になった、という感じです.
というのも、上記以外の方法も試みたのですが、エラーが解消できませんでした.
本編
関連ライブラリのインポート
本セクションでは関連ライブラリをインポートします.
以下はそのためのコードです.
import pandas as pd
import numpy as np
import matplotlib as plt
from sklearn.ensemble import \
RandomForestRegressor, \
GradientBoostingRegressor
from sklearn.tree import DecisionTreeRegressor
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from causalml.inference.meta import \
BaseSRegressor, \
BaseTRegressor, \
BaseXRegressor, \
BaseRRegressor
from causalml.inference.tree import \
UpliftTreeClassifier, \
UpliftRandomForestClassifier
from causalml.dataset.regression import \
synthetic_data
from sklearn.linear_model import LinearRegression
import shap
import matplotlib.pyplot as plt
import time
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
import os
import warnings
warnings.filterwarnings('ignore')
os.environ['KMP_DUPLICATE_LIB_OK'] = 'True' # for lightgbm to work
plt.style.use('fivethirtyeight')
%reload_ext autoreload
%autoreload 2
%matplotlib inline
上記のインポート文で最も気になったのはcausalmlからインポートしているいくつかのモデルです.
何故なら私自身一度も使用したことがないからです.
具体的には以下の4つです.
- BaseSRegressor
- BaseTRegressor
- BaseXRegressor
- BaseRRegressor
個人メモ
- causalmlのBaseSRegressor, BaseTRegressorはそれぞれSingle Treeアルゴリズム, Two Treeアルゴリズムに対応するものと予想※
- XRegressorとRRegressorについては憶測ですが、プロペンシティスコア(傾向スコア)を利用すものと想像
※: 以下リンクの投稿より
以上が関連ライブラリのインポートです.
次は元データの準備です.
元データの準備
本セクションでは元データを準備します。
以下はそのためのコードです.
n_features = 25
n_samples = 10000
y, X, w, tau, b, e = synthetic_data(
mode=1,
n=n_samples,
p=n_features,
sigma=0.5
)
w_multi = np.array(
['treatment_A' if x==1 else 'control' for x in w]
)
e_multi = {'treatment_A': e}
上記から、元データはcausalmlのsynthetic_data()により生成されていることが分かります.
以下は、causalml内のsynthetic_data()に関連するコードです.
def simulate_nuisance_and_easy_treatment(n=1000, p=5, sigma=1.0, adj=0.0):
"""Synthetic data with a difficult nuisance components and an easy treatment effect
From Setup A in Nie X. and Wager S. (2018) 'Quasi-Oracle Estimation of Heterogeneous Treatment Effects'
Args:
n (int, optional): number of observations
p (int optional): number of covariates (>=5)
sigma (float): standard deviation of the error term
adj (float): adjustment term for the distribution of propensity, e. Higher values shift the distribution to 0.
Returns:
(tuple): Synthetically generated samples with the following outputs:
- y ((n,)-array): outcome variable.
- X ((n,p)-ndarray): independent variables.
- w ((n,)-array): treatment flag with value 0 or 1.
- tau ((n,)-array): individual treatment effect.
- b ((n,)-array): expected outcome.
- e ((n,)-array): propensity of receiving treatment.
"""
X = np.random.uniform(size=n * p).reshape((n, -1))
b = (
np.sin(np.pi * X[:, 0] * X[:, 1])
+ 2 * (X[:, 2] - 0.5) ** 2
+ X[:, 3]
+ 0.5 * X[:, 4]
)
eta = 0.1
e = np.maximum(
np.repeat(eta, n),
np.minimum(np.sin(np.pi * X[:, 0] * X[:, 1]), np.repeat(1 - eta, n)),
)
e = expit(logit(e) - adj)
tau = (X[:, 0] + X[:, 1]) / 2
w = np.random.binomial(1, e, size=n)
y = b + (w - 0.5) * tau + sigma * np.random.normal(size=n)
return y, X, w, tau, b, e
上記から以下のことが分かります.
- 特徴量Xは一様分布から生成される
- プロペンシティスコアを表す変数eは定数値ではなく、特徴量Xの一部を用いた関数の出力として生成される
個人メモ
上記コード中のexpitとlogitについてBingに確認したところ以下のとおりでした
expitとlogitは、Scipyのspecialモジュールに含まれる関数です。これらの関数は、ロジスティック関数とその逆関数を表します。
expitは、ロジスティック関数を表します。ロジスティック関数は、次の式で定義されます。
$$
f(x) = \frac{1}{1 + e^{-x}}
$$
この関数は、$x$が大きな正の値になるにつれて、$f(x)$が1に近づきます。一方、$x$が大きな負の値になるにつれて、$f(x)$が0に近づきます。
logitは、ロジスティック関数の逆関数を表します。ロジスティック関数の逆関数は、次の式で定義されます。
$$
f^{-1}(y) = \ln\left(\frac{y}{1-y}\right)
$$
この関数は、$y$が0に近づくにつれて、$f^{-1}(y)$が負の無限大に近づきます。一方、$y$が1に近づくにつれて、$f^{-1}(y)$が正の無限大に近づきます。
これらの関数は、主にロジスティック回帰などの機械学習アルゴリズムで使用されます。
続いて、以下のコードが続きます.
feature_names = ['stars', 'tiger', 'merciful', 'quixotic', 'fireman', 'dependent',
'shelf', 'touch', 'barbarous', 'clammy', 'playground', 'rain', 'offer',
'cute', 'future', 'damp', 'nonchalant', 'change', 'rigid', 'sweltering',
'eight', 'wrap', 'lethal', 'adhesive', 'lip'] # specify feature names
model_tau = LGBMRegressor(importance_type='gain') # specify model for model_tau
上記より、model_tauとしてLightGBM(Gradient Boosting Machine)がインスタンス化されていることが分かります. このモデルを今後どのタイミングで使うことになるのか、気になりますが先に進みたいと思います.
個人メモ
上記コード中のLGBMRegressorの引数についてBingに確認したところ以下のとおりでした.
LGBMRegressorのドキュメントによると、importance_typeという引数があります。この引数は、特徴量の重要度を計算する方法を指定するために使用されます。デフォルト値は"split"ですが、他にも"gain"があります。"split"は、特徴量の分割回数に基づいて重要度を計算します。一方、"gain"は、特徴量の分割によって得られる平均的な情報利得に基づいて重要度を計算します。
importance_typeをgainに設定すると、特徴量の分割によって得られる平均的な情報利得に基づいて重要度を計算します。この方法は、特徴量の重要度をより正確に評価することができます。一方、importance_typeをsplitに設定すると、特徴量の分割回数に基づいて重要度を計算します。この方法は、特徴量の重要度を評価するための単純な方法ですが、情報利得を考慮しないため、より正確な評価ができない場合があります。
以上が元データの準備です.
次はモデルの学習と各種の特徴量重要度の可視化です.
モデル学習と特徴量重要度の可視化
本セクションではアップリフトを予測する複数のモデルを作成し、複数の特徴量重要度についてそれぞれ算出・可視化します.
SLearner
最初のモデルはチュートリアル内でSLearnerと呼ばれるものです.
以下はそのためのコードです.
base_algo = LGBMRegressor()
# base_algo = XGBRegressor()
# base_algo = RandomForestRegressor()
# base_algo = LinearRegression()
slearner = BaseSRegressor(
base_algo,
control_name='control'
)
slearner.estimate_ate(X, w_multi, y)
上記コードから、以下のことが分かります.
- base_algoなる変数にLightGBM(Light Gradient Boosting Machine)が採用されている
- 他のbase_algoはコメントアウトされている
- base_algoには、線形回帰に加えて決定木系のアルゴリズムがラインアップされている
- チュートリアルでSLearnerと謳われているのは、causalmlのBaseSRegressorである
感想めいた話ですが、causalmlでモデルをインスタンス化する際、介入有無は'control'のように文字列で指定させられるのが特徴的だと感じました.
さて、最終行のコードの返り値は以下のとおりでした.
array([0.57322393])
メソッド名estimate_ate()からして、ATE(Average Treatment Effect)の推定値でしょうか. Bingに確認します.
個人メモ
- causalmlのBaseRRegressor().estimate_ate()について教えてください
BaseSRegressor().estimate_ate()は、causalmlライブラリのBaseSRegressorクラスのメソッドであり、因果推論における平均処置効果（ATE）を推定するために使用されます。このメソッドは、2つの回帰モデルを使用してATEを推定します。このメソッドは、以下のように呼び出すことができます。
from causalml.inference.tree import BaseSRegressor # Instantiate the model model = BaseSRegressor() # Fit the model model.fit(Y, T, X) # Estimate the ATE ate = model.estimate_ate(X, p)ここで、Yはアウトカム変数、Tは介入変数、Xは共変量です。このメソッドは、介入変数が二値である場合にのみ使用することができます。また、このメソッドは、2つの回帰モデルを使用してATEを推定するため、介入変数とアウトカム変数の関係が線形である場合に最も適しています。
ここで元データをピボット集計します.
以下はピボットテーブルのコードで、チュートリアル内にはこのような手順はありません.
tmp_df = pd.DataFrame()
tmp_df[feature_names] = X
tmp_df['outcome'] = y
tmp_df['treatment'] = w_multi
tmp_df['tau'] = tau
pivs = []
for target in ['outcome', 'tau']:
piv = tmp_df.pivot_table(
values = target,
index = 'treatment',
aggfunc = [len, np.sum, np.average],
margins=True
)
pivs.append(piv)
下表は上記コードで作成したピボット集計結果の1つです.
| treatment | ('len', 'outcome') | ('sum', 'outcome') | ('average', 'outcome') |
|---|---|---|---|
| control | 4796 | 4908.67 | 1.02349 |
| treatment_A | 5204 | 10126.5 | 1.9459 |
| All | 10000 | 15035.1 | 1.50351 |
上表から次のようなことが分かります.
- 介入群とコントロール群のoutcomeの平均の差が約0.9
- 上記は先ほどの推定値0.57と大きく異なる
今回はRCTではないので、上表のような2群で、outcomeの平均の差をとり、それを母集団における平均的なアップリフトだとするのは誤った考え方です.
そこで、真のアップリフト(synthetic_data上のtau)について、先ほどと同様の体裁で集計したものが下表となります.
| treatment | ('len', 'tau') | ('sum', 'tau') | ('average', 'tau') |
|---|---|---|---|
| control | 4796 | 1872.78 | 0.390489 |
| treatment_A | 5204 | 3114 | 0.598387 |
| All | 10000 | 4986.79 | 0.498679 |
上表から以下のことが分かります.
- 今回のデータにおけるATE(Average Treatment Effect)約0.49
- 上記は先ほどの推定値0.57よりも小さい
したがって、SLearnerによるATEの推定を試みると、単に2群のoutcomeの平均の差をとった時よりも、バイアスが緩和された推定が出来ていることが分かりました.
ここでチュートリアルの手順に戻り、fit_predict()なるメソッドを行います.
以下はそのためのコードです.
slearner_tau = slearner.fit_predict(
X=X,
treatment=w_multi,
y=y
)
以上がSLearnerのモデル学習でした.
次のセクションでは、学習済のSLearnerを使用して、各種の特徴量重要度を見ていきます.
Feature Importance(SLearner)
method = 'auto'
BaseSRegressorのget_importanceメソッドで、特徴量重要度を取得できます. 標題の'auto'とは、このメソッドの引数に与える文字列です. 以下はそのためのコードです.
slearner.get_importance(
X=X,
tau=slearner_tau,
normalize=True,
method='auto',
features=feature_names
)
"""
{'treatment_A': tiger 0.471209
stars 0.354741
quixotic 0.067203
fireman 0.049105
merciful 0.048931
lethal 0.001361
dependent 0.000902
cute 0.000886
offer 0.000623
playground 0.000550
lip 0.000487
barbarous 0.000457
eight 0.000395
change 0.000384
touch 0.000381
future 0.000310
adhesive 0.000308
rigid 0.000301
damp 0.000269
sweltering 0.000251
wrap 0.000227
nonchalant 0.000226
shelf 0.000223
clammy 0.000157
rain 0.000113
dtype: float64}
"""
上記から以下のことが分かります.
- 引数のslearner_tauは学習済モデルに学習データを与えた結果としてのCATE(Conditional Average Treatment Effect)の推定値となっており、正解データではない
個人メモ
一般的な特徴量重要度が、モデルに予測させた推定値をベース算出するものなのかよく分かりませんでした. ただしBing氏によると正解データを用いて算出するようです.
get_importance()メソッドの他に、plot_importance()メソッドがありこちらは描画ができます.
以下はそのためのコードです.
slearner.plot_importance(
X=X,
tau=slearner_tau,
normalize=True,
method='auto',
features=feature_names
)
上記コードの実行結果は下図のとおりです.
method = 'permutation'
続いて、引数をpermutationとする場合です.
permutation importanceとは、任意の特徴量をレコード間で無作為にシャッフルした場合の精度悪化の度合いをもとに算出された特徴量重要度を指すようです. 以下の投稿では、sklearnを利用したpermutation importanceの算出を実施しましたが、今回のBaseSRegressorでは学習器にメソッドが付属しているようです.
さて、以下はそのためのコードです.
slearner.get_importance(
X=X,
tau=slearner_tau,
method='permutation',
features=feature_names,
random_state=42
)
ランダムシードを明示しているのは前述の"シャッフル"が起こるからだと理解しました. 返り値は割愛します. 先ほどと同様に、plot_importance()を実行します. 以下はそのためのコードと描画です.
slearner.plot_importance(
X=X,
tau=slearner_tau,
method='permutation',
features=feature_names,
random_state=42
)
sklaern.inspection.permutation_importance
続いてsklearnを用いてpermutation importanceを計算します. 以下はそのためのコードです.
start_time = time.time()
X_train, X_test, y_train, y_test = train_test_split(
X,
slearner_tau,
test_size=0.3,
random_state=42
)
model_tau_fit = model_tau.fit(X_train, y_train)
perm_imp_test = permutation_importance(
estimator=model_tau_fit,
X=X_test,
y=y_test,
random_state=42
).importances_mean
pd.Series(perm_imp_test, feature_names).sort_values(ascending=False)
print("Elapsed time: %s seconds" % (time.time() - start_time))
# Elapsed time: 0.8685302734375 seconds
上記のコードを見ると、突如時間を計測していることが分かります. 遅いぞ、ということが言いたいのでしょうか.
また、以下のことが分かります.
- 所謂SLearnerではなく、データ準備のセクションでインスタンス化したLGBMRegressorが使用されている
- トレーニングデータとテストデータに分割している
- 学習済のSLearnerに学習データを与えて予測させたCATEの推定値を教師データとしてモデルを学習させています.
何も説明がないものの、このようにして計算したpermutation importanceはBaseSRegressorに付属しているメソッドを使用して計算したpermutation importanceと一致するようです.
以下は、描画のためのコードと、実行結果の図です.
pd.Series(perm_imp_test, feature_names).sort_values(). \
plot(kind='barh', figsize=(12, 8))
トレーニングデータを使用したpermutation importanceについても同様に計算でき、チュートリアルにはコードがありましたが、ここでは割愛します.
Shapley Values
このセクションではshap値を計算します. shap値の計算はBaseSRegressorのメソッドから可能です.
以下の投稿では、shapを利用したshap値の算出を実施しましたが、今回のBaseSLearnerでは学習器にメソッドが付属しているようです.
以下はそのためのコードです.
shap_slearner = slearner.get_shap_values(
X=X,
tau=slearner_tau
)
shap_slearner
以下は上記のコードの実行結果です.
{'treatment_A': array([[ 1.91859822e-01, 2.04385529e-01, -2.12118899e-02, ...,
1.24018015e-04, 1.34670086e-03, 3.12601522e-04],
[-3.66003741e-02, 1.32411912e-01, -5.08025213e-02, ...,
2.83559627e-03, -6.34298699e-04, 1.59553163e-04],
[-2.09073781e-02, -1.33604863e-01, 5.39808331e-02, ...,
-4.66034901e-04, 1.03995845e-03, 2.08040828e-03],
...,
[-8.74254782e-03, 1.95864633e-01, -1.91824122e-02, ...,
8.63983668e-04, -2.68699273e-04, -1.12722375e-03],
[ 8.63892313e-02, 1.32205971e-01, 4.70598575e-02, ...,
1.92042942e-04, -4.32239596e-04, -2.96629240e-04],
[-7.71091385e-02, -2.00123296e-01, 2.23572794e-02, ...,
1.83305175e-04, -6.13922523e-03, 5.71021128e-05]])}
こちらも学習データすべてを使用しているように見えますね.
ndarrayの形状を調べてみます.
以下はそのためのコードで、チュートリアルにこのような手順はありません.
shap_slearner['treatment_A'].shape
# (10000, 25)
チュートリアルに戻ります.
後続で利用しないものの、shap値の平均は次のように計算されます.
np.mean(np.abs(shap_slearner['treatment_A']), axis=0)
"""
array([0.11526197, 0.14080281, 0.02436171, 0.03437657, 0.02983041,
0.00167865, 0.00018388, 0.00054616, 0.00089338, 0.0001909 ,
0.0018401 , 0.00032663, 0.00176994, 0.00369364, 0.00030464,
0.00055004, 0.00036366, 0.00120931, 0.00032772, 0.00066475,
0.00041256, 0.00037257, 0.00214374, 0.00081676, 0.0011478 ])
"""
更に、BaseSRegressorのメソッドplot_shap_values()を使用して、shap値の描画が可能です.
以下はそのためのコードです.
# Plot shap values without specifying shap_dict
slearner.plot_shap_values(
X=X,
tau=slearner_tau,
features=feature_names
)
こちらは先ほどのfeature importanceのようにget系のメソッドとplot系のメソッドで引数の与え方が類似していることが分かります.
下図は上記コードの実行結果です.
以下のように、先ほどget_shap_values()の返り値である、shap_slearnerを引数に与えることもできるようです.
以下はそのためのコードです.
# Plot shap values WITH specifying shap_dict
slearner.plot_shap_values(
X=X,
shap_dict=shap_slearner,
features=feature_names
)
以下は結果です.
介入の種類が複数あるときは複数のグラフが出力されるのかもしれません.
続いて、shapにもあったdependenceプロットです.
以下はそのためのコードです.
# interaction_idx set to None (no color coding for interaction effects)
slearner.plot_shap_dependence(
treatment_group='treatment_A',
feature_idx=1,
X=X,
tau=slearner_tau,
interaction_idx=None,
shap_dict=shap_slearner,
features=feature_names
)
下図は上記コードの実行結果です.
今回使用したsynthetic_dataでは特徴量が一様分布から生成されていたのと平仄があっています.
interaction_idxは'auto'にも設定できるようで、以下はそのためのコードです.
# interaction_idx set to 'auto'
# (searches for feature with greatest approximate interaction)
# specify feature names
slearner.plot_shap_dependence(
treatment_group='treatment_A',
feature_idx='tiger',
X=X,
tau=slearner_tau,
interaction_idx='auto',
shap_dict=shap_slearner,
features=feature_names
)
下図は上記コードの実行結果です.
上図から、tigerが0.4から0.6の間では、starsが大きいとtigerのshap値が高くなり、starsが小さいとtigerのshap値が小さくなっていることが分かります. 真のモデルではX[:, 0]とX[:, 1]の掛け算が含まれていたのでそれが原因でしょうか.
なお、interaction_idxは明示することも可能で、以下はそのためのコードです.
# interaction_idx set to specific index
slearner.plot_shap_dependence(
treatment_group='treatment_A',
feature_idx=1,
X=X,
tau=slearner_tau,
interaction_idx=10,
shap_dict=shap_slearner,
features=feature_names
)
下図は上記コードの実行結果です.
真のモデルに上図の2変数が相乗効果を生む項はないので、このように均一な色付けになるのだと理解しました.
以上が、学習済のSLearnerを用いた各種特徴量重要度の計算でした.
余談ですが、以下の投稿ではshapを用いて特徴量重要度を計算するチュートリアルを学びました.
個人的には不審な点が多く使えないなと判断していたのですが、今回のチュートリアルではうまく計算できているようでスッキリしました.
TLearner
このセクションではTLearnerを実施します.
以下はそのためのコードです.
tlearner = BaseTRegressor(
LGBMRegressor(), control_name='control'
)
tlearner.estimate_ate(X, w_multi, y)
# (array([0.56621952]), array([0.55155255]), array([0.58088649]))
上記から、以下のことが分かります.
- estimate_ate()の返り値が3つのndarray
この3つの数値が何なのかを明確にするため、causalmlのestimate_ate()のコードの説明を確認しました.
以下はその抜粋です.
"""
Args:
X (np.matrix or np.array or pd.Dataframe): a feature matrix
treatment (np.array or pd.Series): a treatment vector
y (np.array or pd.Series): an outcome vector
bootstrap_ci (bool): whether to return confidence intervals
n_bootstraps (int): number of bootstrap iterations
bootstrap_size (int): number of samples per bootstrap
Returns:
The mean and confidence interval (LB, UB) of the ATE estimate.
pretrain (bool): whether a model has been fit, default False.
"""
ご覧のとおり、ATE推定値の平均、信頼区間のLB, UBがそれぞれ対応するようです.
ATEの推定値の数値自体はSLearnerとあまり変わり映えしませんね.
個人的にはATEの信頼区間とはどのように計算するのか気になります(ブートストラップ法などを使うのでしょうか？)
ですが、当初の疑問が解消されたので、先に進みます.
先ほどと同様にfit_predict()なるメソッドを行います.
以下はそのためのコードです.
tlearner_tau = tlearner.fit_predict(X, w_multi, y)
以上がTLearnerのモデル学習でした.
次のセクションでは、学習済のTLearnerを使用して、各種の特徴量重要度を見ていきます.
Feature Importance(TLearner)
method = 'auto'
BaseTRegressorのget_importanceメソッドで、特徴量重要度を取得できます. 標題の'auto'とは、このメソッドの引数に与える文字列です. 以下はそのためのコードです.
tlearner.get_importance(
X=X,
tau=tlearner_tau,
normalize=True,
method='auto',
features=feature_names
)
最初に触れたBaseSRegressorの場合と同様なので、出力は割愛します.
get_importance()メソッドの他に、plot_importance()メソッドがありこちらは描画ができます.
以下はそのためのコードと実行結果です.
tlearner.plot_importance(
X=X,
tau=tlearner_tau,
normalize=True,
method='auto',
features=feature_names
)
method = 'permutation'
続いて、引数をpermutationとする場合です.
以下はそのためのコードです.
tlearner.get_importance(
X=X,
tau=tlearner_tau,
normalize=True,
method='permutation',
features=feature_names
)
最初に触れたBaseSRegressorの場合と同様なので、出力は割愛します.
先ほどと同様に、plot_importance()を実行します. 以下はそのためのコードと描画です.
tlearner.plot_importance(
X=X,
tau=tlearner_tau,
normalize=True,
method='permutation',
features=feature_names
)
sklaern.inspection.permutation_importance
続いてsklearnを用いてpermutation importanceを計算します.
以下はそのためのコードです.
start_time = time.time()
X_train, X_test, y_train, y_test = train_test_split(
X,
tlearner_tau,
test_size=0.3,
random_state=42
)
model_tau_fit = model_tau.fit(X_train, y_train)
perm_imp_test = permutation_importance(
estimator=model_tau_fit,
X=X_test,
y=y_test,
random_state=42
).importances_mean
pd.Series(perm_imp_test, feature_names).sort_values(ascending=False)
print("Elapsed time: %s seconds" % (time.time() - start_time))
先ほどのSLearnerの場合と異なり、sklearnのpermutation importanceを上記の手順で実施した結果と、BaseTRegressorに付属しているメソッドで算出したpermutation importanceは異なる結果になるようです.
以下は、描画のためのコードと、実行結果の図です.
pd.Series(perm_imp_test, feature_names). \
sort_values().plot(kind='barh', figsize=(12, 8))
plt.title("Test Set Permutation Importances")
Shapley Values
このセクションではshap値を計算します. shap値の計算はBaseTRegressorのメソッドから可能です.
以下はそのためのコードです.
shap_tlearner = tlearner.get_shap_values(X=X, tau=tlearner_tau)
shap_tlearner
"""
{'treatment_A': array([[ 0.26873947, 0.27288131, -0.03400783, ..., -0.01144844,
-0.00097293, 0.00075564],
[-0.04374247, 0.12915531, -0.0570392 , ..., 0.03030283,
-0.00238386, -0.00163477],
[-0.04404788, -0.11898599, 0.04792128, ..., -0.01021204,
0.00036785, -0.00298231],
...,
[-0.00631418, 0.20273592, -0.02579956, ..., 0.02101015,
-0.00223981, -0.00325014],
[ 0.08798084, 0.15789791, 0.08645094, ..., -0.00206632,
0.00259849, 0.00188204],
[-0.06466238, -0.22534402, 0.00804794, ..., -0.00801698,
0.0106461 , -0.00172914]])}
"""
更に、BaseTRegressorのメソッドplot_shap_values()を使用して、shap値の描画が可能です.
以下はそのためのコードとその描画です.
# Plot shap values without specifying shap_dict
tlearner.plot_shap_values(
X=X,
tau=tlearner_tau,
features=feature_names
)
チュートリアルではplot_shap_values()への引数の与え方が複数紹介されているものの、先ほどのSLearnerと同様なので割愛します.
また、dependenceプロットについてはそもそもチュートリアル中で触れられていないのでこれもまた割愛します.
以上が、学習済のTLearnerを用いた各種特徴量重要度の計算でした.
XLearner
このセクションではXLearnerを実施します.
以下はそのためのコードです.
xlearner = BaseXRegressor(
learner=LGBMRegressor(),
control_name='control'
)
xlearner.estimate_ate(
X=X,
treatment=w_multi,
y=y,
p=e_multi
)
"""
(array([0.51458419]), array([0.5007696]), array([0.52839878]))
"""
上記から以下のことが分かります.
- 引数pがあり、プロペンシティスコア(傾向スコア)を利用するモデルである
- 今回は引数pに真のプロペンシティスコアが与えられている
- ATEの推定値は約0.5となっており、SLearner, TLearnerによる推定値と比較して、真のATE(約0.47)に肉薄している
ユースケースとしては以下が思い浮かびました.
- 特徴量から介入確率(ベルヌーイ試行の母数)を機械的に計算し、介入割当を確率的に実施するようなケース
なお、個人的にはプロペンシティスコアが未知の場合の予測方法が気になっています. 残念ながら今回のチュートリアルでプロペンシティスコアの推定は取り扱われませんでした.
ここで、XLearnerとは何なのかは以下リンクに記載があります.
XLearnerでは介入群、コントロール群からそれぞれモデルを作成し、傾向スコアで重みづけしているようです. モデルを作成する際の正解データの作り方がXの形に見えなくもないことからXLearnerと呼ぶのだと理解しました.
以上がXLearnerのモデル学習でした.
Feature Importances(XLearner)
XLearnerについても各種の特徴量重要度が計算できるものの、SLearnerと同様の手順で、かつ結果に特筆すべきことがないので本稿では割愛します.
RLearner
このセクションではRLearnerを実施します.
以下はそのためのコードです.
rlearner = BaseRRegressor(
learner=LGBMRegressor(),
control_name='control'
)
rlearner_tau = rlearner.fit_predict(
X=X,
treatment=w_multi,
y=y,
p=e_multi
)
上記から、XLearnerと同様、プロペンシティスコアを利用していることが分かります. チュートリアルでは触れられていないものの、estimate_ate()メソッドも用意されているようです.
以下はそのためのコードです.
rlearner.estimate_ate(
X=X,
treatment=w_multi,
y=y,
p=e_multi
)
# (array([0.5249357]), array([0.52471149]), array([0.5251599]))
こちらも、XLearnerと同様比較的高精度であるようです.
RLearnerについては以下の文献で説明されていることまでは分かりましたが私の理解が及ばず、本稿では参考文献の掲示に留めます.
以上がRLearnerのモデル学習でした.
Feature Importances(RLearner)
RLearnerについても各種の特徴量重要度が計算できるものの、XLearnerの際と同様の理由で割愛します.
まとめ
- causalmlのBaseSRegressor, BaseTRegressor, BaseXRegressor, BaseRRegressorにベースアルゴリズムとしてLightGBMを与えることでアップリフトの推定及び特徴量重要度の算出、図示が簡便に実施できる
- BaseXRegressor, BaseRRegressorについてはいずれもプロペンシティスコア(傾向スコア)を利用するものの、本稿のチュートリアルでは、合成データに付随する真のプロペンシティスコアを用いて学習を実施した
- これらのアルゴリズムはCausalRandomForestやCausalUpliftRandomForestなどと異なり、ベースアルゴリズムに柔軟性がある