【岩波データサイエンス Vol.3】相関関係と因果関係の違いについて説明、更に層別解析をPythonで実装してみた

今日は珍しく、競技プログラミングについてではなく、データサイエンスについてです。

本業なのでこっちをメインにしたいんですが、競プロが楽しいのが問題ですね。
また、こういう記事って結構時間かかるので、時間的余裕がないと情報がたくさんある良い記事を書けないというのもあります。

まぁそんな前置きはさておき、本記事では因果と相関の違い、層別解析について説明していきます。
参考にした本は「岩波データサイエンス Vol.3」です。

www.iwanami.co.jp

データ分析の中でも、因果推論に焦点を当てた本となっています。

因果推論という分野は、もちろん需要はかなりあります。

ですが、データから変数間の因果を推論するというのは、難しいことであり、データ間の因果を完全に推測する手法は見つかっていません。
そのため、そんな手法をめざして今なお活発に研究が進んでいます。
最近だと、Deepmindの研究者らがメタ強化学習による因果推論(Causal reasoning from meta-reinforcement learning)という論文を出していました。
arxiv.org

そんな因果推論の分野ですが、僕も非常に興味があるので勉強をしているところです。
そのまとめを皆様に共有できればと思います。

今日はそのなかでも、岩波データサイエンス Vol.3の1章にある層別解析を実際にpythonで実装・確認していきます。

必要なパッケージ

変数は基本的にnumpyで処理していくことにします。
また、図の作成にはmatplotlibを使用します。

線形回帰をした図の作成用にscikit-learnを、回帰係数の標準誤差を求めるために、scipyのstatsを使っていきます。

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

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

plt.style.use('ggplot')

ちなみに僕はggplot好きです。
そのため、matplotlibのスタイルを変更していますが、変更なしでももちろん問題ありません。

そもそも因果と相関の違いとはなにか

f:id:mashkun:20191022202856p:plain

２つの変数 $X, Y$ があるとします。
それらの関係が、上図のような散布図だとした場合、あたかも $X$ を大きくすると $Y$ も大きくなるように見えますよね。

はたして、それは正しい解釈なのでしょうか。
実は散布図からだけでは、
$X$ の値が大きいときに、 $Y$ の値が大きいという単純な直線関係があるだけなのか、
$X$ の値を大きくしたときに、 $Y$ の値が大きくなるのか、
はわかりません。

ここで、前者が相関関係、後者が因果関係に対応します。

つまり、因果関係というのは、要因 $X$ を変化させたときに要因 $Y$ も変化する、という関係を表します。
また、上記のような原因をしめす変数を原因変数、結果を表す変数を結果変数と言います。

なぜ相関関係が必ずしも因果関係を意味しないのか

では、なぜ相関関係が必ずしも因果関係を意味しないのでしょうか。

その理由について、因果ダイヤグラム(因果グラフ)と呼ばれる図を使って説明していきます。

f:id:mashkun:20191022204204p:plain

上の図はXとYに相関が観察される際の因果グラフです。
矢印は、矢印の根っこの変数を原因、矢印の先の変数を結果とする因果関係を表しています。

図(a)のような場合は $X$ と $Y$ の間に、 $X$ を原因、 $Y$ を結果とする因果関係があることがわかります。
ですが、図(b)の場合はどうでしょうか。

$X$ と $Y$ の間に相関はありますが、因果関係はありません。
というのも、 $X$ と $Y$ に共通の要因 $C$ の影響により、相関があるだけだからです。
このような関係を考慮する必要があるため、相関関係が必ずしも因果関係を意味しないこととなります。

この構造を交絡と呼び、 $C$ のような要因を交絡因子と呼びます。
そして、交絡によって生じる本当は関連のない変数間の相関を疑似相関と呼びます。

実際のデータを分析する際にはこの交絡に対処しなければなりません。
というのも、観察された相関関係が因果関係によって生じたものなのか、交絡によって生じたものなのかを見分けなければ、間違った意思決定をしてしまうからです。
では、それらを見分ける手法としてどのようなものがあるのでしょうか。

その一つに層別解析と呼ばれるものがあります。

層別解析とは

そもそも交絡因子が起こしてしまう問題とは、結果変数が変化した際に、
①原因変数が変化したためなのか
②交絡因子が変化したため、原因変数と結果変数が変化したのか
の区別がつかないことです。

そこで、交絡因子の値を固定してしまおう、という考えが層別解析の考え方です。

