Titanicデータでベイジアンネットワークを実装

先日、DeepLearningを用いたネットワーク分析手法である**SAM(Structural Agnostic Modeling)**をTitanicデータで実装した記事を記載しました。
今回はネットワーク分析の最もポピュラーな手法である**ベイジアンネットワーク**をTitanicデータで実装しましたので紹介します。

SAMの記事は以下をご参照ください。
[SAMの記事](https://qiita.com/YuyaOmori/items/a9c5f9321f9387fe74e1)

記事で紹介していないものも含めて、コードは全てGitHubに置いてあります。
[https://github.com/yuomori0127/Bayesian_Network](https://github.com/yuomori0127/Bayesian_Network)


#目次
1. [ゴールと概要](#1-ゴールと概要)
2. [ベイジアンネットワークとは](#2-ベイジアンネットワークとは)
    1. [独立性の検定](#2-1-独立性の検定)
      1. [0次の独立性の検定](#2-1-1-0次の独立性の検定)
      2. [1次、2次の独立性の検定](#2-1-2-1次2次の独立性の検定)
    2. [V-structure,オリエンテーションルール](#2-2-v-structureオリエンテーションルール)
      1. [V-structure](#2-2-1-v-structure)
      2. [オリエンテーションルール](#2-2-2-オリエンテーションルール)
    3. [スコアベースアプローチ(BIC)](#2-3-スコアベースアプローチbic)
3. [ベイジアンネットワークを関数で簡単に実装](#3-ベイジアンネットワークを関数で簡単に実装)
4. [最後に](#4-最後に)

#1. ゴールと概要
ネットワーク分析のゴールについては以前の記事と重複するので割愛します。
以前の記事をご参照ください。
[ネットワーク分析のゴール](https://qiita.com/YuyaOmori/items/a9c5f9321f9387fe74e1#1-%E3%82%B4%E3%83%BC%E3%83%AB)

簡単に言うと以下のように変数間のグラフ構造を明らかにすることです。
![BIC.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/1f865020-b449-936c-5024-798d7062125c.png)

#2. ベイジアンネットワークとは
ベイジアンネットワークは**条件付き確率をベースに、変数間のグラフ構造を表現したもの**です。

ベイジアンネットワークの中でも、依存関係(矢印)を算出する方法はいくつかあるのですが、
今回は以下の手順で実施しました。

```
1. 独立性の検定により、変数間のスケルトン(方向のない依存関係)を特定
2. V-structure,オリエンテーションルールにより一部のスケルトンに方向を付与
3. BICを用いて残りのスケルトンに方向を付与
```
1,2,は**PCアルゴリズム**[1]と呼ばれる方法で、
3はスコアベースアプローチと呼ばれる方法です。

#2-1. 独立性の検定
まずは下図のようなスケルトン(方向のない依存関係)を特定します。
![スケルトン.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/9467edfd-a367-30fc-d0e1-58cd495338dd.png)

スケルトンを特定するためには**独立性の検定**を行います。
これは**カイ二乗検定**と呼ばれ、**変数間に関連があるかどうかを検定**するものです。
カイ二乗検定 参考：https://bellcurve.jp/statistics/course/9496.html

また、単純な変数間の独立性の検定に加え、**条件付き独立の検定**も実施します。
これはある変数Aを条件としたとき(変数Aを固定した時)の変数B・変数C間の独立性の検定で、これを2次(条件が2変数)まで行います。

つまり、以下3つの独立性の検定を行います。

```
1. 0次の独立性の検定(単純な2変数間の独立性の検定)
2. 1次の独立性の検定(1変数を条件としたときの2変数間の独立性の検定)
3. 2次の独立性の検定(2変数を条件としたときの2変数間の独立性の検定)
```

##2-1-1. 0次の独立性の検定
言語はPythonで、ライブラリは[pgmpy](https://pgmpy.org/index.html)を使用しました。
99%有意で2変数間の検定を行います。

```
#0次の独立性を検定する関数定義
def independence_0(df):
  independence_list=[]
  est = ConstraintBasedEstimator(df)
  for col1,col2 in list(itertools.combinations(df.columns, 2)):
    result = est.test_conditional_independence(col1, col2, method="chi_square", tol=0.01)
    print("{}-{} -> {}".format(col1,col2,result))
    if not result:
      independence_list.append([col1,col2]) 
  return independence_list
```

0次の独立性検定の実行し、全変数間の検定結果を表示させています。
また、検定結果が**False**であるもの(独立ではないと判定されたもの)の一覧を`independence_list_0`に格納しています。

```
#0次の独立性検定の実行
independence_list_0 = independence_0(df)
# 実行結果
# Survived-Pclass -> False
# Survived-Sex -> False
# Survived-Age -> True
# Survived-SibSp -> False
# 以下略
```

検定結果が**False**であるもの(=独立ではない変数)同士を繋げて可視化すると以下のようになりました。
![0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/2119a6b0-39e7-6d30-9a55-07ba30e61192.png)
相当ごちゃごちゃしています。。

ここから1次、2次の独立性の検定を行い、変数間の関係性をさらに見ていきます。

##2-1-2. 1次、2次の独立性の検定
1次の独立性の検定とは、ある1つの変数を条件付きにした際(1つの変数の値を固定した際)の2つの変数間の独立性の検定です。
コードは長くなるので[GitHub](https://github.com/yuomori0127/Bayesian_Network)に置いてあります。

1次の独立性の検定の結果、以下のようになりました。
![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/cd4268ac-6317-d585-75e3-e01fb19b9985.png)
相当スッキリしました。

同様に2次の独立性の検定も実施し、以下のようになりました。
![スケルトン.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/3cf24440-bf58-4a39-8743-f4e6a61df6fc.png)
右下の「Fare」が消えました。
これでスケルトンは完成です。
次はこの図に方向を付与していきます。


#2-2. V-structure,オリエンテーションルール
スケルトンに方向を付与する手段として、**V-structure,オリエンテーションルール**の2つを実施します。

#2-2-1. V-structure
V-structureとは文字通り、V字構造を持った変数間の関係性を指します。
![V.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/923ca7d8-4bc8-3e2f-5f66-315537ca8bce.png)
この時、V字の中心(例では「Survived」)を条件に置き、残りの2つの変数(例では「Pclass」と「Sex」)間で独立性が**無ければ**、以下の矢印を付与することができます。
![V2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/dce0df79-92a9-99f3-a467-32183271d8a9.png)
これはつまり、「Pclass」と「Sex」が決まれば「Survived」も決まりやすい関係にあることを意味します。

このように、V字構造を持つ全ての組み合わせに対して条件付き独立性の検定を行います。
こちらも99%有意で検定しています。

```
# V-structureの探索(条件付き非独立の場合、連結できる)
#①Pclass-Survived-Sex
print(est.test_conditional_independence(
    'Pclass', 'Sex', ['Survived'], method="chi_square", tol=0.01))

#②Age-Pclass-Embarked
print(est.test_conditional_independence(
    'Age', 'Embarked', ['Pclass'], method="chi_square", tol=0.01))

#③Survived-Sex-Parch
print(est.test_conditional_independence(
    'Survived', 'Parch', ['Sex'], method="chi_square", tol=0.01))

#④Sex-Parch-SibSp
print(est.test_conditional_independence(
    'Sex', 'SibSp', ['Parch'], method="chi_square", tol=0.01))

# 実行結果
# False
# True
# True
# True
```

検定の結果、`①Pclass-Survived-Sex`間のみ**False**となり、**独立ではない(=関係性がある)**という結果となりました。
そのため`①Pclass-Survived-Sex`間に矢印が付与され、図は以下のようになります。
![DAG.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/4a9fed21-bcbe-3dea-f168-e373db552dac.png)

#2-2-2. オリエンテーションルール
V-structureを含む部分的な矢印が得られている場合、付加的に矢印の向きを決定できることがあります。
このルールを**オリエンテーションルール**[3]と呼びます。

以下のパターンの場合、付加的に矢印の向きを決定することができます。[3]
![オリエンテーションルール.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/5182b463-a80e-1e91-bbad-bd2cc29eadd5.png)


今回の例ではオリエンテーションルールにあてはまる構造がないため、図は変わりません。
![DAG.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/862d4808-0289-5816-c677-f78fdca2775e.png)


#2-3. スコアベースアプローチ(BIC)
残りのスケルトンに対して矢印を決定する方法として、
**BICによるスコアベースアプローチ**を考えます。

BICの詳細説明は省略します。以下の文献をご参考ください。

 - [4]「作りながら学ぶ！Pythonによる因果分析」(小川 雄太郎 著 2020)
 - [WikiPedia](https://ja.wikipedia.org/wiki/%E3%83%99%E3%82%A4%E3%82%BA%E6%83%85%E5%A0%B1%E9%87%8F%E8%A6%8F%E6%BA%96)

要は**グラフィカルモデルの評価指標**です。

下図でスケルトンとなっている線は4本ありますが、
そのすべてに対して2方向の矢印を仮定した場合、 $ 2^4=16 $ 通りのグラフィカルモデルが考えられます。
その**16パターン全てをBICで評価**し、最も評価値が高いものを採用するという方法で残りのスケルトンに対して矢印を決定していきます。

![DAG.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/a241da94-81a9-34c7-9bf9-f059e014fcd6.png)

16パターン全てをここには記載しませんが、
最も評価が高いもので約-4936.9となりました。(BICはマイナスが大きい方が良いモデル)

```
model = BayesianModel([('Pclass', 'Survived'), ('Sex', 'Survived'),('Age', 'Pclass'), ('Embarked', 'Pclass'),('Sex', 'Parch'), ('SibSp', 'Parch')])  
print(BicScore(df).score(model))
# 実行結果
# -4936.911652064051
```

最もBICでの評価が高いモデルは下図となりました。
![BIC.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/5812463e-99be-4be8-8a46-81c1f7ad8c6e.png)

#3. ベイジアンネットワークを関数で簡単に実装
今まで色々実施してきましたが、実は似たようなことを簡単に実施する関数が存在します。

```
est = ConstraintBasedEstimator(df)

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:  [('Survived', 'Pclass'), ('Survived', 'Sex'), ('Pclass', 'Age'), ('Pclass', 'Embarked'), ('Sex', 'Parch'), ('SibSp', 'Parch')]
# PDAG edges:        [('Age', 'Pclass'), ('SibSp', 'Parch'), ('Parch', 'Sex'), ('Parch', 'SibSp'), ('Embarked', 'Pclass')]
# DAG edges:         [('Age', 'Pclass'), ('Parch', 'Sex'), ('Parch', 'SibSp'), ('Embarked', 'Pclass')]
```

 - 下から3つ目の`Undirected edges`は方向性なし
 - 下から2つ目の`PDAG edges`は部分的な方向性あり
 - 一番下の`DAG edges`は方向性あり

を意味します。

これらをもとに図を描くと以下のようになります。
![pgmpy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/173913/7a16de56-a726-7ace-1878-2e73b6a4cc0b.png)

「Survived」のみ方向性がないですが、それ以外は**全く同じ結果**となりました。
このことからどちらの方法でも**ある程度同じ結果が得られる**ことがわかります。
試しにBICを算出してみます。

```
model = BayesianModel([('Age', 'Pclass'), ('Embarked', 'Pclass'),('Sex', 'Parch'), ('SibSp', 'Parch')])  
print(BicScore(df).score(model))

# 実行結果
# -4517.5517518749075
```
先程のモデルは約-4936.9だったので、今回の方がBICの絶対値が小さいです。
なので、**「Survived」への矢印はあった方が良い**ことがわかります。

#4. 最後に
今回、Titanicデータでベイジアンネットワークを実装しました。
相当簡単に実装できる関数も存在しますし、BICによる評価も簡単にできるため、
**関数により作成したモデルに業務観点で手動で依存関係(矢印)を入れてみて評価してみる**、みたいなこともできます。
是非実施してみてください。


#参考文献
[1] Spirtes, P., Glymour, C., & Scheines, R. (2000). Causation, Prediction, and Search
[2] 「統計WEB 25-5. 独立性の検定」 https://bellcurve.jp/statistics/course/9496.html
[3] Tom S. Verma, Judea Pearl."An Algorithm for Deciding if a Set of Observed Independencies Has a Causal Explanation" https://arxiv.org/abs/1303.5435 (2013)
[4] 「作りながら学ぶ！Pythonによる因果分析」(小川 雄太郎 著 2020)
[5] 「効果検証入門」(安井翔太 著 2020)