bnlearn - Bayesian networks and cross-validation

Bayesian networks and cross-validation

Choosing a Bayesian network learning strategy

Cross-validation is a standard way to obtain unbiased estimates of a model's goodness of fit. By comparing such estimates for different learning strategies (different combinations of learning algorithms, fitting techniques and the respective parameters) we can choose the optimal one for the data at hand in a principled way.

bnlearn implements three cross-validation methods in the bn.cv() function (documented here):

k-fold cross-validation (the default): the data are randomly partitioned into k subsets. Each subset is used in turn to validate the model fitted on the remaining k - 1 subsets.
Custom folds cross-validation: the data are manually partitioned into subsets by the user into subsets, which are then used as in k-fold cross-validation. Subsets are not constrained to have the same size.
Hold-out cross-validation: the data are repeatedly (randomly) partitioned in a training and test subset of given size m and n - m, with observations assigned randomly to each subset at each repetition. Each test subset is used to validate the model fitted on the corresponding training subset.

bn.cv() requires the following parameters for all methods:

method: the label of the cross-validation method, either "k-fold" (the default), "custom-fold" or "hold-out";
data: the data set the Bayesian network will be estimated from;
bn: the structure learning algorithm to be used, or a fixed network structure (to cross-validate just parameter learning);
loss:: the loss function used to evaluate the goodness of fit of the learned network;
fit: the fitting method (such as mle or bayes) to be used to fit the parameters.

Additional parameters for all these components can be specified as well, providing a fine-grained control of the learning strategy, using the following arguments of bn.cv():

algorithm.args for additional arguments to be passed to the learning algorithm;
loss.args for additional arguments to be passed to the loss function;
fit.args for additional arguments to be passed to bn.fit() for fitting the parameters of the Bayesian network.

k-fold cross-validation

For example, let's compare the performance of three different scores in a hill climbing search on the ALARM data set included in bnlearn. The loss funtion we use is the so-called log-likelihood loss, which is the negated expected log-likelihood; so lower values are better. Hill cimbing is implemented in the hc() function.

> library(bnlearn)
> data(alarm)
> bn.cv(alarm, bn = "hc", algorithm.args = list(score = "bde", iss = 1))


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.82857

In this case the bn argument specifies the label of a structure learning algorithm, and algorithm.args is an optional list of arguments for the algorithm specified by bn.

> bn.cv(alarm, bn = "hc", algorithm.args = list(score = "bde", iss = 10))


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.82172

> bn.cv(alarm, bn = "hc", algorithm.args = list(score = "bic"))


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.86586

Partially directed graphs obtained from learning are automatically extended to directed acyclic graphs with cextend, so that the parameters can be fit and the loss function computed.

> bn.cv(alarm, bn = "iamb")


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             IAMB
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         13.70798

Tipically, 10-fold cross-validation is performed 10 times for each learning strategy (golden standard originally introduced by Hastie, Tibshirani and Friedman in «The Elements of Statistical Learning», link), and the resulting sets of loss values compared by plotting them as boxplots. Below we do so that for the BIC and BDe scores using hc().

> cv.bic = bn.cv(alarm, bn = "hc", runs = 10,
+            algorithm.args = list(score = "bic"))
> cv.bde = bn.cv(alarm, bn = "hc", runs = 10,
+            algorithm.args = list(score = "bde", iss = 1))
> plot(cv.bic, cv.bde, xlab = c("BIC", "BDe"))

Loading required namespace: lattice

Note that, while it is customary to use k = 10, a different number of folds can be specified with the k argument.

> bn.cv(alarm, bn = "hc", k = 5)


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       5
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.87307

Custom folds in cross-validation

The custom-folds method works in the same way as k-fold cross-validation, but folds are specified manually with a folds argument.

> folds = list(1:5000, 5001:15000, 15001:20000)
> cv.custom = bn.cv(alarm, bn = "hc", algorithm.args = list(score = "bic"),
+               method = "custom-folds", folds = folds)
> cv.custom


  custom-folds cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.85651

The two methods are otherwise identical, as shown below: the R objects returned by bn.cv() are identical (with the exception of the method's label) if the same folds are used.

> cv.random = bn.cv(alarm, bn = "hc", method = "k-fold")
> folds = lapply(cv.random, `[[`, "test")
> cv.custom = bn.cv(alarm, bn = "hc", method = "custom-folds",
+                 folds = folds)
> all.equal(cv.random, cv.custom, check.attributes = FALSE)

[1] TRUE

Note that:

Each observation can only be included in a single fold.

> folds = list(1:5000, 4000:10000, 10001:20000)
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)

## Error: some observations are included in more than one fold.

Each fold must contain at least one observation.

> folds = list(numeric(0), 1:2000, 2001:20000)
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)

## Error: some folds contain no observations.

All observations must be assigned to the folds.

> folds = list(1:1000, 2000:20000)
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)

## Error: not all observations are assigned to a fold.

