DYNOTEARS: Structure Learning from Time-Series Data

Introduction

Graphical models are a popular approach to understanding large datasets, and provide convenient, interpretable output that is needed in today’s high stakes applications of machine learning and artificial intelligence. In particular, with the growing need for interpretable models and causal insights about an underlying process, directed acyclic graphs (DAGs) have shown promise in many applications. The edges in a DAG provide users with important clues about the relationship between variables in a system. When these edges are not known based on prior knowledge, it is necessary to resort to structure learning, namely, the problem of learning the edges in a graphical model from data. Broadly speaking, structure learning can be divided into static (i.e. equilibrium) and dynamic models, the latter of which explicitly model temporal dependencies. Static models make sense for independent and identically distributed data. Many applications, however, exhibit strong temporal fluctuations that we are interested in modeling explicitly. The problem of learning graphical structures from temporal data collected from dynamic systems has received significant attention from the machine learning (Koller and Friedman, 2009), econometrics (Lütkepohl, 2005), and neuroscience (Rajapakse and Zhou, 2007) communities.

In this paper, we revisit the problem of learning dynamic Bayesian networks (DBNs) (Dean and Kanazawa, 1989; Murphy, 2002) from data. DBNs have been used successfully in a variety of domains such as clinical disease prognosis (Van Gerven et al., 2008; Zandonà et al., 2019), gene regulatory network (Linzner et al., 2019), facial and speech recognition (Meng et al., 2019; Nefian et al., 2002), neuroscience (Rajapakse and Zhou, 2007), among others. DBNs are the standard approach to modeling discrete-time temporal dynamics in directed graphical models. In econometrics, they are also known as structural vector autoregressive (SVAR) models (Demiralp and Hoover, 2003; Swanson and Granger, 1997).

We propose a simple, score-based approach for learning these models that scales gracefully to high-dimensional datasets. To accomplish this, we cast the problem as an optimization problem (i.e. score-based learning), and use standard second-order optimization schemes to solve the resulting program. Our approach is based on the recent algebraic characterization of acyclicity in directed graphs from Zheng et al., 2018, which makes the formulation simple and amenable to different modeling choices.

The main contributions of this paper are the following:

We develop a score-based approach to learning DBNs and use standard optimization routines to optimize the resulting program. The resulting method, which we call DYNOTEARS, can be used to learn time series of arbitrary order, without any implicit assumptions on the underlying graph topologies such as bounded in-degree or treewidth.

We validate our approach with extensive simulation experiments, exhibiting the accuracy of our approach in learning both intra-slice and inter-slice relationships in dynamic models.

We apply our method to two real datasets: A financial dataset consisting of daily stock returns ($d=97$) and the DREAM4 dataset ($d=100$) (Marbach et al., 2009). These examples illustrate the importance of modeling both temporal trends as well as steady state relationships and achieve competitive accuracy among other DBN methods.

The resulting method simultaneously achieves three important goals: 1) Accuracy on high-dimensional data with $d>n$, which allows for application to real world data 2) Robustness to complex graph topologies, and 3) A simple, plug-n-play algorithm for learning based on black-box optimization.

Related work

There are many methods for learning DBNs in the literature. Some approaches ignore contemporaneous dependencies and recover only time-lagged relationships (Haufe et al., 2010; Song et al., 2009). Others learn both types of relationships independently (Haufe et al., 2010; Song et al., 2009). Many methods follow a two-step approach of first learning inter-slice weights and then estimating intra-slice weights from the residuals from the first step (Chen and Chihying, 2007; Hyvärinen et al., 2010; Moneta et al., 2011). There are also hybrid algorithms that combine conditional-independence tests and local search to improve the score (Malinsky and Spirtes, 2018; Malinsky and Spirtes, 2019). While all of these methods can achieve good structure recovery on small graphs, they suffer from the curse of dimensionality. More discussion and comparison can be found in section 2.3. There is also an extensive literature on learning SVAR models in the econometrics and statistics literature (Demiralp and Hoover, 2003; Lanne et al., 2017; Reale and Wilson, 2001; Reale and Wilson, 2002; Swanson and Granger, 1997; Tank et al., 2019).

