因果推論入門：基本概念とPythonによる実践

はじめに

「相関は因果関係を意味しない」という言葉はよく知られていますが、実際に因果関係をどのように推論すればよいのでしょうか。本記事では、因果推論の基本概念をわかりやすく解説し、Pythonを用いた実践的な例を通じて、その理解を深めることを目指します。

この記事を書いたひと

デジタルリアクタ合同会社 代表
望月 優輝 / Yuki Mochizuki   

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

対象読者:

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

記事のポイント:

• 因果推論の基本的な考え方と重要な概念を理解する
• 選択バイアスがなぜ問題となるのかを理解する
• 傾向スコアマッチングの基本的な仕組みを理解する
• Pythonを使って実際に因果推論を体験する

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

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

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

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

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

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

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

しかし、実際に観察できるのは、Y_i(1) と Y_i(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)
  • 処置群と対照群をランダムに割り当てることで、選択バイアスを排除する
  • 最も信頼性の高い方法だが、実施が困難な場合も多い
• 観察データを用いた方法
  • 傾向スコアマッチング: 処置群と対照群で、共変量の分布が似ている個体をマッチングする
  • 操作変数法: 処置に影響を与えるが結果には直接影響しない変数（操作変数）を利用する
  • 差分の差分法: 処置前後の変化の差を、処置群と対照群で比較する

本記事では、これらのうち、傾向スコアマッチング に焦点を当てて解説します。

Pythonによる実践例

ここでは、教育プログラムの効果を検証する例を通じて、因果推論をPythonで実践してみましょう。

Step 1: データの生成

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

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
• 選択バイアス
  • 年齢が高いほどプログラムに参加しやすい
  • 過去の成績が良いほどプログラムに参加しやすい
• 異質な処置効果
  • 過去の成績が低い生徒ほど、プログラムの効果が大きい

Step 2: 選択バイアスの確認

生成したデータを用いて、選択バイアスの存在を視覚的に確認してみましょう。

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()

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

Step 3: ナイーブな効果推定

まず、単純にプログラム参加者と非参加者の平均結果を比較してみましょう。

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

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

Step 4: 傾向スコアマッチングによる推定

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

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度線に近いほど、良いマッチングが実現できていることを意味します。

Step 5: 結果の比較

結果を比較すると、ナイーブな推定では選択バイアスの影響で効果を過大評価していることがわかります。一方、傾向スコアマッチングを用いることで、真のATEに近い推定値を得ることができています。

結果の解釈

今回の分析から、以下のことがわかりました。

1. 選択バイアスの影響: 単純な比較では、教育プログラムの効果を約1.65ポイント過大評価していた。これは、成績の良い生徒がプログラムに参加しやすいという選択バイアスが原因である。
2. 傾向スコアマッチングの有効性: 傾向スコアマッチングを用いることで、選択バイアスを軽減し、真の効果に近い推定値を得ることができた。
3. 注意点と限界:
   • 今回の分析では、観察された共変量（年齢と過去の成績）のみを考慮している。もし、未観測の交絡因子（例えば、学習意欲）が存在する場合、推定結果にバイアスが残る可能性がある。
   • 傾向スコアマッチングは、「共通サポート」の仮定を満たす必要がある。つまり、処置群と対照群で、傾向スコアが重なっている必要がある。
   • ※ ただし、マッチングの手法には様々な種類があり、結果が異なる場合がある。

まとめ

因果推論は、単なる相関関係ではなく、真の因果関係を明らかにするための強力なツールです。しかし、選択バイアスなどの問題に対処するためには、適切な統計的手法を選択し、適用することが重要です。今回の記事では、傾向スコアマッチングという一つの手法を紹介しましたが、他にも様々な手法が存在します。関心があれば、是非調べてみることをお勧めします。

コード

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) 

関連投稿:

Graph Neural Networksを用いた異常検知入門 確率的最適化アルゴリズムの実装と比較：遺伝的アルゴリズムと勾配降下法 次元削減手法の比較検証：定量的・定性的評価によるPCA、t-SNE、UMAPの特徴分析 反実仮想サンプルの生成と使い方