🐕

ERM(Expected Response Metric)によるアップリフトモデルのオフライン評価検証

に公開
因果推論
オフライン評価
tech

Pythonライブラリによる因果推論・因果探索[概念と実践] 因果機械学習の鍵を解くのP.237で紹介しているERM(Expected Response Metric)というアップリフトモデルをオフラインで評価する手法が実務に役立ちそうだったので、これを検証してみることにした。

なお厳密には間違っているかもしれないが、この記事ではUplift ModelingをHTE(異質処置効果)、CATE(条件付き処置効果)とほぼ同義として扱っている。

(さらに言うなら"処置割り当てのパーソナライズ"も同義と考えている)

ERMとは

ERMはYan Zhaoらによる論文「Uplift Modeling with Multiple Treatments and General Response Types」[Zhao et al. 2017]において紹介されたアップリフトモデルの性能評価指標になる。

以下がアップリフトモデルのスコアを計算するためのERMの計算式になる。

Z = \sum_{t=0}^K \frac{1}{P_t}Y \mathbb{I}_{h(x)=t}\mathbb{I}_{T=t}

ここで

  • K: 処置のパターン数
  • P_t: 全データ件数に対して、処置tが割り当てられた件数の比率
  • Y: 各データのアウトカム値
  • \mathbb{I}: 指示関数(条件に一致したら1、そうでない場合0を返す)
  • h(x): 方策モデルによって推薦される処置

となる。

比較的理解しやすい数式ではあるが、後に紹介する数式をコード化したものを見た方が確実に理解出来るかと思われるので、この時点で理解できなくても問題はない。

また実際のサービスのデータを用いて説明した方が分かりやすいので、まずは今回の検証に使用するMineThatDataデータセットについて紹介する。

MineThatData

Kevin HillstromのMineThatDataはメールマーケティング施策のデータセットである。

64,000人の顧客に対して3種類の処置を一様分布に基づいてランダムに割り当てたRCTデータで、以下のような列や値で構成される。

処置

列名 データ型 値
Segment カテゴリ 顧客に対して行われたメール施策
以下の3パターンからランダムに施策(=処置)が割り当てられる

No E-Mail: メールを送信しない(統制群)
Womens E-Mail: 女性向けプロモーションメールを送信
Mens E-Mail: 男性向けプロモーションメールを送信

特徴量

列名 データ型 値
Recency 数値 前回購入からの経過月数
History_Segment カテゴリ 前年度の購入額のグレード(7段階)
History 数値 前年度の購入額
Mens バイナリ 前年度に男性向け商品を購入しているか？
Womens バイナリ 前年度に女性向け商品を購入しているか？
Zip_Code カテゴリ 地域コード
Urban, Suburban, Ruralの3種類
Newbie バイナリ 直近1年以内に登録された顧客か？
Channel カテゴリ 前年度の顧客の購入チャネル
Phone, Web, Multichannelの3種類

アウトカム

列名 データ型 値
Visit バイナリ メール施策後2週間以内にWebサイト訪問が発生したか？
Conversion バイナリ メール施策後2週間以内に購入が発生したか？
Spend 数値 メール施策後2週間以内の購入金額

MineThatDataデータセットの具体的な値を確認したい場合は、この記事と併せて作成したnotebookの以下のコードの実行結果を確認してほしい。

REPOSITORY_URL = "https://raw.githubusercontent.com/PacktPublishing/Causal-Inference-and-Discovery-in-Python/refs/heads/main/data"

hillstrom = pd.read_csv(f"{REPOSITORY_URL}/hillstrom_original.csv")
print("hillstrom.shape:", hillstrom.shape)
hillstrom.head()

なお以降の説明では、Zip_CodeやChannelのワンホットエンコーディング、Segmentの値を数値に変換といった前処理がすでに完了している以下のデータセットを使用する。

特にsegmentの値は No E-Mail -> 0, Womens E-Mail -> 1, Mens E-Mail -> 2 にそれぞれ変換されていて、列の名前もtreatmentに変更されている点を覚えておいてほしい。

hillstrom_clean = pd.read_csv(f"{REPOSITORY_URL}/hillstrom_clean.csv")
print("hillstrom_clean.shape:", hillstrom_clean.shape)
hillstrom_clean.head()

検証

必要なデータも手に入れたので、ここからは本題のERMの検証を行なっていく。

ERMの計算式(再掲)と対応するPythonのコードは以下のようになる。

Z = \sum_{t=0}^K \frac{1}{P_t}Y \mathbb{I}_{h(x)=t}\mathbb{I}_{T=t}

なおコードは「Pythonライブラリによる因果推論・因果探索[概念と実践] 因果機械学習の鍵を解く」のnotebookを参考にした。

