傾向スコア：備忘録

傾向スコア（Propensity Score）

傾向スコアとは、例えば、年齢、性別、基礎疾患などの背景因子（共変量）に基づいて、「ある個体が特定の介入（例：新薬の投与、特定のプログラムへの参加など）を受ける確率」のこと。0から1までの値をとる。

どういうときに使えるかというと、コホート（解析対象とする群）のマッチングに利用される。

例えば、新しいダイエットサプリの効果を調べたいとする。サプリを飲んだ人と飲まなかった人を比べるとき、もし飲んだ人にたまたま『もともと毎日激しい運動をしている若い人』ばかりが集まっていたら、痩せた原因がサプリなのか運動なのか分からなくなる。この手法は、両方のグループから『年齢や運動の習慣がそっくりな人』をペアにして抜き出すことで、周りの条件（運動など）による影響をなくし、純粋にサプリだけの効果を公平に比べられるようにするためのもの。

具体的には、傾向スコアは、複数の背景因子を「1つの確率値」に要約する。このスコアが同じ・または近い人同士は「背景因子が似ている（介入を受ける確率が同程度である）」とみなすことができる。通常は、介入の有無（あり＝1、なし＝0）を目的変数、背景因子を説明変数としたロジスティック回帰分析などを用いて算出される。

どのようなマッチングの方法があるのか

傾向スコアを用いて介入群と対照群のペアを作る（マッチングする）方法には、いくつかのアプローチがある。

最隣接マッチング (Nearest Neighbor Matching)

介入群の個体に対して、対照群の中から最も傾向スコアが近い個体をペアにする、最も直感的で代表的な方法である。

キャリパー・マッチング (Caliper Matching)

最隣接マッチングに「差の許容範囲（キャリパー）」を設ける方法である。例えば「傾向スコアの差が0.05以内」という制限をかけ、その範囲内で最も近いものを選ぶ。適切なペアが見つからない個体は解析から除外されるため、マッチングの精度が高まる（一般的に、傾向スコアのロジットの標準偏差の0.2倍などがキャリパーとして推奨されている）。

最適マッチング (Optimal Matching)

ペア全体の傾向スコアの差の合計が最小になるように、組み合わせを最適化する方法である。

カーネル・マッチング (Kernel Matching) / 局所線形マッチング

特定の1人とペアにするのではなく、傾向スコアが近い複数の対照群のデータに重み付けをして比較する方法である。

マッチング以外の方法

傾向スコアマッチング以外にも、傾向スコアを利用して偏りを調整する方法として以下のような手法がある。

IPW（逆確率重み付け法：Inverse Probability Weighting）

算出した傾向スコアの逆数を「重み」としてアウトカムに乗じることで、背景因子の偏りがない仮想的なデータセットを作り、このデータセットで比較などを行う方法である。

層別化（Stratification / Subclassification）

算出された傾向スコアの大きさに基づいて、集団全体をいくつか（一般的には5等分など）の層に分割し、背景因子が均質な各層の中で比較を行う方法である。

コホート間で偏りがないことを示す方法

マッチングを行った後は、介入群と対照群の間で背景因子（共変量）のバランスに偏りがないかを確認する必要がある。

SMD (Standardized Mean Difference: 標準化平均差)

現在、バランスの評価に最も推奨されている指標である。サンプルサイズに依存せず、各変数の平均値の差を標準偏差で割って求める。連続変数の場合、SMDは以下の式で表される（$\bar{x}$は平均、$s^2$は分散）。
$$SMD = \frac{|\bar{x}_1 - \bar{x}_2|}{\sqrt{\frac{s_1^2 + s_2^2}{2}}}$$

評価基準: 一般的に、SMDが0.1未満（あるいは10%未満）であれば、その変数に関する群間のバランスは十分に取れている（偏りがない）と判断される。

Love Plot (ラブプロット)

複数の変数ごとのSMDを、マッチング前後で比較して視覚化したグラフである。このグラフを論文やレポートに載せてバランスの改善を示すのが標準的である。

傾向スコアマッチング後の群間バランスの評価になぜp値（t検定など）を使わないのか？

