Bayesian Network Models for Incomplete and Dynamic Data

TL;DR

This paper reviews how Bayesian networks can effectively model complex, dynamic, and incomplete data, emphasizing their interpretability and applicability in real-world scenarios.

Contribution

It provides a comprehensive overview of Bayesian network models tailored for dynamic and incomplete data, extending their traditional static, fully observed data applications.

Findings

Bayesian networks effectively model dynamic data.

They handle incomplete observations well.

Their interpretability benefits real-world applications.

Abstract

Bayesian networks are a versatile and powerful tool to model complex phenomena and the interplay of their components in a probabilistically principled way. Moving beyond the comparatively simple case of completely observed, static data, which has received the most attention in the literature, in this paper we will review how Bayesian networks can model dynamic data and data with incomplete observations. Such data are the norm at the forefront of research and in practical applications, and Bayesian networks are uniquely positioned to model them due to their explainability and interpretability.

Figures

Equations

$$P(\mathbf{X}|\mathcal{G},\Theta )=\prod _{i=1}^{N}P\left(X_i|{\Pi }_{X_i},{\Theta }_{X_i}\right),$$

$$X_i|{\Pi }_{X_i}\sim Mul\left({\pi }_{ik|j}\right),$$

$$X_i={\mu }_{X_i}+{\Pi }_{X_i}{\mathbf{\beta }}_{X_i}+{\epsilon }_{X_i},$$

$$X_i|{\Pi }_{X_i}\sim N\left({\mu }_{X_i,{\delta }_{X_i}}+{\Gamma }_{X_i}{\mathbf{\beta }}_{X_i,{\delta }_{X_i}},{\sigma }_{X_i,{\delta }_{X_i}}^2\right)$$

$$X_i={\mu }_{X_i,{\delta }_{X_i}}+{\Gamma }_{X_i}{\mathbf{\beta }}_{X_i,{\delta }_{X_i}}+{\epsilon }_{X_i,{\delta }_{X_i}},$$

$$\underset{\text{learning}}{\underbrace{P(\mathcal{G},\Theta |\mathcal{D})}}=\underset{\text{structure learning}}{\underbrace{P(\mathcal{G}|\mathcal{D})}}\cdot \underset{\text{parameter learning}}{\underbrace{P(\Theta |\mathcal{G},\mathcal{D})}}.$$

$$P(\Theta |\mathcal{G},\mathcal{D})=\prod _{i=1}^{N}P\left({\Theta }_{X_i}|{\Pi }_{X_i},\mathcal{D}\right).$$

$$P(\mathcal{G}|\mathcal{D})\propto P(\mathcal{G})P(\mathcal{D}|\mathcal{G}),$$

$$P(\mathcal{D}|\mathcal{G})=\int P(\mathcal{D}|\mathcal{G},\Theta )P(\Theta |\mathcal{G})d\Theta =\prod _{i=1}^N\int P\left(X_i|{\Pi }_{X_i},{\Theta }_{X_i}\right)P\left({\Theta }_{X_i}|{\Pi }_{X_i}\right)d{\Theta }_{X_i};$$

$$\mathrm{BIC}(\mathcal{G},\Theta |\mathcal{D})=\sum _{i=1}^N\log P\left(X_i|{\Pi }_{X_i},{\Theta }_{X_i}\right)-\frac{\log (n)}{2}\left|{\Theta }_{X_i}\right|.$$

$$P\left({\mathbf{X}}^{(0)},\dots ,{\mathbf{X}}^{(T)}\right)=P\left({\mathbf{X}}^{(0)}\right)\prod _{t=1}^{T}P\left({\mathbf{X}}^{(t)}|{\mathbf{X}}^{(t-1)}\right).$$

$$P\left({\mathbf{X}}^{(t)}|{\mathbf{X}}^{(t-1)}\right)=\prod _{i=1}^{N}P\left(X_i^{(t)}|{\Pi }_{X_i^{(t)}}\right),$$

$$P\left({\mathbf{X}}^{(t)}|{\mathbf{X}}^{(t-1)},{\mathbf{Z}}^{(t)}\right)=\prod _{i=1}^{N}P\left(X_i^{(t)}|{\Pi }_{X_i^{(t)}},{\mathbf{Z}}^{(t)}\right)$$

$${\mathbf{X}}^{(t)}=A_1{\mathbf{X}}^{(t-1)}+\dots +A_L{\mathbf{X}}^{(t-L)}+{\epsilon }_t,$$

$${\mathbf{X}}^{(t)}|{\mathbf{X}}^{(t-1)},\dots ,{\mathbf{X}}^{(t-L)}\sim N\left(A_1{\mathbf{X}}^{(t-1)}+\dots +A_L{\mathbf{X}}^{(t-L)},{\epsilon }_t\right)$$

$${\mathbf{Z}}^{(t)}=A{\mathbf{Z}}^{(t-1)}+B{\mathbf{U}}^{(t)}+{\zeta }_t,$$

$${\mathbf{X}}^{(t)}=C{\mathbf{Z}}^{(t)}+D{\mathbf{U}}^{(t)}+{\epsilon }_t,$$

$$P\left(\mathcal{D},\mathbf{M}|\mathcal{G},\Theta ,\Xi \right)=P\left({\mathcal{D}}^O,{\mathcal{D}}^M,\mathbf{M}|\mathcal{G},\Theta ,\Xi \right)=\int P({\mathcal{D}}^O,{\mathcal{D}}^M|\mathcal{G},\Theta )P\left(\mathbf{M}|{\mathcal{D}}^O,{\mathcal{D}}^M,\mathcal{G},\Xi \right)d{\mathcal{D}}^M;$$

$$P\left(\mathbf{M}|{\mathcal{D}}^O,{\mathcal{D}}^M,\mathcal{G},\Xi \right)=P\left(\mathbf{M}|{\mathcal{D}}^O,\mathcal{G},\Xi \right);$$

$$P\left(\mathbf{M}|{\mathcal{D}}^O,{\mathcal{D}}^M,\mathcal{G},\Xi \right)=P\left(\mathbf{M}|\mathcal{G},\Xi \right).$$

