グラフィカルモデリングさわってみた（python）

本記事は無向グラフィカルモデリングに関して、以下の参考HP（主に1）の内容をまとめた。個人的メモです
※言葉の定義や内容が厳密でない場合があります。定義が間違っている場合はご指摘いただきたいです。

1. 【偽物の相関を見極めろ！】グラフィカルモデリングで変数の相関関係を把握する
2. グラフィカルモデリングの紹介
3. データから因果を推定する
4. グラフィカルモデリングをまとめてみた
5. グラフィカルモデリングの基本を少しまとめました

グラフィカルモデリングについて

グラフィカルモデリングとは

確率変数間の依存関係をグラフで表したもの。
（確率変数を頂点、それらの間の依存関係を辺としたグラフを用いる）

「ほとんどの統計モデルをこのモデルの元で統一的に表せる」とのこと。
（回帰分析、因子分析、SEM、信号検出理論、隠れマルコフモデル、パス解析 etc...）

大きくは以下の二種類に分けられる。

• 無向グラフによるグラフィカルモデリング：マルコフ確率場
• 有向グラフによるグラフィカルモデリング：ベイジアンネットワーク

（無向グラフによるグラフィカルモデリングの方が単純なモデル。今回扱うのは無向グラフ）

作成されたグラフは独立グラフと呼ばれ、上記それぞれの場合の独立グラフは無向独立グラフ、有向独立グラフと呼ばれる。

（無向グラフによるグラフィカルモデリングの場合の表し方）
確率変数間に依存関係がある場合（条件付き独立でない時）のみ、その頂点を辺で結ぶ。

例）AとBは条件付き独立であり、AとC, BとCが条件付き独立ではない場合。

条件付き独立について

独立関係とは

一方の事象が起こったことが分かっても、他方の事象の確率が変化しないこと。
（2つの事象がいずれも起きる確率が、それぞれの確率の積に等しいこと）
つまり...

$$P(A,B) = P(A)P(B)$$

$$P(A|B) = P(A)$$

$$P(B|A) = P(A)$$

の時、AとBは独立といい。

$$A \perp B$$

と表す。

条件つき独立とは

$A$と$B$が条件付き独立とは、他の全ての変数を固定した時に独立になること。
（交絡因子などの影響を除いた場合に独立になるということ）
逆に他の変数全ての変数を固定したときに$A$と$B$に相関関係が残る場合、条件付き独立は成立しない。

$$P(A|B) = P(A|M,B)$$

の時、

AとBはMの元で条件付き独立といい、

$$A \perp B | M$$

と表す。

またこの時、因果の向きは次の3つのうちどれかに絞られる。

$$A \rightarrow M \rightarrow B$$

$$A \leftarrow M \leftarrow B$$

$$A \leftarrow M \rightarrow B$$

条件付き独立かどうかの判定

$x_1 , \cdots , x_n$ が多変量正規分布に従う時、
$\rho_{i,j,rest}=0 $ ならば、$x_i$ と $x_j$ は条件付き独立である。

（$\rho_{i,j,rest}=0 $ は、$x_i$ と$x_j$ の偏相関係数）

マルコフ性

無向独立グラフにおいて、二つの確率変数をつなぐ確率変数を固定した時、両者の間は独立になる。
言い換えると、固定する変数を独立グラフから消去した時、辺で結ばれていない変数（変数のグループ）間は独立であり、辺で繋がっている変数同士には相関が残る。
これを大域マルコフ性（Global Markov Property）という

無向き独立グラフの作成手順

1. 偏相関係数行列を計算
2. 偏相関が小さい変数ペアを条件付き独立と仮定
3. 2のペアの偏相関係数を0とした制約モデルの偏相関行列を計算
4. 検定により3の制約モデルを採択するか判断する。
5. 4がOKの場合strong>：制約モデルを採択し、別の変数ペアの条件付き独立性を仮定する（3に戻る）
   4がNGの場合：制約モデルを棄却し、別の変数ペアの条件付き独立性を仮定する（3に戻る）

