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

最近、[小川雄太郎さんの因果分析](https://www.amazon.co.jp/つくりながら学ぶ-Pythonによる因果分析-因果推論・因果探索の実践入門-Compass-Science/dp/4839973571/ref=sr_1_1?__mk_ja_JP=カタカナ&dchild=1&keywords=因果分析&qid=1608935190&sr=8-1)を勉強していました。
特に因果探索については、自分が今まで知らないことが多く、非常に面白い内容でした。
因果推論は因果方向が分かっている前提で、因果の効果を推論する。
因果探索は、データを用いて因果方向を探索していく。

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
```

![スクリーンショット 2021-01-01 21.53.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/22cbf970-ffcc-7e6a-17f7-0c01d92e84f6.png)

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

![スクリーンショット 2021-01-01 21.56.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/9682bba4-1e67-7292-0b1d-fd90be486551.png)

#### 各変数の相関を確認
```
# 相関行列
corr_matrix = df.corr()
corr_matrix
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/0884b871-341d-b713-93fa-12377638f2d8.png" width=95%>


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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/9a0aabfe-4e5f-f33b-9810-70b01fca5433.png" width=70%>

```
# 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()
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/e88c5f63-55c0-209c-d960-11d7f58eb691.png" width=30%>

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


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

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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/bbda6d8a-bccb-90f1-0a08-0874afec442b.png" width=50%>

```
# 連続値 →　離散値に変更
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()
```

![スクリーンショット 2021-01-02 7.58.45.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/8075b4b5-39a4-a34b-3af0-1ddc45be8b3b.png)




# 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")
```

![Diagram_Boston.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/89fa0791-eb2b-58d6-ce53-58f1c01218e3.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と無難な方向性を当てはめた。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/5e93d40a-4207-5185-194f-ef336d55e1e7.png" width=80%>

# まとめ
因果探索は、データを用いて因果方向を探索する。
ベイジアンネットワーク(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によって影響が強められている。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/eca71f15-95fe-2fcb-5a49-727404bb4e9a.png" width=50%>

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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/cc37508d-85ce-3549-94f5-f66ef7a1f667.png" width=50%>

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

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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/bb567980-08da-5416-0be5-3b093b5cc6b7.png" width=50%>

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

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

![スクリーンショット 2021-01-11 21.21.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/61820e6c-a56c-68c8-aca9-c52245ef45f9.png)


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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/49b273a2-5f9e-f3c8-e3bd-4342ceb92a7b.png" width=50%>

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

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

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/981627/6bde241f-dd92-637f-d685-b6836840df51.png" width=90%>

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

# 参考文献
[Titanicデータでベイジアンネットワークを実装](https://qiita.com/YuyaOmori/items/e051f0360d1f9562620b)
[「つくりながら学ぶ! Pythonによる因果分析」、小川雄太郎、マイナビ出版](https://www.amazon.co.jp/つくりながら学ぶ-Pythonによる因果分析-因果推論・因果探索の実践入門-Compass-Science/dp/4839973571/ref=sr_1_1?__mk_ja_JP=カタカナ&dchild=1&keywords=因果分析&qid=1608935190&sr=8-1)
[統計WEB/偏相関係数](https://bellcurve.jp/statistics/course/9593.html)
[How to Calculate Partial Correlation in Python](https://www.statology.org/partial-correlation-python/)