ベイジアンネットワークを活用して因果探索を実装してみた #Python

製造業出身のデータサイエンティストがお送りする記事
今回は、ベイジアンネットワークを活用して、因果探索を実装してみました。

はじめに

前回（統計的因果探索のLiNGAM）からの続きで、データから因果関係を推測するための機械学習技術の1つ「ベイジアンネットワーク」について整理しました。
細かい理論的な部分は省略します。

ベイジアンネットワークとは

ベイジアンネットワークは条件付き確率をベースに、変数間のグラフ構造を表現したものです。
ベイジアンネットワークの中でも、依存関係を算出する方法はいくつかありますが、今回は以下の方法を実施しました。

独立性の検定により、変数間のスケルトン(方向のない依存関係)を特定
V-structure,オリエンテーションルールにより一部のスケルトンに方向を付与

上記はPCアルゴリズムと呼ばれる手法で、簡単に実装することができました。

今回もUCI Machine Learning Repositoryで公開されているボストン住宅の価格データを用いてベイジアンネットワークを構築します。

pythonのコードは下記の通りです。

# 必要なライブラリーのインポート
import pandas as pd
import numpy as np
import networkx as nx
import os
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston

from pgmpy.models import BayesianModel
from pgmpy.estimators import ConstraintBasedEstimator
from pgmpy.estimators import BicScore

次にデータを読み込みます。

# データセットの読込み
boston = load_boston()

# データフレームの作成
# 説明変数の格納
df = pd.DataFrame(boston.data, columns = boston.feature_names)

# 目的変数の追加
df['MEDV'] = boston.target

# データの中身を確認
df.head()

次に、PCアルゴリズムは連続値は扱えないため、離散値へ変更します。

# 離散値に変更
df_discrete = df.copy()

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

df_discrete.head()

ここから、ライブラリーpgmpyを用いてベイジアンネットワークを構築していきます。

est = ConstraintBasedEstimator(df_discrete)

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は、方向性なし
PDAG edgesは、部分的な方向性あり
DAG edgesは、方向性あり

得られた結果は下記です。

Undirected edges:  [('ZN', 'DIS'), ('INDUS', 'NOX'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('NOX', 'DIS'), ('NOX', 'PTRATIO'), ('RM', 'MEDV'), ('AGE', 'LSTAT'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('TAX', 'PTRATIO'), ('LSTAT', 'MEDV')]
PDAG edges: [('ZN', 'DIS'), ('INDUS', 'NOX'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('NOX', 'INDUS'), ('RM', 'MEDV'), ('AGE', 'LSTAT'), ('DIS', 'ZN'), ('DIS', 'NOX'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('TAX', 'INDUS'), ('TAX', 'RAD'), ('TAX', 'PTRATIO'), ('PTRATIO', 'INDUS'), ('PTRATIO', 'RAD'), ('LSTAT', 'AGE'), ('LSTAT', 'MEDV'), ('MEDV', 'RM'), ('MEDV', 'LSTAT')]
DAG edges: [('INDUS', 'NOX'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('DIS', 'NOX'), ('DIS', 'ZN'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('TAX', 'PTRATIO'), ('LSTAT', 'AGE'), ('MEDV', 'RM'), ('MEDV', 'LSTAT')]

上記で基本的な処理は終了ですが、可視化してみないと分かりにくいので可視化しました。
ただ、LiNGAMみたいに簡単に可視化してくれる方法がなかったため（あったら教えて欲しいです）、Undirected edges方向性なしと、DAG edges方向性ありだけ簡単に可視化しました。

# Undirected edges（方向性なし）可視化
G = nx.Graph()
G.add_nodes_from(df.columns) # ノード
G.add_edges_from(skel.edges()) # 辺
nx.draw_circular(G, with_labels=True, arrowsize=30, node_size=800, alpha=0.3, font_weight='bold')
plt.show()

]:
# DAG edges(方向性あり)可視化
DAG_model = BayesianModel(model.edges())
nx.draw_circular(DAG_model, with_labels=True, arrowsize=30, node_size=800, alpha=0.3, font_weight='bold')
plt.show()

さいごに

最後まで読んで頂き、ありがとうございました。
ベイジアンネットワークも結構簡単に利用できることが分かりました。可視化の部分は、いまいち良い感じに実装できていないので、今後の課題です。

訂正要望がありましたら、ご連絡頂けますと幸いです。