はじめに

本稿の目的は、以下リンク(以下チュートリアル)先チュートリアルを学習し、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 multiprocessing as mp
np.random.seed(42)
from sklearn.model_selection import train_test_split
from sklearn.inspection import permutation_importance
import shap
import causalml 
from causalml.metrics import plot_gain, plot_qini, qini_score 
from causalml.dataset import synthetic_data 
from causalml.inference.tree import plot_dist_tree_leaves_values, get_tree_leaves_mask 
from causalml.inference.tree import CausalRandomForestRegressor, CausalTreeRegressor 
from causalml.inference.tree.utils import timeit 
import matplotlib.pyplot as plt 
import seaborn as sns 
%config InlineBackend.figure_format = 'retina'

今回使用しないライブラリがインポートされていないとも限らないので、深入りせず、先へ進みます.

合成データの生成

以下のコードで合成データを作成します. 変数は全てndarrayで、Xは2次元配列, それ以外は1次元配列です.

# Simulate randomized trial: mode=2 
y, X, w, tau, b, e = synthetic_data(
    mode=2, 
    n=2000, 
    p=10, 
    sigma=3.0
)

説明書きに記載のとおり、今回もRCT(Randomized Controlled Trial)による介入割当を前提としたデータであるようです.

参考: Bingとのチャット

causalmlのsynthetic_data()の引数について教えてください.

synthetic_data() という関数は、人工的なデータを生成するためのものです。この関数には以下の引数があります。

mode: 人工データの生成方法を指定するパラメータです。1から4までの整数値を指定できます。デフォルトは1です。

n: 生成するデータのサンプル数を指定するパラメータです。デフォルトは10000です。

p: 生成するデータの特徴量の数を指定するパラメータです。デフォルトは8です。

sigma: 生成するデータのノイズの標準偏差を指定するパラメータです。デフォルトは1.0です。
random_state: 生成するデータの乱数シードを指定するパラメータです。デフォルトはNoneです。

ノイズは正規分布に従うのでしょうか？

はい、ノイズは正規分布に従います。synthetic_data() 関数のソースコードを見ると、ノイズは np.random.normal(0, sigma, size=n) というように生成されています。ここで、sigma は引数で指定した標準偏差です。正規分布の平均は0です。したがって、ノイズは平均0、標準偏差 sigma の正規分布に従うと言えます。

mode=2の場合どのようなデータが生成されますか？

mode=2 の場合、以下のようなデータが生成されます。
特徴量 X は、標準正規分布から独立にサンプリングされます。
処置変数 treatment は、X[:, 0] + X[:, 3] + np.random.normal(0, 0.01, size=n) に基づいて、ロジスティック回帰で確率を計算し、ベルヌーイ分布からサンプリングされます。
アウトカム変数 y は、treatment * (X[:, 0] + X[:, 9]) + X[:, 1] + np.random.normal(0, sigma, size=n) として計算されます。
因果効果変数 tau は、X[:, 0] + X[:, 9] として計算されます。
つまり、mode=2 の場合、特徴量の一部が処置割り当てやアウトカムに影響を与えることになります。また、因果効果は特徴量に依存することになります。

続いて、これらをpandasデータフレームへ格納します.

df = pd.DataFrame(X) 
feature_names = [f'feature_{i}' for i in range(X.shape[1])]
df.columns = feature_names 
df['outcome'] = y 
df['treatment'] = w 
df['treatment_effect'] = tau

続いて、データを分割します.

# Split data to training and testing samples for model validation
df_train, df_test = train_test_split(df, test_size=0.2, random_state=111) 
n_train, n_test = df_train.shape[0], df_test.shape[0] 

X_train, y_train = df_train[feature_names], df_train['outcome'].values 
X_test, y_test = df_test[feature_names], df_test['outcome'].values 
treatment_train, treatment_test = df_train['treatment'].values, df_test['treatment'].values 
effect_test = df_test['treatment_effect'].values