Since algorithms for learning DBNs typically rely internally on calling methods for learning static BNs, it is worth briefly reviewing this here. These methods can be classified into constraint-based methods and score-based methods, as well as hybrid methods that combine these two approaches. Constraint-based methods use conditional independence tests to recover the Markov equivalence class of DAGs under the assumption of faithfulness (Colombo et al., 2012; Spirtes and Glymour, 1991, e.g.). While this yields fast algorithms, these methods are sensitive to the underlying graph structure and suffer from error propagation (Spirtes, 2010). Score-based methods, on the other hand, use a score function to find the best DAG that fits the given data (Heckerman et al., 1995; Chickering, 2002; Bouckaert, 1993, e.g.). Examples of scores include the BIC, BDe, and BDeu (Spirtes et al., 2000). Score-based methods are computationally expensive due to the acyclicity constraint and the vast number of DAGs to search over (Robinson, 1977). The recent work of Zheng et al., 2018 expresses the acyclicity of a DAG by a smooth equality constraint, which makes it possible to formulate structure learning as a smooth minimization problem subject to this equality constraint.

Dynamic Structure Learning

We model the data using the following standard SVAR model (Demiralp and Hoover, 2003; Kilian, 2011; Swanson and Granger, 1997):

2 Optimization problem

where $F$ is the following smooth augmented objective:

3 Alternative formulations

Experiments

To validate the effectiveness of the proposed method, we performed a series of simulation experiments for which the ground truth is known in advance. We focus here on the main results; more detailed results can be found in Appendix B.

2 Results

In fig. 3, we compare the performance of DYNOTEARS to that of three other algorithms (see section B.1). We consider two distributions for the noise, Gaussian and Exponential, and two choices for the number of samples, $n\in\{50,500\}$. Within each of the four panels, each column corresponds to one choice of intra-slice graph model and mean degree; for instance, ER2 indicates that we used an Erdős–Rényi graph model with a mean degree of 2. Similarly, each row corresponds to one choice of inter-slice graph model and mean degree. For each individual plot, we generate $n$ data samples with autoregression order $p=1$ for five choices of the number of variables, $d\in\{5,10,20,50,100\}$. The vertical axis indicates the performance of each algorithm in terms of the F1 score, which we calculate separately for intra-slice and inter-slice matrices. We discuss the selection of hyperparameter values for the four algorithms in section B.3. The relative ranking of the four algorithms is not especially sensitive to these hyperparameters. In particular, the regularization parameters are largely irrelevant when there is sufficient data ($n\gg dp$).

DYNOTEARS is the best-performing algorithm in fig. 3, with F1 scores close to 1 for $n=500$. DYNOTEARS is also the best-performing algorithm when the number of variables exceeds the number of samples (see panels (c) and (d) for $d=100$). This high-dimensional case is especially difficult, yet common in applications, and we discuss one such example in section 4.2. The second-best algorithm is tsGFCI. However, its performance tends to degrade as we add more edges to the ground-truth graphs; see fig. 7 in the Appendix. The variance in performance is also larger for tsGFCI than for the other algorithms. As expected, learning intra-slice and inter-slice structure separately with NOTEARS + Lasso underperforms DYNOTEARS. In particular, we find that the Lasso step falsely identifies some intra-slice edges as inter-slice. LiNGAM is an algorithm designed for non-Gaussian data, so its poor performance in panels (a) and (c) of fig. 3 is not surprising. However, even on data with exponential noise (panels (b) and (d) of fig. 3), its performance degrades significantly as $d$ increases. section B.7 contains additional figures that compare the performance of the four algorithms using metrics other than the F1 score.

Applications

Our approach allows us to detect whether contemporaneous or time-lagged relationships are more meaningful in a given dataset. We discuss two examples for which each type of interaction is dominant.

2 DREAM4 gene expression data

As in the original DREAM4 challenge, we use mean AUPR and AUROC across the 5 datasets to compare DYNOTEARS to the 23 algorithms presented in Lu et al., 2019. The results are as follows:

