因果推論1 単純な選択バイアス問題

統計的因果推論の勉強

問題

  • 大学の入試の結果y1と入学後の成績y2の関係を調べたい。

設定

  • 真のモデル
    • 入試の平均50、分散40
    • 入学後の成績の平均50、分散40
    • 2つのテストの相関0.8
    • 正規分布に従うと仮定。
    • [0, 100]となるように切り捨てる
  • 選択バイアス
    • 入試の合格点を60とする。
    • 当然、入試の60点以下の人の入学後の成績は欠測となる。
    •  mを欠測のインディケータとする。この場合、合格のとき m=1, 不合格のとき m=0とする。

f:id:nsb248:20170206202538p:plain - 赤い影になっている部分が欠測データ。

回帰分析モデルを仮定する。

  • モデル
    •  y_2 = \theta_1 + \theta_2 y_1 + \epsilon
    •  \epsilon \sim N(0, \sigma ^2)
  • 最尤推定
    •  p(y_2, m | y_1, \theta_1, \theta_2, \sigma)=p(y_2|y_1, \theta_1, \theta_2, \sigma)p(m|y_1)
    • つまり、 (y_1, y_2)を用いた最小二乗法を行えばよい。
  • このときの y_1, y_2の相関は、
    •  corr(y_1, y_2) = \frac{\theta_2 \times V(y_1)}{\sqrt{V(y_1)}\times \sqrt{\theta_2 ^2 \times V(y_1) + V(\epsilon)}}
parameters: [  0.63938213  16.73959436]
correlation: 0.820560498514

2変量正規分布に従うと仮定する。

  • likelihood
    •  p(y_1, y_2, m|\theta) = p(y_1, y_2|\theta) p(m | y_1)
  • observed data likelihood
    • [tex: p(y_1, y_2|\theta) = \prod{i, m_i=1}{p(y(i1), y(i2) | \theta)} \prod{i, m_i=0} p(y(i1) | \theta)]
  • つまり、observed data likelihoodを最大化すれば良い。
mu1: 49.73650529
mu2: 48.54009208
sigma1: 10.18582171
sigma2: 7.93686404
rho: 0.820564 

Python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
from scipy.stats import norm, multivariate_normal
from scipy.optimize import minimize


def negative_log_likelihood(y1, y2, m, m1, m2, s1, s2, rho):
    mean = np.array([m1, m2])
    s12 = rho * s1 * s2
    cov = np.array([[s1 * s1, s12], [s12, s2 * s2]])
    v = 0.0
    for yi1, yi2, mi in zip(y1, y2, m):
        if mi == 1:
            v += multivariate_normal.logpdf(np.array([yi1, yi2]), mean=mean, cov=cov)
        else:
            v += norm.logpdf(yi1, m1, s1)
    return -v


if __name__ == '__main__':
    n = 1000
    u1 = 50
    u2 = 50
    v1 = 100
    v2 = 50
    rho = 0.8
    v12 = rho * math.sqrt(v1 * v2)
    mean = np.array([u1, u2])
    cov = np.array([[v1, v12], [v12, v2]])
    y = np.random.multivariate_normal(mean, cov, n)

    data = pd.DataFrame({'y1': y[:, 0], 'y2': y[:, 1]})
    data['m'] = data.apply(lambda x: 1 if x['y1'] > 60 else 0, axis=1)
    data['y2'] = data.apply(lambda x: None if x['m'] == 0 else x['y2'], axis=1)
    pass_idx = data['m'] == 1

    # using regression model
    coeff = np.polyfit(data[pass_idx]['y1'], data[pass_idx]['y2'], 1)
    c0 = coeff[0]
    var_y1 = np.std(data['y1']) ** 2
    var_epsilon = np.std(data['y2'] - np.poly1d(coeff)(data['y1'])) ** 2
    corr = c0 * var_y1 / (math.sqrt(var_y1) * math.sqrt(c0 ** 2 * var_y1 + var_epsilon))
    print('parameters:', coeff)
    print('correlation:', corr)

    # using multivariate normal model.
    cons = (
        {'type': 'ineq', 'fun': lambda theta: 0.9999 - theta[4]},
        {'type': 'ineq', 'fun': lambda theta: theta[4] + 0.9999},
        {'type': 'ineq', 'fun': lambda theta: theta[0:4] - 1},
        {'type': 'ineq', 'fun': lambda theta: 100 - theta[0:4]},
    )
    theta0 = np.array([50, 60, 10, 30, 0.1])
    theta = minimize(lambda theta: negative_log_likelihood(data['y1'], data['y2'], data['m'],
                                                           theta[0], theta[1], theta[2], theta[3], theta[4]),
                     theta0,
                     constraints=cons,
                     method="SLSQP")
    print(theta)

    axis_x = np.linspace(0, 100, 101)
    plt.scatter(y[:, 0], y[:, 1])
    plt.fill_betweenx((-10, 110), 60, facecolor='r', alpha=0.2)
    plt.xlim(0, 100)
    plt.ylim(0, 100)
    plt.plot(axis_x, np.poly1d(coeff)(axis_x))
    plt.xlabel('y1')
    plt.ylabel('y2')
    plt.show()