相関から因果へ進むには何を仮定するか：傾向スコアマッチングで選択バイアスを測る

はじめに

「相関は因果関係を意味しない」とはよく言われますが、では手元のデータから因果関係をどう推し量ればよいのでしょうか。例えば教育プログラムの参加者が非参加者より成績が良かったとします。それはプログラムの効果なのか、もともと成績の良い生徒が参加しただけなのか。観察データを眺めるだけでは区別がつきません。以下では、この区別をつけるための基本的な考え方を整理したうえで、傾向スコアマッチングをPythonで実装し、選択バイアスが推定をどれだけ歪めるかを数値で確かめます。

対象読者:

• 因果推論に初めて触れる方
• 統計学の基本的な知識はあるが、因果推論は初めてという方
• Pythonを使ってデータ分析を行っている方

記事のポイント:

• 因果推論の基本的な考え方と用語を整理する
• 選択バイアスがなぜ推定を歪めるのかを理解する
• 傾向スコアマッチングの基本的な仕組みを知る
• Pythonの実験で、バイアスの大きさとその補正を数値で確かめる

因果推論の基本的な考え方

因果推論の目的は、「もし別の選択をしていたら、結果はどうなっていたか」という反事実（counterfactual）の問いに答えることです。しかし、同じ人に対して異なる処置（例えば、薬の投与と非投与）を同時に試すことはできません。これが因果推論の出発点にある問題です。

潜在的結果フレームワーク

この問題に対処するために、因果推論では「潜在的結果（potential outcomes）フレームワーク」という考え方を導入します。例えば、ある教育プログラムの効果を検証する場合を考えてみましょう。

• $Y_i(1)$: 個人 i がプログラムを受講した場合の結果（例：テストの点数）
• $Y_i(0)$: 個人 i がプログラムを受講しなかった場合の結果

このとき、個人 i に対する因果効果（処置効果）は、次のように定義されます。

$\tau_i = Y_i(1) - Y_i(0)$

しかし、実際に観察できるのは、$Y_i(1)$ と $Y_i(0)$ のどちらか一方だけです。これを「因果推論の根本問題」と呼びます。

用語の整理

以降で使う用語を整理しておきます。

用語	説明
処置 (Treatment)	因果効果を調べたい介入のこと（例：薬の投与、教育プログラムの実施）
結果 (Outcome)	処置によって影響を受ける変数（例：病気の回復、テストの点数）
平均処置効果 (ATE)	集団全体における平均的な因果効果。 $ATE = E[Y(1) - Y(0)]$ で表されます。

因果効果推定の課題：選択バイアス

現実のデータでは、処置を受けるかどうか（例えば、プログラムに参加するかどうか）は、ランダムに決まるわけではありません。例えば、教育プログラムの効果を調べたい場合、次のようなことが起こりえます。

• もともと成績の良い生徒がプログラムに参加しやすい
• 学習意欲の高い生徒がプログラムに参加しやすい

このような場合、プログラム参加者と非参加者の間で、単純に結果（テストの点数）を比較しても、プログラムの真の効果は測れません。もともとの能力や意欲の差が結果に影響している可能性があるためです。これを「選択バイアス」と呼びます。

選択バイアスが存在する場合、単純な平均の差は以下のように分解できます。

$E[Y|T=1] - E[Y|T=0] = \underbrace{E[Y(1) - Y(0)]}_{\text{真の効果}} + \underbrace{E[Y(0)|T=1] - E[Y(0)|T=0]}_{\text{選択バイアス}}$

因果効果の推定方法

選択バイアスへの対処として、まず挙がるのがランダム化比較試験（RCT）です。処置群と対照群をランダムに割り当てれば選択バイアスは生じず、最も信頼性の高い方法になりますが、コストや倫理面の制約で実施できない場面も多くあります。

RCTが使えない場合は、観察データから推定するしかありません。代表的な手法は、共変量の分布が似ている個体同士を対応づける傾向スコアマッチング、処置には影響するが結果には直接影響しない変数（操作変数）を使う操作変数法、処置前後の変化の差を処置群と対照群で比べる差分の差分法です。以下ではこのうち傾向スコアマッチングを取り上げます。

