Bayesian Networkを用いてボストン住宅価格データセットを因果探索

最近、小川雄太郎さんの因果分析を勉強していました。
特に因果探索については、自分が今まで知らないことが多く、非常に面白い内容でした。
因果推論は因果方向が分かっている前提で、因果の効果を推論する。
因果探索は、データを用いて因果方向を探索していく。

Bayesian Network(PCアルゴリズム)を用いてBoston Housingデータセットを因果探索します。
目的変数で使われる"MEDV(住宅価格)"の因果構造を可視化します。

作成したコードはGoogle colab上で実装し、GitHubにあります。
https://github.com/JunJun100/Bayesian_Network

事前準備

GoogleDriveをマウント

from google.colab import drive
drive.mount('/content/drive')

import os
os.chdir("./drive/My Drive/Causal_Analysys/Boston")

設定

# 乱数のシードを設定
import random
import numpy as np

np.random.seed(1234)
random.seed(1234)

#組み合わせを取得する関数
import itertools

# その他
import pandas as pd
from IPython.display import Image
from sklearn.datasets import load_boston
from sklearn.preprocessing import StandardScaler
import seaborn as sns
%matplotlib inline

# PCアルゴリズム
!pip install pgmpy==0.1.11

# 独立性の検定に使用
from pgmpy.estimators import ConstraintBasedEstimator

# 偏相関
!pip install pingouin
import pingouin as pg

データ作成

boston = load_boston()

df = pd.DataFrame(boston.data, columns=boston.feature_names) # 説明変数
df['MEDV'] = boston.target # 目的変数を追加
df

# 標準化
scaler = StandardScaler()
df_scaled = scaler.fit_transform(df)
df =pd.DataFrame(df_scaled,columns=df.columns)
df.head()

各変数の相関を確認

# 相関行列
corr_matrix = df.corr()
corr_matrix

# ヒートマップ
sns.heatmap(corr_matrix,
            square=True,
            cmap='coolwarm',
            xticklabels=corr_matrix.columns.values,
            yticklabels=corr_matrix.columns.values)

# yとの相関係数を表示
corr_y = pd.DataFrame({"features":df.columns,"corr_y":corr_matrix["MEDV"]},index=None)
corr_y = corr_y.reset_index(drop=True)
corr_y.style.background_gradient()

corr_y = pd.DataFrame({"features":df.columns,"corr_y":corr_matrix["MEDV"]},index=None)
corr_y = corr_y.reset_index(drop=True)
corr_y.style.background_gradient()

MEDV(住宅価格)とLSTAT(低所得者の割合)は強い負の相関(-0.74)、MEDVとRM(平均部屋数)は強い正の相関(0.70)。

データの離散化

PCアルゴリズムは連続値のデータが扱えないので、離散値に変更する必要がある。

# データタイプを確認
df.info()

# 連続値 →　離散値に変更
df_new = df.copy()

df_new['CRIM'] = pd.cut(df_new['CRIM'], 4)
df_new['ZN'] = pd.cut(df_new['ZN'], 4)
df_new['INDUS'] = pd.cut(df_new['INDUS'], 3)
df_new['CHAS'] = pd.cut(df_new['CHAS'], 3)
df_new['NOX'] = pd.cut(df_new['NOX'], 3)
df_new['RM'] = pd.cut(df_new['RM'], 3)
df_new['AGE'] = pd.cut(df_new['AGE'], 4)
df_new['DIS'] = pd.cut(df_new['DIS'], 3)
df_new['RAD'] = pd.cut(df_new['RAD'], 4)
df_new['TAX'] = pd.cut(df_new['TAX'], 4)
df_new['PTRATIO'] = pd.cut(df_new['PTRATIO'], 3)
df_new['B'] = pd.cut(df_new['B'], 4)
df_new['LSTAT'] = pd.cut(df_new['LSTAT'], 3)
df_new['MEDV'] = pd.cut(df_new['MEDV'], 4)

上記の分割数は大きいほど、より精度の高い因果構造を推定できる。
しかし、その分たくさんのデータが必要になる。

離散化されたデータを確認

df_new.head()

PCアルゴリズム

est = ConstraintBasedEstimator(df_new)

skel, seperating_sets = est.estimate_skeleton(significance_level=0.01)
print("Undirected edges: ", skel.edges())

