## Score functions: computing & comparing

First, we load the `learning.test` and `gaussian.test` data sets shipped with bnlearn.

```> library(bnlearn)
```
```Loading required package: methods
```
```> data(learning.test)
> data(gaussian.test)
```

The true network structure for the `learning.test` data, described in the manual page, is the following:

```> learn.net = empty.graph(names(learning.test))
> modelstring(learn.net) = "[A][C][F][B|A][D|A:C][E|B:F]"
> learn.net
```
```
Random/Generated Bayesian network

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

generation algorithm:                  Empty
```

See the relevant examples for an overview on how `bn` objects are created using `modelstring()` and other functions. Similarly, we create a `bn` object for the true structure of `gaussian.test` (described here).

```> gauss.net = empty.graph(names(gaussian.test))
> modelstring(gauss.net) = "[A][B][E][G][C|A:B][D|B][F|A:D:E:G]"
> gauss.net
```
```
Random/Generated Bayesian network

model:
[A][B][E][G][C|A:B][D|B][F|A:D:E:G]
nodes:                                 7
arcs:                                  7
undirected arcs:                     0
directed arcs:                       7
average markov blanket size:           4.00
average neighbourhood size:            2.00
average branching factor:              1.00

generation algorithm:                  Empty
```

Both networks can be correctly learned by all the learning algorithms implemented in bnlearn, and provide one discrete and one continuous test case.

### Computing a network score

We can compute the network score of a particular graph for a particular data set with the `score()` function (manual); if the score function is not specified, the BIC score is returned for both continuous and discrete data.

```> score(learn.net, learning.test)
```
```[1] -24006.73
```
```> score(gauss.net, gaussian.test)
```
```[1] -53191.54
```

Other score functions can be used by changing the `type`, as documented in the manual page of the function.

```> score(learn.net, learning.test, type = "bic")
```
```[1] -24006.73
```
```> score(learn.net, learning.test, type = "aic")
```
```[1] -23873.13
```
```> score(learn.net, learning.test, type = "bde")
```
```[1] -24028.09
```

Many of these score have additional, optional arguments whose defaults are chosen based on literature golden standards and/or the data at hand. The most commonly used is the imaginary sample size (`iss`) of BDe and BGe, which is typically set to a very small value to reduce the relative weight of the prior distribution.

```> score(learn.net, learning.test, type = "bde", iss = 1)
```
```[1] -24028.09
```
```> score(gauss.net, gaussian.test, type = "bge", iss = 3)
```
```[1] -53472.64
```

Several check are run to ensure the network, the data and the test make sense with respect to each other:

• the number of variables in the network and in the data must be the same, although the order is not important.
```> score(learn.net, learning.test[, -1])
```
```## Error: the network and the data have different numbers of variables.
```
```> score(learn.net, learning.test[, sample(1:6)])
```
```[1] -24006.73
```
• the names of the variables must match as well.
```> renamed.test = learning.test
> names(renamed.test) = LETTERS[11:16]
> score(learn.net, renamed.test)
```
```## Error: the variables in the data and in the network do not match.
```
• the data must satisfy the assumptions of the network score (e.g. discrete data and discrete score, continuous data and Gaussian score, etc.).
```> score(learn.net, learning.test, type = "loglik-g")
```
```## Error: score 'loglik-g' may be used with continuous data only.
```
• the label of the score function must be valid; if it is not, valid scores are listed in the error message.
```> score(learn.net, learning.test, type = "???")
```
```## Error: valid score(s) are "loglik" (Log-Likelihood (disc.)), "aic" (AIC (disc.)), "bic" (BIC (disc.)), "bde" (Bayesian Dirichlet (BDe)), "bds" (Bayesian Dirichlet Sparse (BDs)), "bdj" (Bayesian Dirichlet, Jeffrey's prior), "k2" (Cooper & Herskovits' K2), "mbde" (Bayesian Dirichlet (interventional data)), "bdla" (Bayesian Dirichlet, Locally Averaged), "loglik-g" (Log-Likelihood (Gauss.)), "aic-g" (AIC (Gauss.)), "bic-g" (BIC (Gauss.)), "bge" (Bayesian Gaussian (BGe)), "loglik-cg" (Log-Likelihood (cond. Gauss.)), "aic-cg" (AIC (cond. Gauss.)), "bic-cg" (BIC (cond. Gauss.)). See ?bnlearn-package for details.
```

### Testing score equivalence

Arcs whose direction does not influence the v-structures present in the network structure are said to be score equivalent, because their reversal does not alter the score of the network (with the notable exceptions of K2 and BDe/BGe with prior other than the uniform). Usually these arcs are not oriented in the networks learned with constraint-based algorithms.

In `gaussian.test` the only score equivalent arc is `D``B`, as shown below.

```> eq.net = set.arc(gauss.net, "D", "B")
> score(gauss.net, gaussian.test, type = "bic-g")
```
```[1] -53191.54
```
```> score(eq.net, gaussian.test, type = "bic-g")
```
```[1] -53191.54
```

Equivalently, we can check that `gauss.net` and `eq.net` belong to the same equivalence class (represented by a CPDAG, see the manual).

```> all.equal(cpdag(gauss.net), cpdag(eq.net))
```
```[1] TRUE
```

Changing the direction of any other arc alters the score of the network; for example reversing `B``C` lowers the BIC by `6788.81`.

```> noneq1.net = set.arc(gauss.net, from = "B", to = "C")
> score(noneq1.net, gaussian.test, type = "bic-g")
```
```[1] -53191.54
```
```> noneq2.net = set.arc(gauss.net, from = "C", to = "B")
> score(noneq2.net, gaussian.test, type = "bic-g")
```
```[1] -59980.35
```

The `choose.direction()` function (manual) reaches the same conclusion; in the former case the network score is not altered, while in the latter the direction `B``C` results in a higher network score. We can read hoe the comparisons are made by enabling debugging output with `debug = TRUE`.

```> choose.direction(gauss.net, data = gaussian.test, c("B", "D"),
+   criterion = "bic-g", debug= TRUE)
```
```* testing B - D for direction.
> initial score for node D is -12762.57 .
> score delta for arc B -> D is 11242.4 .
> initial score for node B is -7034.42 .
> score delta for arc D -> B is 7636.76 .
@ arc B -> D is better .
```
```> choose.direction(gauss.net, data = gaussian.test, c("B", "C"),
+    criterion = "bic-g", debug = TRUE)
```
```* testing B - C for direction.
> initial score for node C is -16131.65 .
> score delta for arc B -> C is 12402.44 .
> initial score for node B is 496.3417 .
> score delta for arc C -> B is 105.9992 .
@ arc B -> C is better .
```
Last updated on `Tue Nov 14 16:02:10 2017` with bnlearn `4.3-20171001` and `R version 3.0.2 (2013-09-25)`.