（4がNGだった場合）適当なタイミングで1つ前の条件で採択して完了とする。

検定ついて

ある変数ペアを条件付き独立とみなすことで、偏相関行列の計算結果が大きく変わっていないか逸脱度を用いてチェックする。

標本数を$n$、制約なしモデルと制約ありモデル（偏相関係数を0とみなしたモデル）の分散共分散行列をそれぞれ$S_0$、$S_1$とすると逸脱度$d$は、

$$d = n\ln\frac{|S_1|}{S_0}$$

と計算される。
この時、モデル間の逸脱度$d$は、制約の数を自由度をとしたカイ二乗分布に従う。

例） $d=0.052$となるp値が0.819なら、
制約を加えたモデル（ある変数ペアの偏相関係数が0）が正しい時、$d$が0.052以上になるのは、0.819 以上という意味なので、制約モデルを採用する。
（このように、偏相関係数の値が小さい変数間について、その変数間が条件付き独立であると仮定した関連構造モデルを採用する方法を共分散選択という）

コード

AIciaSolidProjectのコード内のデータ、コードを使わせていただいてます
（参考動画内でも使用している、AlciaSolidProjectチャンネルにおける動画の情報（動画時間、視聴回数など）のデータです）

コード内でのモデル選択手順

1. 相関行列の逆行列の標準化したものを計算し、（絶対）値最小のものを探す。
   　（→「偏相関係数が最も小さい変数ペアを探す」と同義）
2. その2つの変数は、条件付き独立であると仮定して、相関行列を推定する。
3. 推定した相関行列を用いて逸脱度とp値を算出。
4. p値が大きければ、上記の手順を続行
5. p値が小さければ終了。

（ある変数とその他変数が全て条件付き独立になるのはあまり良くなく、そうならないように変数ペアを探すほうが良いらしい）

import numpy as np
import os
import pandas as pd
from scipy.stats import chi2

pd.options.display.float_format = '{:.4f}'.format

df = pd.read_csv('AIcia_videos.csv')
df['動画時間_s'] = pd.to_timedelta(df['動画時間']).apply(lambda x: x.seconds)
df = df[['動画時間_s', '視聴回数', 'コメント', '高評価件数', '低評価件数']]

df.head()

一手順ごとに実行してみる

1. 相関行列の逆行列の標準化したものを計算し（絶対）値最小のものを探す。

# 相関行列を計算
sample_cor = df.corr()

# 相関行列の逆行列を標準化したものを計算
def inverse_cor(cor):
    cor_inv = np.linalg.inv(cor)    
    diag = np.diag(1/np.sqrt(np.diag(cor_inv)))
    return diag.dot(cor_inv).dot(diag)

inv_sample_cor = pd.DataFrame(inverse_cor(estimated_cors[-1]), 
                              index=sample_cor.index, 
                              columns=sample_cor.columns)

inv_sample_cor

2. 選んだ2つの変数を条件付き独立であると仮定して、相関行列を推定する。

動画時間_s と 視聴回数 の相関（相関行列の逆行列を標準化したもの）が絶対値最小（0.2586）なので、これらが条件付き独立だと仮定して、相関行列を再計算する。

cond_ind_pairs = [(0, 1)]
estimated_cor = np.array(sample_cor).copy()

error_torelance=1e-4

# 偏相関係数が0になるように相関係数を計算する。
error=1
while error > error_torelance:

    estimated_cor_before = estimated_cor.copy()
    for i, j in cond_ind_pairs:
        # 指定した変数ペアに関して偏相関係数が0になるように相関係数を計算する。
        estimated_cor_inv = np.linalg.inv(estimated_cor)

        new_ij = estimated_cor[i][j] + estimated_cor_inv[i][j]/(estimated_cor_inv[i][i] * estimated_cor_inv[j][j] - estimated_cor_inv[i][j]**2)
        estimated_cor[i][j] = new_ij
        estimated_cor[j][i] = new_ij

    error = np.abs(estimated_cor - estimated_cor_before).max().max()

new_sample_cor = pd.DataFrame(estimated_cor, index=sample_cor.index,columns=sample_cor.columns)
new_sample_cor

