## Bootstrap-based inference

### The general case

A general-purpose bootstrap implementation, similar in scope to the `boot()`

function
in package **boot**, is provided by the `bn.boot()`

function (documented
here). `bn.boot()`

takes a data set, a
structure learning algorithm and an arbitrary function (whose first argument must be an object of
class `bn`

) and returns a list of the values returned by said function for the network
structures learned from the bootstrapped samples.

For example, we may want to know how many arcs we can expect in a network learned with hill
climbing from the `learning.test`

data set (documented
here). We can do it as follows:

> library(bnlearn) > unlist(bn.boot(learning.test, statistic = narcs, + algorithm = "hc", R = 10))

[1] 5 5 5 5 5 5 5 5 5 5

The `R`

argument controls how many bootstrap replicates are performed. Or maybe we
want to compare the computational complexity (measured with the numbers of test/score comparisons)
between hill climbing and Grow-Shrink for a sample size of 500:

> unlist(bn.boot(learning.test, statistic = ntests, + algorithm = "hc", R = 10))

[1] 40 40 40 40 40 40 40 40 40 40

> unlist(bn.boot(learning.test, statistic = ntests, + algorithm = "gs", R = 10))

[1] 371 432 375 374 382 317 448 422 401 335

Many other questions can be answered with this approach; essentially any function of the network
structure can be used for the `statistic`

argument. We can also return the structures
themselves using a dummy function as follows.

> bn.boot(learning.test, statistic = function(x) x, + algorithm = "hc", R = 2)

[[1]] Bayesian network learned via Score-based methods model: [A][C][F][B|A][D|A:C][E|B:F] nodes: 6 arcs: 5 undirected arcs: 0 directed arcs: 5 average markov blanket size: 2.33 average neighbourhood size: 1.67 average branching factor: 0.83 learning algorithm: Hill-Climbing score: BIC (disc.) penalization coefficient: 4.258597 tests used in the learning procedure: 40 optimized: TRUE [[2]] Bayesian network learned via Score-based methods model: [A][C][F][B|A][D|A:C][E|B:F] nodes: 6 arcs: 5 undirected arcs: 0 directed arcs: 5 average markov blanket size: 2.33 average neighbourhood size: 1.67 average branching factor: 0.83 learning algorithm: Hill-Climbing score: BIC (disc.) penalization coefficient: 4.258597 tests used in the learning procedure: 40 optimized: TRUE

In addition, we can control how many observations are included in each bootstrap sample with the
`m`

argument.

> unlist(bn.boot(learning.test, statistic = narcs, + algorithm = "hc", R = 10, m = 50))

[1] 3 4 3 3 4 4 4 4 3 1

### Measuring arc strength

Measuring the degree of *confidence* in a particular graphical feature of a Bayesian
network is a key problem in the inference on the network structure. In the case of single arcs this
quantity is called *arc strength*.

Friedman, Goldszmidt and Wyner (1999)
introduced a very simple way of quantifying such a confidence: generating multiple network structures
by applying nonparametric bootstrap to the data and estimating the relative frequency of the feature
of interest. `boot.strength()`

uses this approach to compute the strength of every possible
arc, and has a syntax similar to that of `bn.boot()`

.

> library(bnlearn) > boot.strength(learning.test, algorithm = "hc")

from to strength direction 1 A B 1.000 0.5000 2 A C 0.000 0.0000 3 A D 1.000 1.0000 4 A E 0.000 0.0000 5 A F 0.005 0.5000 6 B A 1.000 0.5000 7 B C 0.010 0.5000 8 B D 0.000 0.0000 9 B E 1.000 0.9875 10 B F 0.025 0.5000 11 C A 0.000 0.0000 12 C B 0.010 0.5000 13 C D 1.000 1.0000 14 C E 0.000 0.0000 15 C F 0.005 0.5000 16 D A 1.000 0.0000 17 D B 0.000 0.0000 18 D C 1.000 0.0000 19 D E 0.000 0.0000 20 D F 0.000 0.0000 21 E A 0.000 0.0000 22 E B 1.000 0.0125 23 E C 0.000 0.0000 24 E D 0.000 0.0000 25 E F 1.000 0.0125 26 F A 0.005 0.5000 27 F B 0.025 0.5000 28 F C 0.005 0.5000 29 F D 0.000 0.0000 30 F E 1.000 0.9875

Note that this approach computes the joint strength of all the possible arcs; the estimates will
not be independent. For each pair of nodes, the probability that there is an arc between them
regardless of its direction is stored in the `strength`

column, and the probability of
each direction is in the `direction`

column. This parameterization follows
Imoto et al. (2002).

`Tue Nov 8 15:59:15 2022`

with **bnlearn**

`4.9-20221107`

and `R version 4.2.2 (2022-10-31)`

.