BDAGAR: Bivariate DAG Autoregression

• BDAGAR is a framework that fuses directed acyclic graphs with autoregressive models to capture directional and autocorrelated dependencies.
• It is applied in spatial epidemiology and time series analysis, enabling clear separation of endogenous clustering from cross-variable associations.
• The model leverages Bayesian posterior inference with efficient MCMC techniques and sparse precision matrix estimation for robust parameter recovery.

Bivariate Directed Acyclic Graphical Autoregression (BDAGAR) refers to a statistical modeling framework that encodes both directional dependence between two variables and domain-specific autocorrelation—spatial or temporal—through a directed acyclic graphical (DAG) structure combined with autoregressive processes. BDAGAR models have been foundational in both spatial epidemiology, for jointly mapping correlated diseases, and in time series causal inference, for learning contemporaneous influence between variables with stationary error processes. In both domains, BDAGAR leverages the interpretability and sparsity of DAGs while retaining the flexibility of autoregressive latent processes and facilitating Bayesian posterior inference with computationally tractable precision structures (Gao et al., 2019, Roy et al., 30 Mar 2025).

1. Mathematical Foundation and DAG Construction

In BDAGAR, variables or spatial units are ordered such that the dependence structure among random effects or contemporaneous variables is encoded by a DAG. In spatial settings, consider $k$ regions arranged as $\mathcal{V} = \{1, 2, ..., k\}$ with an adjacency graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$, where $(i, j) \in \mathcal{E}$ if $i$ and $j$ are neighboring regions. Assigning a lexicographic ordering, each node $i>1$ has $N(i) = \{j < i : (j,i) \in \mathcal{E}\}$ as its set of “past” neighbors, generating a DAG in which information flows from lower- to higher-indexed units (Gao et al., 2019).

In directed time series formulations, $Y_t = (y_{1,t}, y_{2,t})^\top$ at time $t$ is described by the structural equation model

$\mathcal{V} = \{1, 2, ..., k\}$0

where $\mathcal{V} = \{1, 2, ..., k\}$1 is a $\mathcal{V} = \{1, 2, ..., k\}$2 matrix with either $\mathcal{V} = \{1, 2, ..., k\}$3 or $\mathcal{V} = \{1, 2, ..., k\}$4 nonzero—encoding the DAG—and $\mathcal{V} = \{1, 2, ..., k\}$5 is a diagonal innovation variance matrix. The DAG constraint permits interpretation in terms of Granger-noncausality at the contemporaneous level (Roy et al., 30 Mar 2025).

2. Full Bayesian Hierarchical Model Specification

The spatial BDAGAR is formulated hierarchically:

• Data likelihood: $\mathcal{V} = \{1, 2, ..., k\}$6 with $\mathcal{V} = \{1, 2, ..., k\}$7 independently for observed outcome $\mathcal{V} = \{1, 2, ..., k\}$8 in region $\mathcal{V} = \{1, 2, ..., k\}$9.
• Latent spatial processes: $\mathcal{G} = (\mathcal{V}, \mathcal{E})$0 for the first outcome, and $\mathcal{G} = (\mathcal{V}, \mathcal{E})$1 for the second, where $\mathcal{G} = (\mathcal{V}, \mathcal{E})$2 is the univariate DAGAR precision and $\mathcal{G} = (\mathcal{V}, \mathcal{E})$3 parameterizes endemic and adjacency-mediated cross-association via $\mathcal{G} = (\mathcal{V}, \mathcal{E})$4.
• Priors: $\mathcal{G} = (\mathcal{V}, \mathcal{E})$5, $\mathcal{G} = (\mathcal{V}, \mathcal{E})$6, $\mathcal{G} = (\mathcal{V}, \mathcal{E})$7, $\mathcal{G} = (\mathcal{V}, \mathcal{E})$8, and $\mathcal{G} = (\mathcal{V}, \mathcal{E})$9 (Gao et al., 2019).

In the bivariate DAGAR for time series, the Bayesian hierarchy comprises:

• Horseshoe prior on $(i, j) \in \mathcal{E}$0, promoting sparsity in $(i, j) \in \mathcal{E}$1,
• Stick-breaking prior on $(i, j) \in \mathcal{E}$2 to accommodate heterogeneous variances,
• Nonparametric B-spline prior (with simplex constraints) on each spectral density $(i, j) \in \mathcal{E}$3 of the stationary error processes (Roy et al., 30 Mar 2025).

3. Precision Matrix, Covariance Structure, and Parameter Interpretation

The BDAGAR construction yields sparse precision matrices:

