因果推論に入門( 処置効果が線形なデータで )

因果分析に興味があるということでこちらの書籍を勉強中。

book.mynavi.jp

こちらの書籍で紹介されている因果推論の手法を整理していきます。今回は処置効果が線形な場合に、適用可能な手法の整理となります。

因果推論とは
データの説明
データの可視化
回帰分析
傾向スコアによるIPTW(Inverse Probability of Treatment Weighting : 逆確率重み付け法 )
DR(Doubly Robust法)
各手法の結果まとめ
おわりに

因果推論とは

処置の効果測定をする手法のことです。まず因果関係というのは、ある変数を変化させたときに着目している変数が変化する関係にあります。(ある変数が大きいときに、着目している変数が大きかったり小さかったりしているのは相関関係です。) そしてこのときにどれくらい変化するのか=どれくらい効果があるのかを推定するのが、因果推論になります。

実際に因果推論する場合は因果探索やドメイン知識を活用し、因果ダイアグラムなどを作成して変数間の因果関係を示す必要があります。ただ今回は因果推論の手法の整理ということで以下のような擬似データを作成し、分析していきます。

データの説明

以下の構造を持つデータで因果推論を行います。
テレビCMが商品の購入量に影響するという場面を想定しています。
"商品の購入量: Y"が"テレビCMを見た: Z"という処置の影響を受けるとし、このときどれくらいの効果があるかを推定します。また商品をどのくらい購入するかとテレビCMを見るかどうかは年齢と性別にも影響を受けるとします。

このとき、
購入量を結果変数(処置の結果となる変数)、
CMを見たかどうかを原因変数(処置の原因となる変数)、
年齢、性別を共変量(因果関係に影響を及ぼす変数)
と呼びます。

データは以下で生成します。 処置効果が線形という点ですが、"商品の購入量: Y"を"テレビCMを見た: Z"の一次式(今回は係数を10として)で生成している点になります。

import random
import numpy as np

# 乱数のシードを固定
random.seed(9999)
np.random.seed(9999)

# シグモイド関数をimport
from scipy.special import expit

# 標準正規分布の生成用
from numpy.random import *

# グラフの描画用
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns

# その他
import pandas as pd

# 乱数のシードを固定
random.seed(9999)
np.random.seed(9999)

# データ数
num_data = 500

# 原因変数以外の説明変数( 共変量 )
# 年齢
x_1 = randint(15, 76, num_data) 

# 性別
x_2 = randint(0, 2, num_data)

# 原因変数
# シグモイド関数に入れる部分。テレビCMを見る可能性を確率変数として生成。(大きいほど見る確率は高くなる)
e_z = randn(num_data) # ノイズ(バイアス項)の生成
z_base = x_1 + (1-x_2)*10 - 40 + 5*e_z

# テレビCMを見る確率を計算 = シグモイド関数を計算
z_prob = expit(0.1*z_base)

# テレビCMを見たかどうかの変数（0は見ていない、1は見た）
Z = np.array([])
for i in range(num_data):
  Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0] # 0または1になる確率に対して、ランダムに選択している
  Z = np.append(Z, Z_i)

# 目的変数
# 各効果も指定する。年齢が高いほど減少、男性の方が増加、テレビCMを見ていると増加。
e_y = randn(num_data) # ノイズ(バイアス項)の生成
Y = -x_1 + 30*x_2 + 10*Z + 80 + 10*e_y

# テーブルを作成
df_A = pd.DataFrame({'age': x_1,
                   'gender': x_2,
                   'CM': Z,
                   'pur_vol': Y,
                   })

df_A.head(5)

[output]
# データの先頭行
age  gender  CM   pur_vol
37     0     1.0  56.691968
64     1     1.0  66.464076
62     0     1.0  14.130339
60     1     1.0  73.538865
69     1     1.0  40.870089

データの可視化

データを可視化します。

まずseabornのpairplotで簡単に見てみます。
購入量に対して、年齢、性別は傾向があります。しかしCMの効果については、CMを見た(1)がCMを見ていない(０)に対して、上にも下にもプロットされており、効果があるのかどうかわかりません。

sns.pairplot(df_B)

またCMを見たかどうかで購入量のヒストグラムを比較した図が以下です。
CMを見ていないほうが購入量が右側に寄っており、可視化だけではCMが逆効果になっていそうと判断してしまいそうです。
また各カラムの平均を見ると年齢と性別の平均に乖離があります。

# CMを見たかどうかで購入量のヒストグラムを比較
df_A[df_A["CM"]==1]["pur_vol"].hist(color="r", alpha=0.3, label="CM=1")
df_A[df_A["CM"]==0]["pur_vol"].hist(color="b", alpha=0.3, label="CM=0")
plt.legend()
plt.show()

print(df_A[df_A["CM"]==1].mean())
print("======")
print(df_A[df_A["CM"]==0].mean())

