今日も因果推論についての記事を書いていきます。
参考にした本は、もちろん「岩波データサイエンス Vol.3」です。
本記事ではその中でも、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')
そもそも交絡因子とはなにか
交絡因子については過去のこの記事を参考にしてください。
上記事では、交絡因子についてだけではなく、
・因果と相関の違い
・層別解析の実装
についても書いています。
本記事とかなり近い内容ですので、因果推論に興味がある方は是非読んでいただけると嬉しいです。
因果効果の求めるために回帰モデルはどのように使う？
早速ですが、2つの変数があり、それぞれ原因変数と結果変数であるとします。
ここで、結果変数を従属変数(目的変数)、原因変数
を説明変数として、以下のような回帰モデルを構成したとします。
もし、変数に影響を与える交絡因子がなければパラメータ
が結果変数と原因変数の関係の強さを表すことは言うまでもありません。
しかし、ここでは、変数に影響を与える交絡因子
がある場合を考えていきます。
つまり、以下のような重回帰モデルを考えてきます。
改めて説明すると、交絡因子が存在する問題点は、結果変数が変化した際に、
①原因変数が変化したためなのか
②交絡因子が変化したため、原因変数と結果変数が変化したのか
の区別がつかないことです。
したがって、どうにかして交絡因子の値を固定してしまえば良いわけです。
これが、交絡因子がある中での回帰モデルでの因果効果推定の考え方の肝です。
重回帰モデルでは、各変数の傾きパラメータは、他の変数の影響を考えず、その変数の値を変化させたときの、
の変化量を表しています。
つまり、原因変数が結果変数
に与える因果効果の大きさは、上記の重回帰モデルにおける原因変数
の傾きパラメータ
であると解釈できます。
これが、交絡因子がある中での回帰モデルでの因果効果推定の方法です。
ここでは重回帰モデルと仮定しましたが、ポアソン回帰などの一般化線形モデルでも同様に因果効果の大きさを推定できます。
回帰モデルで因果効果を分析
ここでは、連続量の交絡因子がある場合における、
に対する
の因果効果を推定していきます。
データは、参考書と同じように、
・と
の相関係数を約
・は
を正規乱数として、次のように生成。
としました。
pythonでのデータの生成方法は以下のようにやりました。
の生成方法については、記事上部に掲載した層別解析の記事を御覧ください。
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
理想通り、相関係数の値が約となりました。
乱数を利用しているので、ご自身で試される場合、値は少しブレると思います。ご了承ください。
では、も生成し、
との関係を図示してみます。
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()
きれいな相関が見られますね。
では、交絡因子を考慮せず、このまま線形回帰した場合の因果効果を計算してみます。
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]]
設定したよりもかなり大きく因果効果を推定してしまっていることがわかります。
ここまでは層別解析のときと同じですね。
では、交絡因子を交えて、
から
への因果効果を計算していきます。
重回帰分析のため、と
を同じアレイに統合します。
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]]
ここで、の傾きパラメータは出力の左側なので、
と推定されたことがわかります。
設定したとかなり近い値が推定されました。
したがって、このように回帰モデルを用いることで、交絡因子の影響を排除して因果効果を正しく推定できることがわかります。
ただし、この例で、ここまで正しく因果効果を推定することができたのは、データを生成したモデルが重回帰モデルであることがわかっていたためです。
仮に、用いたモデルがデータを生成したモデルと異なっている場合、因果効果は正しく推定することができません。
気をつけましょう。
そうした場合には、層別解析等、他の手法も試すとうまくいくこともあります。
最後に
今回は、岩波データサイエンス Vol.3を参考に回帰モデルを利用して因果効果を分析してみました。
うまく結果が出て一安心です。
今日はここまでです。