$$\begin{gathered} P\left(\mathcal{D}|\mathcal{G}\right)=P\left({\mathcal{D}}^O,{\mathcal{D}}^M|\mathcal{G}\right)=\int P\left({\mathcal{D}}^O,{\mathcal{D}}^M|\mathcal{G},\Theta \right)P\left(\Theta |\mathcal{G}\right)d\Theta = \\ =\int \underset{\text{missing data}}{\underbrace{P\left({\mathcal{D}}^M|{\mathcal{D}}^O,\mathcal{G},\Theta \right)}}\underset{\text{observed data}}{\underbrace{P\left({\mathcal{D}}^O|\mathcal{G},\Theta \right)}}\underset{\text{averaging over parameters}}{\underbrace{P\left(\Theta |\mathcal{G}\right)d\Theta }}. \end{gathered}$$

$$P\left(\mathcal{D}|\mathcal{G}\right)=\iint P\left({\mathcal{D}}^M|{\mathcal{D}}^O,\mathcal{G},\Theta \right)P\left({\mathcal{D}}^O|\mathcal{G},\Theta \right)P\left(\Theta |\mathcal{G}\right)d\Theta d{\mathcal{D}}^M.$$

Full text

Keywords — Bayesian networks, dynamic data, incomplete data, structure learning, inference

1 Introduction

Bayesian networks (BNs) are a probabilistic graphical model that is ideally suited to study complex phenomena through the interactions of their components. By representing the latter as nodes in a graph, and connecting them with arcs that encode how those components interact with each other, they provide a high-level qualitative view of the phenomenon under investigation built on strong theoretical foundations.

BNs were initially developed in the 1980s in the context of expert systems, with the seminal work summarised in Pearl [87] and Neapolitan [81]. Later monographs by Castillo et al. [16], Neapolitan [82], Cowell et al. [24], Korb and Nicholson [59] provide extensive coverage of the properties of BNs and their use; with Koller and Friedman [58] being the most complete reference book up to date. In that context, “expert systems” were conceived as the combination of a formal representation of domain-specific knowledge gathered from subject matter experts and an “inference engine” that could use that representation to answer arbitrary queries. BNs can represent complex phenomena in a modular way due the conditional independence assumptions they encode, allowing them to scale beyond trivial problems and to be used to develop efficient inference algorithms.

As noted in Cowell et al. [24], by the end of the 1990s the limitations of constructing expert systems had become more and more apparent and were pushing research towards learning the graphical structure and the parameters of BNs from data. For many phenomena that were and have since been of interest in the scientific community, expert systems are too complex for the modeller to be able to elicit the required domain knowledge from experts, even assuming expert knowledge is available in the first place (or with sufficient consensus to be trusted). However, at the same time the availability of increasingly large data sets made it more attractive to use them as an alternative source of information to construct BNs.

Given the limited computational power available at the time, this transition from expert elicitation to learning from data was made possible by making a number of simplifying assumptions on the nature of the data [52]. Two assumptions in particular have been made in the vast majority of the literature since then. Firstly, that data are complete; that is, that all variables are completely observed for all samples. And secondly, that samples are independent from each other; which specifically excludes dynamic data such as time series in which samples are associated with each other.

However, problems at the forefront of research and in practical applications have only been increasing in complexity, making these assumptions too strong in many settings. In this paper, we will review how BNs can be used to model dynamic data and data with incomplete observations, providing several examples in which they have been used successfully to develop cutting-edge applications. The comparatively simple case of complete, static data has certainly received the most attention in the literature, for historical reasons and for its probabilistic and computational simplicity; but BNs are by no means limited to modelling only such data.

The contents of this review are organised as follows. In Section 2 we introduce to BNs in a general setting, including their underlying distributional assumptions and available learning approaches. In Section 3 we discuss BNs targeted at dynamic (that is, temporal) data, starting from their definition and main properties (Section 3.1) and then showing them to be a general approach that subsumes several other models in the literature (Section 3.2). In Section 4 we discuss incomplete data and how they can be correctly used for parameter (Section 4.1) and structure learning. Finally, in Section 5 we provide an overview of available software implementations (Section 5.1) and real-world applications (Section 5.2) of the BNs we discussed.

2 Background and Notation

BNs are a class of graphical models in which the nodes of a directed acyclic graph (DAG) $\mathcal{G}$ represent a set $\mathbf{X}=\{X_{1},\ldots,X_{N}\}$ of random variables describing some quantities of interest. The arcs connecting those nodes express direct dependence relationships, with graphical separation in $\mathcal{G}$ implying conditional independence in probability. As a result, $\mathcal{G}$ induces the factorisation

[TABLE]

in which the joint probability distribution of $\mathbf{X}$ (with parameters $\Theta$) decomposes into one local distribution for each $X_{i}$ (with parameters $\Theta_{X_{i}}$, $\bigcup_{X_{i}\in\mathbf{X}}\Theta_{X_{i}}=\Theta$) conditional on its parents $\Pi_{X_{i}}$. Assuming $\mathcal{G}$ is sparse111There is no universally accepted threshold on the number of arcs for a DAG to be called “sparse”; typically it is taken to have $O(cN)$ arcs, with $c$ between $1$ and $5$., BNs provide a compact representation of both low- and high-dimensional probability distributions.

BNs are also very flexible in terms of distributional assumptions; but while in principle we could choose any probability distribution for $\mathbf{X}$, the literature has mostly focused on three cases for analytical and computational reasons. Discrete BNs [52] assume that both $\mathbf{X}$ and the $X_{i}$ are multinomial random variables. Local distributions take the form

[TABLE]

their parameters $\pi_{ik\operatorname{|}j}$ are the conditional probabilities of $X_{i}$ given each configuration of the values of its parents, usually represented as a conditional probability table for each $X_{i}$. Gaussian BNs [GBNs; 42] model $\mathbf{X}$ with a multivariate normal random variable and assume that the $X_{i}$ are univariate normals linked by linear dependencies. The parameters of the local distributions can be equivalently written [119] as the partial correlations $\rho_{X_{i},X_{j}\operatorname{|}\Pi_{X_{i}}\setminus X_{j}}$ between $X_{i}$ and each parent $X_{j}$ given the other parents; or as the coefficients $\bm{\beta}_{X_{i}}$ of the linear regression model

[TABLE]