def get_expected_response(
    y_true: npt.NDArray[np.float64],
    t_true: npt.NDArray[np.int64],
    t_pred: npt.NDArray[np.int64]
) -> float:
    """アップリフトモデルのaverage expected responseを計算して返却する。
        proposed by:
        Zhao, Y., Fang, X., & Simchi-Levi, D. (2017). Uplift Modeling with Multiple Treatments and General Response Types.
        Proceedings of the 2017 SIAM International Conference on Data Mining, 588-596.
        Society for Industrial and Applied Mathematics.

        Args:
            y_true: 性能評価に使用するデータセットのアウトカムをあらわす列
            t_true: 性能評価に使用するデータセットの処置をあらわす列
            t_pred: モデルが予測した処置

        Returns:
            average expected response

        References:
            https://github.com/PacktPublishing/Causal-Inference-and-Discovery-in-Python/blob/main/Chapter_10.ipynb
    """
    proba_t = pd.Series(t_true).value_counts() / t_true.shape[0]
    treatments = proba_t.index.values

    z_vals = 0
    for treatment in treatments:
        h_indicator = t_pred == treatment
        t_indicator = t_true == treatment
        t_proba_local = proba_t[treatment]

        z_vals += (1 / t_proba_local) * y_true * h_indicator * t_indicator

    return z_vals.mean()

順を追って説明する。

まずproba_tは\frac{1}{P_t}に対応する処理で、各処置の割り当て比率を求めている。

このコードを抜き出して実行してみると以下のようになる。

t_true = hillstrom_clean["treatment"]
pd.Series(t_true).value_counts() / t_true.shape[0]

実行結果:

treatment
1    0.334172
2    0.332922
0    0.332906
Name: count, dtype: float64

一様分布に従ってランダムに処置を割り当てたデータであるため、当然ではあるが全ての処置で割り当て比率が約33%となる。

よって\frac{1}{P_t}は処置の割り当てが偏っている場合にサンプルが少ない処置の影響を補正するための係数だが、RCTデータではほとんど意味のないものであるといえる。

次にh_indicatorは\mathbb{I}_{h(x)=t}の実装で、t_indicatorは\mathbb{I}_{T=t}の実装となる。

この2つの条件を組み合わせることで、MinThatDataの各顧客データにおいて顧客に実際に割り当てられた処置と顧客の属性を元にモデルが選択した処置が一致した場合にのみ、実際のアウトカムがスコアとして加算されるという計算に繋がる。

以上で計算式の説明は完了となるが、これだとまだイメージがしづらいと思われる。

よって試しに以下のような条件のMineThatData仕様のダミーデータを作って、ERMの計算式の挙動を追ってみることにする。

  • 各処置5人ずつ、合計15人の顧客を想定
  • 処置:0では1/5人がコンバージョン、処置:1では2/5人がコンバージョン、処置:2では4/5人がコンバージョンが発生

コードで表すと以下のようになる。

def calc_expected_response(t_pred: npt.NDArray[np.int64]) -> None:
    # 各処置5人ずつ、合計15人の顧客を想定
    t_true = np.array([0, 0, 0, 0, 0] + [1, 1, 1, 1, 1] + [2, 2, 2, 2, 2])
    # 処置:0では1/5人がコンバージョン、処置:1では2/5人がコンバージョン、処置:2では4/5人がコンバージョンが発生
    y_true = np.array([1, 0, 0, 0, 0] + [1, 1, 0, 0, 0] + [1, 1, 1, 1, 0])
    # MineThatData仕様のダミーデータと引数で渡された任意の処置リストでERMスコアを計算
    expected_response = get_expected_response(y_true, t_true, t_pred)
    print(f"Expected Response: {expected_response:.4f}")

まず初めにMineThatData仕様のダミーデータの処置(t_true)と、モデルが選択した処置リスト(t_pred)が完全に一致した場合のERMのスコアを確認してみる。

calc_expected_response(t_pred=np.array([0] * 5 + [1] * 5 + [2] * 5))

Expected Response: 1.4000

これは以下の計算が内部的に行われた結果になる。

処置:0 のスコア = 1 / (5 / 15) * 1 = 3.0
処置:1 のスコア = 1 / (5 / 15) * 1 * 2 = 6.0
処置:2 のスコア = 1 / (5 / 15) * 1 * 4 = 12.0
ERMスコア = (3.0 + 6.0 + 12.0) / 15 = 1.4

一方でMineThatData仕様のダミーデータの処置とモデルが選択した処置リストが全く一致しなかった場合もちろんスコアはゼロになる。

calc_expected_response(t_pred=np.array([2] * 5 + [0] * 5 + [1] * 5))