runs is not a valid argument for this cross-validation method, because the ability of performing multiple runs relies on splitting the data randomly into folds.

To perform multiple runs, each set of manually specified folds (one for each run) should be saved in a list.

> folds = list(
+   list(1:5000, 5001:10000, 10001:20000),
+   list(1:10000, 10001:15000, 15001:20000),
+   list(1:5000, 5001:15000, 15001:20000)
+ )
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)


  custom-folds cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  loss function:                         Log-Likelihood Loss (disc.)
  number of runs:                        3
  average loss over the runs:            10.86871
  standard deviation of the loss:        0.01063191

The following function can be used to generate the folds randomly,

> random.folds = function(data, k) split(sample(nrow(data)), seq(k))

for a single run or cross-validation:

> folds = random.folds(data = alarm, k = 10)
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)


  custom-folds cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.86904

or for multiple runs.

> folds = replicate(10, random.folds(data = alarm, k = 10), simplify = FALSE)
> bn.cv(alarm, bn = "hc", method = "custom-folds", folds = folds)


  custom-folds cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  loss function:                         Log-Likelihood Loss (disc.)
  number of runs:                        10
  average loss over the runs:            10.86851
  standard deviation of the loss:        0.002829784

Hold-out cross-validation

The syntax for hold-out cross-validation is very similar to that for k-fold cross-validation, but:

method must be set to "hold-out";
k denotes the number of times the sample will be split into training and test subsamples (the default is 10);
and optionally m denotes the number of observations to be sampled for the test subsample (the default is 1/10th of the sample size).

So, for instance:

> bn.cv(alarm, bn = "hc", method = "hold-out", k = 25, m = 200,
+   algorithm.args = list(score = "bde", iss = 10))


  hold-out cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of splits:                      25
  size of the test subset:               200
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.7858

Comparing different network structures

bn.cv() can also be used to compare the goodness of fit of different network structures suggested by expert knowledge; using all the data to fit the parameters of those networks and to evaluate their goodness of fit produces optimistically biased results.

Let's start with the most obvious comparison: the true network structure against the empty one.

> res = empty.graph(names(alarm))
> modelstring(res) = paste("[HIST|LVF][CVP|LVV][PCWP|LVV][HYP][LVV|HYP:LVF]",
+   "[LVF][STKV|HYP:LVF][ERLO][HRBP|ERLO:HR][HREK|ERCA:HR][ERCA]",
+   "[HRSA|ERCA:HR][ANES][APL][TPR|APL][ECO2|ACO2:VLNG][KINK]",
+   "[MINV|INT:VLNG][FIO2][PVS|FIO2:VALV][SAO2|PVS:SHNT][PAP|PMB][PMB]",
+   "[SHNT|INT:PMB][INT][PRSS|INT:KINK:VTUB][DISC][MVS][VMCH|MVS]",
+   "[VTUB|DISC:VMCH][VLNG|INT:KINK:VTUB][VALV|INT:VLNG][ACO2|VALV]",
+   "[CCHL|ACO2:ANES:SAO2:TPR][HR|CCHL][CO|HR:STKV][BP|CO:TPR]", sep = "")
> bn.cv(alarm, res, algorithm.args = list(score = "aic"))


  k-fold cross-validation for Bayesian networks

  target network structure:
   [FIO2][MVS][HYP][LVF][APL][ANES][PMB][INT][KINK][DISC][ERLO][ERCA][HIST|LVF]
   [TPR|APL][PAP|PMB][LVV|HYP:LVF][STKV|HYP:LVF][SHNT|PMB:INT][VMCH|MVS]
   [CVP|LVV][PCWP|LVV][VTUB|DISC:VMCH][PRSS|INT:KINK:VTUB][VLNG|INT:KINK:VTUB]
   [MINV|INT:VLNG][VALV|INT:VLNG][PVS|FIO2:VALV][ACO2|VALV][SAO2|SHNT:PVS]
   [ECO2|ACO2:VLNG][CCHL|TPR:SAO2:ANES:ACO2][HR|CCHL][CO|STKV:HR][HRBP|ERLO:HR]
   [HREK|HR:ERCA][HRSA|HR:ERCA][BP|TPR:CO]
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         10.81904

> bn.cv(alarm, empty.graph(names(alarm)))


  k-fold cross-validation for Bayesian networks

  target network structure:
   [CVP][PCWP][HIST][TPR][BP][CO][HRBP][HREK][HRSA][PAP][SAO2][FIO2][PRSS][ECO2]
   [MINV][MVS][HYP][LVF][APL][ANES][PMB][INT][KINK][DISC][LVV][STKV][CCHL][ERLO]
   [HR][ERCA][SHNT][PVS][ACO2][VALV][VLNG][VTUB][VMCH]
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         21.46833

As expected, the first (true) network is a better fit than the second (empty network). The same is true for the following structure, which generated at random and has no relationship with the data.

