ロジスティック回帰と傾向スコアについてグラフで確認する

はじめに　

本稿では筆者が現在取り組んでいる著書である「効果検証入門〜正しい比較のための因果推論／計量経済学の基礎」で取り扱われているデータを用いて、傾向スコアについて確認してみたいと思います.

ロジスティック回帰によって得られたモデルを利用することで、傾向スコアを得ることができます。ここでは予測値として得られる値が重要であって、モデルのパラメータ自体は特に分析に使うことはありません。また、推定されたパラメータの値が直感に即しているかなどの解釈も、効果の分析においては特に質の保証になりません。よって、傾向スコアの推定を行なったモデルに関しては、特に解釈を行う必要はありません。安井 翔太. 効果検証入門〜正しい比較のための因果推論／計量経済学の基礎 (Japanese Edition) (p.151). Kindle 版.

ロジスティック回帰モデルのパラメータとは、ロジスティック回帰の入力となる特徴量の線形結合ということでしょうか？いずれにせよ比較的強いトーンで予測値が重要ということでしたので、気になりました.

メモ

著書の設例とは別のデータを使ってみました.
同一著書内の別の章で取り扱われている販促メールのデータを使用しました.
https://blog.minethatdata.com/2008/03/minethatdata-e-mail-analytics-and-data.html
上記リンク内で紹介されているデータは、キャンペーンメール送信の効果を検証するというものです.
メール送信の有無で消費額（spend）が有意に変化したかどうかを検証するものです.
キャンペーンメールの送信対象者はランダムに選ばれているのですが、著書内の方法に従って消費額にプラスの影響が期待できるようなバイアスをかけます.

# 関連ライブラリをインポートします
import pandas as pd
import numpy as np

# 分析用のデータを読み込みます
filepath = "http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv"
df = pd.read_csv(filepath)

# separating the original data
# to controlled/treatment groups
# that are both biased
np.random.seed(1)
obs_rate_c = 0.5
obs_rate_t = 0.5
male_df = df[df['segment'] != 'Womens E-Mail']
male_df = male_df.assign(treatment = lambda x: 
                         np.where(x['segment'] == 'Mens E-Mail', 1, 0))
biased_data = male_df.assign(obs_rate_c=np.where(
                                (male_df['history'] > 300) 
                                | (male_df['recency'] < 6) 
                                | (male_df['channel'] == "Multichannel"), 
                                obs_rate_c, 
                                1),
                             obs_rate_t=np.where(
                                  (male_df['history'] > 300) 
                                  | (male_df['recency'] < 6) 
                                  | (male_df['channel'] == "Multichannel"), 
                                  1, obs_rate_t),
                             random_number=np.random.uniform(size=len(male_df)))
biased_data = biased_data.query(
    "(treatment == 0 & random_number < obs_rate_c) "
    "| (treatment == 1 & random_number < obs_rate_t)"
    )

上図の内容については本稿の趣旨からやや逸れるので説明は割愛したいと思います.
次に、傾向スコアを算出したいと思います.

はじめに、前提条件の確認となりますが、今回使用したデータのレコードは顧客別です.
さらに、販促メールの送信有無を示す"treatment"というフィールドが(上図にて)追加されています.

ここで傾向スコアとは、"treatment" = 1, or 0 を予測するためのロジスティック関数の出力結果であり、ロジスティック関数の引数は、回帰係数で重みづけられた特徴量の線形結合ということになろうかと思います.

なお、著書内では介入の有無と、ターゲットの両方に影響を与えるという意味で、「共変量」と名付けられており、ロジスティック回帰にあたり、複数の共変量を使用していました.

本稿では単一の共変量で回帰してみたいと思います.
今回は"history"を使用してみました.
このフィールドは、介入対象の選定にバイアスをかけるときに使用したフィールド名で、なおかつターゲットである"spend"にも影響を与えると考えても不自然ではないという考えです.

