bnlearn解析結果のサブグラフを描画する方法

# はじめに
ベイジアンネットワーク分析の代表的なライブラリに R の bnlearn がありますが、その描画部分の不明点を製作者の Marco Scutari 先生に教えてもらいました。とても有難かったのでシェアするために記事を書いています。

# TL;DR
bnlearn で作成したネットワーク図のうち、一部分を描画したいときの tips です。大きなグラフのサブグラフでも可能ですし、他のノードと接続されていない孤立ノードの削除にも使えます。
結論としていいメソッドは存在しないので手動で行う必要があり、簡単な関数化してみました。

# 準備
(自分の備忘録のために残していますが、この記事を読むような人には不要なセクションです)
```R:ライブラリの準備
install.packages("bnlearn");
#devtools::install_github("vspinu/bnlearn") #こちらでもOK
install.packages("BiocManager");
BiocManager::install("Rgraphviz");
library(bnlearn);
library(Rgraphviz);
```
```R:データの準備
set.seed(0);
N <- 100;
x1 <- round(runif(N), 0);
x2 <- round(x1 + rnorm(N) * 0.2, 0);
x3 <- round(runif(N), 0);
x4 <- round(x1 / 2 + x3 / 2 + rnorm(N) * 0.2, 0);
x5 <- round(runif(N), 0);
dat <- data.frame(x1, x2, x3, x4, x5);
dat[dat > 1] <- 1;
dat[dat < 0] <- 0;
dat <- data.frame(lapply(dat, as.factor));
```
データは x1 -> x2、{x1, x3} -> x4 で x5 が独立しています。
今回は離散なので factor 型に変換。

# 解析
ここも本題ではないので必要に応じて読み飛ばしてください。
※[@hrkz_szk](https://qiita.com/hrkz_szk) 様の [bnlearnを使ってベイジアンネットワーク分析をやってみた](https://qiita.com/hrkz_szk/items/a213c2c4ba823cbf78f6) ほとんどそのままです。そちらの方が勉強になります。
```R:基本的なグラフ構造 (hill-climbing)
g <- hc(dat, score = "bic");
graphviz.plot(g, shape = "ellipse");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/539151/a87d5dcb-b898-44d3-44d1-b36169acd168.png)
元の構造を正確に捉えています。
このデータなら必要ないのですが、より複雑なデータではブートストラップを行うと構造推定の精度が向上します。

```R:ブートストラップ
str_g <- boot.strength(dat, R = 200, algorithm = "hc", 
                       algorithm.args = list(score = "bic"));
plot(str_g);
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/539151/5b162746-b716-269c-3d59-384466570ddc.png)

```R:ベストなパラメータでの描画
avg_g <- averaged.network(str_g);
strength.plot(avg_g, str_g, shape = "ellipse");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/539151/78d8b5e7-ca59-37a8-891e-2892f39fd5cc.png)
今回は変わらないですね。
ちなみに閾値を変えると別の結果になります
```R:閾値を下げた場合
g_lower <- averaged.network(str_g);
strength.plot(g_lower, str_g, shape = "ellipse");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/539151/cd2c73b5-0b96-bfab-7827-189fceb3ece7.png)
無関係なはずだった x3とx5 にエッジが引かれてしまっています。

# 本題：部分的なグラフ描画
### 少しだけ中身に触れます
strength.plot 関数のうち、第1引数は bn クラス、第2引数は bn.strength クラスです。このうち bn クラスは subgraph 関数を用いて操作できるのですが、bn.strength クラスにはいいメソッドがないのが問題でした。Marco Scutari 先生曰く

>サブセットを作ろうとすると R が attribute を消してしまうのが問題です

ということです。解決策としては、手動で attribute を付け加えるという形になります。
### コード
```R:bn.strengthクラスのサブセットを作る関数
set_bn_attr <- function(bn_strength, subset_names){
    subset_s <- bn_strength[(bn_strength$from %in% subset_names) &
                            (bn_strength$to %in% subset_names), ];
    attr(subset_s, "threshold") <- attr(bn_strength, "threshold");
    attr(subset_s, "method") <- attr(bn_strength, "method");
    attr(subset_s, "nodes") <- subset_names;
    attr(subset_s, "class") <- attr(bn_strength, "class");
    return(subset_s)
}
```
```R:関数を使って新たなbn.strengthオブジェクトを作り描画
ss <- set_bn_attr(str_g, c("x1", "x2", "x3", "x4"));
strength.plot(subgraph(avg_g, nnames), ss, shape = "ellipse");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/539151/d264da35-7997-a22b-e3bd-365a127c93ef.png)
という訳であっさり解決しました。Marco Scutari 先生に重ねて感謝です。

なおネットワークの中身に関しては、bn オブジェクトを見るといろいろ分かりそうです。

```R:ネットワークの中身
> names(avg_g)
[1] "learning" "nodes"    "arcs"

> avg_g$nodes
$x1
$x1$mb
[1] "x2" "x3" "x4"

$x1$nbr
[1] "x2" "x4"

$x1$parents
character(0)

$x1$children
[1] "x2" "x4"


$x2
$x2$mb
[1] "x1"

$x2$nbr
[1] "x1"

$x2$parents
[1] "x1"

$x2$children
character(0)


$x3
$x3$mb
[1] "x1" "x4"

$x3$nbr
[1] "x4"

$x3$parents
character(0)

$x3$children
[1] "x4"


$x4
$x4$mb
[1] "x1" "x3"

$x4$nbr
[1] "x1" "x3"

$x4$parents
[1] "x1" "x3"

$x4$children
character(0)


$x5
$x5$mb
character(0)

$x5$nbr
character(0)

$x5$parents
character(0)

$x5$children
character(0)


```