ベイジアンネットワークをRのbnlearnパッケージで推定して予測してみる

ベイジアンネットワークの知見が無かったので、調べた情報をまとめています。一応、載せているスクリプトでRを用いて予測するということができます。

【目次】
・ベイジアンネットワークとは
・ベイジアンネットワークの用途
・ベイジアンネットワークの推定のステップ
・Rでの実行例
・おまけ：Webサービスへの応用例
・参考文献

ベイジアンネットワークとは

・複数の確率変数の間の定性的な依存関係をグラフ構造によって表し、個々の変数の間の定量的な関係を条件付確率で表した確率モデル。前提として、有向非循環（グルグルと回らないグラフ）となっているグラフの構造を持つものに限定している。
・入力となる変数と出力となる変数はモデルの中では区別されない。
・時間という明確な因果関係などをモデルに組み込みやすいので、系列データなどを扱うケースが多い。モデル設計者がデータが生成されるプロセスを考慮しながらモデルを組んでいける。
・非循環性とd分離の仮定のみによって導かれる現在考えられる最も自然な離散モデルであり、現在の様々なモデルの中でも、最も表現力と予測力を持つモデルとされる。

ベイジアンネットワークの用途

・物体追跡
・ジェスチャ認識
・Webサイトなどのレコメンデーションサービス
・広告配信
・メルマガなどの配信最適化
など画像処理問題から消費者行動問題までいろいろな分野で活用されているようです。

ベイジアンネットワークの推定のステップ

モデルの選択や推定は学習 (learning)と呼ばれ、以下の2段階のステップを踏みます。

Step1：構造学習(structure learning)…データからネットワーク構造を学習する。
Step2：パラメータ学習(parameter learning)…Step1で学習した構造によって意味付けした分布のパラメータを学習する。

ベイジアンネットワークの推定は条件付き確率やMAP(Maximum a posteriori)推定量などを用います。構造学習の際に、事前知識を導入することができます。

Rでの実行例

・データ
何千もの主要な免疫系細胞から取り出した11種類のリン酸タンパク質とリン脂質の測定値からなるデータです。以下のサイトよりダウンロードできます。
Supporting Online Material Causal Protein-Signaling Networks Derived from Multiparameter Single-Cell Data

・パッケージ
ベイジアンネットワークを実行できるパッケージはいろいろありますが、参考文献に従いbnlearnを用います。加えて、ネットワーク図の可視化のためにRgraphvizを用います。Rgraphvizがcranでは対応したものがなかったので、こちらを参考にして、インストールしました。
Provides plotting capabilities for R graph objects

・ネットワークの図示

#ネットワークの図示
library(Rgraphviz)
library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",
             "[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",
             "[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)
class(net)

graphviz.plot(net, shape = "ellipse")

#ネットワークの図示

library(Rgraphviz)

library(bnlearn)

spec = paste("[PKC][PKA|PKC][praf|PKC:PKA][pmek|PKC:PKA:praf]",

"[p44.42|pmek:PKA][pakts473|p44.42:PKA][P38|PKC:PKA]",

"[pjnk|PKC:PKA][plcg][PIP3|plcg][PIP2|plcg:PIP3]")

net = model2network(spec)

class(net)

graphviz.plot(net, shape = "ellipse")

・データ確認

#データ確認
> protein_dataset <-  read.csv(file = "protain_dataset.csv",as.is = TRUE)
> head(protein_dataset,5)
   praf pmek  plcg PIP2  PIP3 p44.42 pakts473  PKA   PKC  P38 pjnk
1 39.20  8.9  4.26 9.22  1.88   25.9     46.6  223  3.02 31.6 65.5
2 42.60 20.2 12.10 1.45 40.00   30.2     51.9  264  5.94 14.1 16.1
3  5.88 17.9  3.43 2.19 18.10  105.0    300.0 1928  6.98 20.5 12.9
4 34.30 18.8  3.65 1.70 16.50  352.0    519.0 1263 18.30 59.4 35.9
5 24.60 16.0  6.92 9.39  8.06   15.7     53.8  346  9.31 18.1 12.5

#データ確認

> protein_dataset <- read.csv(file = "protain_dataset.csv",as.is = TRUE)

> head(protein_dataset,5)

praf pmek plcg PIP2 PIP3 p44.42 pakts473 PKA PKC P38 pjnk

1 39.20 8.9 4.26 9.22 1.88 25.9 46.6 223 3.02 31.6 65.5

2 42.60 20.2 12.10 1.45 40.00 30.2 51.9 264 5.94 14.1 16.1

3 5.88 17.9 3.43 2.19 18.10 105.0 300.0 1928 6.98 20.5 12.9

4 34.30 18.8 3.65 1.70 16.50 352.0 519.0 1263 18.30 59.4 35.9

5 24.60 16.0 6.92 9.39 8.06 15.7 53.8 346 9.31 18.1 12.5

・構造学習
データセットから、ネットワーク構造を推定します。推定するための手法は様々あるようですが、今回はヒルクライムアルゴリズムを用います。

#構造学習の実行
> print(hc(protein_dataset))

  Bayesian network learned via Score-based methods

  model:
   [PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]
   [PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]
  nodes:                                 11 
  arcs:                                  11 
    undirected arcs:                     0 
    directed arcs:                       11 
  average markov blanket size:           2.18 
  average neighbourhood size:            2.00 
  average branching factor:              1.00 

  learning algorithm:                    Hill-Climbing 
  score:                                 BIC (Gauss.) 
  penalization coefficient:              3.348517 
  tests used in the learning procedure:  185 
  optimized:                             TRUE

