Rでベイジアンネットワークメモ

#データとしては、MASS::birthwt データを使います。
#これは、妊娠時の喫煙有無などの状況が低体重児につながるかどうかの調査データみたいです。
library(MASS)
library(data.table)
library(dplyr)
## 
## Attaching package: 'dplyr'
## 
##  以下のオブジェクトは 'package:data.table' からマスクされています: 
## 
##      between, last 
## 
##  以下のオブジェクトは 'package:MASS' からマスクされています: 
## 
##      select 
## 
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag 
## 
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
data <- birthwt %>% data.table()
data
##      low age lwt race smoke ptl ht ui ftv  bwt
##   1:   0  19 182    2     0   0  0  1   0 2523
##   2:   0  33 155    3     0   0  0  0   3 2551
##   3:   0  20 105    1     1   0  0  0   1 2557
##   4:   0  21 108    1     1   0  0  1   2 2594
##   5:   0  18 107    1     1   0  0  1   0 2600
##  ---                                          
## 185:   1  28  95    1     1   0  0  0   2 2466
## 186:   1  14 100    3     0   0  0  0   2 2495
## 187:   1  23  94    3     1   0  0  0   0 2495
## 188:   1  17 142    2     0   0  1  0   0 2495
## 189:   1  21 130    1     1   0  1  0   3 2495
#このデータから、体重が低いかどうか(low)、人種(race)、妊娠時の喫煙の有無(smoke)、高血圧の病歴の有無(ht)、子宮の過敏性の有無(ui)の5つの変数を取り出し、ベイジアンネットワークを作成します。
#ちなみに人種は、1=白人、2=黒人、3=その他です。
data <- data %>% select(low, race, smoke, ht, ui)
colnames(data) <- c("体重","人種","喫煙","高血圧","過敏性")
## Warning in `names<-.data.table`(`*tmp*`, value = value): The colnames(x)<-
## value syntax copies the whole table. This is due to <- in R itself. Please
## change to setnames(x,old,new) which does not copy and is faster. See
## help('setnames'). You can safely ignore this warning if it is inconvenient
## to change right now. Setting options(warn=2) turns this warning into an
## error, so you can then use traceback() to find and change your colnames<-
## calls.
data
##      体重 人種 喫煙 高血圧 過敏性
##   1:    0    2    0      0      1
##   2:    0    3    0      0      0
##   3:    0    1    1      0      0
##   4:    0    1    1      0      1
##   5:    0    1    1      0      1
##  ---                             
## 185:    1    1    1      0      0
## 186:    1    3    0      0      0
## 187:    1    3    1      0      0
## 188:    1    2    0      1      0
## 189:    1    1    1      1      0
data$体重 <- ifelse(data$体重 == 0,"低体重","正常")
data$人種 <- ifelse(data$人種 == 1,"白人",
                  ifelse(data$人種 == 2,"黒人","その他"))
data$喫煙 <- ifelse(data$喫煙 == 0,"喫煙なし","喫煙あり")
data$高血圧 <- ifelse(data$高血圧 == 0,"高血圧なし","高血圧あり")
data$過敏性<- ifelse(data$過敏性 == 0,"過敏性なし","過敏性あり")
data
##        体重   人種     喫煙     高血圧     過敏性
##   1: 低体重   黒人 喫煙なし 高血圧なし 過敏性あり
##   2: 低体重 その他 喫煙なし 高血圧なし 過敏性なし
##   3: 低体重   白人 喫煙あり 高血圧なし 過敏性なし
##   4: 低体重   白人 喫煙あり 高血圧なし 過敏性あり
##   5: 低体重   白人 喫煙あり 高血圧なし 過敏性あり
##  ---                                             
## 185:   正常   白人 喫煙あり 高血圧なし 過敏性なし
## 186:   正常 その他 喫煙なし 高血圧なし 過敏性なし
## 187:   正常 その他 喫煙あり 高血圧なし 過敏性なし
## 188:   正常   黒人 喫煙なし 高血圧あり 過敏性なし
## 189:   正常   白人 喫煙あり 高血圧あり 過敏性なし
#さて、ベイジアンネットワークを作成するためには deal パッケージを使って次のように行います。
library(deal)
data[] <- lapply(data, as.factor)
pre.network <- network(data)
prior.dist <- jointprior(pre.network)
## Imaginary sample size: 96
update <- learn(pre.network, data, prior.dist)
post.network <- autosearch(getnetwork(update), data, 
                           prior.dist, trace=FALSE)