effect_testという変数は何に使うのか気になりますね.
sklearnのtrain_test_split()について、よく分からない点があったのでBingに確認してみました.

参考: Bingとのチャット

sklearn.model_selection.train_test_split()の返り値のデータ型はなんですか？

返り値のデータ型は、入力配列と同じです。入力がスパースの場合は、scipy.sparse.csr_matrixになります。それ以外の場合は、numpyの配列やpandasのdataframeなど、入力と同じタイプになります。

ここで、データの中身を確認しておきたいと思います.
まずはdf_trainです.

	1696
feature_0	-1.59201
feature_1	0.452005
feature_2	0
feature_3	-0.466855
feature_4	-0.363785
feature_5	-0.770119
feature_6	-0.348852
feature_7	1.03256
feature_8	-0.713953
feature_9	0.545141
outcome	0.831146
treatment	1
treatment_effect	-0.647534

続いて、dfをoutcomeでピボット集計したものです.

treatment	('size', 'outcome')	('sum', 'outcome')	('mean', 'outcome')
0	1009	830.695	0.823286
1	991	1360.44	1.3728
All	2000	2191.14	1.09557

上表のとおり、outcomeの平均はおよそ0.5異なっているようです. 確認のために、dfをtreatment_effectでピボット集計します.

treatment	('size', 'treatment_effect')	('sum', 'treatment_effect')	('mean', 'treatment_effect')
0	1009	895.044	0.88706
1	991	781.618	0.788716
All	2000	1676.66	0.838331

CausalTreeRegressorの実行

ここでは以下のとおり、causalml.inference.tree.CausalTreeRegressorで学習を行います.

ctree = CausalTreeRegressor()
ctree.fit(X=X_train.values, y=y_train, treatment=treatment_train)

公式サイトのコードを見ると、CausalTreeRegressor()の引数control_nameはデフォルトが0になっているようでした. 今回はデフォルト値で問題なさそうです.

CausalRandomForestClassifierの実行

続いてCausal Random Forestを実施します.

crforest = CausalRandomForestRegressor(
    criterion="causal_mse", 
    min_samples_leaf=200, 
    control_name=0, 
    n_estimators=50, 
    n_jobs=mp.cpu_count() - 1
)
crforest.fit(X=X_train, y=y_train, treatment= treatment_train)

不純度ベースの特徴量重要度の計算

以下のように特徴量の重要度をツリーモデル、フォレストモデルそれぞれについて対応する列に格納し、pandasデータフレームを作成します.

df_importances = pd.DataFrame(
    {'tree': ctree.feature_importances_, 
     'forest': crforest.feature_importances_, 
     'feature': feature_names
     })

ツリーモデルやランダムフォレストでは、上記コードのように特徴量の重要度が.feature_importances_に格納されているようです.

続いてランダムフォレストについて標準偏差を計算します.

forest_std = np.std(
    [tree.feature_importances_ for tree in crforest.estimators_], 
    axis=0
)

なお、crforest.estimators_には50個の要素が格納されており、ランダムフォレスト分類器をインスタンス化する際の引数n_estimatorsに紐づいているようです.

以下は参考です.

print(type(crforest.estimators_))
print(len(crforest.estimators_))
print(type(crforest.estimators_[0])) 

"""
<class 'list'>
50
<class 'causalml.inference.tree.causal.causaltree.CausalTreeRegressor'>
"""

続いて、Causal Treeの特徴量重要度を棒グラフで表示します. 以下はそのコードです.

fig, ax = plt.subplots() 
df_importances['tree'].plot.bar(ax=ax) 
ax.set_title("Causal Tree feature importances") 
ax.set_ylabel("Mean decrease in impurity") 
ax.set_xticklabels(feature_names, rotation=45) 
plt.show()

上記のコードの実行結果は下図のとおりです.

同様に、Causal Random Forest についても描画します. 以下はコードです.

