観察研究で因果効果を推定する際、交絡バイアスを制御するために共変量を条件付ける必要があります。しかし、共変量が多い場合、すべての組み合わせに対して処置群と対照群のサンプルを確保することは困難です。10個の二値共変量があれば $2^{10} = 1024$ 通りのセルに分割されるため、各セルのサンプルサイズが極端に小さくなります。
この次元の呪いに対する画期的な解決法が、ローゼンバウムとルービン（Rosenbaum and Rubin, 1983）が提案した傾向スコア（propensity score）です。傾向スコアは多次元の共変量を1次元の値に圧縮しつつ、交絡を制御する能力を保持するという、理論的に美しい性質を持ちます。
傾向スコア法は以下の分野で広く応用されています。
- 疫学: 薬の処方と健康アウトカムの因果関係の分析
- 経済学: 教育プログラムや雇用政策の効果評価
- 社会科学: 社会プログラムの因果効果の推定
- テクノロジー: A/Bテストが不可能な場面での因果推定
本記事では、傾向スコアの理論的基盤（バランシング性とローゼンバウム・ルービンの定理）、傾向スコアマッチングとIPW推定の具体的な手順、そしてPythonでの実装を解説します。
本記事の内容
- 傾向スコアの定義とバランシング性
- ローゼンバウム・ルービンの傾向スコア定理の証明
- 傾向スコアの推定方法
- 傾向スコアマッチング
- 逆確率重み付け（IPW）推定
- Pythonによる傾向スコア法の完全な実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 交絡バイアスとその制御法 — 交絡の数理と制御法の概観
- 平均処置効果（ATE） — ATE、IPW推定量の基礎
- ロジスティック回帰 — 傾向スコアの推定に使用する回帰モデル
傾向スコアの定義と基本定理
傾向スコアの定義
傾向スコア（propensity score）$e(\bm{x})$ は、共変量 $\bm{X} = \bm{x}$ が与えられたときに処置を受ける確率として定義されます。
$$ \begin{equation} e(\bm{x}) = P(T = 1 \mid \bm{X} = \bm{x}) \end{equation} $$
名前の由来は「処置を受ける傾向（propensity）」を表す値だからです。$e(\bm{x})$ は共変量 $\bm{x}$ を持つ個体が処置群に入る確率そのものであり、0から1の値を取ります。
RCTでは傾向スコアは設計で既知です。例えば均等割り当ての場合、$e(\bm{x}) = 0.5$ です。しかし観察研究では傾向スコアは未知であり、データから推定する必要があります。
バランシング性
傾向スコアの最も重要な性質がバランシング性（balancing property）です。
$$ \begin{equation} \bm{X} \perp\!\!\!\perp T \mid e(\bm{X}) \end{equation} $$
これは「傾向スコアの値で条件付けると、共変量 $\bm{X}$ と処置 $T$ が独立になる」ことを意味します。
直感的には、傾向スコアが同じ値（例えば $e = 0.3$）を持つ個体同士では、処置を受けるか否かはランダムに決まっているようなものです。傾向スコアが等しければ、共変量の分布は処置群と対照群で同じになるのです。
バランシング性の証明は直接的です。$P(T = 1 \mid \bm{X}, e(\bm{X}))$ を計算すると、$e(\bm{X}) = P(T = 1 \mid \bm{X})$ は $\bm{X}$ の関数であるため、$\bm{X}$ で条件付ければ $e(\bm{X})$ は定数になります。したがって次のようになります。
$$ P(T = 1 \mid \bm{X}, e(\bm{X})) = P(T = 1 \mid \bm{X}) = e(\bm{X}) $$
これは $T$ が $\bm{X}$ に依存するのは $e(\bm{X})$ を通じてのみであること、すなわち $\bm{X} \perp\!\!\!\perp T \mid e(\bm{X})$ を示しています。
ローゼンバウム・ルービンの傾向スコア定理
傾向スコア定理（Rosenbaum and Rubin, 1983）は、無視可能性が共変量 $\bm{X}$ に対して成り立つならば、傾向スコア $e(\bm{X})$ に対しても成り立つことを主張します。
定理: 強い意味での無視可能性 $(Y(0), Y(1)) \perp\!\!\!\perp T \mid \bm{X}$ が成り立つならば、次が成り立つ。
$$ \begin{equation} (Y(0), Y(1)) \perp\!\!\!\perp T \mid e(\bm{X}) \end{equation} $$
証明: $P(T = 1 \mid Y(t), e(\bm{X}))$ を計算することで示します。反復期待値の法則を使い、$e(\bm{X})$ で条件付けたもとで $\bm{X}$ に対する期待値を取ります。
$$ P(T = 1 \mid Y(t), e(\bm{X})) = E[P(T = 1 \mid Y(t), \bm{X}) \mid Y(t), e(\bm{X})] $$
無視可能性の仮定 $(Y(0), Y(1)) \perp\!\!\!\perp T \mid \bm{X}$ より、$P(T = 1 \mid Y(t), \bm{X}) = P(T = 1 \mid \bm{X}) = e(\bm{X})$ です。右辺を代入すると次のようになります。
$$ P(T = 1 \mid Y(t), e(\bm{X})) = E[e(\bm{X}) \mid Y(t), e(\bm{X})] = e(\bm{X}) $$
最後の等号は $e(\bm{X})$ が $e(\bm{X})$ で条件付ければ定数であることから従います。
$P(T = 1 \mid Y(t), e(\bm{X})) = e(\bm{X}) = P(T = 1 \mid e(\bm{X}))$ が示されたので、$Y(t) \perp\!\!\!\perp T \mid e(\bm{X})$ が成り立ちます。$t = 0, 1$ の両方について同じ議論が成り立つため、$(Y(0), Y(1)) \perp\!\!\!\perp T \mid e(\bm{X})$ が証明されました。$\square$
定理の意味
この定理の実践的な意味は非常に大きいです。$p$ 次元の共変量 $\bm{X}$ で条件付ける代わりに、1次元のスカラー $e(\bm{X})$ で条件付ければ十分なのです。$p = 50$ の共変量があっても、傾向スコアという1つの数値にすべての情報が集約されます。
これにより、層別化やマッチングにおける次元の呪いが劇的に緩和されます。50次元空間での最近傍探索は困難ですが、1次元の傾向スコア上での最近傍探索は容易です。
ただし、傾向スコアは通常未知であり推定が必要です。
傾向スコアの幾何学的直感
傾向スコア定理の意味を幾何学的に理解しましょう。共変量空間 $\mathbb{R}^p$ において、各個体は点 $\bm{x}_i = (x_{i1}, x_{i2}, \ldots, x_{ip})$ として表されます。交絡の制御とは、この $p$ 次元空間の「同じ場所」にいる処置群と対照群を比較することです。
しかし、$p$ が大きい場合、同じ場所にいるペアを見つけるのは次元の呪いにより困難です。傾向スコアは、この $p$ 次元空間を $e(\bm{x})$ という1次元の直線に射影する操作であり、この射影が交絡の制御に必要な情報をすべて保持していることを保証するのがローゼンバウム・ルービンの定理です。
この射影の方向は $\bm{X}$ と $T$ の関係によって決まります。傾向スコアが等しい等値面（$e(\bm{x}) = c$ の超曲面）上では、処置群と対照群の共変量分布が同一であることがバランシング性により保証されています。したがって、同じ傾向スコアを持つ処置群と対照群を比較すれば、あたかもRCT内の比較を行っているのと同等になるのです。
正値性の仮定
傾向スコア法が機能するためには、無視可能性に加えて正値性（positivity / overlap）の仮定が必要です。
$$ 0 < P(T = 1 \mid \bm{X} = \bm{x}) < 1 \quad \text{for all } \bm{x} $$
この仮定は「すべての共変量の値において、処置を受ける確率も受けない確率もゼロでない」ことを意味します。もしある共変量パターンでは処置を受ける確率がゼロ（または1）であるなら、その部分集団では処置群と対照群の比較ができず、因果効果が識別不可能になります。
正値性の違反は、実務上しばしば見られる問題です。例えば、重症患者には必ず手術が行われ（傾向スコアが1に近い）、軽症患者には手術が行われない（傾向スコアが0に近い）場合、手術の効果を推定するための「比較可能な」サブグループが存在しないかもしれません。
次に傾向スコアの推定方法を見ていきましょう。
傾向スコアの推定
ロジスティック回帰による推定
最も標準的な傾向スコアの推定方法はロジスティック回帰です。
$$ \log\frac{e(\bm{x})}{1 – e(\bm{x})} = \bm{\beta}^\top \bm{x} $$
これはすなわち次の式になります。
$$ \hat{e}(\bm{x}) = \frac{1}{1 + \exp(-\hat{\bm{\beta}}^\top \bm{x})} $$
ロジスティック回帰は、傾向スコアが $[0, 1]$ の範囲に収まることを自然に保証し、最尤推定による良好な統計的性質を持ちます。
モデル選択と適合度
傾向スコアモデルに含める共変量の選択は、因果推定の質に直結します。一般的な指針は以下のとおりです。
- アウトカムの予測因子を含める: アウトカム $Y$ と相関のある変数を傾向スコアモデルに含めると、推定の効率性が向上します
- 処置の予測因子を含める: 処置 $T$ と相関のある変数は交絡の除去に必要です
- 操作変数的な変数を除外する: $T$ のみに影響し $Y$ には影響しない変数は、含めると推定の分散が増加します
- 合流点を除外する: $T$ と $Y$ の共通の結果を含めると合流点バイアスが生じます
バランシング診断
傾向スコアの推定後、バランシング性が実際に達成されているかを確認することが重要です。標準的な診断方法には以下のものがあります。
標準化平均差（Standardized Mean Difference, SMD）: 各共変量について、処置群と対照群の平均の差を標準偏差で割った値です。
$$ \text{SMD}_j = \frac{\bar{X}_{j,1} – \bar{X}_{j,0}}{\sqrt{(s_{j,1}^2 + s_{j,0}^2)/2}} $$
慣例として $|\text{SMD}| < 0.1$ であればバランスが良好と判断されます。傾向スコアによる調整後のSMDが調整前と比べて小さくなっていることを確認します。
分散比（Variance Ratio）: 各共変量の分散が処置群と対照群で等しいかも重要な診断指標です。分散比が0.5未満または2を超える場合、分布の形状が大きく異なっていることを示唆し、平均の一致だけでは不十分なバランシングしか達成されていない可能性があります。
傾向スコア分布の重なりの確認: 処置群と対照群の傾向スコアの密度関数を重ねてプロットし、十分な重なり（overlap）があることを視覚的に確認します。重なりが不十分な領域では因果効果の推定が外挿に依存するため、信頼性が低下します。Crump, Hotz, Imbens, and Mitnik（2009）は、傾向スコアが0.1から0.9の範囲内にある個体のみに分析を限定することを提案しました。
バランシングが達成されていれば、傾向スコアマッチングやIPW推定に進みます。
傾向スコアの層別化
マッチングやIPWの前に、傾向スコアの層別化（stratification / subclassification）について触れておきましょう。層別化は傾向スコアの最も単純な利用法であり、傾向スコアの値に基づいてデータを複数の層に分割し、各層内で処置効果を推定してから加重平均します。
Cochran（1968）の古典的な結果は、5つの層に分割するだけで交絡バイアスの約90パーセントを除去できることを示しています。実際には、5から10程度の層が推奨されます。
層別化の利点は実装の単純さと解釈の容易さです。各層内の処置効果を調べることで、効果の異質性（傾向スコアの値による処置効果の変化）を検討することもできます。ただし、層の数を増やしすぎると各層のサンプルサイズが減少し、推定が不安定になります。
傾向スコアマッチング
基本的なアイデア
傾向スコアマッチングは、処置群の各個体に対して、傾向スコアの値が近い対照群の個体をマッチさせ、マッチしたペア間でアウトカムを比較する方法です。
1対1最近傍マッチングの手順は以下のとおりです。
- すべての個体の傾向スコアを推定する
- 処置群の個体 $i$ に対して、$|\hat{e}(\bm{X}_i) – \hat{e}(\bm{X}_j)|$ が最小となる対照群の個体 $j$ を見つける
- マッチした各ペア $(i, j)$ について効果 $Y_i – Y_j$ を計算する
- 全ペアの効果を平均してATTを推定する
$$ \begin{equation} \hat{\tau}_{\text{PSM}} = \frac{1}{n_1}\sum_{i: T_i = 1}(Y_i – Y_{j(i)}) \end{equation} $$
マッチングの変種
キャリパーマッチング: 傾向スコアの差が閾値（キャリパー）$c$ 以内の対照群のみをマッチの候補とします。$|\hat{e}_i – \hat{e}_j| > c$ の場合はマッチなしとして除外します。慣例では $c = 0.2 \times \text{SD}(\text{logit}(\hat{e}))$ が推奨されます。
1対$k$マッチング: 処置群の各個体に対して$k$個の最近傍対照群をマッチさせます。バイアスとのトレードオフで$k$を選択します。$k$が大きいほどバイアスは増加しますが、分散は減少します。
復元/非復元マッチング: 同じ対照群の個体を複数回マッチに使用できる（復元）か、1回のみ使用できる（非復元）かの選択です。復元マッチングはバイアスが小さいですが、独立でないペアが生じるため分散の推定が複雑になります。
マッチングの品質評価
マッチングの品質を評価するために、以下の診断手法が用いられます。
標準化平均差（SMD）: マッチ後のSMDがすべての共変量で0.1以下であることが望ましいとされます。SMDが大きい共変量が残っている場合、傾向スコアモデルの修正（交互作用項や非線形項の追加）が必要です。
分散比（Variance Ratio）: 各共変量の分散が処置群と対照群で等しいかを確認します。分散比が0.5未満または2を超える場合、バランスが不十分と判断されます。
傾向スコアの重なり: マッチ後に処置群と対照群の傾向スコアの分布が十分に重なっていることを確認します。重なりが不十分な場合、マッチングの品質が悪く、外挿に依存した推定が行われている可能性があります。
マッチングの標準誤差
傾向スコアマッチングの標準誤差の計算は、傾向スコアが推定量であることの影響を考慮する必要があり、注意が必要です。Abadie and Imbens（2006）は、1対1マッチング推定量の分散の一致推定量を導出しました。重要な知見として、傾向スコアが既知（RCTの場合）でなく推定されている場合、推定の不確実性は減少するという結果が知られています。これは、傾向スコアの推定が「暗黙的にバランシングを改善する」ため、推定量の分散を減少させる効果があるためです。
ブートストラップ法も標準誤差の推定に使われますが、Abadie and Imbens（2008）はマッチング推定量に対するナイーブなブートストラップが一致性を持たないことを示しています。サブサンプリング法やワイルドブートストラップなどの修正された方法を使うことが推奨されます。
逆確率重み付け（IPW）の詳細
ハジェック推定量の実装
平均処置効果（ATE）で導出したIPW推定量を傾向スコアの推定値を使って実装します。ハジェック推定量は次の式で与えられます。
$$ \begin{equation} \hat{\tau}_{\text{IPW}} = \frac{\sum_{i=1}^{n}\frac{T_i Y_i}{\hat{e}_i}}{\sum_{i=1}^{n}\frac{T_i}{\hat{e}_i}} – \frac{\sum_{i=1}^{n}\frac{(1-T_i) Y_i}{1-\hat{e}_i}}{\sum_{i=1}^{n}\frac{(1-T_i)}{1-\hat{e}_i}} \end{equation} $$
極端な重みへの対処
IPW推定量の実践上の最大の課題は、傾向スコアが0や1に近い個体の重みが非常に大きくなることです。例えば $\hat{e}_i = 0.01$ の個体の重みは $1/0.01 = 100$ となり、この1つの個体の影響がATEの推定値を大きく左右してしまいます。
対処法として以下のものがあります。
トリミング（trimming）: 傾向スコアが極端な個体をサンプルから除外します。例えば $\hat{e}_i < 0.1$ または $\hat{e}_i > 0.9$ の個体を除外します。
クリッピング（clipping / truncation）: 傾向スコアの値を $[\epsilon, 1-\epsilon]$ の範囲に制限します。例えば $\hat{e}_i = \max(0.05, \min(0.95, \hat{e}_i))$ とします。
正規化（normalization）: ハジェック推定量は重みの和で正規化しているため、ホルヴィッツ・トンプソン型よりも極端な重みに対して安定しています。
IPWの漸近的性質
IPW推定量の漸近分布を整理しましょう。傾向スコアの真の値 $e(\bm{X})$ が既知の場合、ホルヴィッツ・トンプソン型IPW推定量は次の漸近分布を持ちます。
$$ \sqrt{n}(\hat{\tau}_{\text{HT}} – \tau) \xrightarrow{d} N(0, V_{\text{HT}}) $$
$$ V_{\text{HT}} = E\left[\frac{\sigma_1^2(\bm{X})}{e(\bm{X})} + \frac{\sigma_0^2(\bm{X})}{1 – e(\bm{X})} + (\tau(\bm{X}) – \tau)^2\right] $$
ここで $\sigma_t^2(\bm{X}) = \text{Var}(Y(t) \mid \bm{X})$ は条件付き分散、$\tau(\bm{X}) = E[Y(1) – Y(0) \mid \bm{X}]$ は条件付き処置効果です。
興味深いことに、傾向スコアが推定された場合、IPW推定量の漸近分散は $V_{\text{HT}}$ 以下になることが知られています（Hirano, Imbens, and Ridder, 2003）。これは傾向スコアの推定が追加情報を利用するためであり、傾向スコアを推定することが効率損失ではなく効率改善につながるという直感的に意外な結果です。
Pythonで傾向スコア法を完全に実装してみましょう。
Pythonによる傾向スコア法の実装
データ生成と傾向スコア推定
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.neighbors import NearestNeighbors
np.random.seed(42)
n = 2000
# --- データ生成 ---
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X3 = np.random.binomial(1, 0.4, n) # 二値共変量
X = np.column_stack([X1, X2, X3])
# 傾向スコア（真値）
logit_e = -0.5 + 0.8 * X1 - 0.6 * X2 + 0.4 * X3
e_true = 1 / (1 + np.exp(-logit_e))
T = np.random.binomial(1, e_true, n)
# 潜在結果
true_ATE = 3.0
Y0 = 2 + 1.5 * X1 + X2 + 0.5 * X3 + np.random.normal(0, 1, n)
Y1 = Y0 + true_ATE
Y = T * Y1 + (1 - T) * Y0
# --- 傾向スコア推定 ---
ps_model = LogisticRegression(max_iter=1000)
ps_model.fit(X, T)
e_hat = ps_model.predict_proba(X)[:, 1]
print(f"真のATE: {true_ATE}")
print(f"ナイーブ推定: {Y[T==1].mean() - Y[T==0].mean():.4f}")
print(f"傾向スコアの範囲: [{e_hat.min():.4f}, {e_hat.max():.4f}]")
推定された傾向スコアの分布と処置群・対照群のオーバーラップを確認します。
バランシング診断
import numpy as np
import matplotlib.pyplot as plt
# 前のブロックの変数を引き続き使用
# --- バランシング診断 ---
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# (a) 傾向スコアの分布
ax = axes[0]
ax.hist(e_hat[T==1], bins=30, alpha=0.5, color="red", density=True,
label="Treated")
ax.hist(e_hat[T==0], bins=30, alpha=0.5, color="blue", density=True,
label="Control")
ax.set_xlabel("Estimated propensity score", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("(a) Propensity score distribution", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) 標準化平均差（調整前後）
ax = axes[1]
var_names = ["X1", "X2", "X3"]
smd_before = []
smd_after = []
# IPWの重み
e_clip = np.clip(e_hat, 0.05, 0.95)
w1 = T / e_clip
w0 = (1 - T) / (1 - e_clip)
w1 = w1 / w1.sum()
w0 = w0 / w0.sum()
for j, name in enumerate(var_names):
xj = X[:, j]
# 調整前のSMD
s_pooled = np.sqrt((xj[T==1].var() + xj[T==0].var()) / 2)
smd_b = (xj[T==1].mean() - xj[T==0].mean()) / s_pooled
smd_before.append(smd_b)
# IPW調整後のSMD
mean_w1 = np.sum(w1 * xj)
mean_w0 = np.sum(w0 * xj)
smd_a = (mean_w1 - mean_w0) / s_pooled
smd_after.append(smd_a)
y_pos = np.arange(len(var_names))
ax.barh(y_pos - 0.2, np.abs(smd_before), 0.35, color="salmon", alpha=0.7,
label="Before adjustment")
ax.barh(y_pos + 0.2, np.abs(smd_after), 0.35, color="steelblue", alpha=0.7,
label="After IPW adjustment")
ax.axvline(0.1, color="gray", linestyle="--", linewidth=1, label="SMD = 0.1 threshold")
ax.set_yticks(y_pos)
ax.set_yticklabels(var_names, fontsize=11)
ax.set_xlabel("|Standardized Mean Difference|", fontsize=12)
ax.set_title("(b) Covariate balance", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis="x")
plt.tight_layout()
plt.savefig("propensity_score_balance.png", dpi=150, bbox_inches="tight")
plt.show()
バランシング診断の結果から、以下のことが確認できます。
-
傾向スコアの分布（左図）: 処置群と対照群の傾向スコアの分布にオーバーラップがあります。これは正値性の仮定がおおむね成り立っていることを示しています。ただし分布にずれがあり、処置群は高い傾向スコア、対照群は低い傾向スコアに偏っています。これが交絡の反映です
-
標準化平均差（右図）: 調整前（赤）はSMDが0.1を超える共変量がありますが、IPW調整後（青）はすべての共変量でSMDが0.1以下になっています。傾向スコアによるIPWが共変量のバランスを回復していることが確認できます
因果効果の推定
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.neighbors import NearestNeighbors
# 前のブロックの変数を引き続き使用
# --- 1. ナイーブ推定 ---
naive = Y[T==1].mean() - Y[T==0].mean()
# --- 2. 回帰推定 ---
reg = LinearRegression().fit(np.column_stack([T, X]), Y)
reg_est = reg.coef_[0]
# --- 3. 傾向スコアマッチング ---
treated_idx = np.where(T == 1)[0]
control_idx = np.where(T == 0)[0]
nn = NearestNeighbors(n_neighbors=1)
nn.fit(e_hat[control_idx].reshape(-1, 1))
distances, indices = nn.kneighbors(e_hat[treated_idx].reshape(-1, 1))
matched_control = control_idx[indices.flatten()]
psm_est = np.mean(Y[treated_idx] - Y[matched_control])
# --- 4. IPW (Hajek) ---
e_clip = np.clip(e_hat, 0.05, 0.95)
ipw_num = np.sum(T * Y / e_clip) / np.sum(T / e_clip)
ipw_den = np.sum((1-T) * Y / (1-e_clip)) / np.sum((1-T) / (1-e_clip))
ipw_est = ipw_num - ipw_den
# --- 5. 二重ロバスト ---
reg1 = LinearRegression().fit(X[T==1], Y[T==1])
reg0 = LinearRegression().fit(X[T==0], Y[T==0])
mu1_hat = reg1.predict(X)
mu0_hat = reg0.predict(X)
dr_est = np.mean(
(mu1_hat - mu0_hat)
+ T * (Y - mu1_hat) / e_clip
- (1-T) * (Y - mu0_hat) / (1-e_clip)
)
print(f"\n=== 推定結果 ===")
print(f"真のATE: {true_ATE:.4f}")
print(f"ナイーブ: {naive:.4f} (バイアス: {naive-true_ATE:+.4f})")
print(f"回帰: {reg_est:.4f} (バイアス: {reg_est-true_ATE:+.4f})")
print(f"PSマッチング: {psm_est:.4f} (バイアス: {psm_est-true_ATE:+.4f})")
print(f"IPW: {ipw_est:.4f} (バイアス: {ipw_est-true_ATE:+.4f})")
print(f"二重ロバスト: {dr_est:.4f} (バイアス: {dr_est-true_ATE:+.4f})")
推定結果を確認すると、ナイーブ推定が大きなバイアスを持つのに対し、傾向スコアを用いた各手法（マッチング、IPW、二重ロバスト）はいずれも真のATEに近い値を返しています。
モンテカルロシミュレーションによる比較
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.neighbors import NearestNeighbors
np.random.seed(0)
n_sims = 500
n = 500
true_ATE = 3.0
results = {"Naive": [], "Regression": [], "PS Matching": [],
"IPW": [], "DR": []}
for _ in range(n_sims):
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X = np.column_stack([X1, X2])
logit_e = -0.3 + 0.7 * X1 - 0.5 * X2
e_true = 1 / (1 + np.exp(-logit_e))
T = np.random.binomial(1, e_true, n)
if T.sum() < 10 or (1-T).sum() < 10:
continue
Y0 = 2 + 1.5 * X1 + X2 + np.random.normal(0, 1, n)
Y1 = Y0 + true_ATE
Y = T * Y1 + (1-T) * Y0
# PS推定
ps = LogisticRegression(max_iter=500).fit(X, T)
e_hat = np.clip(ps.predict_proba(X)[:, 1], 0.05, 0.95)
# ナイーブ
results["Naive"].append(Y[T==1].mean() - Y[T==0].mean())
# 回帰
reg = LinearRegression().fit(np.column_stack([T, X]), Y)
results["Regression"].append(reg.coef_[0])
# PSマッチング
t_idx = np.where(T==1)[0]
c_idx = np.where(T==0)[0]
nn = NearestNeighbors(n_neighbors=1).fit(e_hat[c_idx].reshape(-1,1))
_, idx = nn.kneighbors(e_hat[t_idx].reshape(-1,1))
results["PS Matching"].append(np.mean(Y[t_idx] - Y[c_idx[idx.flatten()]]))
# IPW
ipw = (np.sum(T*Y/e_hat)/np.sum(T/e_hat)
- np.sum((1-T)*Y/(1-e_hat))/np.sum((1-T)/(1-e_hat)))
results["IPW"].append(ipw)
# DR
r1 = LinearRegression().fit(X[T==1], Y[T==1])
r0 = LinearRegression().fit(X[T==0], Y[T==0])
m1 = r1.predict(X)
m0 = r0.predict(X)
dr = np.mean((m1-m0) + T*(Y-m1)/e_hat - (1-T)*(Y-m0)/(1-e_hat))
results["DR"].append(dr)
fig, ax = plt.subplots(figsize=(10, 5.5))
positions = []
labels = []
data = []
for i, (name, vals) in enumerate(results.items()):
arr = np.array(vals)
positions.append(i)
labels.append(name)
data.append(arr)
bp = ax.boxplot(data, positions=positions, widths=0.5, patch_artist=True,
showfliers=False)
colors = ["gray", "blue", "orange", "green", "red"]
for patch, color in zip(bp["boxes"], colors):
patch.set_facecolor(color)
patch.set_alpha(0.5)
ax.axhline(true_ATE, color="black", linewidth=2, linestyle="--",
label=f"True ATE = {true_ATE}")
ax.set_xticks(positions)
ax.set_xticklabels(labels, fontsize=10)
ax.set_ylabel("Estimated ATE", fontsize=12)
ax.set_title(f"Monte Carlo comparison ({n_sims} simulations, n={n})", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("propensity_score_monte_carlo.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"\n=== モンテカルロ結果 ===")
for name, vals in results.items():
arr = np.array(vals)
print(f"{name:15s}: Bias={arr.mean()-true_ATE:+.4f}, RMSE={np.sqrt(np.mean((arr-true_ATE)**2)):.4f}")
モンテカルロシミュレーションの結果から、以下の重要な知見が確認できます。
-
ナイーブ推定は系統的にバイアスがかかっている。箱ひげ図の中央値が真のATEから大きく外れています
-
傾向スコアベースの手法はいずれもバイアスが小さい。PSマッチング、IPW、二重ロバスト推定のすべてが真のATEの周辺に集中しています
-
DR推定量が最も効率的（箱の幅が最も狭い）。アウトカムモデルと傾向スコアモデルの情報を組み合わせることで、推定の分散が減少しています
-
IPW推定量は分散がやや大きい傾向があります。特に傾向スコアが極端な個体がある場合に、重みが大きくなって分散が増加します
傾向スコア法の発展的トピック
傾向スコアの推定方法の比較
傾向スコアの推定にはロジスティック回帰が最も一般的ですが、他の手法も提案されています。
ロジスティック回帰: 最も広く使われている方法です。パラメトリックモデルの仮定を置くため、モデルの指定が正しくなければバイアスが生じます。しかし、傾向スコアの「正しさ」はバランシング性で検証できるため、バランスが達成されていれば問題ありません。ここで重要なのは、傾向スコアモデルの目的は「予測精度を最大化すること」ではなく「共変量のバランスを達成すること」であるという点です。モデルの当てはまりが良くなくても、バランシング診断で良好な結果が得られていればそれで十分です。
プロビットモデル: ロジスティック回帰と類似の結果を与えますが、リンク関数が正規CDF（$\Phi$）です。実践上はロジスティック回帰とほぼ同等の結果を与えます。
一般化ブーステッドモデル（GBM）: McCaffrey et al.（2004）はGBMによる傾向スコア推定を提案しました。GBMはノンパラメトリックな手法であり、共変量間の非線形関係や交互作用を自動的に捉えられます。特に、バランシング性を直接最適化する目的関数を使うことで、より良いバランスが得られることがあります。
コバリアトバランシング傾向スコア（CBPS）: Imai and Ratkovic（2014）のCBPSは、傾向スコアの推定とバランシング条件を同時に最適化する方法です。通常のロジスティック回帰がモデルの当てはまりを最適化するのに対し、CBPSは共変量のバランスを直接の最適化目標とします。
極端な傾向スコアへの対処
傾向スコアが0に近い、あるいは1に近い個体は、正値性の仮定（positivity）の違反に近い状況にあります。これらの個体に対するIPWの重みは非常に大きくなり、推定量の分散を爆発させます。
トリミング（クリッピング）: 傾向スコアの値を $[\epsilon, 1-\epsilon]$（例: $[0.05, 0.95]$）にクリッピングする方法です。分散を減少させますが、バイアスを導入します。
スタビライズドウェイト: IPWの重み $1/e(X)$ の代わりに $P(T=1)/e(X)$ を使う安定化重みです。重みの分散を減少させつつ、推定の一致性を維持します。
オーバーラップ重み（OW）: Li, Morgan, and Zaslavsky（2018）が提案した方法で、重みを $h(e(X)) = e(X)(1-e(X))$ とすることで、正値性の仮定に敏感な極端な個体の影響を自動的に下げます。推定されるのはオーバーラップ母集団（処置の不確実性が高い個体）に対する効果です。
高次元での傾向スコア推定
共変量が非常に多い場合（高次元データ）、通常のロジスティック回帰では過適合が生じ、傾向スコアの推定が不安定になります。
LASSO正則化ロジスティック回帰: L1正則化により変数選択を行いつつ傾向スコアを推定します。ただし、Belloni et al.（2014）は、因果推論における変数選択ではアウトカムモデルの変数も考慮する「二重選択」（double selection）が重要であることを示しました。
高次元推定における注意点: 傾向スコアモデルに「結果に影響するがスコアに影響しない変数」を含めると推定効率が向上しますが、「スコアに影響するが結果に影響しない変数」を含めると推定の分散が増加します。この非対称性は、変数選択の指針として重要です。
一般化傾向スコア
ここまでの議論は二値処置（$T \in \{0, 1\}$）を前提としていましたが、傾向スコアの概念は多値処置（multiple treatments）や連続処置（continuous treatment）に一般化できます。
多値処置の傾向スコア: $T \in \{0, 1, 2, \ldots, K\}$ の場合、一般化傾向スコア（Generalized Propensity Score, GPS）は $e_k(\bm{x}) = P(T = k \mid \bm{X} = \bm{x})$ と定義されます。多項ロジスティック回帰で推定し、IPWでは各個体を $1/e_{T_i}(\bm{X}_i)$ で重み付けします。
連続処置の傾向スコア: $T$ が連続変数の場合、GPS は条件付き密度 $f(t \mid \bm{X} = \bm{x})$ として定義されます。Hirano and Imbens（2004）は、GPSで条件付ければ処置の割り当てが無視可能になることを示し、用量反応関数（dose-response function）の推定法を提案しました。
傾向スコア法の実践的なワークフロー
傾向スコア法を実際の研究で適用する際の標準的なワークフローを整理しましょう。
ステップ1: 研究デザインの設計。因果関係を推定したい問いを明確にし、DAGを描いて必要な共変量を特定します。処置の結果として変化する変数（悪い制御変数）は含めないように注意します。
ステップ2: 傾向スコアモデルの構築。ロジスティック回帰で傾向スコアを推定します。交互作用項や非線形項の追加を検討し、バランシング診断の結果に基づいてモデルを調整します。
ステップ3: バランシング診断。調整後のSMDがすべての共変量で0.1以下であることを確認します。傾向スコアの分布のオーバーラップを視覚的に確認し、正値性の問題がないか点検します。
ステップ4: 因果効果の推定。マッチング、IPW、または二重ロバスト推定量を用いてATEまたはATTを推定します。複数の手法で推定を行い、結果のロバスト性を確認することが推奨されます。
ステップ5: 感度分析。未観測交絡に対するロバスト性を評価します。E-valueやローゼンバウムの感度分析を報告することで、推定結果の信頼性について建設的な議論が可能になります。
感度分析との連携
傾向スコア法は強い無視可能性（$Y(0), Y(1) \perp\!\!\!\perp T \mid \bm{X}$）を前提としていますが、未観測の交絡変数が存在する場合にはこの仮定が破れます。Rosenbaum（2002）の感度分析は、未観測交絡の影響がどの程度まで結果を変えうるかを定量的に評価します。詳細は因果推論の感度分析を参照してください。
まとめ
本記事では、傾向スコア法の理論と実装を解説しました。
- 傾向スコア $e(\bm{x}) = P(T=1 \mid \bm{X}=\bm{x})$ は多次元の共変量を1次元に圧縮する量であり、バランシング性 $\bm{X} \perp\!\!\!\perp T \mid e(\bm{X})$ を持つ
- ローゼンバウム・ルービンの定理により、共変量に対する無視可能性は傾向スコアに対しても成り立つ
- 傾向スコアマッチングは傾向スコアの近い処置群・対照群をペアにして効果を推定する
- IPW推定量は傾向スコアの逆数で重み付けして母集団全体の効果を推定する
- バランシング診断（SMD）により傾向スコアの品質を検証することが不可欠
- Pythonシミュレーションにより傾向スコア法の有効性を確認した
次のステップとして、以下の記事も参考にしてください。