早速ロジスティック回帰をしてみましょう.
Bing, Bardに聞いてみたところ、scikit-learnを使えということでした.
以下がそのコードですが、本日から急にBingチャットの歯切れが悪くなったようで、一部筆者の方で手打ちした部分があります.

import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style="darkgrid")

# データの読み込み
data = biased_data

# 説明変数と目的変数を分離する
X = data[['history']]
y = data['treatment']

# ロジスティック回帰モデルのインスタンス化
model = LogisticRegression()

# モデルのトレーニング
model.fit(X, y)

# 予測を行う
predictions = model.predict(X)

# 線形結合
coefficients = model.coef_
intercept = model.intercept_

作成したオブジェクトにデータを与えて回帰モデルができたら、coefficients, interceptという変数名で取り出せるようです.
predict()メソッドの出力は1, or 0のようでしたので、本当は他に方法があるのだと思いますが、手製のロジスティック関数を作成し、グラフを描画します.
こちらも、Bingチャットが直接的なコードの修正に応じてくれなくなったので、グラフ描画の段でかなりの手作業が発生してしまいました.

# ロジスティック関数
def logistic_function(linear_comb):
    return 1 / (1 + np.exp(-1 * linear_comb))


# グラフ用のデータをロードする
data = biased_data


# 説明変数と正解値を取得する
linear_comb_after_reg = coefficients.item() * X + intercept.item()
probabilities = logistic_function(linear_comb_after_reg)
predictions = model.predict(X)

data = data.assign(
    predictions = predictions, 
    probabilities_predicted= probabilities, 
    linear_comb = linear_comb_after_reg
)

cols = [
    'history', 
    'treatment', 
    'predictions', 
    'probabilities_predicted',
    'linear_comb'

]

data = data[cols]

fig, axes = plt.subplots(2, 2)

sns.scatterplot(ax = axes[0][0], data=data, 
            x ="history", 
            y = "probabilities_predicted", 
            hue = "predictions")
sns.scatterplot(ax = axes[0][1], data=data, 
            x ="history", 
            y = "probabilities_predicted", 
            hue = "treatment")
sns.histplot(ax = axes[1][0], data=data, 
             x = "history", hue = "predictions", 
             element="step", stat="probability", common_norm=False)
sns.histplot(ax = axes[1][1], data=data, 
             x = "history", hue = "treatment", 
             element="step", stat="probability", common_norm=False)
plt.setp(axes[0][0],xlabel = "", ylabel = "probability_predicted")
plt.setp(axes[0][1],xlabel = "", ylabel="")
plt.setp(axes[1][0],xlabel = "history", ylabel="probability")
plt.setp(axes[1][1],xlabel = "history", ylabel="")

肝心のグラフは以下のとおりです.
左上の図が傾向スコアということでしょうか？
元データの"history"の一次元上では、介入と非介入が混在しているので、当然ですが、ロジスティック関数に、1変数の線形結合(プラス切片)を入力してもうまく分離できていないようにみえます.

ここでもう一度筆者の言葉を見てみたいと思います.
予測値として得られた値が重要で、推定されたパラメータが直感に即しているかどうかについては解釈する必要がないようです.
それでは、予測値として得られた値の予測精度がイマイチで、大半の人の直感にも即していない場合はどうなのでしょうか？

ロジスティック回帰によって得られたモデルを利用することで、傾向スコアを得ることができます。ここでは予測値として得られる値が重要であって、モデルのパラメータ自体は特に分析に使うことはありません。また、推定されたパラメータの値が直感に即しているかなどの解釈も、効果の分析においては特に質の保証になりません。よって、傾向スコアの推定を行なったモデルに関しては、特に解釈を行う必要はありません。安井 翔太. 効果検証入門〜正しい比較のための因果推論／計量経済学の基礎 (Japanese Edition) (p.151). Kindle 版.

今回の"history"は、バイアスがかかっているとは言え、介入/非介入で色分けしたときに、分布のオーバーラップ部分が少なくありません.
本題の傾向スコアマッチングはまだ未学習ですが、このあたりがどう整理されるのか期待して読み進めていこうと思います.