従来、2群間の背景因子の違いを比較する際には、t検定などによるp値（有意差検定）がよく用いられてきた。しかし、傾向スコアマッチング後のバランス評価においてp値を用いることは推奨されていない。その理由は、p値がサンプルサイズに大きく依存する指標だからである。マッチングの過程ではペアになれなかった個体が除外されるため、マッチング後は必然的にサンプルサイズが小さくなる。その結果、マッチング後には背景因子に実質的な偏りが残っていたとしても、サンプルサイズが減少したことによって計算上「有意差なし（p > 0.05）」と判定されやすくなってしまう。したがって、マッチングが適切に行われたかどうかの確認には、p値による検定ではなく、前述のサンプルサイズに依存しない指標であるSMDを用いる。

実践例：2群（コントロール群と介入群）のコホートにおける傾向スコアマッチング

取り組む課題のストーリー設定

ある製薬会社が、高血圧に対する「新薬（介入群）」を開発しました。この新薬が「既存薬（コントロール群）」と比べて、半年後の血圧低下量（アウトカム）において優れているかを観察研究のデータから検証したいと考えているとする。
しかし、データを見ると「若い人」や「喫煙歴のない人」ほど新薬を処方されやすい傾向があった。そのまま比較すると新薬の効果を正しく測れないため、傾向スコアマッチングを用いて背景因子を揃えることにした。

このような状況が、標準的な傾向スコアマッチングのユースケースである。では、「新薬を投与した群（介入群）」と「既存薬を投与した群（コントロール群）」を比較したいとして、実際に模擬データでやってみる。

条件設定

• 割り当て (Treatment): 1 = 新薬（介入群）、0 = 既存薬（コントロール群）
• アウトカム (Outcome): 半年後の血圧低下量（数値が大きいほど効果あり）
• 交絡因子 (Covariates): 年齢、性別（1=男性/0=女性）、BMI、高血圧の家族歴（1=あり/0=なし）、喫煙歴（1=あり/0=なし）

マッチングの実行（Pythonによる実践コード）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import NearestNeighbors
from scipy import stats

# ==========================================
# 1. 模擬データの作成と交絡因子の特定
# ==========================================
np.random.seed(42)
# マッチングしない例を増やす場合はNの数を減らす
# N = 1000
N = 250

# 交絡因子の生成
age = np.random.normal(60, 10, N)
gender = np.random.binomial(1, 0.5, N)
bmi = np.random.normal(25, 4, N)
family_history = np.random.binomial(1, 0.3, N)
smoking = np.random.binomial(1, 0.2, N)

# 新薬が割り当てられる確率を計算する数式（模擬）
logit_p = 0.0 - 0.05*(age-60) - 0.1*(bmi-25) - 0.5*family_history - 0.8*smoking

prob = 1 / (1 + np.exp(-logit_p))

# 1000人のうち、新薬群（介入群）に割り当てられた人数「約115人」、残りの「約885人」が既存薬群（コントロール群）
treatment = np.random.binomial(1, prob, N) # 新薬=1, 既存薬=0

# アウトカム（血圧低下量）の生成
outcome = 10 + 5*treatment - 0.2*(age-60) - 0.3*(bmi-25) + np.random.normal(0, 3, N)

df = pd.DataFrame({
    'Treatment': treatment,
    'Outcome': outcome,
    'Age': age,
    'Gender': gender,
    'BMI': bmi,
    'Family_History': family_history,
    'Smoking': smoking
})

# ==========================================
# 2. 傾向スコアの算出
# ==========================================
covariates = ['Age', 'Gender', 'BMI', 'Family_History', 'Smoking']

scaler = StandardScaler()
scaled_covariates = scaler.fit_transform(df[covariates])

# ロジスティック回帰で傾向スコアを計算
model = LogisticRegression(solver='liblinear')
model.fit(scaled_covariates, df['Treatment'])

df['Propensity_Score'] = model.predict_proba(scaled_covariates)[:, 1]

# ==========================================
# 3. マッチングの実行（1:1 最隣接キャリパーマッチング）
# ==========================================
treatment_df = df[df['Treatment'] == 1].reset_index()
control_df = df[df['Treatment'] == 0].reset_index()

