因果推論3 操作変数法
統計的因果推論の勉強
- 理論の詳しい話は下記の本を見てください。
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
- 作者: 星野崇宏
- 出版社/メーカー: 岩波書店
- 発売日: 2009/07/29
- メディア: 単行本
- 購入: 29人 クリック: 285回
- この商品を含むブログ (25件) を見る
問題設定
- 目的
- バウチャー制度の効果を評価したい。
- 背景
- くじ引きで当たると、「私立高校へ進学したら授業料が半額になる権利」がもらえる。
- くじ引きは希望者全員が引ける。当たりは40%とする。
- 私立高校、公立高校ともに入学試験はないとする。
- しかし、私立高校の方が質の高い授業を行っているとする。
- 卒業時に、私立・公立高校で統一テストを実施している。
- 状況整理
- 全部で4つのパターンがある
- くじ引き当たり、私立へ進学する
- くじ引き当たり、公立へ進学する
- くじハズれ、私立へ進学する
- くじハズれ、公立へ進学する
- バウチャー制度の評価ポイント
- くじ引きが当たれば、私立へ進学し、外れたら公立へ進学する生徒の卒業時の成績。
- 全部で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
操作変数法
- 割当
- くじ引きが当たる:
- くじ引きが外れる:
- 操作変数
- のとき、私立へ進学する:
- のとき、公立へ進学する:
- のとき、私立へ進学する:
- のとき、公立へ進学する:
- 結果に影響しないので、仮定する(単調性)
- 私立へ進学するかどうか
- テストの成績
- 独立性
- とは独立
- とは独立
- 推定したいもの
- くじ引きが当たれば私立へ進学するが、外れれば公立へ進学する生徒の成績の差
- ]
- くじ引きが当たれば私立へ進学するが、外れれば公立へ進学する生徒の成績の差
- 計算する上での問題
- 条件を満たす生徒が誰なのかがわからない。(当たっても私立へ進学しない生徒の存在など。)
- 当然、かのどちらかしか観測できていない。
- どう計算するか(詳細な計算は本を見てください)。
- 準備
- -E[y|z=0]=E[y_1-y_0|d_1-d_0=0]p(d_1-d_0=1)]
- 以上より
- 最終的にはデータから計算できるものばかりで表される!!
- 準備
- 結果は期待通り80に近い。
81.2776422538
今後の課題
- 結果の検定をどうするのか?
- 入学試験がある場合はどうすればよいか?(選択バイアスが出るよね?)
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
- 作者: 星野崇宏
- 出版社/メーカー: 岩波書店
- 発売日: 2009/07/29
- メディア: 単行本
- 購入: 29人 クリック: 285回
- この商品を含むブログ (25件) を見る
使ったプログラムコード
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')```