はじめに
本稿の目的は、以下リンク(以下チュートリアル)先チュートリアルを学習し、Pythonのライブラリであるcausalmlを用いた所謂ツリーベースのアップリフトモデリングに関する理解を深めることです.
私個人としては、ゆくゆくは以下のリンクのコンテンツに挑戦したく、本稿の位置付けはそのための基礎トレーニングです.
環境など
- 私用Macパソコン(M2チップ)
- conda仮想環境
- VisualStudioCode
ライブラリのインストール
公式サイト(以下リンク)を参考にターミナル上で以下のコマンドを実行しました.
$ conda create --name envforcausalml python=3.8
$ conda install -c conda-forge causalml
参考
消去法的に上記手法になった、という感じです.
というのも、上記以外の方法も試みたのですが、エラーが解消できませんでした.
本編
関連ライブラリのインポート
関連ライブラリをインポートするためのコードです.
import numpy as np
import pandas as pd
from causalml.dataset import make_uplift_classification
from causalml.inference.tree import UpliftRandomForestClassifier
from causalml.metrics import plot_gain
from sklearn.model_selection import train_test_split
参考: Bingとのチャット
- causalml.dataset.make_uplift_classificationについて教えてください
causalml.dataset.make_uplift_classificationは、Pythonのライブラリであるcausalmlに含まれる関数の一つです。この関数は、Uplift Modelingのためのシミュレーションデータを生成するために使用されます。Uplift Modelingは、マーケティングや広告などの分野で使用される、効果的な施策を提案するための手法です。この関数は、make_classificationと同様に、分類問題のシミュレーションデータを生成することができますが、Uplift Modelingに必要な情報を含むように拡張されています
合成データの生成
以下のコードで合成データを作成します.
なお、x_namesはカラム名のリストでした.
df, x_names = make_uplift_classification()
df.head()
treatment_group_key x1_informative x2_informative x3_informative x4_informative x5_informative x6_irrelevant x7_irrelevant x8_irrelevant x9_irrelevant x10_irrelevant x11_uplift_increase x12_uplift_increase x13_increase_mix x14_uplift_increase x15_uplift_increase x16_increase_mix x17_uplift_increase x18_uplift_increase x19_increase_mix conversion treatment_effect
0 control -0.542888 1.976361 -0.531359 -2.354211 -0.380629 -2.614321 -0.128893 0.448689 -2.275192 -0.934969 0.656869 -1.315304 0.742654 1.891699 -2.428395 1.541875 -0.817705 -0.610194 -0.591581 0 0
1 treatment3 0.258654 0.552412 1.434239 -1.422311 0.089131 0.790293 1.159513 1.578868 0.166540 -0.356051 1.050526 -1.391878 -0.623243 2.443972 -2.889253 2.018585 -1.109296 -0.380362 -1.667606 0 0
2 treatment1 1.697012 -2.762600 -0.662874 -1.682340 1.217443 0.837982 1.042981 0.177398 -0.112409 1.374859 1.072329 -1.132497 1.050179 1.573054 -1.788427 1.341609 -0.749227 -2.091521 -0.471386 0 0
3 treatment2 -1.441644 1.823648 0.789423 -0.295398 0.718509 -0.492993 0.947824 -1.307887 0.123340 -1.957101 1.398966 -2.084619 0.058481 1.369439 0.422538 1.087176 -0.966666 -1.785592 -1.268379 1 1
4 control -0.625074 3.002388 -0.096288 1.938235 3.392424 -0.465860 -0.919897 -1.072592 -1.331181 -0.824105 1.398327 -1.403984 0.760430 1.917635 -2.347675 1.560946 -0.833067 -1.407884 -0.781343 0 0
メモ
みたところ、特徴量が少なくとも10種類あり、うち5個は"irrelevant"というカラム名称になっていることから、介入有無, アウトプットのいずれにも影響がないということかもしれません.
データの概観も確認しておきましょう
df.shape
(4000, 22)
df.describe(include='all').T
count unique top freq mean std min 25% 50% 75% max
treatment_group_key 4000 4 control 1000 NaN NaN NaN NaN NaN NaN NaN
x1_informative 4000.0 NaN NaN NaN 0.046091 0.991839 -3.15541 -0.620203 0.049795 0.683521 3.419016
x2_informative 4000.0 NaN NaN NaN -0.042249 1.98955 -6.973572 -1.387315 0.009056 1.330225 7.579787
x3_informative 4000.0 NaN NaN NaN -0.00527 0.988179 -3.140122 -0.662052 0.014226 0.656202 4.132909
x4_informative 4000.0 NaN NaN NaN 0.037522 1.714961 -6.092744 -1.157709 0.129414 1.277874 5.217351
x5_informative 4000.0 NaN NaN NaN 0.99504 1.311799 -3.724535 0.082171 0.960444 1.900392 6.210711
x6_irrelevant 4000.0 NaN NaN NaN -0.010057 1.009712 -3.564201 -0.692486 0.007682 0.675581 3.895835
x7_irrelevant 4000.0 NaN NaN NaN -0.021449 0.988327 -3.748029 -0.692884 -0.035329 0.652768 3.420657
x8_irrelevant 4000.0 NaN NaN NaN -0.018311 0.999329 -3.390215 -0.702189 -0.012898 0.647542 4.319756
x9_irrelevant 4000.0 NaN NaN NaN 0.005795 1.01828 -3.462369 -0.686555 0.00905 0.69578 3.404471
x10_irrelevant 4000.0 NaN NaN NaN -0.001874 0.995056 -3.352763 -0.678022 0.002378 0.657677 3.771162
x11_uplift_increase 4000.0 NaN NaN NaN 0.998046 0.388101 -1.629526 0.754488 0.999239 1.242388 3.1989
x12_uplift_increase 4000.0 NaN NaN NaN -0.905721 1.124399 -4.933662 -1.685564 -0.93732 -0.151863 2.439626
x13_increase_mix 4000.0 NaN NaN NaN 0.492439 0.802361 -2.946123 -0.057305 0.486911 1.022722 3.107777
x14_uplift_increase 4000.0 NaN NaN NaN 0.983616 0.777571 -2.164286 0.469672 0.974192 1.496734 3.620483
x15_uplift_increase 4000.0 NaN NaN NaN -0.780984 1.328718 -5.399204 -1.690261 -0.820131 0.088473 4.299328
x16_increase_mix 4000.0 NaN NaN NaN 0.810807 0.639944 -1.727029 0.387201 0.80143 1.227189 3.068957
x17_uplift_increase 4000.0 NaN NaN NaN -0.999966 0.390873 -3.401478 -1.159402 -0.996991 -0.838424 1.610457
x18_uplift_increase 4000.0 NaN NaN NaN -0.597995 0.948685 -2.677312 -1.22769 -0.849096 -0.282686 2.621344
x19_increase_mix 4000.0 NaN NaN NaN -0.982159 0.542929 -3.520206 -1.314383 -0.981452 -0.642806 2.108493
conversion 4000.0 NaN NaN NaN 0.546 0.497942 0.0 0.0 1.0 1.0 1.0
treatment_effect 4000.0 NaN NaN NaN 0.04475 0.20678 0.0 0.0 0.0 0.0 1.0
アウトカムは"conversion"ですね.
介入グループごとのアウトカムの集計値も確認してみます.
## Look at the conversion rate and sample size in each group
df.pivot_table(values = 'conversion',
index = 'treatment_group_key',
aggfunc=[np.mean, np.size],
margins=True)
| treatment_group_key | ('mean', 'conversion') | ('size', 'conversion') |
|---|---|---|
| control | 0.511 | 1000 |
| treatment1 | 0.514 | 1000 |
| treatment2 | 0.559 | 1000 |
| treatment3 | 0.6 | 1000 |
| All | 0.546 | 4000 |
メモ
上表を見ると、3種類の介入とコントロール群にそれぞれ0.25の等確率で割り当てるRCT(Randomized Controlled Trial)を想定しているように見えます. conversionは0.5近傍の数字ですが、グループ間で最大0.1弱の変化がありますね. 何のconversionなのかにも寄りますが、ある程度顧客の購買度が高まったところを分母としたconversionなどが現実の事例に近いものと理解しました.
アップリフトランダムフォレスト分類器の実行
このセクションでは、まずトレーニングデータを使用してアップリフトランダムフォレスト分類器をfit()させます。次に、適合したモデルを使用してテストデータを予測します。予測はndarrayを返し、各列には、対応する処置群に割り当てられた場合の予測アップリフトが含まれます。
はじめに、トレーニングデータとテストデータに分割します.
# 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
)
ここで以下のモジュールをインポートします。
from causalml.inference.tree import UpliftTreeClassifier
メモ
冒頭でインポートした"UpliftRandomForestClassifier"との関連性が不明ですが先に進みます.
次に、モデルを学習します.
clf = UpliftTreeClassifier(control_name='control')
clf.fit(df_train[x_names].values,
treatment=df_train['treatment_group_key'].values,
y=df_train['conversion'].values)
p = clf.predict(df_test[x_names].values)
メモ
pは(800, 4)のndarrayです.
ちなみに、x_namesはx1からx19までのリストで、"uplift_increase"や"increase_mix"が何を意味するのか分かりませんが、全て学習に利用されます.
今回はtreatment_group_keyというフィールドにcontrolというフィールドがあり、それをUpliftTreeClassifierオブジェクトを生成するときに指定するようです.
さて上述のコード内の「p」を確認してみたいと思います.
分類器のpredict()メソッドで吐き出したndarrayに分類器のclasses_という属性を付与してデータフレームを作成しています.
df_res = pd.DataFrame(p, columns=clf.classes_)
df_res.head()
| control | treatment1 | treatment2 | treatment3 | |
|---|---|---|---|---|
| 0 | 0.506394 | 0.511811 | 0.573935 | 0.503778 |
| 1 | 0.506394 | 0.511811 | 0.573935 | 0.503778 |
| 2 | 0.580838 | 0.458824 | 0.508982 | 0.452381 |
| 3 | 0.482558 | 0.572327 | 0.556757 | 0.961538 |
| 4 | 0.482558 | 0.572327 | 0.556757 | 0.961538 |
メモ
所謂SingleTreeアルゴリズム(STアルゴリズム)に該当するのか、TwoTreeアルゴリズム(TTアルゴリズム)に該当するのかは理解が及びませんでした.
ここまでで冒頭の文言の内容(以下再掲)が完了したことになりますね.
このセクションでは、まずトレーニングデータを使用してアップリフトランダムフォレスト分類器をfit()させます。次に、適合したモデルを使用してテストデータを予測します。予測はndarrayを返し、各列には、対応する処置群に割り当てられた場合の予測アップリフトが含まれます。
次に、再び同様の手順を今度はUpliftRandomForestClassifierを使用して実行します.
uplift_model = UpliftRandomForestClassifier(control_name='control')
uplift_model.fit(df_train[x_names].values,
treatment=df_train['treatment_group_key'].values,
y=df_train['conversion'].values)
df_res = uplift_model.predict(df_test[x_names].values, full_output=True)
print(df_res.shape)
df_res.head()
# (800, 9)
出力結果の表は以下のとおりです. Pandas DataFrameが直接出力されました. さらに、さらにCATE(Conditional Averate Treatment Effect)のようなフィールド(delta_xxx)が存在するのが分かります.
| control | treatment1 | treatment2 | treatment3 | recommended_treatment | delta_treatment1 | delta_treatment2 | delta_treatment3 | max_delta | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.533659 | 0.499542 | 0.556991 | 0.518805 | 2 | -0.0341178 | 0.0233312 | -0.0148546 | 0.0233312 |
次に、引数full_outputを使わずにy_predを作成します. こちらの方は、先ほどと同様ndarrayとなりましたが、shapeが(800, 4)ではなく(800, 3)に減少しています.
y_pred = uplift_model.predict(df_test[x_names].values)
そのため、後続のPandasデータフレーム生成手順も先ほどのツリーモデルとは少し異なります.
# Put the predictions to a DataFrame for a neater presentation
# The output of `predict()` is a numpy array with the shape of [n_sample, n_treatment] excluding the
# predictions for the control group.
result = pd.DataFrame(y_pred,
columns=uplift_model.classes_[1:])
result.head()
上記コードを見ると丁寧にその旨の注記がなされていますね.
以下は上記コードの出力結果です.
| treatment1 | treatment2 | treatment3 | |
|---|---|---|---|
| 0 | -0.0341178 | 0.0233312 | -0.0148546 |
| 1 | 0.0137888 | 0.109975 | 0.0281498 |
| 2 | 0.0372195 | 0.0446331 | 0.00921715 |
| 3 | 0.0560183 | 0.0963217 | 0.703216 |
| 4 | -0.0831038 | -0.0372902 | 0.197616 |
ツリーモデルの場合とは異なり、CATE(Conditional Average Treatment Effect)の予測結果が返されているように見えます.
アップリフトカーブの描画
合成された母集団の生成
はじめにbest_treatmentというnumpy.ndarray作成します.
# If all deltas are negative, assign to control
# ; otherwise assign to the treatment
# with the highest delta
best_treatment = np.where((result < 0).all(axis=1),
'control',
result.idxmax(axis=1))
numpy.where() の使い方に疎く、調べてみたところ、.all(axis=1) によって、result < 0 の列数が縮減されて1列になるようです. 各行のすべての列においてTrueであればTrue, そうでなければFalseのベクトルが返されるような理解を持ちました. また、.idxmax(axis=1)を用いて、各行について最大値を持つ列のインデックスを返しているようです. データを確認してみます.
pd.DataFrame(best_treatment).head(5)
| 0 | |
|---|---|
| 0 | treatment2 |
| 1 | treatment2 |
| 2 | treatment2 |
| 3 | treatment3 |
| 4 | treatment3 |
次に、acutual_is_best, actual_is_controlというnumpy.ndarrayを作成します.
# Create indicator variables for whether a unit happened to have the
# recommended treatment or was in the control group
actual_is_best = np.where(df_test['treatment_group_key'] == best_treatment, 1, 0)
actual_is_control = np.where(df_test['treatment_group_key'] == 'control', 1, 0)
続いて、syntheticというnumpy.ndarray, synthというpandasデータフレームを作成します.
synthetic = (actual_is_best == 1) | (actual_is_control == 1)
synth = result[synthetic]
テストデータの予測結果を表すデータフレームのうち、実際の割当がコントロール群あるいは、モデルが予測した最良の種類の介入群となっているレコードだけに絞り込まれたものがsynthということになります. データ数を確認してみましょう.
len(df_test), len(synthetic), len(synth)
(800, 800, 339)
synthの中身も見てみましょう.
synth.head(5)
| treatment1 | treatment2 | treatment3 | |
|---|---|---|---|
| 3 | 0.0560183 | 0.0963217 | 0.703216 |
| 4 | -0.0831038 | -0.0372902 | 0.197616 |
| 8 | -0.122634 | -0.0594984 | -0.0822964 |
| 10 | 0.0441345 | 0.0744583 | 0.161097 |
| 11 | 0.117743 | 0.224831 | 0.131922 |
上表の各レコードは、数値が最大の列のインデックス名の介入群、あるいはコントロール群に割り当てられているということになりますが、このデータからだといずれなのかよく分かりませんね.
アップリフトカーブを作成する
はじめに(直訳感を修正するのが面倒なので)リンクの文章の和訳を引用します.
観察された処置効果を使用して、アップリフトカーブを計算します。これは、次のような問いに答えます。「予測されたアップリフトに従って、最も高いものから最も低いものまで並べ替えられた人口のサブセットをターゲットにすることで、累積アップリフトのどれだけをキャプチャできたか？」CausalMLには、処置割り当て、観測結果、および予測された処置効果を含むDataFrameが与えられた場合にアップリフトカーブを計算するplot_gain()関数があります。この関数を使用すると、アップリフトカーブを視覚化し、マーケティングキャンペーンの最適化方法について洞察を得ることができます。
ここではauuc_metricsというpandasデータフレームを作成します.
auuc_metrics = (synth.assign(is_treated = 1 - actual_is_control[synthetic],
conversion = df_test.loc[synthetic, 'conversion'].values,
uplift_tree = synth.max(axis=1))
.drop(columns=list(uplift_model.classes_[1:]))
syntheticはresultからsynthを抜き出すための行指定のような役割を持つnumpy.ndarrayで、ここでも活躍していることが分かりますね.
auuc_metricsの中身を見てみましょう.
auuc_metrics.head(5)
| is_treated | conversion | uplift_tree | |
|---|---|---|---|
| 3 | 0 | 0 | 0.703216 |
| 4 | 0 | 1 | 0.197616 |
| 8 | 0 | 1 | -0.0594984 |
| 10 | 1 | 0 | 0.161097 |
| 11 | 0 | 0 | 0.224831 |
メモ
フィールド"is_treated"が0の場合、フィールド"uplift_tree"にはモデルが予測した介入種類ごとのCATE(Conditional Average Treatment Effect)のうち最大のCATEが表示されていることになります. 一方で"is_treated"が
1の場合には、その介入種類について予測されたCATEがフィールド"uplift_tree"に表示されていることになります.
なお、実際の介入種類と、最も効果的と予測された介入種類が異なるテストデータ中のレコードはそもそもこのデータフレームには含まれていません.
最後にアップリフトカーブを描画します.
plot_gain(auuc_metrics,
outcome_col='conversion',
treatment_col='is_treated'
)
上図でリンク先のコード例は終了となりますが、上図のy軸の"Gain"が何なのかよく分からなかたので、ピボットテーブルで確認してみます.
auuc_metrics.pivot_table(values = ['conversion', 'uplift_tree'],
index = 'is_treated',
aggfunc=[np.size, np.sum, np.mean],
margins=True
)
上図を見ると、テストーデータに所定の絞り込みを行ったデータではコントロール群と介入群でコンバージョンが約0.16異なることとが分かります. さらに、sum列のuplifet_treeを見ると、このデータの339人全体ではコンバージョンの合計が55.8増加するとモデルは予測していることが分かります. この339, 55.8という数値がアップリフトカーブ上のある1点におけるx軸, y軸座標にそれぞれ対応しているということでしょう.
とはいえ、この55.8という数値は予測値なので個人的にはやや分かりづらいです.
調べてみるとplot_gain()の他に、plot_lift()という関数もあったので試してみます.
from causalml.metrics import plot_lift
plot_lift(auuc_metrics,
outcome_col='conversion',
treatment_col='is_treated'
)
初心者としては上図のようにLiftがyになっている方が分かりやすいです.
以下サイトに関数の簡単な解説がありました.
Plot the lift chart of model estimates in cumulative population.
If the true treatment effect is provided (e.g. in synthetic data), it's calculated as the mean of the true treatment effect in each of cumulative population.
Otherwise, it's calculated as the difference between the mean outcomes of the treatment and control groups in each of cumulative population.
For details, see Section 4.1 of Gutierrez and G{'e}rardy (2016),Causal Inference and Uplift Modeling: A review of the literature.
For the former,treatment_effect_colshould be provided. For the latter, bothoutcome_colandtreatment_colshould be provided.
これをみると、ランダムのプロットはy = 定数 の直線になりそうですが、おそらく原点(0, 0)を表現するレコードがデータに挿入されているので、プロットは原点からスタートするようです. さらにy軸の数値は例えばランダムプロット(上図赤線のプロット)の場合、データをランダムに並び替えたときの最初のx個のデータについて介入群とコントロール群のアウトカムの差を平均したものなので、x軸が0に近いときはばらつきが大きくなるものと理解しました.
単にコンバージョンが1のデータと0のデータを降順で並び替えても、良いカーブが描けそうな気がするのでグラフで確認します.
import matplotlib.pyplot as plt
tmp_df = auuc_metrics.copy()
tmp_df = tmp_df.sort_values(by='uplift_tree', ascending=False)
x = tmp_df['uplift_tree'] * (-1)
tmp_df["cumsum_is_treated"] = tmp_df.cumsum()['is_treated']
y = tmp_df["cumsum_is_treated"]
plt.plot(x, y)
plt.xlabel('uplift_tree_sign_inverted')
plt.ylabel('cumsum_is_treated')
カーブを見やすくするためx軸の符号を反転させました. そのためxが小さいほど予測されたCATEが高いということになります. これを見ると、予測されたCATEが高い順に並び替えたとき、前半のデータに介入群が集中しているというわけではなさそうなので、疑問解消です.
最後に
元データがRCT(Random Controlled Trial)を想定しているのか確認してみたいと思います.
RCTであれば、特徴量の分布が割り当てによらないことが期待されるので、ヒストグラムで確認します.
import seaborn as sns
for x in x_names:
sns.displot(df, x=x, hue="treatment_group_key", element="step")
全てのヒストグラムは掲載しませんが、RCTを想定した元データであるようです.
ほとんどの特徴量の分布は釣鐘型を描いていましたが、そうでないものもあり、上記コード実行結果の1例として以下に示します.
まとめ
- causalmlというライブラリを用いてアップリフトランダムフォレストを実施することができる
- causalmlというライブラリを用いてアップリフトーカーブ(gain, lift)を簡便に描画することができる
- 上記のコード例では元データとしてRCTデータが使用されていた