# キャリパーの設定（傾向スコアの標準偏差の0.2倍）
caliper = df['Propensity_Score'].std() * 0.2

# 最近傍探索アルゴリズム
nn = NearestNeighbors(n_neighbors=1)
nn.fit(control_df[['Propensity_Score']])

distances, indices = nn.kneighbors(treatment_df[['Propensity_Score']])

matched_treatment_idx = []
matched_control_idx = []

# キャリパーの範囲内のペアだけを採用
for i in range(len(distances)):
    if distances[i][0] <= caliper:
        matched_treatment_idx.append(treatment_df.loc[i, 'index'])
        matched_control_idx.append(control_df.loc[indices[i][0], 'index'])

matched_df = pd.concat([df.loc[matched_treatment_idx], df.loc[matched_control_idx]])

# --- ドロップした数と割合の計算と出力 ---
'''
ここでいうドロップは、介入群からドロップされた人数が対象。
'''
original_treatment_count = len(treatment_df)
matched_treatment_count = len(matched_treatment_idx)
dropped_count = original_treatment_count - matched_treatment_count
dropped_ratio = (dropped_count / original_treatment_count) * 100

print("="*50)
print("【マッチング結果のサマリー】")
print(f"介入群（新薬）の初期データ数: {original_treatment_count} 例")
print(f"マッチング成功数: {matched_treatment_count} ペア ({matched_treatment_count * 2} 例)")
print(f"ドロップした介入群の数: {dropped_count} 例 (ドロップ率: {dropped_ratio:.1f}%)")
print("="*50)

# ==========================================
# 4. グラフの描画（傾向スコアの分布 ＆ Love Plot）
# ==========================================
sns.set_theme(style="whitegrid")
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# ① マッチング前の傾向スコア分布
sns.histplot(data=df, x='Propensity_Score', hue='Treatment', bins=30, kde=True,
             palette=['blue', 'red'], ax=axes[0], alpha=0.5)
axes[0].set_title('Propensity Score Distribution (Before Matching)')
axes[0].legend(labels=['Treatment (New Drug)', 'Control (Existing Drug)'])

# ② マッチング後の傾向スコア分布
sns.histplot(data=matched_df, x='Propensity_Score', hue='Treatment', bins=30, kde=True,
             palette=['blue', 'red'], ax=axes[1], alpha=0.5)
axes[1].set_title('Propensity Score Distribution (After Matching)')
axes[1].legend(labels=['Treatment (New Drug)', 'Control (Existing Drug)'])

# ③ Love Plot の描画
def calculate_smd(data, covar, treatment_col):
    t_data = data[data[treatment_col] == 1][covar]
    c_data = data[data[treatment_col] == 0][covar]
    smd = np.abs(t_data.mean() - c_data.mean()) / np.sqrt((t_data.var() + c_data.var()) / 2)
    return smd

smd_before = [calculate_smd(df, c, 'Treatment') for c in covariates]
smd_after = [calculate_smd(matched_df, c, 'Treatment') for c in covariates]

axes[2].scatter(smd_before, covariates, color='red', label='Before Matching', marker='o')
axes[2].scatter(smd_after, covariates, color='blue', label='After Matching', marker='x')
axes[2].axvline(x=0.1, color='gray', linestyle='--', label='Threshold (0.1)')
axes[2].set_xlabel('Standardized Mean Difference (SMD)')
axes[2].set_title('Love Plot for Covariate Balance')
axes[2].legend()

plt.tight_layout()
plt.show()

# ==========================================
# 5. アウトカムの比較　(t-test)
# ==========================================
mean_diff_before = df[df['Treatment']==1]['Outcome'].mean() - df[df['Treatment']==0]['Outcome'].mean()

matched_treatment = matched_df[matched_df['Treatment']==1]['Outcome']
matched_control = matched_df[matched_df['Treatment']==0]['Outcome']
mean_diff_after = matched_treatment.mean() - matched_control.mean()

t_stat, p_val = stats.ttest_ind(matched_treatment, matched_control)

