因果推論2 傾向スコア:母集団に偏りがあるときの施策の効果を測る

統計的因果推論の勉強

問題

  • AとBの2つの学校があり、Aには何か数学の点数が上がるような施策をし、Bには何もしない。
  • AとBには、親の学歴の分布が異なる。
  • 2つのグループの設定が異なるので、単純には比較できない。
  • 親の学歴とかを除いて、施策が効果があったのかを比較したい。

単純に無作為サンプリングできる場合。

  • 上記のように異なる特徴をもつ集団があったとしても、AとBの生徒をグチャ混ぜにして、無細工に選んで施策を行った場合、単純に比較できる。
  • この場合、施策以外の特徴量を揃えることができるから。
  • 例えば、数学の点数は施策に関係なく、親の学歴に依存している場合を考える。 f:id:nsb248:20170207200206p:plain
  • 赤が施策あり。青が施策なし。
  • 縦軸が数学の点数、横軸が親の学歴。
  • 赤の平均が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に施策を行った。 f:id:nsb248:20170207200905p:plain
  • 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推定を使ったときの分散はどう求める? どうやって検定を行う?
  • 傾向スコアのモデルが適切でないと最終結果がよくない。今回  y = \alpha x + \betaでなく y = \alpha x としていたときはかなり悪かった。
  • 潜在的な共変量があるときはどうすれば?