[output]
==CMを見た集団の各カラムの平均==
age        51.702454
gender      0.463190
CM          1.000000
pur_vol    52.210321
dtype: float64
==CMを見ていない集団の各カラムの平均==
age        31.781609
gender      0.557471
CM          0.000000
pur_vol    64.855318

回帰分析

まずは回帰分析です。今回は結果変数が連続値なので、線形回帰を用います。
回帰モデルによる因果推論は比較的分かりやすいです。
目的変数として結果変数を、説明変数として考量すべき変数(原因変数、共変量)を用いたモデルを構築し、そのときの係数が効果の量になります。以下で実行します。

# 説明変数
X = df_A[["age", "gender", "CM"]]

# 被説明変数（目的変数）
y = df_A["pur_vol"]

# 学習(回帰の実施)
model = LinearRegression()
model.fit(X, y)

# 回帰した結果の係数を出力
print("係数：", model.coef_)

[output]
係数： [-1.02228972 29.62534334 10.51298688]

各係数は、年齢: -1.0、性別: 29.6、CM: 10.5となりました。
データ生成時が、年齢: -1、性別: 30、CM: 10なのでそれなりに良い結果かと。
また各変数の影響が一括で出せるのもお手軽で良いかなと思いました。

傾向スコアによるIPTW(Inverse Probability of Treatment Weighting : 逆確率重み付け法 )

次は傾向スコアを用いたIPTWになります。こちらは傾向スコア(ある属性情報に応じた処置を受ける確率)を用いて、ATE(Average Treatment effect, 平均処置効果)と呼ばれる推定した処置効果を算出する手法です。
手順としては、
1. 原因変数を共変量でモデリング
2. 共変量の元での傾向スコアを算出
3. ATEを調整化公式で変形した数式に傾向スコアを代入し推定したえ処置効果を算出
となります傾向スコアは、今回はロジスティック回帰で求めます。
また3について、調整化公式とはdoオペレータで表した数式を、doオペレータ無しで表現する方法であり、
こちらと条件付き確率から同時確率へ変換で、ATEは以下のように変形できます。
【調整化公式】

$P(Y=y | do(Z=z))$

$= \sum_{x=1}^n P(Y=y|Z=z, X=x) P(X=x)$

$= \sum_{x=1}^n \frac{P(Y=y, Z=z, X-x)}{P(Z=z, X=x)}$

【ATEの式変形】
$ATE$

$= E(Y1|do(Z=1)) - E(Y0|do(Z=0))$

$= \frac{1}{N} \sum_{x=1}^n\frac{yi}{P(Z=1|X=xi)}Zi - \frac{1}{N} \sum_{x=1}^n\frac{yi}{P(Z=0|X=xi)}(1-Zi)$

IPTWを以下で実行します。

# 傾向スコアを求めるために、原因変数を目的変数としてモデリング
# 説明変数
X = df_A[["age", "gender"]]

# 目的変数（原因変数）
Z = df_A["CM"]

# 学習(回帰の実施)
model = LogisticRegression()
model.fit(X,Z)

# 係数を出力
print("--係数--")
print("係数beta：", model.coef_)
print("係数alpha：", model.intercept_)

# 各サンプルの傾向スコア
Z_pred = model.predict_proba(X)
print("--傾向スコアの推論値--")
print(Z_pred[0:5])  # 推定結果
print("--正解--")
print(Z[0:5])  # 正解

# 処置の効果を推定
Y = df_A["pur_vol"]
ATE_i = Y / Z_pred[:, 1] * Z - Y / Z_pred[:, 0] * (1-Z) # 式変形後
ATE = 1 / len(Y) * ATE_i.sum() # 平均処置効果 ATEを計算
print("推定したATE: ", ATE)

[output]
--係数--
係数beta： [[ 0.09254961 -0.67589802]]
係数alpha： [-2.84351622]
--傾向スコアの推論値--
[[0.35874405 0.64125595]
 [0.08288635 0.91711365]
 [0.05242315 0.94757685]
 [0.11572386 0.88427614]
 [0.05383419 0.94616581]]
--正解--
0    1.0
1    1.0
2    1.0
3    1.0
4    1.0

推定したATE 13.269226060428515

まずシグモイド関数の係数について正解: 0.1, -1, -3に対して0.09、-0.68、-2.84という結果でした。
傾向スコアは出力した分に関しては、1に対してのスコアが高くそこそこ妥当な結果かと思います。処置効果の推定値は回帰分析よりは、正解から離れた結果となりました。　　

以下に傾向スコアとのヒストグラムを描画します。

# 傾向スコアのヒストグラムを確認
temp_df = pd.merge(pd.DataFrame(Z_pred).reset_index(), df_A["CM"].reset_index())　# サンプル内で処置に対応する傾向スコアのテーブル作成
pd.concat([temp_df[temp_df["CM"]==0][0], temp_df[temp_df["CM"]==1][1]], ignore_index=True).hist() # 傾向スコアのヒストグラムを描画

