グラフィカルモデリングさわってみた(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()
動画時間_s視聴回数コメント高評価件数低評価件数
0108271418690
11025153314980
214891219231110
35302118281642
412831462181121

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

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
動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.0000-0.0003-0.25860.00010.0000
視聴回数-0.00031.0000-0.3939-0.7903-0.6753
コメント-0.2586-0.39391.00000.00000.0000
高評価件数0.0001-0.79030.00001.00000.4176
低評価件数0.0000-0.67530.00000.41761.0000

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
動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.34140.45210.41010.1251
視聴回数0.34141.00000.84860.91040.8311
コメント0.45210.84861.00000.83220.6137
高評価件数0.41010.91040.83221.00000.6606
低評価件数0.12510.83110.61370.66061.0000

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

new_inv_cor = pd.DataFrame(inverse_cor(new_sample_cor),
                                       index=sample_cor.index,
                                       columns=sample_cor.columns)
new_inv_cor
動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.0000-0.2375-0.11120.1729
視聴回数-0.00001.0000-0.4131-0.6980-0.7369
コメント-0.2375-0.41311.0000-0.13300.1751
高評価件数-0.1112-0.6980-0.13301.00000.3369
低評価件数0.1729-0.73690.17510.33691.0000

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)
結果

========== cycle 1 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.0577-0.2602-0.15090.1315
視聴回数0.05771.0000-0.4250-0.7019-0.7329
コメント-0.2602-0.42501.0000-0.11840.1821
高評価件数-0.1509-0.7019-0.11841.00000.3367
低評価件数0.1315-0.73290.18210.33671.0000

独立を仮定する変数ペア: 視聴回数 | 動画時間_s
偏相関係数の絶対値: 0.05765421774048792
逸脱度(vs 元の相関行列): 0.01665
p値(vs 元の相関行列): 0.89734
逸脱度(vs 1つ前の相関行列): 0.01665
p値(vs 1つ前の相関行列): 0.89734
=> 制約モデルを採用
========== cycle 2 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.0000-0.2375-0.11120.1729
視聴回数-0.00001.0000-0.4131-0.6980-0.7369
コメント-0.2375-0.41311.0000-0.13300.1751
高評価件数-0.1112-0.6980-0.13301.00000.3369
低評価件数0.1729-0.73690.17510.33691.0000

独立を仮定する変数ペア: 高評価件数 | 動画時間_s
偏相関係数の絶対値: 0.11123492706884651
逸脱度(vs 元の相関行列): 0.13878
p値(vs 元の相関行列): 0.93296
逸脱度(vs 1つ前の相関行列): 0.12213
p値(vs 1つ前の相関行列): 0.72673
=> 制約モデルを採用
========== cycle 3 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.0000-0.0001-0.32420.00000.1370
視聴回数-0.00011.0000-0.4023-0.7024-0.7412
コメント-0.3242-0.40231.0000-0.15620.1674
高評価件数-0.0000-0.7024-0.15621.00000.3604
低評価件数0.1370-0.74120.16740.36041.0000

独立を仮定する変数ペア: 低評価件数 | 動画時間_s
偏相関係数の絶対値: 0.13698168568594546
逸脱度(vs 元の相関行列): 0.3785
p値(vs 元の相関行列): 0.94465
逸脱度(vs 1つ前の相関行列): 0.23972
p値(vs 1つ前の相関行列): 0.62441
=> 制約モデルを採用
========== cycle 4 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.0002-0.2440-0.00010.0000
視聴回数0.00021.0000-0.4125-0.7024-0.7482
コメント-0.2440-0.41251.0000-0.16020.2192
高評価件数-0.0001-0.7024-0.16021.00000.3639
低評価件数-0.0000-0.74820.21920.36391.0000

独立を仮定する変数ペア: コメント | 高評価件数
偏相関係数の絶対値: 0.16015368761151882
逸脱度(vs 元の相関行列): 0.51679
p値(vs 元の相関行列): 0.97185
逸脱度(vs 1つ前の相関行列): 0.13829
p値(vs 1つ前の相関行列): 0.70999
=> 制約モデルを採用
========== cycle 5 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.0000-0.0003-0.24660.00010.0000
視聴回数-0.00031.0000-0.5025-0.7401-0.7541
コメント-0.2466-0.50251.00000.00000.2766
高評価件数0.0001-0.74010.00001.00000.4002
低評価件数0.0000-0.75410.27660.40021.0000

独立を仮定する変数ペア: 動画時間_s | コメント
偏相関係数の絶対値: 0.24656561863024512
ある変数とその他の全ての変数が条件付き独立になりそうなので、仮定を取り消す
========== cycle 6 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.0000-0.0003-0.24660.00010.0000
視聴回数-0.00031.0000-0.5025-0.7401-0.7541
コメント-0.2466-0.50251.00000.00000.2766
高評価件数0.0001-0.74010.00001.00000.4002
低評価件数0.0000-0.75410.27660.40021.0000

独立を仮定する変数ペア: コメント | 低評価件数
偏相関係数の絶対値: 0.2765835272303214
逸脱度(vs 元の相関行列): 1.02683
p値(vs 元の相関行列): 0.96037
逸脱度(vs 1つ前の相関行列): 0.51004
p値(vs 1つ前の相関行列): 0.47512
=> 制約モデルを採用
========== cycle 7 ==========

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.0000-0.0003-0.25860.00010.0000
視聴回数-0.00031.0000-0.3939-0.7903-0.6753
コメント-0.2586-0.39391.00000.00000.0000
高評価件数0.0001-0.79030.00001.00000.4176
低評価件数0.0000-0.67530.00000.41761.0000

最終的な偏相関係数行列

動画時間_s視聴回数コメント高評価件数低評価件数
動画時間_s1.00000.00030.2586-0.0001-0.0000
視聴回数0.00031.00000.39390.79030.6753
コメント0.25860.39391.0000-0.0000-0.0000
高評価件数-0.00010.7903-0.00001.0000-0.4176
低評価件数-0.00000.6753-0.0000-0.41761.0000

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

参考: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