教育プログラムの効果をシミュレーションで確かめる

ここからは、教育プログラムの効果検証を例に、選択バイアスのあるデータで何が起きるかをPythonで確かめます。真の効果が分かっている合成データを使うので、推定値がどれだけずれたかを答え合わせできます。

選択バイアスを仕込んだデータを作る

まず、選択バイアスが存在する状況をシミュレーションしたデータを作成します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

def generate_synthetic_data(n_samples=1000):
    """教育プログラムの効果を検証するための合成データを生成"""
    # 共変量（年齢、過去の成績）
    age = np.random.normal(20, 2, n_samples)
    previous_score = np.random.normal(70, 10, n_samples)

    # 処置への割り当て（選択バイアスあり）
    propensity = 1 / (1 + np.exp(-(age - 20) * 0.1 - (previous_score - 70) * 0.03))
    treatment = np.random.binomial(1, propensity)

    # 潜在的結果の生成
    y0 = 60 + 0.5 * previous_score + np.random.normal(0, 5, n_samples)
    treatment_effect = 15 - 0.1 * previous_score  # 成績が低い学生ほど効果が大きい
    y1 = y0 + treatment_effect + np.random.normal(0, 2, n_samples)

    # 観察される結果
    y_obs = treatment * y1 + (1 - treatment) * y0

    return pd.DataFrame({
        'age': age,
        'previous_score': previous_score,
        'treatment': treatment,
        'outcome': y_obs,
        'true_effect': y1 - y0
    })

このコードでは、以下の状況をシミュレートしています。

• 共変量
  • 年齢: 平均20歳、標準偏差2
  • 過去の成績: 平均70点、標準偏差10
• 選択バイアス
  • 年齢が高いほどプログラムに参加しやすい
  • 過去の成績が良いほどプログラムに参加しやすい
• 異質な処置効果
  • 過去の成績が低い生徒ほど、プログラムの効果が大きい

参加者と非参加者で共変量がどう違うか

生成したデータで、選択バイアスがどう見えるかを確認します。

def plot_selection_bias(df):
    """選択バイアスの可視化"""
    plt.figure(figsize=(10, 6))
    sns.scatterplot(data=df, x='previous_score', y='treatment', alpha=0.5)
    plt.title('過去の成績とプログラム参加の関係（選択バイアス）')
    plt.xlabel('過去の成績')
    plt.ylabel('プログラム参加（1=参加, 0=非参加）')
    plt.savefig('selection_bias.png')
    plt.close()

この図を見ると、過去の成績が高い生徒ほど、プログラムに参加している傾向がはっきりとわかります。

平均の差をそのまま取るとどうなるか

まず、プログラム参加者と非参加者の平均をそのまま比べます。

def naive_ate(df):
    """単純な平均処置効果の推定"""
    treated = df[df['treatment'] == 1]['outcome'].mean()
    control = df[df['treatment'] == 0]['outcome'].mean()
    return treated - control

結果は、傾向スコアマッチングの結果と合わせて後で並べて確認します。

傾向スコアでペアを作って推定する

次に、選択バイアスに対処するため、傾向スコアマッチングを行います。

from sklearn.linear_model import LogisticRegression

def estimate_propensity_scores(df):
    """傾向スコア（処置を受ける確率）を推定"""
    X = df[['age', 'previous_score']]
    treatment = df['treatment']

    model = LogisticRegression()
    model.fit(X, treatment)
    return model.predict_proba(X)[:, 1]

def match_and_estimate_ate(df):
    """傾向スコアマッチングによる処置効果の推定"""
    df['propensity_score'] = estimate_propensity_scores(df)

    treated = df[df['treatment'] == 1]
    control = df[df['treatment'] == 0]

    # 最近傍マッチング
    matched_pairs = []
    for _, treated_unit in treated.iterrows():
        distances = abs(control['propensity_score'] - treated_unit['propensity_score'])
        best_match_idx = distances.idxmin()
        matched_pairs.append({
            'treated': treated_unit['outcome'],
            'control': control.loc[best_match_idx, 'outcome'],
            'treated_ps': treated_unit['propensity_score'],
            'control_ps': control.loc[best_match_idx, 'propensity_score']
        })

    return np.mean([p['treated'] - p['control'] for p in matched_pairs]), matched_pairs