# 各サンプルのTE(ITE)のヒストグラム
pd.DataFrame(ATE_i).hist()

IPTW法は傾向スコアで割っているので、傾向スコアが0に近いと計算が不安定になります。
ヒストグラムを見ると小さい値が確認されるので、このあたりでそこそこの推定結果になったことが伺えます。実際に各サンプルに対する処置効果の推定値のヒストグラムは-600~1000と広範な値をとっていました。

DR(Doubly Robust法)

最後はIPTW法と回帰分析を組み合わせたDoubly Robust(DR法)になります。
こちらはIPTW法では考慮できていなかった各サンプルの反実仮想を、回帰分析などで推定し計算に加える手法です。

手順としては、
1. 反実仮想の推定値を計算: 目的変数として結果変数を、説明変数として原因変数と共変量を用いてモデル構築し、原因変数が各場合における(今回は0または1)推定値を出力する
2. 傾向スコアを求める
3. 下記の式に従ってサンプルごとの推定処置効果を計算し、ATEを求める

1については、購入量に対して回帰モデルを構築して、推論時にCMを見たかどうかについて0、1それぞれにおいて予測して反実仮想の値としてやります。
3については、1で求めた値を使ってIPTWで用いた前半の項と後半の項をそれぞれ加えてやります。
処置の推定効果をそれぞれ
前半の項を、

$\frac{yi}{P(Z=1|X=xi)}Zi → \frac{yi}{P(Z=1|X=xi)}Zi+(1-\frac{Zi}{P(Z=1|X=xi)})\hat{Y1}^i$

後半の項を、

$\frac{yi}{P(Z=0|X=xi)}(1-Zi) → \frac{yi}{P(Z=0|X=xi)}(1-Zi)+(1-\frac{(1-Zi)}{P(Z=0|X=xi)})\hat{Y0}^i$

として差を求めます。
DR法を以下で実行します。

# 反実仮想の推定値(購入量)を求める

# 説明変数
X = df_A[["age", "gender", "CM"]]

# 目的変数
y = df_A["pur_vol"]

# モデリング
model = LinearRegression()
model.fit(X, y)

# 介入Z=0(テレビCMを見ていない)の場合の購入量を予測
X_0 = X.copy()
X_0["CM"] = 0
Y_0 = model.predict(X_0)

# 介入Z=1(テレビCMを見た)の場合の購入量を予測
X_1 = X.copy()
X_1["CM"] = 1
Y_1 = model.predict(X_1)

# 傾向スコアを求める
# 説明変数
X = df_A[["age", "gender"]]

# 被説明変数（目的変数）
Z = df_A["CM"]

# 回帰の実施
model = LogisticRegression()
model.fit(X, Z)

# 傾向スコアを求める
Z_pre = model.predict_proba(X)

# ATEを求める
Y = df_A["pur_vol"]

ATE_1_i = Y/Z_pre[:, 1]*Z + (1-Z/Z_pre[:, 1])*Y_1
ATE_0_i = Y/Z_pre[:, 0]*(1-Z) + (1-(1-Z)/Z_pre[:, 0])*Y_0
pred_eff = ATE_1_i-ATE_0_i
ATE = 1/len(Y)*pred_eff.sum()
print("推定したATE", ATE)

[output]
推定したATE 11.578048448712533

推定した処置効果は、11.6とIPTWの結果(13.2)より10に近く良好となりましたが、
線形回帰の方が良好な結果となりましたこちらも各サンプルの処置効果の推定値を確認してみましょう。

IPTW法と比べてとる値の範囲は狭くなりましたが、-200や500くらいの絶対値が大きな値が見られます。

各手法の結果まとめ

以下に今回の結果とデータ量を変えてみた(500 -> 250)結果をまとめました。

データ量	正解	回帰分析	IPTW 法	DR 法
250	10	7.3	11.6	9.1
500	10	10.5	13.3	11.6

データ量が少なると線形回帰は結果が悪化し、傾向スコアを用いた手法は結果が良くなりました。とはいってもこのあたりはよく分かりません。(線形回帰の結果が悪化するのはイメージつきますが)

おわりに

今回は「つくりながら学ぶ! Pythonによる因果分析」を参考に因果推論の各手法を簡単に試しましたが、これらの手法でもまだ色々バリエーションはありそうです。例えば傾向スコアはロジスティック回帰ではなく、ニューラルネットワークや樹木系モデル(ランダムフォレスト、勾配ブースティング系)でも求めることができるよう(できそう)です (ここまで書いて思いましたが、シグモイド関数のモデル係数はあまり良くなかった)。また今回は処置効果が線形なデータを使用しましたが、現実の課題は交互作用があったり非線形であったりすることがほとんどです。
どんなものが実務で使えるものかわかりませんがとりあえず勉強続けます。