so that $X_{i}\operatorname{|}\Pi_{X_{i}}\sim N\left(\mu_{X_{i}}+\Pi_{X_{i}}\bm{\beta}_{X_{i}},\sigma^{2}_{X_{i}}\right)$. Finally, conditional linear Gaussian BNs [CLGBNs; 63] combine discrete and continuous random variables in a mixture model:

• •

discrete $X_{i}$ are only allowed to have discrete parents (denoted $\Delta_{X_{i}}$), and are assumed to follow a multinomial distribution as in (2);

• •

continuous $X_{i}$ are allowed to have both discrete and continuous parents (denoted $\Gamma_{X_{i}}$, $\Delta_{X_{i}}\cup\Gamma_{X_{i}}=\Pi_{X_{i}}$), and their local distributions are

[TABLE]

which can be written as a mixture of linear regressions

[TABLE]

against the continuous parents with one component for each configuration $\delta_{X_{i}}$ of the discrete parents $\Delta_{X_{i}}$. If $X_{i}$ has no discrete parents, the mixture reverts to a single linear regression like that in (3).

The task of learning a BN from a data set $\mathcal{D}$ containing $n$ observations is performed in two steps:

[TABLE]

Structure learning consists in finding the DAG $\mathcal{G}$ that encodes the dependence structure of the data, thus maximising $\operatorname{P}(\mathcal{G}\operatorname{|}\mathcal{D})$ or some alternative goodness-of-fit measure; parameter learning consists in estimating the parameters $\Theta$ given the $\mathcal{G}$ obtained from structure learning. Both steps can integrate data with expert knowledge through the use of suitable prior distributions on $\mathcal{G}$ and $\Theta$ [see for example 15, 77, 30]. If we assume that parameters in different local distributions are independent and that the data contain no missing values [52], we can perform parameter learning independently for each $X_{i}$ because following (1)

[TABLE]

Furthermore, assuming $\mathcal{G}$ is sparse, each local distribution $X_{i}\operatorname{|}\Pi_{X_{i}}$ will involve only a few variables and thus will have a low-dimensional parameter space, making parameter learning computationally efficient.

On the other hand, structure learning is well known to be both NP-hard [19] and NP-complete [17], even under unrealistically favourable conditions such as the availability of an independence and inference oracle [21]222Interestingly, some relaxations of BN structure learning are not NP-hard; see for example Dojer [29] on learning the structure of causal networks.. This is despite the fact that if we take

[TABLE]

following (1) we can decompose the marginal likelihood $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ into one component for each local distribution

[TABLE]

and despite the fact that each component can be written in closed form for discrete BNs [52], GBNs [42] and CLGBNs [11]. The same is true if we replace $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ with frequentist goodness-of-fit scores such as BIC [105], which is commonly used in structure learning because of its simple expression:

[TABLE]

Compared to marginal likelihoods, BIC also has the advantage that it does not depend on any hyperparameter, while converging to $\log\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ as $n\to\infty$. These score functions have two important properties:

• •

they allow local computations because, following (1), they decompose into one component for each local distribution;

• •

they take the same value for all the DAGs that encode the same probability distribution (score equivalence), which can then be grouped in equivalence classes [20].333All DAGs in the same equivalence class have the same underlying undirected graph and v-structures (patterns of arcs like $X_{i}\rightarrow X_{j}\leftarrow X_{k}$, with no arcs between $X_{i}$ and $X_{k}$).

Structure learning via score maximisation is usually based on general-purpose heuristic optimisation algorithms, adapted to take advantage of these two properties to increase the speed of structure learning [109]. The most common are greedy search strategies such as hill-climbing and tabu search [99] that employ local moves designed to affect only one or two local distributions in each iteration; other options explored in the literature include genetic algorithms [61] and ant colony optimisation [27]. Learning equivalence classes directly (as opposed to DAGs) can be done along the same lines with the Greedy Equivalence Search [GES; 18] algorithm. Exact maximisation of $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ and BIC has also become feasible in recent years thanks to increasingly efficient pruning of the space of DAGs and tight bounds on the scores [25, 114, 102].

Another option for structure learning is using conditional independence tests to learn conditional independence constraints from $\mathcal{D}$, and thus to identify which arcs should be included in $\mathcal{G}$. The resulting algorithms are called constraint-based algorithms, as opposed to the score-based algorithms we introduced in the previous paragraph; for an overview and a comparison of these two approaches see Scutari et al. [108]. Chickering et al. [21] proved that constraint-based algorithms are also NP-hard for unrestricted DAGs; and they are in fact equivalent to score-based algorithms given a fixed topological ordering of the nodes in $\mathcal{G}$ when independence constraints are assessed with statistical tests related to cross-entropy [23].

Finally, once both $\mathcal{G}$ and $\Theta$ have been learned, we can answer queries about our quantities of interest using the resulting BN as our model of the world. Common types are conditional probability queries, in which we compute the posterior probability of some variables given evidence on others; and most probable explanation queries, in which we identify the configuration of values of some variables that has the highest posterior probability given the values of some other variables. The latter is especially suited to implement both prediction and imputation of missing data. These queries can be automated, for any given BN, using either exact or approximate inference algorithms that work directly on the BN without the need for any manual calculation; for an overview of such algorithms see Scutari and Denis [107].

Given their ability to concisely but meaningfully represent the world, to automatically answer arbitrary queries, and to combine data and expert knowledge in the learning process, BNs can model a wide variety of phenomena effectively. However, their applicability is not always apparent to practitioners in other fields due to the strong focus of the literature on the simple scenario in which data are static (as opposed to dynamic, that is, with a time dimension) and complete (as in, completely observed). Dynamic data are central to a number of cutting-edge applications and research in fields as different as genetics and robotics; and incomplete data are a fact of life in almost any real-world data analysis. BNs can handle both in a rigorous way, as we will see in the following.

3 Dynamic Bayesian Networks

Dynamic BNs (DBNs444Confusingly, discrete BNs are sometimes called DBNs to be consistent with Gaussian BNs being called GBNs.) combine classic (static) BNs and Markov processes to model dynamic data in which each individual is measured repeatedly over time, such as longitudinal or panel data. An approachable introduction is provided for instance in Murphy [78]. They have major applications in engineering [86, 39], medicine [118], genetics and systems biology [89]. The term “dynamic” in this context implies we are modelling a dynamic system, not necessarily that the network changes over time.