print("\n【解析結果】")
print(f"マッチング前 の血圧低下量の差: {mean_diff_before:.2f} mmHg (交絡の影響あり)")
print(f"マッチング後 の血圧低下量の差: {mean_diff_after:.2f} mmHg (RCTに近い比較)")
print(f"マッチング後の群間比較 p値: {p_val:.2e}")

• 「1:1の最隣接キャリパーマッチング」を適用し、介入群の被験者1人につき、傾向スコアが最も近いコントロール群の患者1人をペアにする。ペアになれなかったデータは解析から除外される。
• バランスの確認 (SMD): 全ての変数についてマッチング後のSMDを計算し、0.1未満になっていることをLove plotで確認する。
• アウトカムの比較: バランスが取れた（＝背景が均質な）2群間で、目的とする結果を比較する。ここでようやく、RCTに近い形での因果推論が可能になる。

実践例：交絡因子の影響を加味したアウトカムで比較したい

マッチングが「ペアを作れないデータを捨てる」手法であるのに対し、IPWは「データを1件も捨てず、一人ひとりのデータに『重み』をつけることで、仮想的に偏りのない集団（擬似母集団）を作り出す」という手法である。
本当は、1人の人間に対して「新薬を飲んだ場合のアウトカム」と「飲まなかった場合のアウトカム」の両方を同時に観察できれば、それが一番確実な効果測定であるが、現実にはどちらか一方しか観察できない。例えば、「もともと新薬を飲む確率が10%しかない若者」が実際に新薬を飲んでいた場合、その人に「10倍」の重みをつける。すると、データ上では「その若者と同じような背景を持つ人が10人いて、全員が新薬を飲んだ」という仮想のコピー（擬似母集団）が生まれる。このように全員を重み付けして膨らませることで、「偏りのない仮想的な全員の集団」を作り、その仮想集団どうしのアウトカムを比較する。

取り組む課題のストーリー設定

あるECサイト（オンラインショップ）が、売上アップのために「特別割引クーポン付きのメール（介入）」を配信した。このメールによって「翌月の購入金額（アウトカム）」がどれくらい増えたのか（メールの純粋な効果）を検証したいとする。
しかし、マーケティング担当者は「過去の購入金額が高い人」や「ログイン頻度が高い人」など、もともと優良な顧客を優先してメールを配信していた。そのまま比較すると、「メールの効果」なのか「もともと優良顧客だから買った」のかが分からない（交絡要因）。
データを捨てずに、会員全体（ATE：平均処置効果）に対するメールの効果を測るため、IPWを用いて背景因子の偏りを調整する。

条件設定

• 割り当て (Treatment): 1 = メール配信あり、0 = メール配信なし
• アウトカム (Outcome): 翌月の購入金額（円）
• 交絡因子 (Covariates): 年齢、過去の累計購入金額、月間ログイン日数

IPWの実行（Pythonによる実践コード）

IPWの「重み」を計算し、重み付け後のバランス（SMD）の確認と、効果の推定を行う。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

# ==========================================
# 1. 模擬データの作成と交絡因子の特定
# ==========================================
np.random.seed(42)
N = 2000

# 交絡因子の生成
age = np.random.normal(35, 10, N)
past_purchase = np.random.normal(50000, 20000, N)
login_days = np.random.poisson(10, N)

# メールが配信されやすい条件（傾向スコアのロジット）
# （ログインが多く、過去の購入額が大きいほどメールが届きやすい）
# logit_p = -4.0 + 0.05*(age-35) + 0.0001*(past_purchase-50000) + 0.2*(login_days-10)
logit_p = -1.0 + 0.02*(age-35) + 0.00002*(past_purchase-50000) + 0.05*(login_days-10)

prob = 1 / (1 + np.exp(-logit_p))
treatment = np.random.binomial(1, prob, N) # メールあり=1, なし=0

# アウトカム（翌月の購入金額）の生成
# メールの真の効果は「+1500円」だが、交絡因子の影響も強く受ける
outcome = 2000 + 1500*treatment + 20*(age-35) + 0.05*(past_purchase-50000) + 100*(login_days-10) + np.random.normal(0, 1000, N)

