## Score-based structure learning from data with missing values

Score-based algorithms in the literature are typically defined to use a generic score function to compare different network structures. However, for the most part network scores assume that data are complete.

### The Structural Expectation-Maximization (Structural EM) algorithm

A possible approach to sidestep this limitation is the *Structural EM* algorithm from Nir Friedman
(link), which scores candidate network structures on completed data by
iterating over:

- an
*expectation step*(E): in which we complete the data by imputing missing values from a fitted Bayesian network; - the
*maximization step*(M): in which we learn a Bayesian network by maximizing a network score over the completed data.

This algorithm is implemented in the `structural.em()`

function in **bnlearn**
(documented here). The arguments of `structural.em()`

reflect its modular nature:

`maximize`

, the label of a score-based structure learning learning algorithm, and`maximize.args`

, a list containing its arguments (other than the data);`fit`

, the label of a parameter estimator in`bn.fit()`

(documented here), and`fit.args`

, a list containing its arguments (other than the data);`impute`

, the label an imputation method in`impute()`

(documented here), and`impute.args`

, a list containing its arguments (other than the data).

The number of iterations of the E and M steps is controlled by the `max.iter`

argument, which defaults to
5 iterations.

#### With partially observed variables

Consider some simple MCAR data in which 5% of the values are missing for each variable.

> incomplete.data = learning.test > for (col in seq(ncol(incomplete.data))) + incomplete.data[sample(nrow(incomplete.data), 100), col] = NA

With the default arguments, `structural.em()`

uses hill-climbing as the structure learning algorithm,
maximum likelihood for estimating the parameters of the Bayesian network, and likelihood weighting to impute the
missing values.

> dag = structural.em(incomplete.data) > dag

Bayesian network learned from Missing Data 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: Structural EM score-based method: Hill-Climbing parameter learning method: Maximum Likelihood (disc.) imputation method: Posterior Expectation (Likelihood Weighting) penalization coefficient: 4.258597 tests used in the learning procedure: 148 optimized: TRUE

We can change that using the arguments listed above.

> dag = structural.em(incomplete.data, + maximize = "tabu", maximize.args = list(tabu = 50, max.tabu = 50), + fit = "bayes", fit.args = list(iss = 1), + impute = "exact", max.iter = 3)

> dag

Bayesian network learned from Missing Data 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: Structural EM score-based method: Tabu Search parameter learning method: Bayesian Dirichlet imputation method: Exact Inference penalization coefficient: 4.258597 tests used in the learning procedure: 1627 optimized: TRUE

In particular, changing the default `impute = "bayes-lw"`

into `impute = "exact"`

may be
useful because the convergence of the Structural EM is not guaranteed if the imputation is performed using approximate
Monte Carlo inference. However, it is usually much slower.

In addition, we can set the argument `return.all`

to `TRUE`

to have
`structural.em()`

return its complete status at the last iteration: the network structure it has learned,
the completed data it was learned from and the fitted Bayesian network used to perform the imputation.

> info = structural.em(incomplete.data, return.all = TRUE, + maximize = "tabu", maximize.args = list(tabu = 50, max.tabu = 50), + fit = "bayes", fit.args = list(iss = 1), + impute = "exact", max.iter = 3) > names(info)

[1] "dag" "imputed" "fitted"

The network structure is the same as that returned when `return.all = FALSE`

, which is the default.

> info$dag

Bayesian network learned from Missing Data 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: Structural EM score-based method: Tabu Search parameter learning method: Bayesian Dirichlet imputation method: Exact Inference penalization coefficient: 4.258597 tests used in the learning procedure: 1627 optimized: TRUE

The completed data are stored in a data frame with the same structure as the original data.

> head(info$imputed)

A B C D E F 1 b c b a b b 2 b a c a b b 3 a a a a a a 4 a a a a b b 5 a a b c a a 6 c c a c c a

The fitted Bayesian network is a `bn.fit`

object.

> info$fitted

Bayesian network parameters Parameters of node A (multinomial distribution) Conditional probability table: a b c 0.3328001 0.3353996 0.3318003 Parameters of node B (multinomial distribution) Conditional probability table: A B a b c a 0.85746712 0.44481982 0.11758404 b 0.02289872 0.21707737 0.09227267 c 0.11963415 0.33810281 0.79014330 Parameters of node C (multinomial distribution) Conditional probability table: a b c 0.74551756 0.20382590 0.05065654 Parameters of node D (multinomial distribution) Conditional probability table: , , C = a A D a b c a 0.80793106 0.08461902 0.10756427 b 0.08708492 0.81320565 0.10274929 c 0.10498403 0.10217533 0.78968643 , , C = b A D a b c a 0.17334596 0.89455451 0.22772841 b 0.12790575 0.06149183 0.51686033 c 0.69874829 0.04395366 0.25541126 , , C = c A D a b c a 0.45766488 0.34145241 0.11391341 b 0.19295900 0.40234551 0.47709121 c 0.34937611 0.25620207 0.40899538 Parameters of node E (multinomial distribution) Conditional probability table: , , F = a B E a b c a 0.81880294 0.20615269 0.10985155 b 0.09354673 0.17238792 0.11082321 c 0.08765033 0.62145939 0.77932524 , , F = b B E a b c a 0.38949250 0.35905681 0.23665660 b 0.50039979 0.33976420 0.51428420 c 0.11010771 0.30117899 0.24905920 Parameters of node F (multinomial distribution) Conditional probability table: a b 0.5023995 0.4976005

#### With completely unobserved (latent) variables