3.1 Definitions and Properties

For simplicity, let’s assume at first that we are operating in discrete time: our system consists of one set $\mathbf{X}^{(t)}$ of random variables for each of $t=1,\ldots,T$ time points. We can model it as a DBN with a Markov process of the form

[TABLE]

where $\operatorname{P}(\mathbf{X}^{(0)})$ gives the initial state of the process and $\operatorname{P}(\mathbf{X}^{(t)}\operatorname{|}\mathbf{X}^{(t-1)})$ defines the transition between times $t-1$ and $t$. We can model this transition with a 2-time BN (2TBN) defined over $(\mathbf{X}^{(t-1)},\mathbf{X}^{(t)})$, in which we naturally assume that any arc between a node in $t-1$ and a node in $t$ must necessarily be directed towards the latter following the arrow of time. When modelling $\mathbf{X}^{(t)}$, the nodes in $\mathbf{X}^{(t-1)}$ only appear in the conditioning; we take them to be essentially fixed and to have no free parameters, so we leave them as root nodes. After all, $\mathbf{X}^{(t-1)}$ will be stochastic in $\operatorname{P}(\mathbf{X}^{(t-1)}\operatorname{|}\mathbf{X}^{(t-2)})$ and it would not be consistent with (6) to treat $\mathbf{X}^{(t-1)}$ as a stochastic quantity twice! Then, following (1) we can write

[TABLE]

and we usually assume that the parameters associated with the local distributions do not change over time to make the process time-homogeneous.

These choices are motivated by computational and statistical simplicity: there is no intrinsic limitation in the construction of DBNs that prevents them from modelling dependencies that stretch further back in time or, for that matter, trends or seasonality. (See Figure 1 for an overview of the models we will discuss in detail below.) Grzegorczyk and Husmeier [47, 48], for instance, constructed non-stationary, non-homogeneous DBNs for modelling continuous data using change-points to capture different regimes in different time intervals; they balanced this increased flexibility by requiring $\mathcal{G}$ to be the same for all regimes, allowing only the parameters of the local distributions to change. This restriction was later relaxed in Kamalabad and Grzegorczyk [56] by allowing both $\mathcal{G}$ and the parameters to vary but coupling them between adjacent time intervals. The hyperparameters controlling the coupling were allowed to vary between segments, following an inverse Gamma hyperprior. Robinson and Hartemink [95] also defined a non-stationary DBN with change-points and associated arc changesets, with truncated geometric priors on the size of the changesets, number of and interval between change-points; and Song et al. [112] constructed DBNs that vary smoothly (not piecewise) over time in both structure and parameters. Kim et al. [57] introduced an even more flexible model that used spline regression with B-splines to identify $\Pi_{X_{i}^{(t)}}$. Augmenting temporal (panel) data with non-temporal (cross-sectional) data for learning DBNs has been explored as well by Lähdesmäki and Shmulevich [60] using an score that approximates the resulting intractable likelihood.

Modelling DBNs as discrete time processes is likewise a choice motivated by mathematical simplicity. Nodelman et al. [83] originally proposed a class of continuous-time DBNs (CTBNs) with independent exponential waiting times and discrete nodes; they are uniquely identifiable since all arcs are non-instantaneous, and they have a closed-form marginal likelihood as well. More recently, this work has been expanded in Liu et al. [68] by replacing exponential waiting times with hypoexponentials to better reflect the behaviour of data in several domains. Discrete-time DBNs are certainly simpler than any of these CTBNs, but that mathematical simplicity comes with important practical consequences. Firstly, in order to work in discrete time we must choose a uniform time step (the length of time between $t-1$ and $t$) for the whole DBN; but in many real-world phenomena different variables can have very different time granularities, and those time granularities may vary as well in the course of data collection, making any single choice for the time step inappropriate. Secondly, the choice of the time step may obscure the dynamics of the phenomenon. The implication of using discrete time is that we aggregate all the state changes in the DBN over the entire course of each time step. On the one hand, if variables evolve at slower pace than the time step we are forced to model the DBN as a higher-order Markov process555A stochastic process is a Markov process of order $L$ if $\mathbf{X}^{(t)}$ depends only on $\mathbf{X}^{(t-1)},\ldots\mathbf{X}^{(t-L)}$ and is independent from $\mathbf{X}^{(t-L-1)},\ldots,\mathbf{X}^{(0)}$. The higher the order, the further back in time the dependencies can reach., resulting in a much more complex model. On the other hand, if variables evolve at a faster pace than time step, this averaging will effectively hide state changes and their interplay into a single summary statistic; hence the DBN will provide a very poor approximation of the underlying phenomenon. Furthermore, if we are unable to correctly identify the $\Pi_{X_{i}^{(t)}}$, $X_{i}^{(t)}$ may end up being recursively linked to the parents of the $\Pi_{X_{i}^{(t)}}$666This issue is often called entanglement., resulting into dense DBNs that are much more complex than the real underlying phenomenon and that are difficult to learn from limited data.

A partial solution to the latter problem is to allow the parents $\Pi_{X_{i}^{(t)}}$ to be either in the same time or in the previous time point, modelling the DBN as a first-order Markov process. When observations represent average or aggregate measurements over a period of time (say, the $t$th time point corresponds to the $t$th week’s worth of data), it makes sense to allow instantaneous dependencies between variables in the same time point, since the instantaneousness of the dependence is just a fiction arising from our model definition. On the other hand, when observations actually correspond to instantaneous measurements (such as from synchronised sensors) only non-instantaneous dependencies are usually allowed in the model, on the grounds that conditioning event ($\Pi_{X_{i}^{(t)}}=\pi_{X_{i}^{(t)}}$) should precede in time the conditioned event ($X_{i}^{(t)}=x_{i}^{(t)}$). From a causal perspective, we can similarly argue that each of the $\Pi_{X_{i}^{(t)}}$ can only cause $X_{i}^{(t)}$ if it precedes $X_{i}^{(t)}$ in time; if that $\Pi_{X_{i}^{(t)}}$ is in the same time point as $X_{i}^{(t)}$ then what we are modelling is co-occurrence and not causation. This the core idea of Granger causality [46], which states that one time series (such as the $\Pi_{X_{i}^{(t)}}$) can be said to have a causal influence on a second time series (such as the $X_{i}^{(t)}$) if and only if incorporating past knowledge about the former improves predictive accuracy for the latter. Therefore, allowing instantaneous dependencies makes causal reasoning on the DBN markedly more difficult. The same is true for learning the structure of the DBN in the first place: learning a general BN is NP-hard [19, 21], while learning a DBN containing only non-instantaneous dependencies is not [29]. Intuitively, the space of the possible DBNs is much smaller if we do not allow instantaneous dependencies because there are fewer candidate arcs that we can include, and because their directions are fixed to follow the arrow of time and Granger causality777Interestingly, when the number of available time points is small, inference based on DBNs is more accurate than methods directly based on Granger causality; but the opposite is true for longer time series [123]..

