## Score functions: computing & comparing

### Loading the reference data sets

First, we load the `learning.test`

and `gaussian.test`

data sets shipped
with **bnlearn**.

> library(bnlearn)

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

`Tue Nov 14 16:02:10 2017`

with **bnlearn**

`4.3-20171001`

and `R version 3.0.2 (2013-09-25)`

.