fig, ax = plt.subplots() 
df_importances['forest'].plot.bar(yerr=forest_std, ax=ax)
ax.set_title("Causal Forest feature importances") 
ax.set_ylabel("Mean decrease in impurity") 
ax.set_xticklabels(feature_names, rotation=45) 
plt.show()

下図は上記コードの実行結果です.

Causal Random Forestの場合棒グラフ上にエラーバーが表示されていることが分かります. これは標準誤差でしょうか. 不純度ベースの特徴量重要度とは一体何を指すのかイマイチ判然としませんが、全ての特徴量について、エラーバーの範囲が0.0を含んでいます. つまり今回のケースですと、サンプル数不足で不純度ベースの特徴量重要度については何も言えないということでしょうか.

Permutation importance　

今度はpermutation importanceと呼ばれる数値を計算し、不純度ベースの特徴量重要度と同様に、棒グラフへ描画します.

permutation importanceとは、任意の特徴量をレコード間で無作為にシャッフルした場合の精度悪化の度合いをもとに算出された特徴量重要度のことみたいです.

以下はその計算及び描画のためのコードです.

for name, model in zip(('Causal Tree', 'Causal Forest'), (ctree, crforest)): 
    imp = permutation_importance(model, X_test, y_test, 
                                 n_repeats=50, 
                                 random_state=0)
    fig, ax = plt.subplots() 
    ax.set_title(f'{name} feature importances') 
    ax.set_ylabel('Mean decrease in impurity') 
    plt.bar(feature_names, imp['importances_mean'], yerr=imp['importances_std']) 
    ax.set_xticklabels(feature_names, rotation=45) 
    plt.show()

permutation_importance()は冒頭にインポートしたsklearnの関数です. 返り値には標準誤差(標準偏差?)のフィールドも入っているようですね.
下図は上記コードの実行結果です.

今度は標準誤差(偏差?)が0.0から離陸することができたようです.
先ほどBingチャットに教えてもらった真のモデルの生成ロジックと照らし合わせると違和感があるのでもとのコードを確認したところ以下のようになっていました.

