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).

Last updated on `Tue Nov 8 15:59:15 2022` with bnlearn `4.9-20221107` and `R version 4.2.2 (2022-10-31)`.