df = pd.DataFrame({
    'Treatment': treatment,
    'Outcome': outcome,
    'Age': age,
    'Past_Purchase': past_purchase,
    'Login_Days': login_days
})

# ==========================================
# 2. 傾向スコアの算出
# ==========================================
covariates = ['Age', 'Past_Purchase', 'Login_Days']

scaler = StandardScaler()
scaled_covariates = scaler.fit_transform(df[covariates])

# ロジスティック回帰で傾向スコア（メール配信確率）を計算
model = LogisticRegression(solver='liblinear')
model.fit(scaled_covariates, df['Treatment'])

df['Propensity_Score'] = model.predict_proba(scaled_covariates)[:, 1]

# ==========================================
# 3. IPW（逆確率重み付け）の計算
# ==========================================
# 介入群（メールあり）の重み = 1 / 傾向スコア
# 対照群（メールなし）の重み = 1 / (1 - 傾向スコア)
df['IPW_Weight'] = np.where(
    df['Treatment'] == 1,
    1 / df['Propensity_Score'],
    1 / (1 - df['Propensity_Score'])
)

# === 重みのトリミング（1%〜99%タイルでクリップ） ===
lower_bound = df['IPW_Weight'].quantile(0.01)
upper_bound = df['IPW_Weight'].quantile(0.99)
df['IPW_Weight'] = df['IPW_Weight'].clip(lower_bound, upper_bound)

# ==========================================
# 4. バランスの確認 (重み付け前後のSMDとLove Plot)
# ==========================================
def calculate_unweighted_smd(data, covar, treatment_col):
    t1 = data[data[treatment_col] == 1][covar]
    t0 = data[data[treatment_col] == 0][covar]
    smd = np.abs(t1.mean() - t0.mean()) / np.sqrt((t1.var() + t0.var()) / 2)
    return smd

def calculate_weighted_smd(data, covar, treatment_col, weight_col):
    t1 = data[data[treatment_col] == 1]
    t0 = data[data[treatment_col] == 0]

    # 重み付き平均の計算
    m1 = np.average(t1[covar], weights=t1[weight_col])
    m0 = np.average(t0[covar], weights=t0[weight_col])

    # 重み付き分散の計算
    v1 = np.average((t1[covar] - m1)**2, weights=t1[weight_col])
    v0 = np.average((t0[covar] - m0)**2, weights=t0[weight_col])

    smd = np.abs(m1 - m0) / np.sqrt((v1 + v0) / 2)
    return smd

# SMDの計算
smd_before = [calculate_unweighted_smd(df, c, 'Treatment') for c in covariates]
smd_after = [calculate_weighted_smd(df, c, 'Treatment', 'IPW_Weight') for c in covariates]

# Love Plotの描画
plt.figure(figsize=(8, 4))
plt.scatter(smd_before, covariates, color='red', label='Before Weighting', marker='o')
plt.scatter(smd_after, covariates, color='blue', label='After Weighting (IPW)', marker='x')
plt.axvline(x=0.1, color='gray', linestyle='--', label='Threshold (0.1)')
plt.xlabel('Standardized Mean Difference (SMD)')
plt.title('Love Plot for Covariate Balance (IPW)')
plt.legend()
plt.tight_layout()
plt.show()

# ==========================================
# 5. アウトカムの比較（効果の推定）
# ==========================================
# ① 単純な平均の差（交絡の影響を含んだ不正確な結果）
naive_effect = df[df['Treatment']==1]['Outcome'].mean() - df[df['Treatment']==0]['Outcome'].mean()

# ② IPWによる重み付け平均の差（交絡を調整した正確な結果）
weighted_mean_treated = np.average(df[df['Treatment']==1]['Outcome'], weights=df[df['Treatment']==1]['IPW_Weight'])
weighted_mean_control = np.average(df[df['Treatment']==0]['Outcome'], weights=df[df['Treatment']==0]['IPW_Weight'])
ipw_effect = weighted_mean_treated - weighted_mean_control