> rand = random.graph(names(alarm))
> rand


  Random/Generated Bayesian network

  model:
   [CVP][PCWP][TPR][BP][CO][HRBP][HREK][PAP][PRSS][ECO2][MINV][LVF][APL][ANES]
   [KINK][HR][ERCA][HIST|PCWP][SAO2|HREK][PMB|PCWP][INT|PAP][LVV|BP:ECO2:KINK]
   [PVS|CO:HREK:LVF][VALV|ANES][VLNG|HRBP][HRSA|PCWP:HIST:TPR][FIO2|HIST:PAP]
   [DISC|HIST:APL][ERLO|SAO2:PMB][VTUB|CVP:KINK:LVV][MVS|TPR:FIO2]
   [STKV|PRSS:KINK:DISC][CCHL|DISC][SHNT|HRSA][VMCH|PCWP:HRSA][HYP|FIO2:MVS]
   [ACO2|HYP]
  nodes:                                 37
  arcs:                                  36
    undirected arcs:                     0
    directed arcs:                       36
  average markov blanket size:           2.86
  average neighbourhood size:            1.95
  average branching factor:              0.97

  generation algorithm:                  Full Ordering
  arc sampling probability:              0.05555556

> bn.cv(alarm, bn = rand)


  k-fold cross-validation for Bayesian networks

  target network structure:
   [CVP][PCWP][TPR][BP][CO][HRBP][HREK][PAP][PRSS][ECO2][MINV][LVF][APL][ANES]
   [KINK][HR][ERCA][HIST|PCWP][SAO2|HREK][PMB|PCWP][INT|PAP][LVV|BP:ECO2:KINK]
   [PVS|CO:HREK:LVF][VALV|ANES][VLNG|HRBP][HRSA|PCWP:HIST:TPR][FIO2|HIST:PAP]
   [DISC|HIST:APL][ERLO|SAO2:PMB][VTUB|CVP:KINK:LVV][MVS|TPR:FIO2]
   [STKV|PRSS:KINK:DISC][CCHL|DISC][SHNT|HRSA][VMCH|PCWP:HRSA][HYP|FIO2:MVS]
   [ACO2|HYP]
  number of folds:                       10
  loss function:                         Log-Likelihood Loss (disc.)
  expected loss:                         21.24115

Cross-validation and predictive error

Cross-validation is also commonly used in classification problems to compute the prediction error of a particular model. We can do that with bn.cv() by specifying

loss = "pred" to make frequentist predictions with predict(..., method = "parents"); or
loss = "pred-lw" to make Bayesian predictions with predict(..., method = "bayes-lw");

and passing the label of the target node to the loss function via the loss.args argument. The first predicts the target node from its parents; the second from all the available nodes or from a subset of nodes specified in loss.args list as element called from.

> xval = bn.cv(alarm, bn = "hc", loss = "pred", loss.args = list(target = "STKV"))
> xval


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Classification Error
  training node:                         STKV
  expected loss:                         0.1745

We may also want to compute other statistics from that run of cross-validation, like sensitivity and specificity, or just to tabulate observed and predicted values. To that end, we can extract the observed and predicted values for each fold from xval as follows.

> OBS = unlist(lapply(xval, `[[`, "observed"))
> PRED = unlist(lapply(xval, `[[`, "predicted"))
> str(OBS)

 Factor w/ 3 levels "HIGH","LOW","NORMAL": 3 2 3 3 3 3 2 3 3 3 ...

> str(PRED)

 Factor w/ 3 levels "HIGH","LOW","NORMAL": 3 3 3 3 3 3 3 3 3 3 ...

> table(OBS, PRED)

        PRED
OBS       HIGH   LOW NORMAL
  HIGH       0     7    826
  LOW        0   950   2624
  NORMAL     0    33  15560

Cross-validation and predictive correlation

Predictive correlation and mean square error are used to evaluate Gaussian Bayesian networks in the same way as prediction error is used to evaluate discrete Bayesian networks. bn.cv() implements:

loss = "cor" to compute the correlation between observed values and frequentist predictions from predict(..., method = "parents");
loss = "cor-lw" to compute the correlation between observed values and frequentist predictions from predict(..., method = "bayes-lw");
loss = "mse" to compute the correlation between observed values and frequentist predictions from predict(..., method = "parents");
loss = "mse-lw" to compute the correlation between observed values and frequentist predictions from predict(..., method = "bayes-lw").

> bn.cv(gaussian.test, bn = "hc", loss = "cor", loss.args = list(target = "C"))


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Predictive Correlation
  training node:                         C
  expected loss:                         0.9967885

> bn.cv(gaussian.test, bn = "hc", loss = "mse-lw", loss.args = list(target = "C"))


  k-fold cross-validation for Bayesian networks

  target learning algorithm:             Hill-Climbing
  number of folds:                       10
  loss function:                         Mean Squared Error (Posterior, Gauss.)
  training node:                         C
  expected loss:                         0.2592348

Last updated on Thu Nov 24 19:16:36 2022 with bnlearn 4.9-20221107 and R version 4.2.2 (2022-10-31).