偏相関係数が0になっているかチェックしてみる

new_inv_cor = pd.DataFrame(inverse_cor(new_sample_cor),
                                       index=sample_cor.index,
                                       columns=sample_cor.columns)
new_inv_cor

3. 推定した相関行列を用いて逸脱度とp-値を見る

new_sample_cor_array = np.array(new_sample_cor)
sample_cor_array = np.array(sample_cor)

# dim => サンプル数
dim = sample_cor_array.shape[0]

# 逸脱度を計算
dev = dim * (np.log(np.linalg.det(new_sample_cor_array)) - np.log(np.linalg.det(sample_cor_array)))

# p値を計算
p = 1 - chi2.cdf(dev, df=1)

print('逸脱度:',dev)
print('p値:',p)

逸脱度: 0.01664772806781567
p値: 0.8973370162960439

4. p-値が大きければ、上記の手順を続行

p-値が充分大きいので、仮定（動画時間_s と 視聴回数が条件付き独立）を採択する。

一気通貫で実行

※本来は、データの中身と1サイクルごと結果を加味しながら正しく、解析できているかチェックするべきとのこと。

関数定義

# 相関行列の（標準化された）逆行列を計算する関数
def inverse_cor(cor):
    cor_inv = np.linalg.inv(cor)    
    diag = np.diag(1/np.sqrt(np.diag(cor_inv)))
    return diag.dot(cor_inv).dot(diag)

# 指定した変数ペアの偏相関係数が0になるような相関行列を算出する関数
def estimate(sample_cor, cond_ind_pairs, error_torelance=1e-4):
    # cond_ind_pairs：相関係数を0と仮定した変数ペアのリスト
    # sample_cor：相関係数
    sample_cor = np.array(sample_cor)    
    estimated_cor = sample_cor.copy()    

    error = 1
    while error > error_torelance:

        estimated_cor_before = estimated_cor.copy()
        for i, j in cond_ind_pairs:
            estimated_cor_inv = np.linalg.inv(estimated_cor)
            new_ij = estimated_cor[i][j] + estimated_cor_inv[i][j]/(estimated_cor_inv[i][i] * estimated_cor_inv[j][j] - estimated_cor_inv[i][j]**2)
            estimated_cor[i][j] = new_ij
            estimated_cor[j][i] = new_ij

        error = np.abs(estimated_cor - estimated_cor_before).max().max()

    return estimated_cor

# 逸脱度とp値を計算する関数
def deviance_and_p(original_cor, estimated_cor, df):
    # original_cor：比較する相関行列
    # estimated_cor：新たに制約を設けた相関行列
    dim = original_cor.shape[0]
    dev = dim * (np.log(np.linalg.det(estimated_cor)) - np.log(np.linalg.det(original_cor)))
    p = 1 - chi2.cdf(dev, df)

    return dev, p

実行

# 相関行列を格納するリスト
estimated_cors = []

# 条件付き独立を仮定する変数ペアを格納するリスト
cond_ind_pairs = []

# 条件付き独立の仮定が不採用になった変数ペアを格納するリスト
not_cond_ind_pairs = []

# 有意水準を設定
alpha = 0.05

# 偏相関係数の閾値を設定（この値より偏相関係数が小さい変数ペアに関して、条件付き独立を検討する）
corr_threshold = 0.3

estimated_cors.append(df.corr())
variable_names = df.columns.tolist()

