早稲田卒のAI研究者のブログ

データサイエンス、機械学習、競技プログラミング、音楽、勉強

【岩波データサイエンス Vol.3】交絡因子がある場合の因果効果を、回帰モデルを利用してPythonで分析してみた

今日も因果推論についての記事を書いていきます。

参考にした本は、もちろん「岩波データサイエンス Vol.3」です。

www.iwanami.co.jp

本記事ではその中でも、1章にある交絡因子がある場合の因果効果を、回帰モデルを利用して分析していきます。

必要なパッケージ


変数は基本的にnumpyで処理し、図の作成にはmatplotlibを使用することにします。
また、重回帰係数を求めるためにscikit-learnのlinear_modelを用います。

そのため、以下の3つのパッケージをインポートしておきます。

matplotlibのスタイルはお好みで変更してください。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import linear_model

plt.style.use('ggplot')

そもそも交絡因子とはなにか

交絡因子については過去のこの記事を参考にしてください。

www.mashkun.com

上記事では、交絡因子についてだけではなく、
・因果と相関の違い
・層別解析の実装
についても書いています。

本記事とかなり近い内容ですので、因果推論に興味がある方は是非読んでいただけると嬉しいです。

因果効果の求めるために回帰モデルはどのように使う?


早速ですが、2つの変数X, Yがあり、それぞれ原因変数と結果変数であるとします。
ここで、結果変数Yを従属変数(目的変数)、原因変数Xを説明変数として、以下のような回帰モデルを構成したとします。
\hat{Y} = a + bX

もし、変数X, Yに影響を与える交絡因子がなければパラメータbが結果変数と原因変数の関係の強さを表すことは言うまでもありません。

しかし、ここでは、変数X, Yに影響を与える交絡因子Zがある場合を考えていきます。
つまり、以下のような重回帰モデルを考えてきます。
\hat{Y} = a + bX + cZ

改めて説明すると、交絡因子が存在する問題点は、結果変数が変化した際に、
原因変数が変化したためなのか
交絡因子が変化したため、原因変数と結果変数が変化したのか
の区別がつかないことです。

したがって、どうにかして交絡因子Zの値を固定してしまえば良いわけです。

これが、交絡因子がある中での回帰モデルでの因果効果推定の考え方の肝です。

重回帰モデルでは、各変数の傾きパラメータは、他の変数の影響を考えず、その変数の値を1変化させたときの、Yの変化量を表しています。

つまり、原因変数Xが結果変数Yに与える因果効果の大きさは、上記の重回帰モデルにおける原因変数Xの傾きパラメータbであると解釈できます。

これが、交絡因子がある中での回帰モデルでの因果効果推定の方法です。



ここでは重回帰モデルと仮定しましたが、ポアソン回帰などの一般化線形モデルでも同様に因果効果の大きさを推定できます。

回帰モデルで因果効果を分析

ここでは、連続量の交絡因子Zがある場合における、Yに対するXの因果効果を推定していきます。

データは、参考書と同じように、
XZの相関係数を約0.8
Yeを正規乱数として、次のように生成。
    Y=1.5X+1.1Z+e
としました。

pythonでのデータの生成方法は以下のようにやりました。
Xの生成方法については、記事上部に掲載した層別解析の記事を御覧ください。

size = 400  # データ数を400とする
r = 0.8  # xとzの相関係数を設定
z = np.random.randn(size)  # 標準正規分布から生成
z_lie = np.random.randn(size)  # x生成のための仮データ
x = r * z + (1 - r ** 2) ** 0.5 * z_lie  # zと相関係数が0.8となるxを生成

r_out = np.corrcoef(x, z)[1, 0]  # x, zの相関係数
print("correlation: {}".format(r_out))  # 相関係数の出力

# correlation: 0.8023020385889894


理想通り、相関係数の値が約0.8となりました。
乱数を利用しているので、ご自身で試される場合、値は少しブレると思います。ご了承ください。

では、Yも生成し、Xとの関係を図示してみます。

y = 1.5 * x + 1.1 * z + np.random.randn(400)  # y

# 散布図
plt.scatter(x, y, s = 8)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

f:id:mashkun:20191023220721p:plain

きれいな相関が見られますね。
では、交絡因子を考慮せず、このまま線形回帰した場合の因果効果を計算してみます。

x_lm = x.reshape(-1, 1)  # xを二次元に
y_lm = y.reshape(-1, 1)  # yを二次元に

# 線形回帰
lm = linear_model.LinearRegression()
lm.fit(x_lm, y_lm)

print(lm.coef_)  # xの傾きパラメータを出力

# [[2.43229932]]

設定した1.5よりもかなり大きく因果効果を推定してしまっていることがわかります。

ここまでは層別解析のときと同じですね。


では、交絡因子Zを交えて、XからYへの因果効果を計算していきます。

重回帰分析のため、XYを同じアレイに統合します。

z_lm = z.reshape(-1, 1)  # zを2次元に
xz = np.concatenate([x_lm, z_lm], axis=1)  # xとzを統合


そして、scikit-learnで重回帰し、各変数の傾きパラメータを出力していきます。

xz_lm = linear_model.LinearRegression()
xz_lm.fit(xz_0, y_0)

print(xz_lm.coef_)

# [[1.42759184 1.34237261]]

ここで、Xの傾きパラメータは出力の左側なので、1.43と推定されたことがわかります。

設定した1.5とかなり近い値が推定されました。
したがって、このように回帰モデルを用いることで、交絡因子の影響を排除して因果効果を正しく推定できることがわかります。



ただし、この例で、ここまで正しく因果効果を推定することができたのは、データを生成したモデルが重回帰モデルであることがわかっていたためです。
仮に、用いたモデルがデータを生成したモデルと異なっている場合、因果効果は正しく推定することができません。
気をつけましょう。

そうした場合には、層別解析等、他の手法も試すとうまくいくこともあります。

最後に

今回は、岩波データサイエンス Vol.3を参考に回帰モデルを利用して因果効果を分析してみました。
うまく結果が出て一安心です。

今日はここまでです。