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:

  1. maximize, the label of a score-based structure learning learning algorithm, and maximize.args, a list containing its arguments (other than the data);
  2. fit, the label of a parameter estimator in bn.fit() (documented here), and fit.args, a list containing its arguments (other than the data);
  3. 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)

Loading required namespace: gRain


Attaching package: 'gRbase'

The following objects are masked from 'package:bnlearn':

    ancestors, children, nodes, parents

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

## 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:  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")

## 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:                     Exact Inference 
  penalization coefficient:              4.258597 
  tests used in the learning procedure:  94 
  optimized:                             TRUE
Last updated on Fri Dec 9 03:12:25 2022 with bnlearn 4.9-20221107 and R version 4.2.2 (2022-10-31).