具体的な手順として、
①交絡因子の値によって対象をいくつかの層に分ける
②層ごとに解析を行う
③それらの結果を統合して、全体における原因変数と結果変数の関係の推定をする
となります。

今は、交絡因子が一つの場合について説明していますが、もちろん交絡因子が複数あっても層別解析を行うことはできます。

では層別解析についてわかったところで、実装していきます。

層別解析の例

ここでは、連続量の交絡因子 $Z$ を離散化して層別解析を行い、 $Y$ に対する $X$ の因果効果を推定していきます。

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

上記のデータについてpythonで分析していきます。
以下のようにデータを生成しました。

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

"""
$X$ の生成に関しては以下の記事を参考にしました。
ありがとうございました。
丁寧に理論も記載されているので、興味ある方はぜひご覧ください。

相関のある2つの擬似乱数の生成（Pythonサンプルつき） - Qiita
"""

きちんとデータが生成されたか見てみましょう。

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

# correlation: 0.8012456796209242

理想通りデータを生成することができました。

では、結果変数である、 $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:20191022213104p:plain

予想通り相関があることがわかります。
これだけ見るとやはり、 $X$ と $Y$ に因果関係があると思ってしまいますね。

ではまず、交絡因子を考慮せず、 $X$ と $Y$ の関係を分析してみます。
つまり、すべてのデータを用いて、 $\hat{Y} = a + bX$ における $X$ の回帰係数 $b$ を計算していきます。

Pythonではscikit-learnを用いて以下のように求めていきます。

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_)  # 回帰係数の出力

# [[2.4127323]]

この値が、 $Y$ に対する $X$ の因果効果の推定値となります。
生成したデータでは、 $X$ の傾きのパラメータを $1.5$ としていたので、因果効果が過大に推定されていることがわかります。

次に交絡因子を考慮して分析していきます。
今回は $Z$ を4分割に離散化し、それぞれを $Z_{0}, Z_{1}, Z_{2}, Z_{3}$ とします。
そして、それぞれの層ごとに、 $\hat{Y_{i}} = a_{i} + b_{i}X_{i}$ のモデルで傾きのパラメータを推定していきます。

Pythonでの実装では以下のように $Z$ をパーセンタイルで分割しました。

xyz = np.concatenate([x.reshape(1, -1), y.reshape(1, -1), z.reshape(1, -1)])  # 3変数を同じアレイに格納

xyz_0 = np.zeros(3).reshape(3,1)  # 3行1列のアレイを用意
xyz_1 = np.zeros(3).reshape(3,1)
xyz_2 = np.zeros(3).reshape(3,1)
xyz_3 = np.zeros(3).reshape(3,1)

for i in range(size):
    if xyz[2,i] <= np.percentile(z, 25):  # 25パーセンタイル以下
        xyz_0 = np.append(xyz_0, xyz[:,i].reshape(-1,1), axis=1)
    elif np.percentile(z, 25) < xyz[2,i] <= np.percentile(z, 50):  # 25パーセンタイル以上50パーセンタイル以下
        xyz_1 = np.append(xyz_1, xyz[:,i].reshape(-1,1), axis=1)
    elif np.percentile(z, 50) < xyz[2,i] <= np.percentile(z, 75):  # 50パーセンタイル以上75パーセンタイル以下
        xyz_2 = np.append(xyz_2, xyz[:,i].reshape(-1,1), axis=1)
    else:  # 75パーセンタイル以上
        xyz_3 = np.append(xyz_3, xyz[:,i].reshape(-1,1), axis=1)

xyz_0 = np.delete(xyz_0, 0, 1)  # 1列目のすべて0の要素をドロップ
xyz_1 = np.delete(xyz_1, 0, 1)
xyz_2 = np.delete(xyz_2, 0, 1)
xyz_3 = np.delete(xyz_3, 0, 1)

これで、 $Z$ の大きさによって $X$ と $Y$ の値を層分けすることができました。

散布図を書いて確認してみましょう。

# 散布図作成
plt.scatter(x_0, y_0, s = 8, label = "i = 0")
plt.scatter(x_1, y_1, s = 8, label = "i = 1")
plt.scatter(x_2, y_2, s = 8, label = "i = 2")
plt.scatter(x_3, y_3, s = 8, label = "i = 3")
plt.title("X - Y")
plt.xlabel("X")
plt.ylabel("Y")
plt.legend()

