bnlearn - Boostrap-based inference

Boostrap-based inference

Measuring arc strength

Measuring the degree of confidence in particular graphical features of a Bayesian network has always been one key problem in the inference on the network structure. In the case of single arcs this quantity has also been called arc strength.

Friedman 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. The simplest application of this approach allows an easy computation of the arc strength for the arcs present in a given network:

> res = hc(learning.test)
> arc.strength(res, learning.test, criterion = "bootstrap")
  from to strength
1    A  B        1
2    A  D        1
3    C  D        1
4    B  E        1
5    F  E        1

Note that this approach does not require discrete data; it works fundamentally in the same way for any type of data supported by the learning algorithms.

Another interesting application is checking how the strength coefficients vary with the sample size:

> a50 = arc.strength(res, learning.test, criterion = "bootstrap", m = 50)
> a100 = arc.strength(res, learning.test, criterion = "bootstrap", m = 100)
> a200 = arc.strength(res, learning.test, criterion = "bootstrap", m = 200)
> a300 = arc.strength(res, learning.test, criterion = "bootstrap", m = 300)
> a400 = arc.strength(res, learning.test, criterion = "bootstrap", m = 400)
> data.frame(arcs(res), strength50 = a50$strength, strength100 = a100$strength, 
+   strength200 = a200$strength, strength300 = a300$strength, strength400 = a400$strength)
  from to strength50 strength100 strength200 strength300 strength400
1    A  B      0.250       0.930        1.00        1.00       1.000
2    A  D      0.240       0.925        1.00        1.00       1.000
3    C  D      0.000       0.000        0.14        0.79       0.995
4    B  E      0.020       0.340        0.93        1.00       1.000
5    F  E      0.125       0.635        0.97        1.00       1.000

Measuring arc strength and confidence in direction

Another, related, problem is to assess arc strength (irrespective of its direction) and confidence in its direction (given its presence in the graph) at the same time. The solution presented in Imoto et al. is essentally the same as in the previous case.

> boot.strength(learning.test, R = 200, m = 5000, "gs")
   from to strength  direction
1     A  B    1.000 0.49250000
2     A  C    0.060 0.25000000
3     A  D    1.000 0.96750000
4     A  E    0.295 0.68644068
5     A  F    0.030 0.58333333
6     B  A    1.000 0.50750000
7     B  C    0.385 0.35064935
8     B  D    0.465 0.93548387
9     B  E    1.000 0.90750000
10    B  F    0.165 0.15151515
11    C  A    0.060 0.75000000
12    C  B    0.385 0.64935065
13    C  D    1.000 0.97250000
14    C  E    0.210 0.91666667
15    C  F    0.035 0.42857143
16    D  A    1.000 0.03250000
17    D  B    0.465 0.06451613
18    D  C    1.000 0.02750000
19    D  E    0.340 0.26470588
20    D  F    0.005 0.00000000
21    E  A    0.295 0.31355932
22    E  B    1.000 0.09250000
23    E  C    0.210 0.08333333
24    E  D    0.340 0.73529412
25    E  F    1.000 0.03000000
26    F  A    0.030 0.41666667
27    F  B    0.165 0.84848485
28    F  C    0.035 0.57142857
29    F  D    0.005 1.00000000
30    F  E    1.000 0.97000000

Interpretation of the output of boot.strength (documented here) is usefult both to identify arcs representing strong relationships and to distinguish compelled arcs from score equivalent ones.

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.

Let's assume for example that we want to know how many arcs we can expect in a network learned with a hill climbing search from the learning.test data set. It can be done in the following way:

> unlist(bn.boot(learning.test, statistic = function(x) { nrow(arcs(x)) },
+   algorithm = "hc", R = 10))
 [1] 5 5 5 5 5 5 5 5 6 5

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 = function(x) { x$learning$ntests },
+   algorithm = "hc", R = 10, m = 50))
 [1] 35 20 30 25 30 35 30 35 20 25
> unlist(bn.boot(learning.test, statistic = function(x) { x$learning$ntests },
+   algorithm = "gs", R = 10, m = 50))
 [1] 26 16 31 26 27 25 17 27 22 22

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