因果分析に興味があるということでこちらの書籍を勉強中。
book.mynavi.jp
こちらの書籍で紹介されているベイジアンネットワークによる因果探索の手法を整理していきます。また以下の記事も参考にさせていただきました。
・統計的因果探索とAI
・グラフィカルモデルに基づく因果探索手法の調査
目次
因果探索とベイジアンネットワーク
因果探索とは、データ項目のつながり(因果グラフ)を推測する手法になります。
因果探索の手法としてベイジアンネットワークがあり、基本的には
- 有向非巡回グラフ(矢印が巡回せず一方通行)である(DAG: Directed Acyclic Graphと呼ばれる)
- 未観測共通原因（潜在的な交絡因子）が無いこと
というように因果モデルに仮定を置いて分析を進めるものになります。
ベイジアンネットワークによる因果探索の手法は以下に示すよう大きく3つあります。
[スコアリングによる構造学習]
さまざまなDAGに対してスコア(※)を求め、最も指標値の良いDAGを選ぶものです。
つまりいろんなパターンのつながりをもったネットワークを考えて(規定して)、それぞれのスコアを算出して比較する手法になります。
(※)観測データに対するネットワークの当てはまりの良さを示すスコア：AIC, BICなど
[条件付き独立性検定による構造学習]
変数間で独立性検定を各条件付きパターンで行い、スケルトン構造を同定します。その後ルールベースによって部分的な方向づけを繰り返して、スケルトン構造->PDAG(Partially DAG)->DAGと構築していく手法になります。
[ハイブリッド型構造型学習]
スコアリングと条件付き独立性検定を組み合わせたもので、代表的な手法にMMHCアルゴリズム(Max-Min Hill Climbing algorithm)があります。
今回は、
- 条件付き独立性検定による構造学習の代表的なアルゴリズムであるPCアルゴリズム
- 条件付き確率表(CPT)とベイズ情報量基準(BIC)を用いたスコアリングによる構造学習
を実装していきます。
また実装していくにあたりpythonライブラリとして主にpgmpyを使っていますが、書籍で紹介されているバージョンが0.1.9でありエラーが出ますので、0.1.21で実装しています(2023.05.13時点最新バージョンは0.1.22)。またv0..で開発段階のため、今後も大きな修正は行われそうなので使用する際は注意となります。
データの説明
それでは、分析対象の擬似データの構造を示します。 今回も以前の因果推論の記事で用いた構造と同じもので行います。
職場での上司の研修受講が部下の面談満足度にどう影響するのか知りたいという場面を想定します。
双方とも上司の部下育成の熱心さと仕事のパフォーマンスにも影響を受けるとします。加えて、研修受講の効果は上司の部下育成の熱心さが高いと効果も高くなるとします。
( そこそこ妥当、そこそこ違和感というような場面設定ですが)
(効果の大きさまでは分析しません。この部分に関しては、他の記事でも紹介している因果推論の出番となります。)
生成コードは以下です。
# データ数 num_data = 2000 # x1 : 上司の部下育成の熱心さ。-1から1の一様乱数。 x1 = np.random.uniform(low=-1, high=1, size=num_data) # x2 : 上司の仕事のパフォーマンス。-1から1の一様乱数。 x2 = np.random.uniform(low=-1, high=1, size=num_data) # x2 = np.random.choice([-1, -0.5, 0, 0.5, 1], num_data, p=[0.1, 0.2, 0.1, 0.3, 0.3]) # Z : 上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか e_z = randn(num_data) # ノイズの生成 # シグモイド関数により受講するかどうかを熱心さに影響を受ける確率で表現(あくまで確率) z_prob = expit(7*x1 + 5*x2 + 5*e_z) Z = np.array([]) for i in range(num_data): Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0] # Z = np.append(Z, Z_i) # Y : 部下の面談の満足度 ## 研修の効果は熱心さで効果が変わってくるとする(介入効果の非線形性) t = np.zeros(num_data) for i in range(num_data): if x1[i] < 0: t[i] = 0.5 elif x1[i] >= 0 and x1[i] < 0.5: t[i] = 0.7 elif x1[i] >= 0.5: t[i] = 1.0 e_y = randn(num_data) Y1 = 2.0 + t*Z + 0.3*x1 + 0.5*x2 + 0.2*e_y # テーブルを作成 df = pd.DataFrame({'x1': x1, 'x2': x2, 'Z': Z, 't': t, 'Y1': Y1 }) df.head() # 先頭を表示
[output]
x1 x2 Z t Y1
0 -0.616961 0.731925 0.0 0.5 2.112008
1 0.244218 0.424114 1.0 0.7 3.067455
2 -0.124545 0.327542 0.0 0.5 2.022895
3 0.570717 -0.230410 1.0 1.0 3.148009
4 0.559952 -0.641392 0.0 1.0 2.032487
分析パターン①(PCアルゴリズム)
条件付き独立性検定による構造学習の代表的なアルゴリズムでもあり、制約ベースの因果探索とも位置付けられるものになります。
2. オリエンテーションルールで部分的に方向づけを繰り返してPDAG(Partially DAG)をつくる
3. 2を繰り返してDAGを推定
独立性検定とは統計的仮説検定の1つであり、背理法により"独立でない" or "独立でないとはいえない"を判定する手法です。 統計WEBの記事などが参考になるかと思います。 https://bellcurve.jp/statistics/course/9309.html https://bellcurve.jp/statistics/course/9496.html
今回は、"独立でない"を"関連があるので、因果関係にある"という様に解釈して分析を進めていく手法となります。 また今回は検定統計量として、以下の数式で表されるカイ二乗統計量を用います。
(: 度数分布表の行数、
: 度数分布表の列数、
:
行目
列目の度数、
:
行目
列目の推定期待度数)
また
(:
行目のデータ数、
:
行目のデータ数、
: 全データ数)
データの離散化(ビン分割)
PCアルゴリズムは独立性の検定を行うため連続値のデータは扱えないので、まず離散値に変換します。 今回は以下のように、pandasのcut関数で閾値での分割を行います。
df_bin = df.copy() del df_bin["t"] # x1:部下育成への熱心さ df_bin["x1"], x_bins = pd.cut(df["x1"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # x2:仕事のパフォーマンス df_bin["x2"], x_bins = pd.cut(df["x2"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # Z:上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか df_bin["Z"], Z_bins = pd.cut(df["Z"], 2, labels=["0", "1"], retbins=True) # Y1:部下の面談の満足度 df_bin["Y1"], Y_bins = pd.cut(df["Y1"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # データ型変更 df_bin["x1"] = df_bin["x1"].astype("int") df_bin["x2"] = df_bin["x2"].astype("int") df_bin["Z"] = df_bin["Z"].astype("int") df_bin["Y1"] = df_bin["Y1"].astype("int") # 確認 df_bin.head()
[output] x1 x2 Z Y1 0 1 5 0 2 1 4 4 1 4 2 3 4 0 2 3 4 2 1 4 4 4 1 0 2
離散化できました。また各ビンのデータ量も確認しておきます。
for col_name in df_bin.columns: print(f"カラム名：{col_name}") display(df_bin[col_name].value_counts()) print("=====")
カラム名：x1 1 452 5 414 3 410 4 381 2 343 Name: x1, dtype: int64 ===== カラム名：x2 4 432 3 409 1 400 5 393 2 366 Name: x2, dtype: int64 ===== カラム名：Z 1 1009 0 991 Name: Z, dtype: int64 ===== カラム名：Y1 2 693 3 545 4 469 1 189 5 104 Name: Y1, dtype: int64
閾値で分割したため、Y1ではデータ数に偏りが出ていますね。x1の場合分けによるtの値とrandn()の正規分布でのノイズの影響が出ている様です。
データの準備もできたので、PCアルゴリズムによる因果探索をしていきます。
ちなみに分析パターン①は以下のコードのみで結果が出ますが、中身を見ていくために今回は条件付き独立性検定とオリエンテーションルールの手順を一つずつ実行していきます。
# 以下のコードでPCアルゴリズムが実行できる from pgmpy.estimators import PC from pgmpy.base import DAG from pgmpy.independencies import Independencies est = PC(df_bin) skel, seperating_sets = est.build_skeleton(ci_test="chi_square", max_cond_vars=6, significance_level=0.05) pdag = est.skeleton_to_pdag(skel, seperating_sets) model = pdag.to_dag() print("DAG edges: ", model.edges())
条件付き独立性検定
1. 0次の独立性検定
まずは0次の独立性の検定です。0次なので、2つの変数間にある変数0個を条件付きにした際の検定ですのでそのまま行います。カイ二乗検定なので、以下のコードで一つづつ実行します。
# こちらは今回のパラメータ設定です。 from pgmpy.estimators.CITests import chi_square # 独立性の検定手法としてカイ二乗検定を選択 # パラメーターの設定 boolean_flg = True # Trueだと検定結果を返し、Falseだと検定統計量とp値を返す level_value = 0.05 # 有意水準
# 4変数なので6通り print(chi_square(X="x1", Y="Z", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print(chi_square(X="x1", Y="x2", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print(chi_square(X="x1", Y="Y1", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print(chi_square(X="Z", Y="x2", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print(chi_square(X="Z", Y="Y1", Z=[], data=df_bin, boolean=boolean_flg, significance_level=level_value)) print(chi_square(X="x2", Y="Y1", Z=[], data=df_bin, boolean=boolean_flg, significance_level=level_value))
[output] False True False False False False
x1 - x2間のみが独立でないとはいえないという結果なので、この部分だけエッジがない構造を以下の様に描くことができます。
2. 1次の独立性検定
引き続いて1次の独立性の検定です。1次なので、2つの変数間にある変数1個を条件付きにした際の検定です。また条件付きにする変数は検定対象の変数と関連があるものから選びます。
例えば変数x1についての1次の独立性の検定は、以下のコードで実行します。
# 変数x1について： 2変数(Z, Y1)と関連があったので、1*2=2通り print(chi_square(X="x1", Y="Z", Z=["Y1"], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print(chi_square(X="x1", Y="Y1", Z=["Z"], data=df_bin, boolean= boolean_flg, significance_level=level_value))
[output] True False
x1とZが独立でないとはいえないという結果なので、x1 - Zのエッジが消去されて以下の様に構造が描けます。
省略しますが他の変数についても検定をしていきます。結果としては、どの検定でも独立とは判定されなかったので1次の独立性検定の結果としては上記になります。
3. 2次の独立性検定
続いて2次の独立性の検定です。1次の独立性検定と同様の説明ですが、2つの変数間にある変数2個を条件付きにした際の検定です。また条件付きにする変数は検定対象の変数と関連があるものから選びます。
2次ということから3変数以上と関連がある変数についてしか検定できないので、検定対象はY1のみとなります。以下のコードで実行します。
# 変数Y1について: 4変数(x1, Z, x2, Y2)なので、3*4=12通り print(chi_square(X="Y1", Y="x1", Z=["Z", "x2"], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print("=====") print(chi_square(X="Y1", Y="Z", Z=["x1", "x2"], data=df_bin, boolean= boolean_flg, significance_level=level_value)) print("=====") print(chi_square(X="Y1", Y="x2", Z=["x1", "Z"], data=df_bin, boolean= boolean_flg, significance_level=level_value))
[output] False ===== False ===== False
いずれも独立でないという結果なので、2次の独立性検定で消去されるエッジはありませんでした。
オリエンテーションルールによる方向づけ
独立性の検定で得られた構造がDAGのベースとなるスケルトンが得られました。ここから方向づけのためのオリエンテーションフェーズです。
1. V字構造に着目してスケルトンからPDAGへ変換
まずはスケルトンをPDAGへ変換していくために、はじめにV字構造の部分に対して方向づけを行います。V字構造では下記の図の様に4パターンが考えられます。因果関係を絞るためにCI(A, C | B)を検定しますが、検定結果が独立でないと判定された場合のみ、パターン1の因果構造とできます。
スケルトンからx1-Y1-Zとx1-Y1-x2がV字構造と判断できるので、これらについて独立性の検定をしていきます。 以下のコードで実行します。
# v字構造での方向づけ # x1-Y1-Z print(chi_square(X="x1", Y="Z", Z=["Y1"], data=df_bin, boolean= boolean_flg, significance_level=level_value)) # x1-Y1-x2 print(chi_square(X="x1", Y="x2", Z=["Y1"], data=df_bin, boolean= boolean_flg, significance_level=level_value))
[output] True False
x1-Y1-Zは独立でないという結果なので、因果関係は以下のようになります。
2. PDAGからの部分的な方向づけ
できたPDAGに対して、以下の図を参考にさらにオリエンテーションルールに従った方向づけを行います。(成り立つ理由は理解できていませんが、PCアルゴリズムの原著に示されている様です。)
上記のパターンを適用すると、Y1->Zの方向をつけることができ、さらにそこからx2->Zも方向付けることができます。
これで因果構造を描くことは出来ました。しかし真の構造とはかなり違ったものとなり、良い分析はできておりません。そもそも今回の場面設定ではZ: 研修の受講がY1: 部下の面談の満足度にどの程度影響しているかを知りたいというのが因果分析を行うモチベーションとなるはずです。ですので、Y1 -> Zという因果構造になったことからも分析結果に対して懸念を持つことになるかと思います。 そこで分析の方法を変更してみます。
分析パターン②(条件付き独立性検定 -> スコアリングによる構造学習)
それでは分析の改善のため、データの離散化の方法を変更します。
分析パターン①では、離散化を閾値で等間隔になる様に分割を行なったため、ビン間にデータ量の偏りがありました。PCアルゴリズムは独立性の検定により、エッジの有無判定、スケルトンからPDAGへの変換を行うので離散化の方法による影響、つまりデータの偏りや量は注目するべき点になります。そこで分析パターン②では、データ量が偏らない様に分割していきます。
データの生成は同じなので、説明は割愛します。
ビンの分割方法の変更
以下のコードで実行します。pandasのqcut関数を用います。
df_bin = df.copy() del df_bin["t"] # x1:部下育成への熱心さ df_bin["x1"], x_bins = pd.qcut(df["x1"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # x2:仕事のパフォーマンス df_bin["x2"], x_bins = pd.qcut(df["x2"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # Z:上司が「上司向け：部下とのキャリア面談のポイント研修」に参加したかどうか df_bin["Z"], Z_bins = pd.cut(df["Z"], 2, labels=["0", "1"], retbins=True) # Y:部下の面談の満足度 df_bin["Y1"], Y_bins = pd.qcut(df["Y1"], 5, labels=[1, 2, 3, 4, 5], retbins=True) # データ型変更。 df_bin["x1"] = df_bin["x1"].astype("int") df_bin["x2"] = df_bin["x2"].astype("int") df_bin["Z"] = df_bin["Z"].astype("int") df_bin["Y1"] = df_bin["Y1"].astype("int") # 確認 df_bin.head()
[output] x1 x2 Z Y1 0 2 5 0 3 1 4 4 1 5 2 3 4 0 2 3 4 2 1 5 4 4 1 0 2
for col_name in df_bin.columns: print(f"カラム名：{col_name}") display(df_bin[col_name].value_counts()) print("=====")
[output] カラム名：x1 2 400 4 400 3 400 5 400 1 400 Name: x1, dtype: int64 ===== カラム名：x2 5 400 4 400 2 400 1 400 3 400 Name: x2, dtype: int64 ===== カラム名：Z 1 1009 0 991 Name: Z, dtype: int64 ===== カラム名：Y1 3 400 5 400 2 400 4 400 1 400 Name: Y1, dtype: int64 =====
偏りなくビン分割できています。それではこちらのデータで同様に条件付き独立性の検定を行なっています。
条件付き独立性検定
こちらも分析パターン①と同様に、0次、1次、2次と条件付き独立性の検定を行なっていくだけなので割愛します。
ただ得られた構造は分析パターン①と異なり、下記のように真の構造と同じエッジが残りました。
結果の流れとしては、0次で得られた構造からエッジが無くなることなく2次までの検定が完了しました。
ただこの構造ではV字構造が無いためPDAGへの変換ができず、オリエンテーションルールが適用できません。そこでスコアリングによる構造学習で方向づけを考えていきます。
スコアリングによる構造学習
スコアリングによる構造学習は、
という手法です。今回はBICを指標とします。
手順としては、
2. 観測データから条件付き確率表(CPT)を推定する
3. 得られた推定したCPTから各変数の各サンプルの確率値とパラメーター数が得られる
4. 各モデルのBICが求められるので比較して、最も良好な値のモデルを採用する
となります。
はじめにモデルを規定して指標を比較するので、変数が多くなるとDAGのパターンが爆発的に増加してしまうのが欠点になります(そのため厳密性を犠牲にして計算量を減らす手法もいくつか提案されている様です)。
BICについて
BICが以下の式で計算されます。
mはモデルを表しており、対象としているベイジアンネットワークのDAGになります。は観測されたデータ、
はモデルのパラメータ数で自由度に対応する様な値、
は観測されたデータ数になります。
またはそのモデルと規定した場合の確率値から求められる対数尤度であり、比例する値として以下の様にも計算されます。
（は変数、
は変数
についてのとある条件パターン、
は変数
のとる値）
負の対数尤度に、過剰な適合を考慮するためのペナルティとしてパラメータ数とデータ数を掛け合わせたものを足しているので、値が小さいほど良いモデルと考える指標になります。
それでは分析を進めるために、まず因果構造を考えます。今回はZ: 研修受講がY1: 部下の面談の満足度にどう影響するかを分析することが目標になります。ですので、Z -> Y1という因果方向を前提としていくつかのモデルを考えます。DAGであることから循環しないモデルを考えると、以下の9パターンが考えられます。
以下のコードで各モデル(因果構造)におけるBICを計算します。pgmpyでのBICの計算値は上記で示したものに-0.5をかけたものになっているので、値が大きいほど良いモデルと考えます。
# スコアリング構造学習で判断する from pgmpy.models import BayesianNetwork from pgmpy.estimators import BicScore model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('x1', 'Y1'), ('Z', 'x2'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('x1', 'Y1'), ('Z', 'x2'), ('Y1', 'x2')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('x2', 'Z'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('Z', 'x2'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('Z', 'x2'), ('Y1', 'x2')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('x2', 'Z'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('Z', 'x2'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('Z', 'x2'), ('Y1', 'x2')]) print(BicScore(df_bin).score(model)) model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('x2', 'Z'), ('x2', 'Y1')]) print(BicScore(df_bin).score(model))
[output] -9652.730917876488 -9468.571910366429 -9468.571910366429 -9468.571910366429 -9468.571910366429 -9468.571910366429 -9652.730917876488 -9468.571910366429 -9675.163065400724
出力結果より、真のモデルは一番値が小さく分析結果としては一番当てはまらないモデルという結果になってしまいました。。。。。
おわりに
今回は因果探索への入門として、ベイジアンネットワークを用いた手法を実装していきました。 手法の原理もある程度とばしながらですが一通りの実装はできたと思います。データの分割方法の変更、PCアルゴリズムとスコアリングによる構造学習との併用も実験的に行えて勉強になりました。ただ4変数しかないにもかかわらず真の構造を推定できず、因果探索自体が難しいタスクという感触でした。データからの結果だけで無く、ドメイン知識による因果関係の推定も重要になってきそうです。また手法などは色々と提案されている様ですが、データの分割方法などの前処理の方法論は特に無いようですので、試しながら肌感覚を養っていけたらと思います。
(追記 注意点)
便利に思われるPCアルゴリズムですが、検定する変数の順番が異なればエッジが変わることもあります。これについては私も実験して確認し、著者の方に問い合わせたところ(以下のリンク)事実の様です。またこれに対して、
どうするのが正解なのか、きちんと解説している文献を見たことがありません。
とのことです。これは非常に大事な注意点かと思います。
https://github.com/YutaroOgawa/causal_book/issues/39