f:id:mashkun:20191022220542p:plain

重なって見にくいですが、4つの層に分割できていることがわかります。

では、各層について線形回帰してみます。
コードは少し長いですが、同じことを4回しているだけですので、単純です。

# データ抽出
x_0 = xyz_0[0, :].reshape(-1, 1)
y_0 = xyz_0[1, :].reshape(-1, 1)
x_1 = xyz_1[0, :].reshape(-1, 1)
y_1 = xyz_1[1, :].reshape(-1, 1)
x_2 = xyz_2[0, :].reshape(-1, 1)
y_2 = xyz_2[1, :].reshape(-1, 1)
x_3 = xyz_3[0, :].reshape(-1, 1)
y_3 = xyz_3[1, :].reshape(-1, 1)

# 線形回帰
lm_0 = linear_model.LinearRegression()
lm_0.fit(x_0, y_0)
lm_1 = linear_model.LinearRegression()
lm_1.fit(x_1, y_1)
lm_2 = linear_model.LinearRegression()
lm_2.fit(x_2, y_2)
lm_3 = linear_model.LinearRegression()
lm_3.fit(x_3, y_3)

# 回帰係数計算
b_0 = lm_0.coef_[0, 0]
b_1 = lm_1.coef_[0, 0]
b_2 = lm_2.coef_[0, 0]
b_3 = lm_3.coef_[0, 0]

# 作図
fig = plt.figure(figsize = (12,8))
ax0 = fig.add_subplot(221)
ax1 = fig.add_subplot(222)
ax2 = fig.add_subplot(223)
ax3 = fig.add_subplot(224)

ax0.scatter(x_0, y_0, s = 8)
ax0.plot(x_0, lm_0.predict(x_0), color = 'red', label = "b = {}".format(round(b_0,2))) 
ax0.set_title("i = 0")
ax0.set_xlabel("X_0")
ax0.set_ylabel("Y")
ax0.legend(loc="upper left")

ax1.scatter(x_1, y_1, s = 8)
ax1.plot(x_1, lm_1.predict(x_1), color = 'red', label = "b = {}".format(round(b_1,2))) 
ax1.set_title("i = 1")
ax1.set_xlabel("X_1")
ax1.set_ylabel("Y")
ax1.legend(loc="upper left")

ax2.scatter(x_2, y_2, s = 8)
ax2.plot(x_2, lm_2.predict(x_2), color = 'red', label = "b = {}".format(round(b_2,2))) 
ax2.set_title("i = 2")
ax2.set_xlabel("X_2")
ax2.set_ylabel("Y")
ax2.legend(loc="upper left")

ax3.scatter(x_3, y_3, s = 8)
ax3.plot(x_3, lm_3.predict(x_3), color = 'red', label = "b = {}".format(round(b_3,2))) 
ax3.set_title("i = 3")
ax3.set_xlabel("X_3")
ax3.set_ylabel("Y")
ax3.legend(loc="upper left")

plt.tight_layout()  # 図の体裁を整えるため

f:id:mashkun:20191022215742p:plain

各層における回帰係数が求まりました。
最後にこれらを統合すれば、 $Y$ に対する $X$ の因果効果の推定値を求めることができあます。
統合の方法は、単なる平均ではなく、回帰係数をその分散の逆数で重みをついた平均です。

分散を求めるためには回帰係数の標準誤差が必要になります。
scilit-learnだと難しいので、scipyのstatsを用いて標準誤差を求めます。

# 標準誤差を入手
_, _, _, _, std_0 = stats.linregress(xyz_0[0, :], xyz_0[1, :])
_, _, _, _, std_1 = stats.linregress(xyz_1[0, :], xyz_1[1, :])
_, _, _, _, std_2 = stats.linregress(xyz_2[0, :], xyz_2[1, :])
_, _, _, _, std_3 = stats.linregress(xyz_3[0, :], xyz_3[1, :])

# 分散を計算
r_var_0 = 1 / (std_0**2)
r_var_1 = 1 / (std_1**2)
r_var_2 = 1 / (std_2**2)
r_var_3 = 1 / (std_3**2)

# 分散の重み付け平均を計算
mean_ = ((r_var_0 * b_0) + (r_var_0 * b_1) + (r_var_0 * b_2) + (r_var_0 * b_3)) \
        / (r_var_0 + r_var_1 + r_var_2 + r_var_3)

print(mean_)

# 1.6984607973650623