print("="*50)
print("【メール配信による翌月購入金額の増加効果】")
print(f"単純比較の差 (重み付け前): {naive_effect:.0f} 円 (優良顧客の偏りで過大評価)")
print(f"IPWを用いた差 (重み付け後): {ipw_effect:.0f} 円 (真の効果[+1500円]に近い)")
print("="*50)

==================================================
【メール配信による翌月購入金額の増加効果】
単純比較の差 (重み付け前): 2000 円 (優良顧客の偏りで過大評価)
IPWを用いた差 (重み付け後): 1582 円 (真の効果[+1500円]に近い)
==================================================

• IPWの計算式 (1 / 傾向スコア):例えば、「もともとメールが届く確率が10%しかなかったのに、たまたまメールが届いた人」がいたとする。この人のデータには 1 ÷ 0.1 = 10 という大きな重みが与えられる（つまり1人で10人分としてカウント）。逆に、「メールが届いて当然（確率90%）の人」の重みは 1 ÷ 0.9 ≒ 1.1 と小さくなる。この重みの調整によって、集団全体のバランスを均等に（くじ引きで決めたかのような状態に）強制補正する。
• 重み付きSMD:マッチングの時はデータを捨てた残りの集団で単純にSMDを計算するが、IPWの場合は計算式に「重み（weights）」を含めて加重平均・加重分散を計算する。Love Plotを見ると、重み付けによってSMDが0.1未満に収まっていることが確認できる。
• データ欠損ゼロ:マッチングの例ではキャリパーによってデータがドロップする可能性があるが、IPWは一人ひとりのデータを「重み」として生かすため、2000人全員のデータを使って全体の平均的な効果（ATE）を算出できる。

実践例：1群のコホートから、均質なグループのサブグループを層別化したい

層別化（Stratification / Subclassification）は「データを捨てずに、全体をいくつかの似た者同士のグループ（層）に切り分けて分析する」手法である。

取り組む課題のストーリー設定

あるIT企業（従業員2000人）が、希望者を対象とした「オンライン健康指導プログラム（介入群）」を導入した。1年後に、プログラムに参加した人と参加しなかった人の「健康スコアの改善度（アウトカム）」を比較して、プログラムの効果を検証したい。

しかし、データを見ると「年齢が高い人」や「事務職（デスクワーク）の人」、「残業時間が少ない人」ほど、健康意識が高くプログラムに参加しやすい（偏りがある）ことが分かった。
全従業員のデータを活かしてプログラムの本当の効果を測るために、「傾向スコアによる5つの層への分割（層別化）」を行い、それぞれの層を均質に含んだデータセットを作成する。

条件設定

• 割り当て (Treatment): 1 = プログラム参加、0 = 不参加
• アウトカム (Outcome): 1年後の健康スコア改善度（点数が高いほど健康になった）
• 交絡因子 (Covariates): 年齢、性別（1=男性/0=女性）、職種（1=事務職/0=その他）、月間残業時間

層別化の実行（Pythonによる実践コード）

以下のコードでは、データを傾向スコアによって5つに分割（5分位）し、各層の中で効果を計算した後、全体の結果として統合する。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression

# ==========================================
# 1. 模擬データの作成と交絡因子の特定
# ==========================================
np.random.seed(42)
N = 2000  # 従業員2000人の1群コホート

# 交絡因子の生成
age = np.random.normal(40, 10, N)
gender = np.random.binomial(1, 0.5, N)
office_worker = np.random.binomial(1, 0.6, N) # 事務職が60%
overtime = np.random.normal(30, 15, N) # 残業時間

# プログラムに参加しやすい条件（傾向スコアのロジット）
# （年齢が高い、事務職、残業が少ないほど参加しやすい）
logit_p = -1.0 + 0.05*(age-40) + 0.8*office_worker - 0.03*(overtime-30)
prob = 1 / (1 + np.exp(-logit_p))
treatment = np.random.binomial(1, prob, N) # 参加=1, 不参加=0

# アウトカム（健康スコアの改善度）の生成
# プログラムの真の効果は「+3点」だが、年齢や職種の影響も受ける
outcome = 2 + 3*treatment + 0.1*(age-40) + 1.5*office_worker - 0.05*(overtime-30) + np.random.normal(0, 2, N)