## [Autosearch (1) -619.4821 [体重][人種|喫煙][喫煙][高血圧][過敏性]
## (2) -616.5 [体重][人種|喫煙][喫煙][高血圧|過敏性][過敏性]
## (3) -613.438 [体重][人種|喫煙][喫煙][高血圧|体重:過敏性][過敏性]
## (4) -610.4685 [体重][人種|喫煙][喫煙][高血圧|体重:過敏性][過敏性|体重]
## (5) -607.8736 [体重][人種|体重:喫煙][喫煙][高血圧|体重:過敏性][過敏性|体重]
## (6) -606.3248 [体重|喫煙][人種|体重:喫煙][喫煙][高血圧|体重:過敏性][過敏性|体重]
## (7) -605.3213 [体重|喫煙][人種|体重:喫煙:高血圧][喫煙][高血圧|体重:過敏性][過敏性|体重]
## ...(8) -605.2876 [体重|喫煙][人種|体重:喫煙:高血圧][喫煙][高血圧|体重:過敏性][過敏性|体重:喫煙]
## ....Total 0.17 add 0.09 rem 0.01 turn 0.07 sort 0 choose 0 rest 0 ]
plot(getnetwork(post.network))

これ見ると、低体重であることが高血圧の原因になっていますが、低体重は結果変数なので、原因にはなりえません。また、喫煙の有無が人種の原因になっていますが、これもおかしいですね。

というわけで、ベイジアンネットワークを作成するときに、ここからここへは辺をつなげないという制約を加えることができます。

ここでは、「低体重はどの変数の原因にもならない」と 「人種はどの変数の結果にもならない」という制約を置いてみましょう。

低体重のノード番号は1、人種のノード番号は2なので、

pre.network <- getnetwork(update)
banlist(pre.network) <- rbind(cbind(1, 2:5),
                              cbind(1:5, 2))
post.network <- autosearch(pre.network, data,
                           prior.dist, trace=FALSE)
## [Autosearch (1) -619.4821 [体重][人種][喫煙|人種][高血圧][過敏性]
## (2) -616.5 [体重][人種][喫煙|人種][高血圧|過敏性][過敏性]
## (3) -613.5305 [体重|過敏性][人種][喫煙|人種][高血圧|過敏性][過敏性]
## (4) -610.4685 [体重|高血圧:過敏性][人種][喫煙|人種][高血圧|過敏性][過敏性]
## (5) -608.9411 [体重|高血圧:過敏性][人種][喫煙|人種][高血圧|人種:過敏性][過敏性]
## (6) -607.8464 [体重|高血圧:過敏性][人種][喫煙|人種:過敏性][高血圧|人種:過敏性][過敏性]
## (7) -606.9597 [体重|喫煙:高血圧:過敏性][人種][喫煙|人種:過敏性][高血圧|人種:過敏性][過敏性]
## (8) -606.0825 [体重|喫煙:高血圧:過敏性][人種][喫煙|人種:過敏性][高血圧|人種][過敏性|高血圧]
## .Total 0.15 add 0.13 rem 0 turn 0.01 sort 0 choose 0.01 rest 0 ]
plot(getnetwork(post.network), showban=FALSE)

localprob(getnetwork(post.network))
## $体重
## , , 高血圧あり, 過敏性あり
## 
##         喫煙あり  喫煙なし
## 正常   0.5000000 0.5000000
## 低体重 0.5000000 0.5000000
## 
## , , 高血圧なし, 過敏性あり
## 
##         喫煙あり  喫煙なし
## 正常   0.5200000 0.4814815
## 低体重 0.4800000 0.5185185
## 
## , , 高血圧あり, 過敏性なし
## 
##         喫煙あり  喫煙なし
## 正常   0.5294118 0.5263158
## 低体重 0.4705882 0.4736842
## 
## , , 高血圧なし, 過敏性なし
## 
##         喫煙あり  喫煙なし
## 正常   0.3823529 0.2285714
## 低体重 0.6176471 0.7714286
## 
## 
## $人種
## 
##    その他      黒人      白人 
## 0.3473684 0.2035088 0.4491228 
## 
## $喫煙
## , , 過敏性あり
## 
##         喫煙あり  喫煙なし
## その他 0.3243243 0.4102564
## 黒人   0.2162162 0.2820513
## 白人   0.4594595 0.3076923
## 
## , , 過敏性なし
## 
##         喫煙あり  喫煙なし
## その他 0.1882353 0.4435484
## 黒人   0.2117647 0.1693548
## 白人   0.6000000 0.3870968
## 
## 
## $高血圧
##        高血圧あり 高血圧なし
## その他  0.3333333  0.3511111
## 黒人    0.3166667  0.1733333
## 白人    0.3500000  0.4755556
## 
## $過敏性
##            過敏性あり 過敏性なし
## 高血圧あり  0.3157895  0.1722488
## 高血圧なし  0.6842105  0.8277512