CausalMLのUpliftRandomForestClassifierを用いたアップリフトモデリング



# はじめに 

本稿の目的は、以下リンク(以下チュートリアル)先チュートリアルを学習し、Pythonのライブラリであるcausalmlを用いた所謂ツリーベースのアップリフトモデリングに関する理解を深めることです. 

https://github.com/uber/causalml/blob/master/examples/uplift_trees_with_synthetic_data.ipynb

私個人としては、ゆくゆくは以下のリンクのコンテンツに挑戦したく、本稿の位置付けはそのための基礎トレーニングです. 
 
https://causal-machine-learning.github.io/kdd2021-tutorial/ 


# 環境など

- 私用Macパソコン(M2チップ)
- conda仮想環境
- VisualStudioCode

# ライブラリのインストール　

公式サイト(以下リンク)を参考にターミナル上で以下のコマンドを実行しました. 


https://github.com/uber/causalml/tree/master


```terminal 
$ conda create --name envforcausalml python=3.8
$ conda install -c conda-forge causalml 
```
::: note 
参考
消去法的に上記手法になった、という感じです. 
というのも、上記以外の方法も試みたのですが、エラーが解消できませんでした. 
:::


# 本編

## 関連ライブラリのインポート　
関連ライブラリをインポートするためのコードです.
```python 
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
``` 


::: note
参考: Bingとのチャット
- causalml.dataset.make_uplift_classificationについて教えてください
> causalml.dataset.make_uplift_classificationは、Pythonのライブラリであるcausalmlに含まれる関数の一つです。この関数は、Uplift Modelingのためのシミュレーションデータを生成するために使用されます。Uplift Modelingは、マーケティングや広告などの分野で使用される、効果的な施策を提案するための手法です。この関数は、make_classificationと同様に、分類問題のシミュレーションデータを生成することができますが、Uplift Modelingに必要な情報を含むように拡張されています
:::


## 合成データの生成 

以下のコードで合成データを作成します. 
なお、x_namesはカラム名のリストでした. 
```python
 df, x_names = make_uplift_classification()
```

```python
 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 
``` 

::: note
メモ
みたところ、特徴量が少なくとも10種類あり、うち5個は"irrelevant"というカラム名称になっていることから、介入有無, アウトプットのいずれにも影響がないということかもしれません. 
:::


データの概観も確認しておきましょう

```python
df.shape
```
``` 
(4000, 22)
```

```python
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"ですね. 

介入グループごとのアウトカムの集計値も確認してみます. 

```python 
## 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 |

::: note 
メモ
上表を見ると、3種類の介入とコントロール群にそれぞれ0.25の等確率で割り当てるRCT(Randomized Controlled Trial)を想定しているように見えます.  conversionは0.5近傍の数字ですが、グループ間で最大0.1弱の変化がありますね. 何のconversionなのかにも寄りますが、ある程度顧客の購買度が高まったところを分母としたconversionなどが現実の事例に近いものと理解しました. 
::: 



# アップリフトランダムフォレスト分類器の実行　
このセクションでは、まずトレーニングデータを使用してアップリフトランダムフォレスト分類器をfit()させます。次に、適合したモデルを使用してテストデータを予測します。予測はndarrayを返し、各列には、対応する処置群に割り当てられた場合の予測アップリフトが含まれます。


はじめに、トレーニングデータとテストデータに分割します. 


```python 
# 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
                                     )
``` 

ここで以下のモジュールをインポートします。

```python 
from causalml.inference.tree import UpliftTreeClassifier
```
::: note
メモ
冒頭でインポートした"UpliftRandomForestClassifier"との関連性が不明ですが先に進みます.  
:::


次に、モデルを学習します. 

```python 
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)
```
::: note
メモ
pは(800, 4)のndarrayです. 
ちなみに、x_namesはx1からx19までのリストで、"uplift_increase"や"increase_mix"が何を意味するのか分かりませんが、全て学習に利用されます. 

今回はtreatment_group_keyというフィールドにcontrolというフィールドがあり、それをUpliftTreeClassifierオブジェクトを生成するときに指定するようです. 
::: 


さて上述のコード内の「p」を確認してみたいと思います. 
分類器のpredict()メソッドで吐き出したndarrayに分類器のclasses_という属性を付与してデータフレームを作成しています. 


```python 
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 |


::: note
メモ
所謂SingleTreeアルゴリズム(STアルゴリズム)に該当するのか、TwoTreeアルゴリズム(TTアルゴリズム)に該当するのかは理解が及びませんでした. 
::: 

ここまでで冒頭の文言の内容(以下再掲)が完了したことになりますね. 
> このセクションでは、まずトレーニングデータを使用してアップリフトランダムフォレスト分類器をfit()させます。次に、適合したモデルを使用してテストデータを予測します。予測はndarrayを返し、各列には、対応する処置群に割り当てられた場合の予測アップリフトが含まれます。



次に、再び同様の手順を今度はUpliftRandomForestClassifierを使用して実行します. 
```python 
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)に減少しています. 


```python 
y_pred = uplift_model.predict(df_test[x_names].values)
```