df = pd.DataFrame({
    'Treatment': treatment,
    'Outcome': outcome,
    'Age': age,
    'Gender': gender,
    'Office_Worker': office_worker,
    'Overtime': overtime
})

# ==========================================
# 2. 傾向スコアの算出
# ==========================================
covariates = ['Age', 'Gender', 'Office_Worker', 'Overtime']

scaler = StandardScaler()
scaled_covariates = scaler.fit_transform(df[covariates])

# ロジスティック回帰で傾向スコア（参加確率）を計算
model = LogisticRegression(solver='liblinear')
model.fit(scaled_covariates, df['Treatment'])

df['Propensity_Score'] = model.predict_proba(scaled_covariates)[:, 1]

# ==========================================
# 3. 傾向スコアによる層別化（5分位：Quintiles）
# ==========================================
# 傾向スコアの大きさに従って、全体を人数の等しい5つのグループ（層）に分ける
df['Stratum'] = pd.qcut(df['Propensity_Score'], q=5, labels=['Q1', 'Q2', 'Q3', 'Q4', 'Q5'])

# ==========================================
# 4. 各層でのバランス確認と効果の計算
# ==========================================
results = []
stratum_effects = []

print("="*55)
print("【各層（サブグループ）の分析結果】")

for stratum in ['Q1', 'Q2', 'Q3', 'Q4', 'Q5']:
    # その層のデータだけを抽出
    stratum_df = df[df['Stratum'] == stratum]

    # 参加者と不参加者の人数
    n_treated = len(stratum_df[stratum_df['Treatment'] == 1])
    n_control = len(stratum_df[stratum_df['Treatment'] == 0])

    # アウトカムの平均値
    mean_treated = stratum_df[stratum_df['Treatment'] == 1]['Outcome'].mean()
    mean_control = stratum_df[stratum_df['Treatment'] == 0]['Outcome'].mean()

    # その層におけるプログラムの効果（差）
    effect = mean_treated - mean_control
    stratum_effects.append(effect)

    print(f"層 {stratum}: 参加者={n_treated}人, 不参加者={n_control}人 | 改善度の差 = {effect:.2f}点")

print("="*55)

# ==========================================
# 5. 全体の効果の統合（プーリング）
# ==========================================
# 単純な全体比較（交絡の影響を含んだ不正確な結果）
naive_effect = df[df['Treatment']==1]['Outcome'].mean() - df[df['Treatment']==0]['Outcome'].mean()

# 層別化によって統合された全体効果（5つの層の効果の平均をとる）
stratified_effect = np.mean(stratum_effects)

print("【最終的な効果の比較】")
print(f"単純比較の差（層別化前）: {naive_effect:.2f} 点 (交絡の影響で真の効果(+3)より過大評価)")
print(f"層別化による効果（統合後）: {stratified_effect:.2f} 点 (真の効果に近い)")

# ==========================================
# 6. グラフの描画（層ごとの傾向スコア分布）
# ==========================================
plt.figure(figsize=(10, 6))
sns.boxplot(data=df, x='Stratum', y='Propensity_Score', hue='Treatment', palette=['blue', 'red'])
plt.title('Propensity Score Distribution by Stratum (Q1 to Q5)')
plt.xlabel('Stratum (Quintiles)')
plt.ylabel('Propensity Score')
plt.legend(title='Treatment (1:Participated, 0:Did not)')
plt.tight_layout()
plt.show()

• qcutによる5等分: pd.qcut という機能を使って、算出した傾向スコアが低い順から全体を20%ずつ、Q1〜Q5の5つの層（Stratum）に切り分けている。
• 各層の中は「均質」になる: Q1の層の中、Q2の層の中…といったそれぞれの箱の中では、プログラムに参加した人（赤）と参加しなかった人（青）の傾向スコアが非常に近い値になっている。つまり、「各層の中では背景因子の偏りがない（均質である）」状態が作られている。
• 効果の統合（プーリング）: 均質になった各層の中で「参加者と不参加者の改善度の差」を計算している（出力の「各層の分析結果」部分）。そして最後に、その5つの層の計算結果の平均をとることで、データを1件も捨てることなく、企業全体としてのプログラムの純粋な効果を算出している。