pdag = est.skeleton_to_pdag(skel, seperating_sets)
print("PDAG edges:       ", pdag.edges())

model = est.pdag_to_dag(pdag)
print("DAG edges:        ", model.edges())

## 以下、実行結果
# Undirected edges:  [('INDUS', 'NOX'), ('INDUS', 'RAD'), ('INDUS', 'PTRATIO'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('RM', 'LSTAT'), ('RM', 'MEDV'), ('AGE', 'DIS'), ('AGE', 'LSTAT'), ('RAD', 'TAX'), ('PTRATIO', 'MEDV'), ('LSTAT', 'MEDV')]
# PDAG edges:        [('NOX', 'INDUS'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('RM', 'LSTAT'), ('RM', 'MEDV'), ('DIS', 'AGE'), ('RAD', 'INDUS'), ('RAD', 'TAX'), ('TAX', 'NOX'), ('TAX', 'RAD'), ('PTRATIO', 'INDUS'), ('LSTAT', 'MEDV'), ('MEDV', 'RM')]
# DAG edges:         [('NOX', 'INDUS'), ('NOX', 'PTRATIO'), ('NOX', 'TAX'), ('RM', 'LSTAT'), ('RM', 'MEDV'), ('DIS', 'AGE'), ('RAD', 'INDUS'), ('TAX', 'RAD'), ('PTRATIO', 'INDUS'), ('LSTAT', 'MEDV')]

# 可視化

# MEDVとPTRATIOの因果方向は決まらなかった
Image("Diagram_Boston.png")

LSTAT→MEDV、RM→MEDVの因果構造が可視化された。
ただし、LSAT→RM→MEDVの因果構造となっており擬似相関を生み出している。
MEDVとPTRATIOの因果方向は決まらなかった。

MEDV(住宅価格)とLSTAT(低所得者の割合)は強い負の相関(-0.74)、MEDVとRM(平均部屋数)は強い正の相関(0.70)

# 可視化

# MEDVとPTRATIOは、無難な方向性を当てはめる
Image("Diagram_2.png")

MEDVとPTRATIO(教師当たりの生徒数)の因果方向は決まらない。
PTRATIO→MEDVと無難な方向性を当てはめた。

まとめ

因果探索は、データを用いて因果方向を探索する。
ベイジアンネットワーク(PCアルゴリズム)を用いてボストン住宅価格データセットを因果探索すると、
LSTAT→MEDV、RM→MEDVの因果構造が可視化された。
ただし、LSAT→RM→MEDVの因果構造となっており擬似相関を生み出している。
MEDVとPTRATIOの因果方向は決まらず、無難な方向性(PTRATIO→MEDV)を当てはめた。

補足

Pythonで偏相関を計算できるライブラリpingouinがある。
https://pingouin-stats.org/generated/pingouin.partial_corr.html

偏相関

RMとMEDVは強い相関があるが、これはLSTATによって影響が強められている。

RMとMEDVの相関から、LSTATの影響を取り除いたものを偏相関という。
その値を、偏相関係数という。

偏相関係数の定義式

「zの影響を除いたx」と「zの影響を除いたy」の偏相関係数は、

\rho_{xy,z} = \frac{r_{xy} - r_{xz} r_{yz}}{ \sqrt{1 - r^2_{xz}} \sqrt{ 1 - r^2_{yz} } }

偏相関計数を計算

RMとMEDVの偏相関計数を計算すると、

# LSTATの影響を取り除き、RMとMEDVの偏相関係数を計算
pg.partial_corr(data=df, x='RM', y='MEDV', covar='LSTAT')

RMとMEDVの偏相関係数は、0.46となり正の相関がある。

ついでに偏相関行列を計算すると、

#　すべての変数について偏相関計数を計算 (小数点以下3桁まで)
df.pcorr().round(3)

RMとMEDVには偏相関がある

RMとMEDVは強い正の相関(0.70)がある。
LSTATの影響を取り除くと、RMとMEDVは偏相関(0.46)がある。

参考文献

Titanicデータでベイジアンネットワークを実装
「つくりながら学ぶ! Pythonによる因果分析」、小川雄太郎、マイナビ出版
統計WEB/偏相関係数
How to Calculate Partial Correlation in Python