以下のコードで、マッチングの質を確認します。

def plot_matching_result(pairs):
    """マッチングの質の可視化"""
    plt.figure(figsize=(10, 6))
    treated_ps = [p['treated_ps'] for p in pairs]
    control_ps = [p['control_ps'] for p in pairs]

    plt.scatter(treated_ps, control_ps, alpha=0.5)
    plt.plot([0, 1], [0, 1], 'r--')
    plt.title('マッチングの質の確認')
    plt.xlabel('処置群の傾向スコア')
    plt.ylabel('対照群の傾向スコア')
    plt.savefig('matching_quality.png')
    plt.close()

この図は、マッチングされたペアの傾向スコアを並べたものです。点が45度線に近いほど、処置群の各個体に対して傾向スコアの近い相手を対照群から見つけられています。

二つの推定値を並べる

推定方法	推定されたATE	真のATE
ナイーブな推定	9.45	7.80
傾向スコアマッチングによる推定	7.65	7.80

ナイーブな推定は選択バイアスの影響で効果を過大評価しています。傾向スコアマッチングでは、真のATEに近い推定値が得られています。

1.65ポイントの過大評価はどこから来たか

単純な比較は、教育プログラムの効果を約1.65ポイント過大評価していました。成績の良い生徒ほどプログラムに参加しやすい設定にしたため、参加者の平均点には「もともとの成績の高さ」が上乗せされていたのが原因です。傾向スコアマッチングは、共変量の近い個体同士を比べることでこの上乗せ分を打ち消し、真の効果に近い推定値を返しました。

ただし、マッチングをすれば選択バイアスが常に消えるわけではありません。今回の分析で考慮した共変量は年齢と過去の成績の2つで、処置への割り当てもこの2つだけで決まるようにデータを作ってあります。実際のデータでは、学習意欲のような未観測の交絡因子が割り当てに効いている可能性があり、その場合は観察された共変量をいくらマッチングしてもバイアスが残ります。傾向スコアマッチングには「共通サポート」の仮定、つまり処置群と対照群で傾向スコアの範囲が重なっているという条件も必要です。マッチングの手法自体にも様々な種類があり、どれを選ぶかで結果が変わる場合があります。

まとめ

因果推論は、相関の背後にある因果関係に迫るための枠組みです。今回の実験では、選択バイアスのあるデータで単純な平均の比較が効果を約2割過大評価し、傾向スコアマッチングで真値の近くまで戻せることを確かめました。ただしこれは、割り当てを決める変数がすべて観測されているという、合成データだからこそ成り立つ前提の上での結果です。実データに適用するときは、処置への割り当てが何で決まっているかを先に洗い出し、観測できていない要因が効いていそうであれば、操作変数法や差分の差分法など仮定の異なる手法を候補に入れて検討することになります。

コード

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

# 再現性のために乱数シードを設定
np.random.seed(42)

def generate_synthetic_data(n_samples=1000):
    """教育プログラムの効果を検証するための合成データを生成"""
    # 共変量（年齢、過去の成績）
    age = np.random.normal(20, 2, n_samples)
    previous_score = np.random.normal(70, 10, n_samples)

    # 処置への割り当て（プログラムへの参加）
    # 年齢と過去の成績が高いほど参加しやすい（選択バイアス）
    propensity = 1 / (1 + np.exp(-(age - 20) * 0.1 - (previous_score - 70) * 0.03))
    treatment = np.random.binomial(1, propensity)

    # 潜在的結果の生成
    # プログラムに参加しない場合の結果
    y0 = 60 + 0.5 * previous_score + np.random.normal(0, 5, n_samples)

    # プログラムに参加した場合の結果
    # 過去の成績が低い人ほどプログラムの恩恵を受けやすい
    treatment_effect = 15 - 0.1 * previous_score
    y1 = y0 + treatment_effect + np.random.normal(0, 2, n_samples)

    # 観察される結果
    y_obs = treatment * y1 + (1 - treatment) * y0

    # データフレームの作成
    df = pd.DataFrame({
        'age': age,
        'previous_score': previous_score,
        'treatment': treatment,
        'outcome': y_obs,
        'true_effect': y1 - y0
    })

    return df