3.2 Models That Can Be Expressed as DBNs

DBNs have a strong expressive power, and they subsume and generalise a variety of classic models that have been studied individually in the literature. Here we will discuss three such models: hidden Markov models, vector auto-regressive models and Kalman filters. Hidden Markov models and Kalman filters can be seen as particular instances of state-space models, and their DBN representations make clear the relationship between these two classes of models. Representing these models as DBNs has advantages beyond making their comparison more convenient; it makes it possible to use the probabilistic machinery we covered in Sections 2 and 3.1 when more convenient or computationally efficient than the alternatives [as in the case of inference in hidden Markov models; 78]. Some examples are shown in Figure 2, and will be illustrated below.

Hidden Markov models [HMMs; 124] are one of the most widespread approaches to model phenomena with hidden state, that is, in which the behaviour of the observed variables $\mathbf{X}$ depends on that of one or more discrete latent variables $\mathbf{Z}$ as well as on other variables in $\mathbf{X}$. This scenario commonly arises when technical, economic or ethical considerations make it unfeasible to completely observe the underlying state of the phenomenon of interest. Notable examples are imputation [73] and phasing [28] in genome-wide association studies, due to limitations in the technology to probe and tag DNA; tracking animals in ecology [85], where we can observe their movements using radio beacons but not their behaviour; and confirming mass migrations through history by combining archaeological artefacts and ancient DNA samples [104]. Until the advent of deep neural networks, HMMs were also the choice model for speech [40] and handwriting recognition [90].

In DBN terms, a typical HMM model with $M$ latent variables can be written as

[TABLE]

with the restriction that the parents of $Z_{j}^{(t)}$ can only be other latent variables. Latent variables are assumed to be discrete; and in the vast majority of the literature observed variables are assumed to be discrete as well. Depending on the choice of $\Pi_{Z_{j}^{(t)}}$, we can obtain various HMM variants such as hierarchical HMMs [34], in which each $Z_{j}^{(t)}$ is defined as an HMM itself to produce a multi-level stochastic model; and factorial HMMs [45], in which the $\mathbf{X}^{(t)}$ are driven by the configuration of a set of mutually independent $Z_{j}^{(t)}$ (shown in Figure 2, top-left panel).

Vector auto-regressive models [VARs; 12] are a straightforward multivariate extension of univariate auto-regressive time series for continuous variables. As such, their major applications are forecasting in finance [7] and more recently in the analysis of fMRI data [41]. They are defined as

[TABLE]

for some fixed Markov order $L$. We can rewrite (9) as

[TABLE]

and then restrict the parents of each $X_{i}^{(t)}$ to those for which the corresponding regression coefficients in $A_{1},\ldots,A_{L}$ are different from zero using the one-to-one correspondence between regression coefficients and partial correlations [119]. Formally, $X_{j}^{(t-l)}\in\Pi_{X_{i}^{(t)}}$ if and only if $A_{l}[i,j]\neq 0$, which makes it possible to write (9) in a similar form to (6) and obtain a Gaussian DBN. The same construction is used by Song et al. [112] for their non-homogeneous DBN models, which are parameterised as VAR processed and estimated using $L_{1}$-penalised regressions. In the special case in which $L=1$, VARs can be graphically represented in two equivalent ways shown in the bottom panels of Figure 2: an “unrolled” DBN in which each node corresponds to a single $X_{i}^{(t)}$; and a more compact “rolled-up” DBN in which each node corresponds to a variable $X_{i}$, an arc from $X_{i}$ to $X_{j}$ implies $X_{i}^{(t-1)}\to X_{j}^{(t)}$. An arc from $X_{i}$ to itself implies $X_{i}^{(t-1)}\to X_{i}^{(t)}$ as a special case for $i=j$.

Kalman filters [KFs; 50] combine traits of both HMMs and VARs, as discussed in depth in Roweis and Ghahramani [97] and Ghahramani [44]: like VARs, they are linear Gaussian DBNs; but they also have latent variables like HMMs. They are widely used for filtering (that is, denoising) and prediction in GPS positioning systems [120]; atmospheric modelling and weather prediction [14]; and seismology [101]. In their simplest form (see Figure 2, top-right panel), KFs include a layer of one or more latent variables that model the unobservable part of the phenomenon,

[TABLE]

feeding into one or more observed variables

[TABLE]

with independent Gaussian noise added in both layers. Both layers often include additional (continuous) explanatory variables $\mathbf{U}$ and can also be augmented with (discrete) switching variables to allow for different regimes as in Grzegorczyk and Husmeier [47, 48]. If we exclude the latter, the assumption is that the system is jointly Gaussian: that makes it possible to frame KFs as a DBN in the same way we did for VARs.

4 Bayesian Networks from Incomplete Data

The vast majority of the literature on learning BNs rests on the assumption that $\mathcal{D}$ is complete, that is, a data set in which every variable has an observed value for each sample. However, in real-world applications we frequently have to deal with incomplete data; some samples will be completely observed while others will contain missing values for some of the variables. While it is tempting to simply impute the missing values as a preprocessing step, it has long been known that even fairly sophisticated techniques like hot-deck imputation are problematic in a multivariate setting [55]. Just deleting incomplete samples can also bias learning, depending on how missing data are missing.

Rubin [98] and Little and Rubin [67] formalised three possible patterns (or mechanisms) of missingness (illustrated in Figure 3):