DYNOTEARS achieves an average AUROC of 0.664 and an average AUPR of 0.173. Among the six DBN methods tested, this ranks 1st and 2nd in AUPR and AUROC, respectively, as shown in Table 1.

Furthermore, DYNOTEARS is within one standard deviation of the best performing method (G1DBN) based on AUPR, and no other method is within one standard deviation of DYNOTEARS based on AUROC.

Overall, this ranks 4th in AUPR and 8th in AUROC (see Table 3 and Table 2, respectively). While the top performing methods were based on nonparametric models such as GPs, DYNOTEARS still outperforms several other nonparametric methods despite its use of the linear model (1).

Discussion

In this paper, we proposed DYNOTEARS, an algorithm for learning dynamic Bayesian networks, inspired by recent work on structure learning for static Bayesian networks using differentiable acyclicity constraints (Zheng et al., 2018). Our algorithm learns both intra-slice and inter-slice dependencies between variables simultaneously, in contrast with some existing methods that perform these estimations in succession.

An important feature of DYNOTEARS is its simplicity, both in terms of formulating an objective function and in terms of optimizing it. It performs well on simulated data across a wide range of parameter choices in the data-generation process.

We also applied DYNOTEARS to two empirical datasets from different application domains, finance and molecular biology. The results reveal insightful patterns in the data. Both of these applications have $d\approx 100$, confirming that DYNOTEARS can be applied to larger datasets than those considered in most existing work on DBNs.

To conclude, we briefly discuss some limitations and possible extensions of DYNOTEARS.

We have assumed that the structure of the DBN is fixed through time and is identical for all time series in the data (i.e., it is the same for all $m\in\{1,\ldots,M\}$). It would be useful to relax these assumptions in various ways, for example by allowing the structure to change smoothly over time (Song et al., 2009) or at discrete change points that we infer from the data (Grzegorczyk and Husmeier, 2011). Another topic for future work is to investigate the behaviour of the algorithm on nonstationary or cointegrated time series (Malinsky and Spirtes, 2019), or in situations with confounders (Huang et al., 2015; Malinsky and Spirtes, 2018). A possible approach is to apply a post-processing step to the output of DYNOTEARS so as to remove spurious relationships between variables (e.g., by using statistical tests).

Undersampling

As pointed out by a reviewer, as with most DBN models, we implicitly assume that the sampling rate of the process is at least as high as the fluctuations in the underlying causal process. See for example Gong et al., 2015; Hyttinen et al., 2016; Plis et al., 2015; Cook et al., 2017. As a check on the sensitivity of DYNOTEARS to this (strong) assumption, we tested a modification of our approach for undersampled data adapted from Cook et al., 2017. In essence, by running DYNOTEARS on the data and then post-processing any intra-slice edges by making them bi-directed, we can obtain a graph which is comparable to the methods suggested in Cook et al., 2017. Although our method was not designed to handle undersampling, it achieves lower false-negative rates compared to other methods. Of course, more careful modifications to handle undersampling is an interesting direction for future work.

Nonlinear dependence

Finally, we emphasize that linear assumption in (1-2) is made purely for simplicity, in order to keep the focus on the most salient dynamic and temporal aspects of this problem. For example, using the general approach outlined in Zheng et al., 2019, it is possible to model complicated nonlinear dependencies via neural networks or orthogonal basis expansions. Furthermore, it is straightforward to replace the least squares loss with the logistic loss (or more generally, any exponential family log-likelihood) to model binary data. It is also possible to go a step further and consider combinations of continuous and discrete data (Andrews et al., 2019), which is important for many real-world applications.

References

Appendices

Appendix A Comparison of one-stage and two-stage algorithms

It is possible to minimize the DYNOTEARS objective using either a one-stage algorithm (see section 2.2) or a two-stage algorithm (see section 2.3). The two formulations give nearly identical results when the number of samples exceeds the number of variables (i.e., when $n\gg dp$). However, the two-stage algorithm runs somewhat faster, so it should be the preferred option in cases where there is sufficient data.