おわりに

傾向スコアは機械学習のデータセットの分割にも利用できる。
「層別化して、それぞれの層から一定の割合で分割（抽出）する」という手順になる。
機械学習の分野において、このアプローチは「傾向スコアを用いた層化抽出（Stratified Sampling）」という概念に近い。

例えば、全データを「Train : Val : Test ＝ 8 : 1 : 1」の割合で分割したいとする。
先ほどの実践例と同じように、まずデータ全体に対して傾向スコアを算出し、スコアの大きさに基づいて複数の層（例：5つの層 Q1〜Q5）に分割する。これにより、各層の中は「背景因子が似た者同士」の均質なグループになる。
そして、全体をまとめて「8:1:1」に割るのではなく、それぞれの層（Q1〜Q5）の中で個別に「8:1:1」の割合でランダムに分割する。

• Q1の中から → Train用(80%)、Val用(10%)、Test用(10%)を抽出
• Q2の中から → Train用(80%)、Val用(10%)、Test用(10%)を抽出
• ...（Q5まで繰り返す）
• 最後に、各層から抽出されたTrain用データ、Val用データ、Test用データをそれぞれ結合して、最終的な3つのデータセットを完成させる。

単純に全体をランダムに分割する（単純無作為抽出）だけでも、比較的バラバラになるのだが、偶然によって生じる偏りは説明できない。これに対して、傾向スコアで層別化した上で分割すると、Train、Val、Testのどのデータセットを見ても、「傾向スコアの分布（＝年齢や性別などの背景因子の構成比）が元の全体データと同じ」になる。

訓練、検証、テストへの分割

from sklearn.model_selection import train_test_split

# 1. 傾向スコアの算出と層の作成（前のステップのおさらい）
# df['Propensity_Score'] に傾向スコアが入っている前提
df['Stratum'] = pd.qcut(df['Propensity_Score'], q=5, labels=['Q1', 'Q2', 'Q3', 'Q4', 'Q5'])

# ==========================================
# 2. scikit-learnによる層化抽出（Stratified Sampling）
# ==========================================

# まず、全体を Train(80%) と 残り(20%) に分ける
# ★ stratifyパラメータに 'Stratum' を指定することで、Q1〜Q5の割合を保ったまま分割される
train_df, temp_df = train_test_split(
    df, 
    test_size=0.20, 
    stratify=df['Stratum'], 
    random_state=42
)

# 次に、残り(20%) を Val(10%) と Test(10%) に半分ずつ分ける
# temp_df の中でも層の割合が保たれているため、再度 stratify を指定して半分にする
val_df, test_df = train_test_split(
    temp_df, 
    test_size=0.50, 
    stratify=temp_df['Stratum'], 
    random_state=42
)

# ==========================================
# 3. 分割結果の確認
# ==========================================
print(f"全体データ数: {len(df)}")
print(f"Trainデータ数: {len(train_df)} ({(len(train_df)/len(df))*100:.0f}%)")
print(f"Valデータ数: {len(val_df)} ({(len(val_df)/len(df))*100:.0f}%)")
print(f"Testデータ数: {len(test_df)} ({(len(test_df)/len(df))*100:.0f}%)")

print("\n【Trainデータ内の層の割合】")
print(train_df['Stratum'].value_counts(normalize=True).sort_index())

print("\n【Testデータ内の層の割合】")
print(test_df['Stratum'].value_counts(normalize=True).sort_index())

交差検証の場合

from sklearn.model_selection import StratifiedKFold

skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# stratifyの代わりに、y（ターゲット）の引数にStratumを渡すことで層別化される
for fold, (train_index, val_index) in enumerate(skf.split(df, df['Stratum'])):
    train_fold = df.iloc[train_index]
    val_fold = df.iloc[val_index]
    # ここでモデルの訓練と評価を行う

傾向スコア、使っていこう。

---

Stay visionary