• •

Missing completely at random (MCAR): when complete samples are indistinguishable from incomplete ones. In other words, the probability that a value will be missing is independent from both observed and missing values. For instance, in the left panel of Figure 3 both $X_{2}$ and $X_{6}$ are only partially observed; hence $X_{2}=(X_{2}^{O},X_{2}^{M})$ and $X_{6}=(X_{6}^{O},X_{6}^{M})$ where $X_{2}^{0}$ and $X_{6}^{O}$ are observed and $X_{2}^{M}$ and $X_{6}^{M}$ are missing. The patterns of missingness are controlled by $M_{2}$ (for $X_{2}$) and $M_{6}$ (for $X_{6}$) and are completely random; say, $M_{2}$ and $M_{6}$ are binary variables that encode the probability of two instruments (independently) breaking down and thus failing to measure $X_{2}$ and $X_{6}$ for some individuals.

• •

Missing at random (MAR): cases with incomplete data differ from cases with complete data, but the pattern of missingness is predictable from other observed variables. In other words, the probability that a value will be missing is a function of the observed values. An example is the central panel of Figure 3: compared to the left panel, we now know that the two instruments are likely to fail when (say) high values of $X_{1}$ (and $X_{3}$, in the case of $X_{6}$) are observed.

• •

Missing not at random (MNAR): the pattern of missingness is not random or it is not predictable from other observed variables; the probability that an entry will be missing depends on both observed and missing values. Common examples are variables that are missing systematically or for which the patterns of missingness depends on the missing values themselves. In the right panel of Figure 3, say that $X_{2}$ is censored (that is, it is never observed if its value is higher than a fixed threshold); but when $X_{2}$ is missing the probability that $X_{6}$ is missing (encoded by $M_{6}$) is also extremely high. Since we never observe $X_{2}$ when $X_{6}$ is missing, we are unable to correctly model this relationship; as far as we know both $X_{2}$ and $X_{6}$ may be missing due to some common external factor, since they appear to be missing together in the data.

MCAR and MAR are ignorable patterns of missingness; the probability that some value is missing may depend on observed values but not on missing values, and thus can be properly modelled. If we denote with $\mathcal{D}^{O}$ and $\mathcal{D}^{M}$ the observed and unobserved portions of $\mathcal{D}$, and we group all the binary missingness indicators $M_{i}$ in $\mathbf{M}$ with parameters $\Xi$, then we can write

[TABLE]

if the missing data are MAR then $\mathbf{M}$ only depends on $\mathcal{D}^{O}$,

[TABLE]

and if the missing data MCAR then $\mathbf{M}$ does not depend on either $\mathcal{D}^{O}$ or $\mathcal{D}^{M}$,

[TABLE]

In both cases it is possible to model $\mathbf{M}$ from the available data. However, this is not the case for MNAR since $\mathbf{M}$ depends on the unobserved $\mathcal{D}^{M}$.

Modelling incomplete data is analytically intractable and computationally prohibitive compared to the complete data scenario: an exact analysis requires the computation of the joint posterior distribution of $\Theta$ considering all possible completions of $\mathcal{D}$. However, completing $\mathcal{D}$ has a computational complexity that grows exponentially with the number of missing entries, since it involves finding the set of missing data completions with the highest joint probability given the observed data. Considering just the most probable completion can also induce over-confidence in the results of the analysis, since completed observations will have lower variability by construction; but full Bayesian inference averaging over all possible $\Theta$ and completed $\mathcal{D}$ is so computationally challenging as to be unfeasible even in simple settings. Furthermore, this joint maximisation breaks the parameter independence assumption and makes it impossible to define decomposable scores for structure learning without resorting to some approximation.

4.1 Parameter Learning

Many approaches have been developed for parameter learning from incomplete data given a fixed, known structure888Although almost all the literature specific to BNs assumes discrete data., ranging from iterative methods like Data Augmentation [DA; 115] and the Expectation-Maximisation algorithm [EM; 62], to methods based on probability intervals such as Bound and Collapse [BC; 91] and uncertain probabilities [26]. (See Figure 4 for an overview of those we will discuss in detail below.) These methods usually assume missing data are MCAR or MAR to work on $\operatorname{P}(\Theta,\mathcal{D}^{M}\operatorname{|}\mathcal{G},\mathcal{D}^{O})$, and their accuracy decreases dramatically if that is not the case [113]; however, BC has been found to be robust also for MNAR data. Oniśko et al. [84] also noted that simple imputation approaches can perform well in learning the parameters of a BN given a fixed, sparse network structure, which further expands available options.

In the context of parameter learning, the EM algorithm retains its classic structure:

• •

the expectation (E) step consists in computing the expected values of the sufficient statistics (such as the counts $n_{ijk}$ in discrete BNs, partial correlations in GBNs), using exact inference along the lines described above to make use of incomplete as well as complete samples;

• •

the maximisation (M) step takes the sufficient statistics from the E-step and estimates the parameters of the BN, either using maximum likelihood or Bayesian posterior estimators.

The parameter estimates are then used in the next E-step to update the expected values of the sufficient statistics; repeated iterations of these two steps will in the limit return the maximum likelihood or maximum a posteriori estimates for the parameters. Using the notation of (10), the E-step is equivalent to computing $\operatorname{E}(\mathcal{D}^{M},\mathcal{D}^{O}\operatorname{|}\mathcal{G},\Theta)$ and the M-step to maximising $\operatorname{P}(\Theta\operatorname{|}\mathcal{G},\mathcal{D}^{O},\mathcal{D}^{M})$.

DA is quite similar, but instead of converging iteratively to a single set of parameter estimates it uses Gibbs Sampling [43] to generate values from the posterior distributions of both $\mathcal{D}^{M}$ and $\Theta$. The two steps are as follows:

• •

in the imputation (I) step the data are completed with values drawn from the predictive distributions of the missing values;

• •

in the parameter (P) estimation step a parameter value is drawn from the posterior distribution of $\Theta$ conditional on the completed data from the I-step.