cycle = 1
while True:
    print(f'========== cycle {cycle} ==========')
    inv_corr_temp = pd.DataFrame(inverse_cor(estimated_cors[-1]), 
                                 index=variable_names, 
                                 columns=variable_names)
    display(inv_corr_temp)

    abs_inv_corr = abs(inv_corr_temp).values
    for pair in cond_ind_pairs + not_cond_ind_pairs:
        abs_inv_corr[pair[0],pair[1]] = 1
        abs_inv_corr[pair[1],pair[0]] = 1

    # corr_thresholdよりも偏相関係数の高い変数ペアしかないのであれば計算終了とする
    if abs_inv_corr.min() > corr_threshold:
        break

    ind_index = np.where(abs_inv_corr == abs_inv_corr.min())
    ind_index = (ind_index[0][0],ind_index[1][0])

    print('独立を仮定する変数ペア:',variable_names[ind_index[0]],'|',variable_names[ind_index[1]])
    print('偏相関係数の絶対値:',abs(inv_corr_temp.values)[ind_index[0]][ind_index[1]])

    # ある変数とその他の全ての変数が条件付き独立にならないようにする
    abs_inv_corr[ind_index[0],ind_index[1]] = 1
    abs_inv_corr[ind_index[1],ind_index[0]] = 1

    if (abs_inv_corr[ind_index[0],:].sum()==abs_inv_corr.shape[0]).any():
        print('ある変数とその他の全ての変数が条件付き独立になりそうなので、仮定を取り消す')
        not_cond_ind_pairs.append(ind_index)
        cycle += 1
        continue

    # 相関行列を計算する
    estimated_cor = pd.DataFrame(estimate(estimated_cors[-1], cond_ind_pairs+[ind_index]),
                                 index=variable_names, 
                                 columns=variable_names)

    # 元々の相関行列と新たに算出した相関行列の逸脱度とp値
    dev_original, p_original = deviance_and_p(estimated_cors[0].values, estimated_cor.values, 
                                              df=len(cond_ind_pairs + [ind_index]))
    print('逸脱度（vs 元の相関行列）:', round(dev_original,5))
    print('p値（vs 元の相関行列）:', round(p_original,5))

    # ひとつ前に算出した相関行列と新たに算出した相関行列の逸脱度とp値
    dev_prev, p_prev = deviance_and_p(estimated_cors[-1].values, estimated_cor.values, df=1)
    print('逸脱度（vs 1つ前の相関行列）:', round(dev_prev,5))
    print('p値（vs 1つ前の相関行列）:', round(p_prev,5))

    p = min(p_prev,p_original)
    if p > alpha:
        cond_ind_pairs.append(ind_index)
        estimated_cors.append(estimated_cor)
        print('=> 制約モデルを採用')
    else:
        not_cond_ind_pairs.append(ind_index)
        print('=> 制約モデルを棄却')

    cycle += 1

# 偏相関係数行列
inv_corr_temp = pd.DataFrame(inverse_cor(estimated_cors[-1]), 
                             index=variable_names, 
                             columns=variable_names)

partial_corr = inv_corr_temp.applymap(lambda x: -x if x < 1-1e-4 else 1)

print('最終的な偏相関係数行列')
display(partial_corr)

最終的な偏相関係数行列

無向独立グラフを書いてみる

参考：https://programgenjin.hatenablog.com/entry/2019/02/26/075121

※graphvizをインストールしておく必要あり。

pip install graphviz

import itertools
from graphviz import Graph

model_graph = Graph(format='png')

# 無向グラフ
for variable in variable_names:
    model_graph.node(variable)

# 全ての変数の組み合わせのリスト
comb_list = list(itertools.combinations([i for i in range(len(df.columns))], 2))

# 変数組み合わせの順番が逆にしたものも追加
cond_ind_pairs += [(pair[1],pair[0]) for pair in cond_ind_pairs]

for pair in cond_ind_pairs:
    try:
        comb_list.remove(pair)
    except:
        pass

for pair in comb_list:
    model_graph.edge(variable_names[pair[0]], variable_names[pair[1]])

model_graph.render('./graph', view=True)
model_graph

参考

参考HP

1. 【偽物の相関を見極めろ！】グラフィカルモデリングで変数の相関関係を把握する
2. グラフィカルモデリングの紹介
3. データから因果を推定する
4. グラフィカルモデリングをまとめてみた
5. グラフィカルモデリングの基本を少しまとめました

関連書籍

Udemy

【PythonとStanで学ぶ】仕組みが分かるベイズ統計学入門

Pythonで実践する統計モデリング入門

PyMC3 で学ぶ統計モデリングの基礎