因果推論2 傾向スコア:母集団に偏りがあるときの施策の効果を測る
統計的因果推論の勉強
- 理論の詳しい話は下記の本を見てください。
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
- 作者: 星野崇宏
- 出版社/メーカー: 岩波書店
- 発売日: 2009/07/29
- メディア: 単行本
- 購入: 29人 クリック: 285回
- この商品を含むブログ (25件) を見る
問題
- AとBの2つの学校があり、Aには何か数学の点数が上がるような施策をし、Bには何もしない。
- AとBには、親の学歴の分布が異なる。
- 2つのグループの設定が異なるので、単純には比較できない。
- 親の学歴とかを除いて、施策が効果があったのかを比較したい。
単純に無作為サンプリングできる場合。
- 上記のように異なる特徴をもつ集団があったとしても、AとBの生徒をグチャ混ぜにして、無細工に選んで施策を行った場合、単純に比較できる。
- この場合、施策以外の特徴量を揃えることができるから。
- 例えば、数学の点数は施策に関係なく、親の学歴に依存している場合を考える。
- 赤が施策あり。青が施策なし。
- 縦軸が数学の点数、横軸が親の学歴。
- 赤の平均が49.56。青の平均が49.67。
- 期待通り、今回の施策には効果がないことがわかる。
n = 1000 parent_school_career = np.random.normal(3, 1, n) treated = np.random.binomial(1, 0.5, n) math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3) data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated}) # plot all sampling data plt.scatter(data['parent'], data['math'], c=data['treated']) plt.show() idx1 = data['treated'] == 1 idx0 = data['treated'] == 0 math0 = data[idx0]['math'] math1 = data[idx1]['math'] print('z=0:', np.average(math0), np.std(math0)) print('z=1:', np.average(math1), np.std(math1))
母集団に偏りがある場合
- Aの親の方がBに比べ、学歴が高いとする。
- Aに施策を行った。
- Aが赤。Bが青。
- 赤の平均が54.74。青の平均が45.49。
- 施策には効果がないのに、平均を単純に比較すると有意に差が出てしまう。
n = 1000 treated = np.random.binomial(1, 0.5, n) parent_school_career = np.random.normal(2.5, 1, n) + treated math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3) data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated}) # plot all sampling data plt.scatter(data['parent'], data['math'], c=data['treated']) plt.ylabel('math score') plt.xlabel("parent's school career") plt.show() idx1 = data['treated'] == 1 idx0 = data['treated'] == 0 math0 = data[idx0]['math'] math1 = data[idx1]['math'] print('z=0:', np.average(math0), np.std(math0)) print('z=1:', np.average(math1), np.std(math1))
傾向スコアを用いて共変量(親の学歴)の影響をなくす。IPW推定を行う。
- 赤の平均が49.81。青の平均が49.52。
- 共変量の影響がなくなり、期待通り平均の差がなくなった。
def calculate_score(x, alpha): a = np.dot(x, np.array([alpha]).T).flatten() s = 1 / (1 + np.exp(-a)) return s def loglikelihood(x, z, alpha): s = calculate_score(x, alpha) p = z * np.log(s) + (1 - z) * np.log(1 - s) return np.sum(p) def non_random_sampling_ipw(): n = 1000 treated = np.random.binomial(1, 0.5, n) parent_school_career = np.random.normal(2.5, 1, n) + treated math_score = np.random.normal(50, 10, n) + 10 * (parent_school_career - 3) data = pd.DataFrame({'math': math_score, 'parent': parent_school_career, 'treated': treated}) data['const'] = np.ones(n) idx1 = data['treated'] == 1 idx0 = data['treated'] == 0 cons = ( {'type': 'ineq', 'fun': lambda theta: 10 - theta}, {'type': 'ineq', 'fun': lambda theta: theta + 10}, ) alpha = np.random.uniform(0, 1, 2) res = minimize(lambda a: -loglikelihood(data[['parent', 'const']], data['treated'].as_matrix().flatten(), a), alpha, constraints=cons, method="SLSQP") score = calculate_score(data[['parent', 'const']], res.x) w0 = (1 - data['treated']) / (1 - score) w1 = data['treated'] / score n0 = np.sum(w0) n1 = np.sum(w1) y0 = np.sum(w0 * data['math']) / n0 y1 = np.sum(w1 * data['math']) / n1 print(y0) print(y1)
今後の課題
- IPW推定を使ったときの分散はどう求める? どうやって検定を行う?
- 傾向スコアのモデルが適切でないと最終結果がよくない。今回 でなく としていたときはかなり悪かった。
- 潜在的な共変量があるときはどうすれば?