因果推論3 操作変数法

統計的因果推論の勉強

問題設定

  • 目的
    • バウチャー制度の効果を評価したい。
  • 背景
    • くじ引きで当たると、「私立高校へ進学したら授業料が半額になる権利」がもらえる。
    • くじ引きは希望者全員が引ける。当たりは40%とする。
    • 私立高校、公立高校ともに入学試験はないとする。
    • しかし、私立高校の方が質の高い授業を行っているとする。
    • 卒業時に、私立・公立高校で統一テストを実施している。
  • 状況整理
    • 全部で4つのパターンがある
      • くじ引き当たり、私立へ進学する
      • くじ引き当たり、公立へ進学する
      • くじハズれ、私立へ進学する
      • くじハズれ、公立へ進学する
    • バウチャー制度の評価ポイント
      • くじ引きが当たれば、私立へ進学し、外れたら公立へ進学する生徒の卒業時の成績。

観察データ作成(ケース1)

  • 全体の10%当たりは結果に関係なく必ず公立へ進学する。
  • 生徒の20%は裕福で、上記の10%に入った生徒以外はくじ引きの結果に関係なく私立へ進学する。
  • 生徒の成績の平均は、y = 600 + 80 * (私立) + 150 * (裕福)
    • 観測できない影響を大きくするために、150 * (裕福)とする。
  • 標準偏差は50
  • 生徒数は1000
  • 生成したデータはこんな感じ。
  • 生徒が裕福かどうかは不明とする。
 lottery  private       score
0        0        0  605.830138
1        0        0  595.145346
2        0        0  660.159643
3        1        1  676.014558
4        0        0  612.031560
5        0        0  592.204991
6        0        1  704.909145
7        0        0  618.559375
8        0        0  626.251403
9        0        0  536.343088
       lottery      private        score
count  1000.00000  1000.000000  1000.000000
mean      0.41700     0.485000   649.887460
std       0.49331     0.500025    72.562223
min       0.00000     0.000000   476.666647
25%       0.00000     0.000000   595.276742
50%       0.00000     0.000000   644.046202
75%       1.00000     1.000000   704.860884
max       1.00000     1.000000   856.068556

簡単な分析(ケース1)

  • 単純に私立と公立で平均をとってみる。
    • データ生成時のモデル上、差は80程度のはず。
    • 他の要素の影響が含まれている。
106.589641791

操作変数法

  • 割当
    • くじ引きが当たる: z=1
    • くじ引きが外れる: z=0
  • 操作変数
    • z=1のとき、私立へ進学する: d_1=1
    • z=1のとき、公立へ進学する: d_1=0
    • z=0のとき、私立へ進学する: d_0=1
    • z=0のとき、公立へ進学する: d_0=0
    • 結果に影響しないので、d_1 \geq d_0仮定する(単調性)
  • 私立へ進学するかどうか
    • d=z d_1 + (1 - z) d_0
  • テストの成績
    • y=d y_1 + (1 - d) y_2
  • 独立性
    • yz|dは独立
    • d_zzは独立
  • 推定したいもの
    • くじ引きが当たれば私立へ進学するが、外れれば公立へ進学する生徒の成績の差
      • E[y_1 - y_0|d_1=1, d_0=0]
  • 計算する上での問題
    • 条件を満たす生徒が誰なのかがわからない。(当たっても私立へ進学しない生徒の存在など。)
    • 当然、y_1y_0のどちらかしか観測できていない。
  • どう計算するか(詳細な計算は本を見てください)。
    • 準備
      • E[y|z=1-E[y|z=0]=E[y_1-y_0|d_1-d_0=0]p(d_1-d_0=1)]
      • p(d_1=1)=p(d_1=1, d_0=1) + p(d_1=1, d_0=0)
      • p(d_0=1)=p(d_0=1,d_1=1)
      • E(d|z=1)=p(d_1=1)
      • E(d|z=0)=p(d_0=1)
    • 以上より
      • E(y_1 - y_0|d_1=1, d_0=0)=\frac{E(y|z=1)-E(y|z=0)}{p(d_1-d_0=1)}=\frac{E(y|z=1)-E(y|z=0)}{E(d|z=1)-E(d|z=0)}
    • 最終的にはデータから計算できるものばかりで表される!!
  • 結果は期待通り80に近い。
81.2776422538

今後の課題

  • 結果の検定をどうするのか?
  • 入学試験がある場合はどうすればよいか?(選択バイアスが出るよね?)

使ったプログラムコード

import numpy as np
import pandas as pd


def sampling():
    n = 1000
    score_base = 600
    score_private = 80
    score_rich = 150
    score_sigma = 50
    rich_ratio = 0.2
    public_ratio = 0.1
    winning_ticket_ratio = 0.4

    data = pd.DataFrame({
        'rich': np.random.binomial(1, rich_ratio, n),
        'lottery': np.random.binomial(1, winning_ticket_ratio, n)
    })
    data['private'] = data.apply(lambda d: 1 if d['rich'] == 1 or d['lottery'] == 1 else 0, axis=1)
    data['private'] = data['private'] * np.random.binomial(1, 1 - public_ratio, n)
    data['score'] = score_base + score_private * data['private'] + score_rich * data['rich'] \
        + np.random.normal(0, score_sigma, n)

    return data.drop(['rich'], axis=1)


def simple_analysis(data):
    d0 = data[data['private'] == 0]
    d1 = data[data.apply(lambda d: d['private'] == 1 and d['lottery'] == 1, axis=1)]
    print(np.average(d1['score']) - np.average(d0['score']))


def instrumental_variable_method(data):
    data0 = data[data['lottery'] == 0]
    data1 = data[data['lottery'] == 1]
    y0 = np.average(data0['score'])
    y1 = np.average(data1['score'])
    d0 = np.average(data0['private'])
    d1 = np.average(data1['private'])
    print((y1 - y0) / (d1 - d0))


if __name__ == '__main__':
    np.random.seed(123)
    data = sampling()
    print(data.head(10))
    print(data.describe())
    simple_analysis(data)
    instrumental_variable_method(data)
    print('hello')```