Expected Response: 0.0000

これだと評価に使用するデータセットに対して忠実なイエスマンであればあるほどスコアが増えるチャンスが多くなる一方で、どれだけ処置の最適化を頑張ってもデータセットと意見が合わないとスコアはもらえないといった世知辛い問題が発生するのではないかと思える。

しかしこれはそもそもテスト設計が間違っているだけで、スコア加算のチャンスはデータセットのサンプルサイズが一定以上あればかなり公平に訪れる上に、ERMスコアはモデルの性能を正しく評価できる指標であると考えて間違いはないと言っておきたい。

このあたりの解釈に関しては、説明がかなりややこしい話である上にかなりの部分で自分の推測を含んだ話となっているため、興味がある人だけ次の節を確認してほしい。

ERMの解釈

(以下の説明は完全に個人的な解釈であるため、信憑性にやや疑問有)

スコア加点のチャンスが一度も訪れないケースとは、式で考えると \mathbb{I}_{h(x)=t}\mathbb{I}_{T=t} が全てゼロを返す状況と言い換えられる。

しかしそもそも評価に使うデータセットがRCTで収集したデータである以上、処置の割り当てはランダムに均等であるため、どんな処置選択を行うモデルであっても常に約3分の1の確率で処置は一致するはずである。

これを確かめるには以下のコードを実行すればよい。

# 試しにproba_tに[0.50, 0.35, 0.15]や[1.00, 0.00, 0.00]などを指定しても、処置が一致する確率は常に約33%になる
proba_t = [0.33, 0.33, 0.34]
# モデルの処置選択結果に見立てた値
t_pred = np.random.default_rng(0).choice([0, 1, 2], size=hillstrom_clean.shape[0], p=proba_t)
# 処置が一致した数 / データ総数
print((hillstrom_clean["treatment"] == t_pred).sum() / len(hillstrom_clean))

実行結果:

0.334125

つまり上記テストのように15回連続処置が一致しないという事態が起きる可能性はけしてゼロではないが、その後必ず揺り戻しが起きるはずというのが答えとなる。

なお特徴量を一切考慮せずに話をしているが、特徴量 x 処置 x アウトカムの関係を上手く捉えたなんかすごいアルゴリズムで作られたモデルで処置を割り当てたとしても、最終成果物は上記コードのt_predの形になるため結局は同じ話になる。

ちなみにRCTではなく処置が偏っているデータセットで評価を行った場合でも、式の係数\frac{1}{P_t}で補正が行われるため、スコア加点のチャンスは偏るかもしれないがスコアの補正は行われるので一応問題にはならないような気がしている。

次にERMスコアはモデルの性能を正しく評価できる指標であるという根拠について考えてみる。

まずシンプルに問題を捉えるために、顧客を「前年度に女性向け商品を購入しているか？」の2パターンだけで分ける、つまりは特徴量としてwomansだけを持つと考えてみる。

ここで特徴量womansと各処置毎のコンバージョン数やCVRを求めると以下のようになる。

hillstrom_clean.groupby(["womens", "treatment"]).agg(
    cv_cnt=("conversion", "sum"),
    total_cnt=("conversion", "count"),
    cvr=("conversion", lambda x: x.sum() / x.count() if x.count() > 0 else 0)
).reset_index()
womens treatment cv_cnt total_cnt cvr
0 0 53 9638 0.00549907
0 1 59 9622 0.00613178
0 2 110 9558 0.0115087
1 0 69 11668 0.00591361
1 1 130 11765 0.0110497
1 2 157 11749 0.0133628

なおこの集計はRCTデータを使って行われたものとなる。

上記のようにRCTデータで処置ごとのCVRを求めるのは、A/Bテストや選択バイアスのない処置効果推定の方法として正しい。

そしてここでのCVRはERMスコアを加算できる確率と同じである。

つまり最終的なERMスコアが大きいほど、トータルでのCVRが高くなるような処置選択を行なっているモデルといえる。

仮に特徴量を増やした場合、例えばmensを増やしたところで上記の理屈は変わらない。

一応特徴量にmensを追加した場合の集計結果を記載しておくが形式はほぼ同じである。

womens mens treatment cv_cnt total_cnt cvr
0 1 0 53 9638 0.00549907
0 1 1 59 9622 0.00613178
0 1 2 110 9558 0.0115087
1 0 0 49 9519 0.0051476
1 0 1 99 9647 0.0102623
1 0 2 104 9568 0.0108696
1 1 0 20 2149 0.00930665
1 1 1 31 2118 0.0146364
1 1 2 53 2181 0.0243008