If the data contain a latent variable which we do not observe for any observation, the E step in the fist iteration
fails because it cannot fit a Bayesian network to impute the missing values. (If all variables are at least partially
observed, `structural.em()`

uses locally complete observations for this purpose. `bn.fit()`

does
the same as illustrated here.)

> incomplete.data[, "A"] = factor(rep(NA, nrow(incomplete.data)), levels = levels(incomplete.data[, "A"])) > structural.em(incomplete.data)

## Warning in check.data(x, allow.levels = TRUE, allow.missing = TRUE, ## warn.if.no.missing = TRUE, : at least one variable has no observed values.

## Error: the data contain latent variables, so the 'start' argument must be a 'bn.fit' object.

As the error message suggests, we can side-step this issue by providing a `bn.fit`

object ourselves via
the `start`

argument: it will be used to perform the initial imputation.

> start.dag = empty.graph(names(incomplete.data)) > cptA = matrix(c(0.3336, 0.3340, 0.3324), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptB = matrix(c(0.4724, 0.1136, 0.4140), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptC = matrix(c(0.7434, 0.2048, 0.0518), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptD = matrix(c(0.351, 0.314, 0.335), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptE = matrix(c(0.3882, 0.2986, 0.3132), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptF = matrix(c(0.5018, 0.4982), ncol = 2, dimnames = list(NULL, c("a", "b"))) > start = custom.fit(start.dag, list(A = cptA, B = cptB, C = cptC, D = cptD, E = cptE, F = cptF)) > dag = structural.em(incomplete.data, start = start, max.iter = 3)

## Warning in check.data(x, allow.levels = TRUE, allow.missing = TRUE, ## warn.if.no.missing = TRUE, : at least one variable has no observed values.

## Warning in check.data(x, allow.missing = TRUE): variable A has levels that are ## not observed in the data.

> dag

Bayesian network learned from Missing Data model: [A][B][C][F][D|B:C][E|B:F] nodes: 6 arcs: 4 undirected arcs: 0 directed arcs: 4 average markov blanket size: 2.00 average neighbourhood size: 1.33 average branching factor: 0.67 learning algorithm: Structural EM score-based method: Hill-Climbing parameter learning method: Maximum Likelihood (disc.) imputation method: Posterior Expectation (Likelihood Weighting) penalization coefficient: 4.258597 tests used in the learning procedure: 83 optimized: TRUE

Unfortunately, the latent variable will almost certainly end up as an isolated nodes unless we connect it to at
least some nodes that are partially observed: the noisiness of Monte Carlo inference can easily overwhelm the
dependence relationships we encode in the network in `start`

argument.

> start.dag = model2network("[A][B|A][C][D][E][F]") > cptA = matrix(c(0.3336, 0.3340, 0.3324), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptB = matrix(c(0.856, 0.025, 0.118, 0.444, 0.221, 0.334, 0.114, 0.094, 0.790), nrow = 3, ncol = 3, + dimnames = list(B = c("a", "b", "c"), A = c("a", "b", "c"))) > start = custom.fit(start.dag, list(A = cptA, B = cptB, C = cptC, D = cptD, E = cptE, F = cptF)) > dag = structural.em(incomplete.data, start = start, max.iter = 3)

## Warning in check.data(x, allow.levels = TRUE, allow.missing = TRUE, ## warn.if.no.missing = TRUE, : at least one variable has no observed values.

> dag

Bayesian network learned from Missing Data model: [A][C][F][B|A][D|A:C][E|A: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: Structural EM score-based method: Hill-Climbing parameter learning method: Maximum Likelihood (disc.) imputation method: Posterior Expectation (Likelihood Weighting) penalization coefficient: 4.258597 tests used in the learning procedure: 94 optimized: TRUE

Passing a whitelist to the structure learning algorithm is the simplest way to do that.

> start.dag = model2network("[A][B|A][C][D][E][F]") > cptA = matrix(c(0.3336, 0.3340, 0.3324), ncol = 3, dimnames = list(NULL, c("a", "b", "c"))) > cptB = matrix(c(0.856, 0.025, 0.118, 0.444, 0.221, 0.334, 0.114, 0.094, 0.790), nrow = 3, ncol = 3, + dimnames = list(B = c("a", "b", "c"), A = c("a", "b", "c"))) > start = custom.fit(start.dag, list(A = cptA, B = cptB, C = cptC, D = cptD, E = cptE, F = cptF)) > dag = structural.em(incomplete.data, + maximize.args = list(whitelist = data.frame(from = "A", to = "B")), + start = start, max.iter = 3)

> dag

Bayesian network learned from Missing Data model: [A][C][F][B|A][D|A:C][E|A: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: Structural EM score-based method: Hill-Climbing parameter learning method: Maximum Likelihood (disc.) imputation method: Posterior Expectation (Likelihood Weighting) penalization coefficient: 4.258597 tests used in the learning procedure: 91 optimized: TRUE

Exact inference does not have this issue because it has no stochastic noise: the imputed values are deterministic given the observed values in each observation.

> dag = structural.em(incomplete.data, start = start, max.iter = 3, impute = "exact")

> dag

Bayesian network learned from Missing Data model: [A][C][F][B|A][D|A:C][E|A: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: Structural EM score-based method: Hill-Climbing parameter learning method: Maximum Likelihood (disc.) imputation method: Exact Inference penalization coefficient: 4.258597 tests used in the learning procedure: 94 optimized: TRUE

`Fri Dec 9 03:12:25 2022`

with **bnlearn**

`4.9-20221107`

and `R version 4.2.2 (2022-10-31)`

.