そのため、後続のPandasデータフレーム生成手順も先ほどのツリーモデルとは少し異なります. 


```python 
# 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作成します. 

```python 
# 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)を用いて、各行について最大値を持つ列のインデックスを返しているようです. データを確認してみます. 

```python 
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を作成します. 

```python 
# 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データフレームを作成します. 

```python 
synthetic = (actual_is_best == 1) | (actual_is_control == 1)
synth = result[synthetic]
```
テストデータの予測結果を表すデータフレームのうち、実際の割当がコントロール群あるいは、モデルが予測した最良の種類の介入群となっているレコードだけに絞り込まれたものがsynthということになります. データ数を確認してみましょう. 


```python 
len(df_test), len(synthetic), len(synth)
```
``` 
(800, 800, 339)
```

synthの中身も見てみましょう. 

```python 
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データフレームを作成します. 


```python
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の中身を見てみましょう. 

```python 
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  |

::: note
メモ
フィールド"is_treated"が0の場合、フィールド"uplift_tree"にはモデルが予測した介入種類ごとのCATE(Conditional Average Treatment Effect)のうち最大のCATEが表示されていることになります. 一方で"is_treated"が
1の場合には、その介入種類について予測されたCATEがフィールド"uplift_tree"に表示されていることになります. 
なお、実際の介入種類と、最も効果的と予測された介入種類が異なるテストデータ中のレコードはそもそもこのデータフレームには含まれていません.  
::: 


最後にアップリフトカーブを描画します. 

```python 
plot_gain(auuc_metrics, 
          outcome_col='conversion', 
          treatment_col='is_treated'
          )
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/0d7e4d96-e359-9aa1-5ce9-711352b9aca5.png)

上図でリンク先のコード例は終了となりますが、上図のy軸の"Gain"が何なのかよく分からなかたので、ピボットテーブルで確認してみます. 

```python 
auuc_metrics.pivot_table(values = ['conversion', 'uplift_tree'], 
                        index = 'is_treated', 
                        aggfunc=[np.size, np.sum, np.mean], 
                        margins=True

)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/ac193295-c000-48fc-e08e-7ee2c5876e26.png)


上図を見ると、テストーデータに所定の絞り込みを行ったデータではコントロール群と介入群でコンバージョンが約0.16異なることとが分かります. さらに、sum列のuplifet_treeを見ると、このデータの339人全体ではコンバージョンの合計が55.8増加するとモデルは予測していることが分かります. この339, 55.8という数値がアップリフトカーブ上のある1点におけるx軸, y軸座標にそれぞれ対応しているということでしょう. 

とはいえ、この55.8という数値は予測値なので個人的にはやや分かりづらいです. 
調べてみるとplot_gain()の他に、plot_lift()という関数もあったので試してみます. 

```python 
from causalml.metrics import plot_lift

plot_lift(auuc_metrics, 
          outcome_col='conversion', 
          treatment_col='is_treated'
          )
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/b7758699-76eb-ff36-6927-ea4b8e661d92.png) 

初心者としては上図のようにLiftがyになっている方が分かりやすいです. 
以下サイトに関数の簡単な解説がありました. 

https://github.com/uber/causalml/blob/master/causalml/metrics/visualize.py

> 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_col` should be provided. For the latter, both `outcome_col` and `treatment_col` should be provided.

これをみると、ランダムのプロットはy = 定数　の直線になりそうですが、おそらく原点(0, 0)を表現するレコードがデータに挿入されているので、プロットは原点からスタートするようです. さらにy軸の数値は例えばランダムプロット(上図赤線のプロット)の場合、データをランダムに並び替えたときの最初のx個のデータについて介入群とコントロール群のアウトカムの差を平均したものなので、x軸が0に近いときはばらつきが大きくなるものと理解しました. 

単にコンバージョンが1のデータと0のデータを降順で並び替えても、良いカーブが描けそうな気がするのでグラフで確認します. 

```python 
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')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/09be9cd7-b562-d942-60ba-205f2aec6168.png)　

カーブを見やすくするためx軸の符号を反転させました. そのためxが小さいほど予測されたCATEが高いということになります. これを見ると、予測されたCATEが高い順に並び替えたとき、前半のデータに介入群が集中しているというわけではなさそうなので、疑問解消です. 

# 最後に 
元データがRCT(Random Controlled Trial)を想定しているのか確認してみたいと思います. 
RCTであれば、特徴量の分布が割り当てによらないことが期待されるので、ヒストグラムで確認します. 

```python 
import seaborn as sns
for x in x_names: 
    sns.displot(df, x=x, hue="treatment_group_key", element="step")
```
全てのヒストグラムは掲載しませんが、RCTを想定した元データであるようです. 
ほとんどの特徴量の分布は釣鐘型を描いていましたが、そうでないものもあり、上記コード実行結果の1例として以下に示します. 

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/2fce8a62-0225-20ee-82d0-16594274f602.png)




# まとめ 
- causalmlというライブラリを用いてアップリフトランダムフォレストを実施することができる 
- causalmlというライブラリを用いてアップリフトーカーブ(gain, lift)を簡便に描画することができる
- 上記のコード例では元データとしてRCTデータが使用されていた