def naive_ate(df):
    """単純な平均処置効果の推定"""
    treated = df[df['treatment'] == 1]['outcome'].mean()
    control = df[df['treatment'] == 0]['outcome'].mean()
    return treated - control

def plot_treatment_effects(df):
    """処置効果の可視化"""
    plt.figure(figsize=(10, 6))

    # 処置群と対照群の成績分布
    sns.kdeplot(data=df[df['treatment'] == 1], x='outcome', 
                label='処置群（プログラム参加）', alpha=0.5)
    sns.kdeplot(data=df[df['treatment'] == 0], x='outcome', 
                label='対照群（非参加）', alpha=0.5)

    plt.title('プログラム参加の有無による成績分布の比較')
    plt.xlabel('最終成績')
    plt.ylabel('密度')
    plt.legend()
    plt.savefig('treatment_effects.png')
    plt.close()

def plot_selection_bias(df):
    """選択バイアスの可視化"""
    plt.figure(figsize=(10, 6))

    sns.scatterplot(data=df, x='previous_score', y='treatment', alpha=0.5)
    plt.title('過去の成績とプログラム参加の関係（選択バイアス）')
    plt.xlabel('過去の成績')
    plt.ylabel('プログラム参加（1=参加, 0=非参加）')
    plt.savefig('selection_bias.png')
    plt.close()

def estimate_propensity_scores(df):
    """傾向スコア（処置を受ける確率）を推定"""
    # 共変量を使って処置の有無を予測するモデルを構築
    X = df[['age', 'previous_score']]
    treatment = df['treatment']

    # ロジスティック回帰で確率を推定
    model = LogisticRegression()
    model.fit(X, treatment)

    # 傾向スコアを計算
    propensity_scores = model.predict_proba(X)[:, 1]
    return propensity_scores

def match_and_estimate_ate(df):
    """傾向スコアマッチングによる処置効果の推定"""
    # 傾向スコアの推定
    df['propensity_score'] = estimate_propensity_scores(df)

    # 処置群と対照群に分割
    treated = df[df['treatment'] == 1]
    control = df[df['treatment'] == 0]

    # マッチング（最近傍法）
    matched_pairs = []
    for _, treated_unit in treated.iterrows():
        # 最も傾向スコアが近い対照群の個体を見つける
        distances = abs(control['propensity_score'] - treated_unit['propensity_score'])
        best_match_idx = distances.idxmin()
        matched_pairs.append({
            'treated': treated_unit['outcome'],
            'control': control.loc[best_match_idx, 'outcome'],
            'treated_ps': treated_unit['propensity_score'],
            'control_ps': control.loc[best_match_idx, 'propensity_score']
        })

    # マッチングされたペアの平均差を計算
    matched_effect = np.mean([pair['treated'] - pair['control'] 
                            for pair in matched_pairs])
    return matched_effect, matched_pairs

def plot_matching_result(pairs):
    """マッチング結果の可視化"""
    plt.figure(figsize=(10, 6))

    # マッチングされたペアの傾向スコアをプロット
    treated_ps = [p['treated_ps'] for p in pairs]
    control_ps = [p['control_ps'] for p in pairs]

    plt.scatter(treated_ps, control_ps, alpha=0.5)
    plt.plot([0, 1], [0, 1], 'r--')  # 45度線

    plt.title('マッチングの質の確認')
    plt.xlabel('処置群の傾向スコア')
    plt.ylabel('対照群の傾向スコア')
    plt.savefig('matching_quality.png')
    plt.close()

# メインの実行部分
if __name__ == "__main__":
    # データの生成
    df = generate_synthetic_data()

    # 結果の表示
    print("ナイーブな平均処置効果推定値:", naive_ate(df))
    print("真の平均処置効果:", df['true_effect'].mean())

    # 傾向スコアマッチングによる推定
    matched_ate, pairs = match_and_estimate_ate(df)
    print("傾向スコアマッチングによる処置効果推定値:", matched_ate)

    # 可視化
    plot_treatment_effects(df)
    plot_selection_bias(df)
    plot_matching_result(pairs)