#構造学習の実行

> print(hc(protein_dataset))

Bayesian network learned via Score-based methods

model:

[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]

[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]

nodes: 11

arcs: 11

undirected arcs: 0

directed arcs: 11

average markov blanket size: 2.18

average neighbourhood size: 2.00

average branching factor: 1.00

learning algorithm: Hill-Climbing

score: BIC (Gauss.)

penalization coefficient: 3.348517

tests used in the learning procedure: 185

optimized: TRUE

構造学習より得られたネットワーク構造を可視化します。因果関係がダメな感じのネットワークが出来ていないかチェックする必要がありますね。機械的に作るだけでなく、背景的知識も考慮するという進め方が推奨されています。

#推定したネットワークを可視化
spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",
                        "[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")
net.estimated = model2network(spec.estimated)
class(net.estimated)
graphviz.plot(net.estimated, shape = "ellipse")

#推定したネットワークを可視化

spec.estimated = paste("[PIP2][p44.42][PKC][PIP3|PIP2][pakts473|p44.42][pjnk|PKC][plcg|PIP3]",

"[PKA|p44.42:pakts473][P38|PKC:pjnk][pmek|plcg:P38][praf|pmek]")

net.estimated = model2network(spec.estimated)

class(net.estimated)

graphviz.plot(net.estimated, shape = "ellipse")

・パラメータ学習
先ほど、データセットから求めた構造について、以下のコードでパラメータ推定します。パラメータの推定や構造推定は試行錯誤するところなので、こんな簡単に済む訳では無いようです。

#パラメータ学習
bn.fit(学習した構造, データセット,method = "mle")

1 2	#パラメータ学習 bn.fit(学習した構造, データセット,method = "mle")

・予測
テストデータを用意して、データセットから求めた構造・パラメータを用いて、任意の変数の予測を行います。予測したい変数の値は欠損していたらエラーで予測してくれないですが、全部0にしておけば回ります。運用上は予測したい変数の値はわからないはずなので。

#ベイジアンネットワークによる予測
library(forecast)
training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ
test.set = protein_dataset[406:810, ]  #テスト用データ
baysian_structure = hc(training.set)   #訓練用データでベイジアンネットワークの構造を学習
fitted = bn.fit(baysian_structure, training.set,method = "mle")     #パラメータの学習
test.set2 <- test.set
test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。
pred = predict(fitted, "PKC", test.set2)  #テストデータが与えられたもとでのPKCの予測
head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出


> head(cbind(pred, test.set[, "PKC"]))      #予測値と実績値の比較
         pred     
[1,] 19.76000 22.3
[2,] 21.66300 26.4
[3,] 11.82080 27.9
[4,] 23.85982 53.3
[5,] 21.82902 34.3
[6,] 17.29319 16.8
> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出
              ME     RMSE      MAE       MPE     MAPE
Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#ベイジアンネットワークによる予測

library(forecast)

training.set = protein_dataset[1:405, ] #パラメータ推定用の訓練用データ

test.set = protein_dataset[406:810, ] #テスト用データ

baysian_structure = hc(training.set) #訓練用データでベイジアンネットワークの構造を学習

fitted = bn.fit(baysian_structure, training.set,method = "mle") #パラメータの学習

test.set2 <- test.set

test.set2$PKC <- 0 #本番想定でこの値を知らないということにしておく。

pred = predict(fitted, "PKC", test.set2) #テストデータが与えられたもとでのPKCの予測

head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

> head(cbind(pred, test.set[, "PKC"])) #予測値と実績値の比較

pred

[1,] 19.76000 22.3

[2,] 21.66300 26.4

[3,] 11.82080 27.9

[4,] 23.85982 53.3

[5,] 21.82902 34.3

[6,] 17.29319 16.8

> accuracy(f = pred, x = test.set[, "PKC"]) #正確度の算出

ME RMSE MAE MPE MAPE

Test set 0.47376 10.02839 7.479115 -90.60837 114.6982

#予測結果の可視化
library(ggplot2)
df <- data.frame(prediction=pred,actual=test.set[, "PKC"])
p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point() 
p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")
p

#予測結果の可視化

library(ggplot2)

df <- data.frame(prediction=pred,actual=test.set[, "PKC"])

p <- ggplot(df,aes(x = prediction,y = actual)) + geom_point()

p <- p + geom_line(data = data.frame(x = c(0,40), y = c(0,40)),aes(x = x, y = y), colour = "red")

全然45度線上に乗っていないので、あまり精度は高くないようです。

cpquery関数と言って、条件付き確率を予測できる関数もあるようです。
Perform conditional probability queries

おまけ：Webサービスへの応用例

参考文献には、Webサイトの閲覧データをもとに、ユーザーの行動を予測するというモデルの紹介がなされていました。ただ、Webサイトの場合、ページ数や商品の数・ユーザーの数も膨大なことから、ネットワーク構造の推定が難しいようです。そこで、ネットワークを作るに際して、変数を確率的潜在意味解析（pLSA）で扱いやすい数に絞るなどの工夫をされていました。自社のサイトに関しても適用する際に、変数の多さは待った無しだと思うので、いざ適用する際は情報圧縮技術を駆使したいですね。