• In the spatial setting, the joint $(i, j) \in \mathcal{E}$4 precision matrix $(i, j) \in \mathcal{E}$5 incorporates $(i, j) \in \mathcal{E}$6 and $(i, j) \in \mathcal{E}$7 (within-disease autocorrelation parameters), and $(i, j) \in \mathcal{E}$8, $(i, j) \in \mathcal{E}$9 (cross-disease effects—endemic and neighborhood-mediated, respectively). Off-diagonal blocks encode the endemic and spatial cross-linkage.
• In the time-series context, the contemporaneous precision is $i$0.

Key parameters:

• $i$1: Within-disease or within-variable spatial/temporal clustering,
• $i$2: Endemic or same-unit association,
• $i$3: Cross-variable association mediated by adjacency,
• $i$4: Directed edge encoding the direction and magnitude of influence between contemporaneous series (Gao et al., 2019, Roy et al., 30 Mar 2025).

4. Posterior Inference and Computation

Inference for BDAGAR models is carried out via MCMC:

• Spatial: Block Gibbs sampling for Gaussian latent effects and regression coefficients, Metropolis-Hastings for autocorrelation and cross-disease parameters. The sparse/structured precision ensures computational tractability even at large $i$5. Model comparison is performed using WAIC.
• Temporal: Posterior sampling proceeds without enforcing the DAG constraint during MCMC; each sample of $i$6 is subsequently projected to the closest DAG-matrix via a quasi-Newton optimization with acyclicity penalty. The resulting "projection posterior" yields consistent inferences with negligible computational overhead for $i$7. Posterior marginal edge probabilities (direction selection) are estimated by thresholding across the sampled and projected $i$8 matrices (Roy et al., 30 Mar 2025).

5. Comparative Advantages over Existing Models

The principal distinction between BDAGAR and standard bivariate (multi)variates spatial autoregressive models such as MCAR/GMCAR is the explicit separation of within-outcome spatial clustering from cross-disease association. Traditional MCAR models require an unwieldy and often non-interpretable cross-covariance matrix, whereas BDAGAR parameterizes and identifies each structural and endemic association individually (Gao et al., 2019). The directed acyclic graphical assumption leads to sparse, factorizable, and efficiently invertible precision matrices.

In both spatial epidemiology and time series contexts, BDAGAR attains superior model fit—as evidenced by lower WAIC and more stable effective parameter counts—compared to MCAR/GMCAR and unrestricted vector-autoregressive alternatives. For instance, in spatial mapping of California cancer rates with eight demographic covariates, the BDAGAR structure with [lung] $i$9 [esophagus|lung] ordering yielded the lowest WAIC (417.05), interpretable cross- and within-disease correlations, and provided clear epidemiological interpretations through separation of clustering and endemic linkage (Gao et al., 2019).

6. Applied Examples and Scientific Implications

BDAGAR has been applied to jointly model county-level incidence rates of esophagus and lung cancer across California (2012–2016). The posterior median and credible intervals were:

• $j$0 (weak esophagus clustering),
• $j$1 (moderate lung clustering adjusted for esophagus),
• $j$2 (strong endemic association),
• $j$3 (adjacency-mediated cross-linkage). Posterior spatial correlations between outcomes ranged from 0.97–1.00 across counties. The model identified demographic covariates, such as age and unemployment, differentially influencing each cancer (Gao et al., 2019).

In the time series context, the bivariate DAGAR (with Whittle pseudo-likelihood) supports recovery of causality direction, stationarity, and identifiability under mild conditions. The projection-posterior methodology performed consistently in simulation and real data settings (Roy et al., 30 Mar 2025).

7. Model Selection, Identifiability, and Implementation

Identifiability for BDAGAR is established under either known innovation variances or linearly independent spectral densities. The mapping from $j$4 to the distribution of $j$5 is one-to-one under the DAG constraint. Model order (directionality) can be selected empirically by comparing posterior mass on candidate DAGs or through WAIC in the spatial context.

BDAGAR implementations for spatial data can be realized in R or Python using sparse linear algebra, leveraging conjugacy for most parameters. For time series, unconstrained samplers utilize efficient Gibbs or Hamiltonian steps, augmented by post hoc projections utilizing modern quasi-Newton routines, with hyperparameters tunable through cross-validation or default weakly-informative priors (Roy et al., 30 Mar 2025). Embarrassingly parallel computation is possible across MCMC samples due to the decoupling of projection and likelihood evaluation.

---

References:

(Gao et al., 2019): "Spatial Modeling for Correlated Cancers Using Bivariate Directed Graphs" (Roy et al., 30 Mar 2025): "Bayesian Inference for High-dimensional Time Series with a Directed Acyclic Graphical Structure"