def simulate_randomized_trial(n=1000, p=5, sigma=1.0, adj=0.0):
    """Synthetic data of a randomized trial
        From Setup B 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): no effect. added for consistency
    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.normal(size=n * p).reshape((n, -1))
    b = np.maximum.reduce([np.repeat(0.0, n), X[:, 0] + X[:, 1], X[:, 2]]) + np.maximum(
        np.repeat(0.0, n), X[:, 3] + X[:, 4]
    )
    e = np.repeat(0.5, n)
    tau = X[:, 0] + np.log1p(np.exp(X[:, 1]))

    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

上記だと、outcomeにはインデックス0〜4番目の特徴量のみが使用されているようです. またプロペンシティスコアは定数0.5です. どうやら生成系AIに出まかせを吹き込まれていたようです.

ともあれ、こちらの棒グラフは、真のモデルと平仄がとれていそうな佇まいであることが分かりました.

SHAP値

続いて、SHAP値です.

以下のコードでは、レコードごとのSHAP値を特徴量ごとに打点します.

# The error "Model type not yet supported by TreeExplainer" will occur 
# if causalml support PR is not closed yet. 
# git clone https://github.com/alexander-pv/shap.git && cd shap
# pip install .
import shap
tree_explainer = shap.TreeExplainer(ctree) 
shap.initjs() 
shap.summary_plot(tree_explainer.shap_values(X_test)[0], X_test, plot_type="dot")

shap_values()[0]のように最初のインデックスを指定する必要があるようで、チュートリアルどおりのコードだとうまくいかなかったので修正しました. shap_values(X_test)は2要素あり、実は2要素目が何なのか、イマイチ理解が出来ていないのですが先に進みます. 下図は上記コードの実行結果です.

これは一見すると何が言いたいのか分かりませんね.
続いてCausal Random Forestについても同様の手順を実施します. 以下はそのためのコードです.

cforest_explainer = shap.TreeExplainer(crforest)
shap.initjs()
shap.summary_plot(cforest_explainer.shap_values(X_test)[0], X_test)

以下は実行結果です.

これはよく分かりませんね. 例えば真のモデルではfeature_0が使われていますが、このグラフを見ると説明力があるのかないのか判然としません. また、そもそも論ですが、チュートリアルの描画と全く異なります. そこで、shap_values(X_test)[1]で再描画してみました. 以下はそのコードです.

shap.summary_plot(cforest_explainer.shap_values(X_test)[1], X_test)

下図は上記コードの結果です.

ここで、shap.TreeExplainer()の挙動を知るために公式サイト（以下リンク）を確認しました.

shap.TreeExplainer().shap_values()の返り値はlistかarrayであるようです. アウトプットがベクトルであるようなモデルの場合返り値はリストになり、リストの各要素はそれぞれのアウトプットに対応するshap値の行列になるそうです. 今回は返り値が要素数2のリストになっているので、モデルのアウトプットがベクトルということになりますね.

以下はpredict()のコードで,返り値は1個, 又は3個になるように見えます.

    def predict(self, X: np.ndarray, with_outcomes: bool = False) -> np.ndarray:
        """Predict individual treatment effects

        Args:
            X (np.matrix): a feature matrix
            with_outcomes (bool), default=False,
                                  include outcomes Y_hat(X|T=0), Y_hat(X|T=1) along with individual treatment effect
        Returns:
           (np.matrix): individual treatment effect (ITE), dim=nx1
                        or ITE with outcomes [Y_hat(X|T=0), Y_hat(X|T=1), ITE], dim=nx3
        """
        if with_outcomes:
            self.n_outputs_ = self.max_outputs_
            for estimator in self.estimators_:
                estimator._with_outcomes = True
        else:
            self.n_outputs_ = 1
        y_pred = super().predict(X)
        return y_pred

残念ながらこれ以上のことは分かりませんでしたが、先に進みたいと思います.
以下では、shap.dependence_plotで、各特徴量のshap値について詳しく確認することができます. 以下はそのためのコードとその実行結果です.

shap.dependence_plot("feature_0", 
                     cforest_explainer.shap_values(X_test)[1], 
                     X_test, 
                     interaction_index="feature_2")

feature_0とそのSHAP valueには正の相関があるといってよさそうです. またfeature_0とfeature_2に相関がありそうなことも見て取れます.ちなみにshap_values(X_test)[0]だとこうなります.

shap.dependence_plot("feature_0", 
                     cforest_explainer.shap_values(X_test)[0], 
                     X_test, 
                     interaction_index="feature_2")

上図を見ると、feature_0は数値が負のときにshap値と負の相関を持っているように見えますね. 数値が小さいほど、予測値の押し上げ要因になるということでしょうか. 真のモデルではそのようになっていなかったのでイマイチ釈然としません.

チュートリアルは以上です.
最後にアップリフトカーブを描画して終了したいと思います. 以下はそのコードと実行結果です.

usecols = [
    'treatment', 
    'outcome'
]
auuc_metrics = df_test[usecols]
pred_uplift = crforest.predict(X=X_test, with_outcomes=True)
auuc_metrics['pred_uplift'] = pred_uplift[:, 0] 

from causalml.metrics import plot_lift, plot_gain 

plot_lift(auuc_metrics, 
          outcome_col='outcome', 
          treatment_col='treatment'
          )

uplift降順に並び替えて最初の50人あたりまでが施策効果が平均以上に効いているようです.

まとめ

causalmlのCausalRandomForestClassifierを用いて学習、予測タスクを行うことができた
causalmlのCausalRandomForestClassifierのfeature_importancesやsklearnのpermutation importance、shapのshap_valuesを用いて特徴量の重要度を可視化することができた
SHAP値の可視化については、その数値の生成過程に不可解な点があり、未解決
UpliftRandomForestClassifierとCausalRandomForestClassifierの棲み分けは不明