More formally, we define the augmented parameter vector $\{\mathcal{D}^{M},\Theta\}$ containing both the missing values and the parameters of the BN. Given an initial set of values, it updates each element of $\{\mathcal{D}^{M},\Theta\}$ by sampling a new value for each missing value from $\operatorname{P}(\mathcal{D}^{M}_{i}\operatorname{|}\mathcal{D}^{M}_{-i},\Theta)$, and by sampling a new value for each parameter from $\operatorname{P}(\Theta_{i}\operatorname{|}\Theta_{-i},\mathcal{D}^{M})$, each in turn. After an initial burn-in phase, this process will converge to its stationary distribution and return parameter realisations from the posterior distribution of $\Theta$ conditional on the observed data. A similar Gibbs sampling approach has been proposed more recently by Riggelsen [94]: it samples from a simpler, approximate predictive distribution and makes use of weights to implement an efficient importance sampling scheme. In addition, the weights make it possible to use samples generated in the burn-in phase as well as those from the stationary distribution of the Gibbs samples because they weight samples according to their estimated predictive accuracy.

Finally, BC and its successor the Robust Bayesian Estimator [RBE; 92] exploit the discrete nature of categorical variables to produce rough interval estimates of the conditional probabilities learned from incomplete variables, which are then reduced to point estimates either via a convex combination of the intervals’ bounds, expert knowledge or both. The first (bound) step, uses the fact that each $\pi_{ik\operatorname{|}j}$ can be bounded below by assuming that none of the missing values for $X_{i}$ are completed with their $k$th value when the $j$th parent configuration is observed; and that $\pi_{ik\operatorname{|}j}$ can be bounded above by assuming that all missing values are completed with their $k$th value. This approach has the merit of not making any assumptions on the distribution of missing data. Furthermore, the width of each interval provides an explicit representation of the reliability of the estimates, which can be taken into account in inference and prediction. The second (collapse) step assumes missing data are MAR or MCAR to be able to compute the expected completions for incomplete samples, which are then used to compute the mean and variance of the $\pi_{ik\operatorname{|}j}$. Interestingly, the intervals from the bound step can be used to augment both EM and Gibbs sampling and obtain more precise inferences, but the predictive accuracy of RBE was shown to be superior to both in Ramoni and Sebastiani [92]. A similar investigation in the context of BN classifiers can be found in Peña et al. [88]. A conceptually similar approach was also proposed by Liao and Ji [65], which first used qualitative expert knowledge on the parameter values to bound them, and then estimated their values using convex optimisation embedded in the EM algorithm.

Note that we can work on latent variables using similar approaches as long as $\mathcal{G}$ is fixed and we just need to learn $\Theta$. A recent example is given in Yamazaki and Motomura [121], who show it is possible to learn the domain of a latent discrete variable as well as the associated parameters as long as it has observed parents and children. Several other examples are discussed in the context of DBNs in Murphy [78].

4.2 Structure Learning

Learning the structure of a BN from incomplete data is, in many respects, an extension of the techniques covered in Section 4.1. (See Figure 4 for an overview of the approaches we will discuss in detail below as well as those for parameter learning.) In its general form it is computationally unfeasible because we need to perform a joint optimisation over the missing values and the parameters to score each candidate network. Starting from (4), we can make this apparent by rewriting $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ as a function of $\mathcal{D}^{O},\mathcal{D}^{M}$:

[TABLE]

From this expression we can see that in order to maximise $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ we should jointly maximise the probability of the observed data $\mathcal{D}^{O}$ and the probability of the missing data $\mathcal{D}^{M}$ given the observed data, for each candidate $G$ and averaging over all possible $\Theta$. This gives us the maximum a posteriori solution to structure learning; a full Bayesian approach would require averaging over all the possible configurations of the missing data as well, leading to

[TABLE]

Compared to (11), (12) contains one extra dimension for each missing value (in addition to one dimension for each parameter in $\Theta$) and thus it is too high-dimensional to compute in practical applications. An additional problem is that, while $\operatorname{P}(\mathcal{D}^{O}\operatorname{|}\mathcal{G},\Theta)$ decomposes as in (5), $\operatorname{P}(\mathcal{D}^{M}\operatorname{|}\mathcal{D}^{O},\mathcal{G},\Theta)$ does not in the general case.

In order to sidestep these computational issues, the literature has pursued two possible approaches: iteratively completing and refining the data, an using standard algorithms and scores for complete data; or using scoring functions that approximate BIC and $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$ but that are decomposable and can be computed efficiently even on incomplete data.

The Structural EM algorithm [SEM999This acronym is another source of confusion, since SEM can also stand for “structural equation models” which are closely related to BNs. See Gupta and Kim [49] for a discussion of their similarities and differences.; 36] is the most famous implementation of the first approach; it has important applications in phylogenetics [38], clinical record [118] and clinical trial analysis [66]. SEM makes structure learning computationally feasible by searching for the best structure inside of EM, instead of embedding EM inside a structure learning algorithm. It consists of two steps like the classic EM:

• •

in the E-step, we complete the data by computing the expected sufficient statistics using the current network structure;

• •

in the M-step, we find the structure that maximises the expected score function for the completed data.

Since the scoring in the M-step uses the completed data , structure learning can be implemented efficiently using standard algorithms. The original proposal by Friedman [36] used BIC and greedy search; Friedman [37] later extended SEM to a fully Bayesian approach based posterior scores, and proved the convergence of the resulting algorithm.

In fact, any combination of structure learning algorithm and score can be used in the M-step; most recently Scanagatta et al. [103] proposed learning BNs with a bounded-treewidth structure using their k-MAX algorithm and a variant of BIC. (This has the two-fold advantage of speeding up the M-step and of yielding BNs for which completing data in the E-step is relatively fast.) Singh [111] proposed a similar approach based on DA, generating sets of completed data sets and averaging the resulting learned networks in each iteration. Myers et al. [80] also chose to iteratively learn both the network structures and the missing data at the same time, but did so using evolutionary algorithms and encoding both as “genes”. Hence, the individuals in the population being evolved comprise both a completed data set and the associated BN. This was combined with Metropolis-Hastings to speed up learning as discussed in Myers et al. [79].

More recently, Adeel and de Campos [1] proposed an exact learning algorithm that explicitly models the patterns of missingness with auxiliary variables, which are included as separate nodes in the BN rather than just being computational devices. Additionally, they showed that its computational complexity is the same as that of other exact learning algorithms for complete data; and they adapted the proposed algorithm into a (faster) heuristic that is then proven to be consistent.

