ベイジアンネットワークで変数間の関連性を見る
実務でデータ分析を行うときにしばしば変数間の因果関係の有無およびその強さについて問われることがあるのですが、その回答としてベイジアンネットワークを使いたいと思っていたので調べてみました。
この記事ではベイジアンネットワークを扱うためのパッケージ （{bnlearn}および{BNSL}）を用いた構造の推定およびそのTipsを示します。
分析環境の準備
ライブラリの読み込み
まずは必要なライブラリを読み込みます。今回はRにおけるベイジアンネットワークのメジャーなライブラリである{bnlearn}に加え、{BNSL}を用います。
library(bnlearn) library(BNSL)
サンプルデータの作成
続いてこちらの記事を参考にサンプルデータを作成します。一部の係数を修正しましたが、記事同様にHeightとSBPがBMIに刺さり、BMIがFBSに刺さる構造となっています。
set.seed(123) norm <- rnorm(4000) Data <- matrix(norm, nrow = 1000, ncol = 4, byrow = T) colnames(Data) <- c("Height", "BMI", "SBP", "FBS") Data2 <- Data <- as.data.frame(Data) Data2$Height <- 170 + Data$Height * 10 Data2$SBP <- 120 + Data$SBP * 10 Data2$BMI <- 22 + Data$Height * 5 + Data$SBP * 5 Data2$FBS <- 90 + (Data2$BMI - 22) * 5 + Data$FBS * 10
以下のようなデータが得られます。
> head(Data2) Height BMI SBP FBS 1 164.3952 26.99116 135.5871 115.66090 2 171.2929 24.95102 124.6092 92.10449 3 163.1315 24.68614 132.2408 107.02886 4 174.0077 21.22465 114.4416 103.99239 5 174.9785 27.99603 127.0136 115.25225 6 159.3218 11.53086 109.7400 30.36538
分析実行その１
構造学習
それではベイジアンネットワークによる関連性の推定に移ります。ベイジアンネットワークはデータから変数間の関連性を推定する構造学習と，推定された構造を元にパラメータを推定するフィッティングの二段階で実施されます。
今回、構造学習のアルゴリズムにはGrow-Shrinkを指定しました。Webで調べた限りだとHill-Climbingが使われる事も多い様子ですが、先ほど参照した記事によるとGrow-Shrinkはグラフ理論による裏付けがあるとのこと。関数はgsです（Hill-Climbingはhc）。
str_01_gs <- gs(Data2) str_01_hc <- hc(Data2) par(mfrow = c(1, 2)) plot(str_01_gs) plot(str_01_hc)
gsでは設定した通りの構造が推定されましたが、hcではBMIとSBPの関係性が逆転しています。これ以降、基本的にgsを使うことにします。
フィッティング
ではこの構造を用いてパラメータを推定しましょう。なおここで得られるパラメータ推定値は個別にlmをかけたものと同じとなり、例えばBMIであればlm(BMI ~ Height + SBP)と等しい値となります。
fit_01 <- bn.fit(str_01_gs, Data2) coef_01 <- coefficients(fit_01) > coef_01 $Height (Intercept) 170.091 $BMI (Intercept) Height SBP -123.0 0.5 0.5 $SBP (Intercept) 120.1408 $FBS (Intercept) BMI -20.240828 5.019595
> coef(lm(BMI ~ Height + SBP, Data2)) (Intercept) Height SBP -123.0 0.5 0.5
同じ値となっていますね。
分析実行その２
構造学習
次にサンプルデータのgaussian.testを用いてパラメータを推定します。元の構造はわかりませんが、gsによって推定した形は以下のようです。
str_02 <- gs(gaussian.test) plot(str_02)
フィッティング
この構造を用いてパラメータの推定を行おうとするとエラーが出ます。どうやらグラフの構造が悪く、フィッティングできないようです。
> fit_02 <- bn.fit(str_02, gaussian.test) bn.fit(str_02, gaussian.test) でエラー: the graph is only partially directed.
そこでstr_02の中身を確認すると、modelが[partially directed graph]となっていることがわかります。またarcsのundirected arcsも1となっています。
> str_02 Bayesian network learned via Constraint-based methods model: [partially directed graph] nodes: 7 arcs: 7 undirected arcs: 1 directed arcs: 6 average markov blanket size: 4.00 average neighbourhood size: 2.00 average branching factor: 0.86 learning algorithm: Grow-Shrink conditional independence test: Pearson's Correlation alpha threshold: 0.05 tests used in the learning procedure: 144 optimized: FALSE
再度str_02の構造を確認するとノードBD間のエッジが無向になっていることがわかります。
そこでこの構造にB -> Dのエッジを設定してみましょう。ベイジアンネットワークではこのように、試行錯誤の結果や分析者の知見を取り込みながら、データ全体の構造を推定していきます。
この修正によりstr_03のundirected arcsが0となり、modelがpartially directed graphではなくノード間の関連性を表したものとなっていることを確認しましょう。
> (str_03 <- set.arc(str_02, "B", "D")) Bayesian network learned via Constraint-based methods model: [A][B][E][G][C|A:B][D|B][F|A:D:E:G] nodes: 7 arcs: 7 undirected arcs: 0 directed arcs: 7 average markov blanket size: 4.00 average neighbourhood size: 2.00 average branching factor: 1.00 learning algorithm: Grow-Shrink conditional independence test: Pearson's Correlation alpha threshold: 0.05 tests used in the learning procedure: 144 optimized: FALSE
ちゃんとB -> D間の矢印が有向になっています。
> plot(str_03)
では、この構造を用いて再度フィッティングしてみましょう。
> fit_03 <- bn.fit(str_03, gaussian.test) > (coef_03 <- coefficients(fit_03)) $A (Intercept) 1.007493 $B (Intercept) 2.039499 $C (Intercept) A B 2.001083 1.995901 1.999108 $D (Intercept) B 5.995036 1.498395 $E (Intercept) 3.493906 $F (Intercept) A D E G -0.006047321 1.994853041 1.005636909 1.002577002 1.494373265 $G (Intercept) 5.028076
今度は上手くいったようですね！
分析実行その３
続いて直接効果と間接効果の分離を試みます。再度gsで構造を推定し、B -> Dのエッジと共にB -> Fのエッジも追加しましょう。これでBはFに対して直接刺さるエッジと、Dを経由してFに刺さるエッジを持つことになります。
str_04 <- set.arc(str_03, "B", "F") par(mfrow = c(1, 3)) plot(str_02) # gsで推定した構造 plot(str_03) # B -> Dに修正 plot(str_04) # B -> Fを追加
フィッティング
続いてフィッティングです。B -> Fのエッジの有無で係数がどのように変化するかを比較してみましょう。必要な部分だけ抜粋します。
fit_04 <- bn.fit(str_04, gaussian.test) > coefficients(fit_04) ~ 省略 ~ $F (Intercept) A B D E G -0.4123157 1.9947919 -0.1025991 1.0737532 1.0023721 1.4943286 ~ 省略 ~
> coefficients(fit_03) ~ 省略 ~ $F (Intercept) A D E G -0.006047321 1.994853041 1.005636909 1.002577002 1.494373265 ~ 省略 ~
切片に多少の変化があっただけで、係数としては大きな変化はないようです。どうやらBの変動による寄与はほとんどDを経由してFに影響しており、BのFへの直接の効果はなさそうです。
分析実行その４
構造学習
今度は構造の推定に{BNSL}を試してみましょう。
str_05 <- bnsl(gaussian.test) plot(str_05)
どうもgsを用いた時と比べて構造がかなり異なるようですね。比較してみましょう。
par(mfrow = c(1, 2)) plot(str_02) plot(str_05)
"A -> C"や"G -> F"など、方向性が逆転する関係が見られる上、Eも独立してしまいました。やはりデータのみから因果関係を推定するのは不安定ということになるのでしょうか。
フィッティング
この構造でフィッティングしてみます。フィッティング自体はbnlearn::bn.fitを用いています（結果は省略）。
fit_05 <- bn.fit(str_05, gaussian.test)
分析実行その５
さらに続いて、データを標準化した場合の結果も確認してみます。scaleでデータを標準化します。
my_gaussian <- as.data.frame(scale(gaussian.test)) str_06 <- bnsl(my_gaussian) par(mfrow = c(1, 2)) plot(str_04) plot(str_06)
何と、標準化すると構造が変わってしまいました。。。フィッティングは省略します。
分析実行その６
最後に、ライブラリ{MASS}に格納されているサンプルデータbirthwtを用いて推定してみます。なおbnlearnではintegerは使えないので、各列をnumericに変換し、簡単のために変数も一部に限定します。
構造学習
最後なので、構造学習のアルゴリズムには色々なものを指定してみます。
vars <- c("low", "race", "smoke", "ht", "ui") my_birthwt <- as.data.frame((do.call("cbind", lapply(birthwt[, vars], as.numeric)))) str_07_01 <- hc(my_birthwt) str_07_02 <- gs(my_birthwt) str_07_03 <- iamb(my_birthwt) str_07_04 <- fast.iamb(my_birthwt) str_07_05 <- inter.iamb(my_birthwt) str_07_06 <- tabu(my_birthwt) str_07_07 <- mmhc(my_birthwt) str_07_08 <- rsmax2(my_birthwt) str_07_09 <- mmpc(my_birthwt) str_07_10 <- si.hiton.pc(my_birthwt) str_07_11 <- chow.liu(my_birthwt) str_07_12 <- aracne(my_birthwt) par(mfrow = c(1, 3)) plot(str_07_01) plot(str_07_02) plot(str_07_03) par(mfrow = c(1, 3)) plot(str_07_04) plot(str_07_05) plot(str_07_06) par(mfrow = c(1, 3)) plot(str_07_07) plot(str_07_08) plot(str_07_09) par(mfrow = c(1, 3)) plot(str_07_10) plot(str_07_11) plot(str_07_12)
アルゴリズムによって少しずつ推定された構造が違うのですが、low -> smokeおよびrace -> smokeのエッジは概ね共通して見られるようです。しかし残念ながらこのデータは「出生児の低体重の有無」と「母体環境」であるため、エッジの向きは逆ですね。。。