Appendix B Numerical experiments

B.2 Interpreting a PAG as a DAG

The tsGFCI algorithm returns a partial ancestral graph (PAG), which one cannot immediately compare to a ground-truth DAG. Thus, we developed a set of rules to convert the PAG output to a DAG, making sure to do so in a manner that favors tsGFCI. The rules are as follows:

If an edge is directed (i.e., A $\rightarrow$ B in the PAG), then we treat it as a directed edge in the DAG.

If an edge in the PAG is either directed or it indicates the presence of a latent factor (i.e., A $\circ\rightarrow$ B), then we check whether the directed edge exists in the ground truth graph and assume that tsGFCI made the correct choice.

If two nodes are related through a latent variable (i.e., A $\longleftrightarrow$ B in the PAG), then we disregard the edge.

If the edge is ambiguous (i.e., A $\circ-\circ$ B), then we assume that tsGFCI made the correct choice; in other words, we check whether A $\rightarrow$ B, B $\leftarrow$ A, or A is not connected to B in the ground-truth DAG and we assume that tsGFCI made the same choice.

Using these rules, we pick the outcomes most favorable for tsGFCI in ambiguous cases. This implies, in particular, that our results slightly overstate the performance of tsGFCI on simulated data.

B.3 Hyperparameter selection

Although we did not attempt to optimize hyperparameters for our experiments on simulated data, our work on the S&P100 and the DREAM4 datasets (see Section 4) indicates ways in which one can estimate these parameters through cross-validation.

B.4 Data generation process

We provide more details about the data generation process that we use in our numerical experiments from Section 3.

To go from an unweighted to a weighted DAG, we follow Zheng et al., 2018 and we sample weights uniformly at random from $[-2.0,\,-0.5]\cup[0.5,\,2.0]$.

Inter-slice model

B.5 Autoregressive order

B.6 Running times

Although tsGFCI and LiNGAM run significantly faster than DYNOTEARS, we believe that the gain in accuracy from using the latter makes it worthwhile even in cases with hundreds of variables. As a reminder, running DYNOTEARS on the S&P100 dataset (which has $d\approx 100$) takes a few minutes on a typical laptop.

B.7 Additional results

The relative performance of different algorithms varies as we change the density of the ground-truth graphs. In fig. 7, we show the F1 scores for simulated data with Gaussian noise, $n=500$, four choices of intra-slice graphs (columns), and four choices of inter-slice graphs (rows). The performance of tsGFCI is especially sensitive to changes in graph densities, with a notable drop in F1 scores when the intra-slice graph is ER4.

Additional performance metrics

In figs. 8, 9, 10 and 11, we plot four additional performance metrics to complement the F1 scores from fig. 3 in the main part of the paper. The metrics are standard and are defined as follows:

True positive rate (TPR): number of correctly-identified edges divided by the number of edges in the ground-truth graph.

False discovery rate (FDR): number of incorrectly-identified edges divided by the number of edges in the estimated graph.

Structural Hamming distance (SHD): number of changes (i.e., edge removals, edge additions, and edge reversals) required to go from one (unweighted) graph to another.

Frobenius norm (FRO) of the difference between two weighted matrices (i.e., between a ground-truth adjacency matrix and an estimated matrix).

Note that the Frobenius norm does not apply to the tsGFCI algorithm, which only returns unweighted edges.

For $n=500$, DYNOTEARS generally outperforms the other algorithms for both Gaussian (fig. 8) and exponential noise (fig. 9); there are some exceptions to this when the number of variables is small, $d\in\{5,10,20\}$. For $n=50$ (see figs. 10 and 11), NOTEARS + Lasso and LiNGAM output estimated graphs that are significantly denser than the ground truth. As a result, while these two algorithms have large TPRs, their overall performance is not competitive due to a large number of false positives.

Appendix C S&P100 application

Appendix D DREAM4 application

D.2 Comparison to other methods

In tables 2 and 3, we compare the performance of DYNOTEARS to that of other methods. We obtain performance metrics for other algorithms from Lu et al., 2019.