はじめに
以下の投稿でロジスティック回帰と傾向スコアについてグラフで確認しました.
その後、傾向スコアマッチングについて著書を読み進めたので、標題のとおり手を動かして実感を得てみたいと思います.
使用するデータ
例のごとく、同一著書内の別の章で取り扱われている販促メールのデータを使用しました.
以下が大元のリンクです.
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)"
)
ロジスティック回帰
今朝（２０２３年５月２５日）は昨日までと異なりBingチャットが回答をpythonコード付きで提供してくれるようなので、Bing氏に色々と教えてもらいました.
ひょっとすると給料日で機嫌が良いのかもしれません.
以下がロジスティック回帰です.
# 必要なライブラリのインポート
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from scipy.stats import norm
# データの読み込み
data = biased_data
# 処置変数と共変量を指定
treatment = "treatment"
covariates = ["history"]
# 傾向スコアをロジスティック回帰で推定
model = LogisticRegression()
model.fit(data[covariates], data[treatment])
propensity_scores = model.predict_proba(data[covariates])[:, 1]
どうやら、scikit-learnのLogisticRegressionというモジュール内で、ロジスティック関数の出力を得るメソッドが用意されているようです.
Bingに教えてもらいましょう.
predict_proba()メソッドは、各クラスに属する確率を予測するために使用されます。このメソッドは、各クラスに属する確率を含む2次元配列を返します。predict_proba()メソッドは、分類器の予測とともに、各クラスに属する確率を返すため、確率と予測した分類結果を返すと言えます。(Bing)
回帰係数を取り出し、ロジスティック関数を自作しそれに投入しなくても便利なメソッドが用意されているということになりますね.
傾向スコアマッチング
Bingチャットによると、傾向スコアマッチングを行うための関数は以下になるようです.
本当は本稿で、マッチングされた傾向スコアを可視化するところまで辿り着きたかったのですが、以下の関数が私にとっては見慣れない佇まいなので、本稿はこの関数について理解するところまでで終了とさせていただきます.
# 傾向スコアマッチングを行う関数の定義
def propensity_score_matching(treated, control, pscores, caliper):
matched = []
unmatched = []
for i in range(len(treated)):
t_score = pscores[treated[i]]
min_dist = np.inf
min_index = -1
for j in range(len(control)):
c_score = pscores[control[j]]
dist = abs(t_score - c_score)
if dist < min_dist:
min_dist = dist
min_index = j
if min_dist <= caliper:
matched.append((treated[i], control[min_index]))
control.remove(control[min_index])
else:
unmatched.append(treated[i])
return matched, unmatched
個人的には、関数定義の最初に変数名が唐突に出てきても、どういったデータ型のものを想定しているのか、経験からすぐに理解できるわけではありません.
ですので、Bing チャットに説明をお願いしました.
-
matched = []は、matchedが空リストであることを示しています。 -
unmatched = []は、unmatchedが空リストであることを示しています。 -
for i in range(len(treated)):は、処置群treatedに含まれる各観測値についてループすることを示しています。 -
t_score = pscores[treated[i]]は、treated[i]番目の観測値に対応する傾向スコアをt_scoreに代入することを示しています。 -
min_dist = np.infは、min_distが無限大で初期化されることを示しています。 -
min_index = -1は、min_indexが-1で初期化されることを示しています。 -
for j in range(len(control)):は、対照群controlに含まれる各観測値についてループすることを示しています。 -
c_score = pscores[control[j]]は、control[j]番目の観測値に対応する傾向スコアをc_scoreに代入することを示しています。 -
dist = abs(t_score - c_score)は、t_scoreとc_scoreの差の絶対値をdistに代入することを示しています。 -
if dist < min_dist:は、distがmin_distより小さい場合に実行されます。 -
min_dist = distは、min_distにdistの値を代入することを示しています。 -
min_index = jは、min_indexにjの値を代入することを示しています。 -
if min_dist <= caliper:は、min_distがcaliper以下の場合に実行されます。 -
matched.append((treated[i], control[min_index]))は、matchedリストに(treated[i], control[min_index])のタプルを追加することを示しています。 -
control.remove(control[min_index])は、controlリストからcontrol[min_index]番目の要素を削除することを示しています。 -
else:は、上記の条件式が偽である場合に実行されます。 -
unmatched.append(treated[i])は、unmatchedリストにtreated[i]番目の要素を追加することを示しています。
なんとなくイメージが湧いたものの、よく分からない箇所もあるので、見通しを良くするためにpythonコード内に追記してもらいました.
以下は説明書きが追加された先ほどのコードです.
def match(treated, control, pscores, caliper):
"""
Propensity score matching.
Parameters:
-----------
treated: list
List of treated observations.
control: list
List of control observations.
pscores: dict
Dictionary of propensity scores for all observations.
caliper: float
Maximum distance between propensity scores for a match.
Returns:
--------
matched: list
List of matched pairs.
unmatched: list
List of unmatched treated observations.
"""
matched = [] # matchedが空リストであることを示しています。
unmatched = [] # unmatchedが空リストであることを示しています。
for i in range(len(treated)):
# 処置群treatedに含まれる各観測値についてループすることを示しています。
t_score = pscores[treated[i]]
# treated[i]番目の観測値に対応する傾向スコアをt_scoreに代入することを示しています。
min_dist = np.inf
# min_distが無限大で初期化されることを示しています。
min_index = -1
# min_indexが-1で初期化されることを示しています。
for j in range(len(control)):
# 対照群controlに含まれる各観測値についてループすることを示しています。
c_score = pscores[control[j]]
# control[j]番目の観測値に対応する傾向スコアをc_scoreに代入することを示しています。
dist = abs(t_score - c_score)
# t_scoreとc_scoreの差の絶対値をdistに代入することを示しています。
if dist < min_dist:
# distがmin_distより小さい場合に実行されます。
min_dist = dist
# min_distにdistの値を代入することを示しています。
min_index = j
# min_indexにjの値を代入することを示しています。
if min_dist <= caliper:
# min_distがcaliper以下の場合に実行されます。
matched.append((treated[i], control[min_index]))
# matchedリストに(treated[i], control[min_index])のタプルを追加することを示しています。
control.remove(control[min_index])
# controlリストからcontrol[min_index]番目の要素を削除することを示しています。
else:
unmatched.append(treated[i])
# unmatchedリストにtreated[i]番目の要素を追加することを示しています。
return matched, unmatched
冒頭は英語で説明がありますが、特にそういったお願いをしたわけではないので、始まり行の直下スペースには、引数と返り値のデータ型と対応する説明をそれぞれ英語で記載する.
実務上そうするのが常識なのだと理解しました.
また、コメントアウト部分は改行扱いとするようハンドで修正しました.
改めて確認すると、しらみつぶしにマッチングをかけ、caliper（日本語でいうとノギスのようです）を基準としてmatched, ummatchedへ配分しているというように見えます.
想像していたより単純な繰り返し作業でした.
まとめ
確かに、著書で説明があったとおりのマッチングが行われているようです.
今後については、実際に引数を投入ししてみたいと思います.