なおwomensとmensが共に0のデータはMineThatDataデータセット上に存在しないため、その特徴量を持つ顧客が新たに大量に現れたデータで処置割り当てモデルを学習した場合については対策が必要になりそうな気はしている。

いずれにせよ、ここからさらに特徴量を増やしていっても特徴量の組み合わせのパターンを網羅するだけの必要十分なサンプルサイズさえあれば拡大可能な理屈であると考えている。

ERMによるモデル評価

最後に簡単なアップリフトモデルを作成したのち、ERMを使用してランダムな処置割り当てと性能を比較してみる。

まずMineThatDataデータセットを前処理したのち、モデル学習用とERMによる性能評価用に分割する。

(以降のコードは先に紹介した本のnotebookのコードを基本的に参考にしている)

# 多重共線性を考慮して冗長な列を削除
hillstrom_clean = hillstrom_clean.drop(["zip_code__urban", "channel__web"], axis=1)

# 特徴量、アウトカム、処置
hillstrom_X = hillstrom_clean.drop(["visit", "conversion", "spend", "treatment"], axis=1)
hillstrom_Y = hillstrom_clean["conversion"]
hillstrom_T = hillstrom_clean["treatment"]

# 学習用と評価用に半分に分割
X_train, X_test, Y_train, Y_test, T_train, T_test = train_test_split(
    hillstrom_X,
    hillstrom_Y,
    hillstrom_T,
    test_size=.5
)

次にT-Learnerを使用して特徴量に対する各処置の処置効果量の推定を行うモデルを学習する。

各処置の処置効果量の予測用モデルにはLightGBMの回帰を使用している。

(なおT-Learner選んだ理由は完全に個人の趣味となる)

def create_model(
    random_state: int = 0,
    n_estimators: int = 100,
    max_depth: int = 10,
    learning_rate: float = 0.01
) -> LGBMRegressor:
    return LGBMRegressor(
        random_state=random_state,
        n_estimators=n_estimators, 
        max_depth=max_depth, 
        learning_rate=learning_rate)

t_learner = TLearner(models=[create_model() for _ in range(3)])
t_learner.fit(Y=Y_train, T=T_train, X=X_train)

これで全顧客の特徴量に対して各処置の処置効果量をT-Learnerのモデルから求める準備が整ったことになる。

最後に各顧客の処置効果量を渡すと最適な処置を選択してくれる関数を以下のように定義する。

def get_effects_argmax(
    effects_arrays: list[npt.NDArray[np.float64]],
    return_matrix: bool = False
) -> npt.NDArray[np.int64] | tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]:
    """各行における処置効果の一番大きい処置のインデックスを返す

    Args:
        effects_arrays: 統制群と各処置群の処置効果量の差をあらわす配列 # shape: (顧客数, 処置数(ただし統制群除く))
        return_matrix: 各行の処置効果量を完全に保持する行列を戻り値として返却するか

    Returns:
        各行の処置効果量最大の処置をあらわすインデックスの配列
    
    References:
        https://github.com/PacktPublishing/Causal-Inference-and-Discovery-in-Python/blob/main/Chapter_10.ipynb
    """
    
    n_rows = effects_arrays[0].shape[0]
    null_effect_array = np.zeros(n_rows)
    stacked = np.stack([null_effect_array] + effects_arrays).T
    
    if return_matrix:
        return np.argmax(stacked, axis=1), stacked
    
    return np.argmax(stacked, axis=1)

あとはT-Learnerモデルの処置割り当てとランダムな処置割り当てのERMスコアを計算して、検証は完了となる。

# 各顧客の特徴量に対する、処置0(No-Email)と処置1(Womens E-Mail)の予測処置効果量の差、
# 処置0と処置2(Mens E-Mail)の予測処置効果量の差を求めている
effects_argmax = get_effects_argmax([t_learner.effect(X_train, T0=0, T1=t) for t in [1, 2]])

# T-Learnerモデルで処置を割り当てた場合のERMスコアを計算
expected_response = get_expected_response(y_true=Y_test, t_true=T_test, t_pred=effects_argmax)
print(f"[T-Learner] Expected Response: {expected_response}")

# ランダムに処置を割り当てた場合のERMスコアを計算
expected_response = get_expected_response(
    y_true=Y_test, 
    t_true=T_test, 
    t_pred=np.random.default_rng(0).choice([0, 1, 2], size=len(T_test))
)
print(f"[Random] Expected Response: {expected_response}")

実行結果:

[T-Learner] Expected Response: 0.011936667524715741
[Random] Expected Response: 0.009084576366055452

一切チューニングなどを行わなくても流石にランダムな処置割り当てよりはマシなモデルであることが確認できた。

おわり

上記について何か間違いがあれば遠慮なく指摘していただきたい。

Discussion