The second group of approaches includes the variational-Bayesian EM from Beal and Ghahramani [10] that maximises a variational approximation of $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$, which in turn is a lower bound to the true marginal likelihood. Balov [5] proved that structure learning with BIC is not consistent even under MCAR, and suggested replacing it with node-average penalised log-likelihoods computed from locally complete observations. An alternative consists in using approximations based on mixtures of truncated exponentials, as was showcased in Fernández et al. [33]: they combined EM and DA to fit regression-like BNs with structures that resemble naive Bayes and tree-augmented naive Bayes classifiers, and approximating explanatory variables. Approachable introductions to this area of research are provided in Chickering and Heckerman [20] and Heckerman [51], which describe the relationship between $\operatorname{P}(\mathcal{D}\operatorname{|}\mathcal{G})$, its Laplace approximation and BIC in mathematical detail.

5 Software and Applications

To conclude this review, we consider the practicalities of using the BN models discussed above to investigate real-world problems: the availability of software implementations and notable examples of applications in various fields.

5.1 Available Software

Among the DBN papers we covered, only Grzegorczyk and Husmeier [47, 48] and [60] provide implementations of the methods they propose as Matlab scripts available upon request from the authors. Additionally, Liu et al. [68] provides links to the custom software they use in the paper. Hugin [4] implements discrete-time, homogeneous DBNs as defined in (6), as does GeNIe [9].

Software implementations for handling incomplete are more readily available, even though typically not from the papers that propose the methods. Both Hugin and GeNIe implement EM for parameter learning, but not Structural EM for parameter learning. The R packages bnlearn [106] and bnstruct [35] implement the structural EM algorithm; another implementation in Matlab is made available as part of the supplementary material of [93]. Adeel and de Campos [1] have later used it to implement their exact learning algorithm along with Gobnilp [8]; and [94] uses a custom C++ implementation of his Gibbs sampling approach. The approaches proposed by Balov [5] are implemented in the catnet R package [6], and those proposed by Fernández et al. [33] are included in Elvira [31].

5.2 Notable Applications

Dynamic and incomplete data are common in many application fields in which BNs are used; the need to work with such data is the main motivation for developing the approaches we covered in Sections 3 and 4. Therefore, many of them have been used in advanced applications especially in the clinical and life sciences.

This is especially true for DBNs: since they subsume a number of statistical models that are themselves very popular (see Section 3.2), it would be fair to say that they are used in most branches of science under different guises. However, their main application fields in the literature remain genetics and systems biology, and in particular the modelling of cell signalling pathways in the form of networks as they evolve over time. This involves probing the state of multiple proteins or other omics simultaneously through time, and studying their interactions and how they evolve in response to external stimuli. In vitro experiments, in which cell lines are probed in a controlled environment, have been modelled as homogeneous DBNs [see, for instance, 53, 89]. On the other hand, more complex and in vivo experiments often involve molecular systems that are not in a steady-state regime and thus require non-homogeneous DBNs. Some examples based on the models presented in Grzegorczyk and Husmeier [47, 48], Kamalabad and Grzegorczyk [56] and other papers from the same authors study embryonic stem cells [74]; sepsis as a cause of acute lung injury [32]; blunt trauma and traumatic spinal cord injury association with inflammation [122] and hypotension [3]. Other omics data including transcriptomics [69] and microbiome data [75, 70] have been studied using DBNs as well.

DBNs have also seen application in ecology, specifically in studying species dynamics. Like the human body, an ecosystem acts as a complex system where species interact with each other (for example, through predation or competition) and with the environment they live in. The resulting dynamics have been modelled in Aderhold et al. [2] with the non-homogeneous DBN from Grzegorczyk and Husmeier [48]. A DBN that corrects for spatial correlation as well through the use of latent variables have been developed in Trifonova et al. [117] to study fish species dynamics in 7 geographically and temporally varied areas within the North Sea. Furthermore, ecosystems have been studied using DBNs from the point of view of environmental sciences as well, considering for example the impact of climate change on groundwater [76] and how to best manage water reservoirs under infrequent rainfalls [96].

Finally, a third class of applications of DBNs is modelling mechanical systems in engineering, reliability and quality control studies. For instance, studying the deterioration of steel structures subject to fatigue [71] including warships [64]; implementing real-time reliability monitoring of subsea pipes [13] and CNC industrial systems [116]; fault detection, identification and recovery in autonomous spacecrafts [22]; and predicting traffic dynamics using probes in San Francisco [54] In each of these applications, DBNs capture the interactions of the components of the mechanical systems they are modelling while at the same time incorporating the effects of the surrounding environment, thus allowing the operators (or the systems itself) to adjust its mode of operation or to stop it completely before causing any damage or safety issue. Software systems have been similarly modelled, for example to study the relationship between security vulnerabilities and attack vectors in network security [39].

As for incomplete data, they are such a widespread issue in science and engineering that it is difficult to identify any particular field in which they are prevalent; hence the BN learning approaches in Section 4 have found very different applications in the literature. We note, however, that the analysis of clinical data has historically been one of the key applications that has driven methodological developments in this area. We mentioned earlier the study of clinical records of patients chronic obstructive pulmonary disease [118] and the recovery process of patients with whiplash-associated disorders [66]; two other examples are the diagnosis of brain diseases such as dementia and Alzheimer [110] and of sleep apnea [100]. Completely different applications have also been explored: one example is the classification of radar images [72].

6 Summary

In this paper we have reviewed the fundamental definitions and properties of BNs, and how BNs can be stretched to encode more complex probabilistic models than what might be apparent from reference material and most of the literature. Both have an overwhelming focus on the straightforward case in which the data being modelled are both static (that is, with no time dimension) and complete (that is, with no missing values). However, dynamic and incomplete data are central to many cutting-edge applications in research fields ranging from genetics to robotics; BNs can play an important role in many of these settings, as has been evidenced by the examples referenced in Sections 3 and 4. Given their expressive power, BNs also subsume several classic probabilistic models and can augment them with automatic reasoning capabilities through various kinds of queries that can be performed algorithmically. Much research has been and is being developed to adapt BNs to these applications and make them a competitive choice for modelling complex data.