Random data generation from Gaussian DAG models

Generating DAGs and parameters: functions rDAG() and rDAGWishart()

Function rDAG() can be used to randomly generate a DAG structure 𝒟 = (V, E), where V = {1, …, q} and E ⊆ V × V is the set of edges. rDAG() has two arguments: the number of nodes (variables) q and the prior probability of edge inclusion w ∈ [0, 1]; the latter can be tuned to control the degree of sparsity in the resulting DAG. By fixing a probability of edge inclusion w = 0.2, a DAG structure with q = 10 nodes can be generated as follows:

set.seed(1)
q <- 10
w <- 0.2
DAG <- rDAG(q,w)

DAG
#>    1 2 3 4 5 6 7 8 9 10
#> 1  0 0 0 0 0 0 0 0 0  0
#> 2  0 0 0 0 0 0 0 0 0  0
#> 3  0 0 0 0 0 0 0 0 0  0
#> 4  0 0 1 0 0 0 0 0 0  0
#> 5  1 0 0 0 0 0 0 0 0  0
#> 6  0 0 0 0 0 0 0 0 0  0
#> 7  1 0 1 0 0 0 0 0 0  0
#> 8  1 0 0 0 0 0 0 0 0  0
#> 9  0 0 0 1 0 0 1 0 0  0
#> 10 0 0 0 0 1 0 0 0 0  0

Output of rDAG() is the 0-1 (q, q) adjacency matrix of the generated DAG, with element 1 at position (u, v) indicating the presence of an edge u → v. Notice that the generated DAG is topologically ordered, meaning that edges are allowed only from high to low nodes (nodes are labeled according to rows/columns indexes); accordingly the DAG adjacency matrix is lower-triangular.

Generating Gaussian DAG parameters

Consider a Gaussian DAG model of the form where (L, D) are model parameters providing the decomposition of the precision (inverse-covariance) matrix Ω = LD⁻¹L^⊤; specifically, L is a (q, q) matrix of coefficients such that for each (u, v)-element L_uv with u ≠ v, we have L_uv ≠ 0 if and only if (u, v) ∈ E, while L_uu = 1 for each u = 1, …, q; also, D is a (q, q) diagonal matrix with (u, u)-element D_uu. The latter decomposition follows from the equivalent Structural Equation Model (SEM) representation of a Gaussian DAG model:

where x = (X₁, …, X_q)^⊤; see also Castelletti & Mascaro (2021).

Function rDAGWishart implements random sampling from (L, D) | 𝒟 ∼ DAG-Wishart(a_c^𝒟, U), where U is the rate parameter (a (q, q) s.p.d. matrix) and a_c^𝒟 (a (q, 1) vector) is the shape parameter of the DAG-Wishart distribution. This class of distributions was introduced by Ben David et al. (2015) as a conjugate prior for Gaussian DAG model-parameters. In its compatible version (Peluso & Consonni, 2020), elements of the vector parameter a_c^𝒟 are uniquely determined from a single common shape parameter a > q − 1.

Inputs of rDAGWishart are: the number of samples n, the underlying DAG 𝒟, the common shape parameter a and the rate parameter U. Given the DAG 𝒟 generated before, the following example implements a single (n = 1) draw from a compatible DAG-Wishart distribution with parameters a = q, U = I_q:

a <- q
U <- diag(1,q)
outDL <- rDAGWishart(n=1, DAG, a, U)
class(outDL)
#> [1] "list"

L <- outDL$L; D <- outDL$D
class(L); class(D)
#> [1] "matrix" "array"
#> [1] "matrix" "array"

The output of rDAGWishart() consists of two elements: a (q, q, n)-dimensional array collecting the n sampled matrices L⁽¹⁾, …, L⁽ⁿ⁾ and a (q, q, n)-dimensional array collecting the n sampled matrices D⁽¹⁾, …, D⁽ⁿ⁾. We refer the reader to Castelletti & Mascaro (2021) and Castelletti & Mascaro (2022+) for more details.

Generating data from a Gaussian DAG model

Data generation from a Gaussian DAG model is then straightforward. Recall that Ω = LD⁻¹L^⊤, where Ω is the inverse-covariance (precision) matrix of a multivariate Gaussian model satisfying the constraints imposed by a DAG. Accordingly, we can recover the precision and covariance matrices as:

# Precision matrix
Omega <- L %*% solve(D) %*% t(L)
# Covariance matrix
Sigma <- solve(Omega)

Next, i.i.d. draws from a Gaussian DAG model can be obtained through the function rmvnorm() provided within the R package mvtnorm:

n <- 1000
X <- mvtnorm::rmvnorm(n = n, sigma = Sigma)

References

• Ben-David E, Li T, Massam H, Rajaratnam B (2015). “High dimensional Bayesian inference for Gaussian directed acyclic graph models.” arXiv pre-print.

• Cao X, Khare K, Ghosh M (2019). “Posterior graph selection and estimation consistency for high-dimensional Bayesian DAG models.” The Annals of Statistics, 47(1), 319–348.

• Castelletti F, Mascaro A (2021). “Structural learning and estimation of joint causal effects among network-dependent variables.” Statistical Methods & Applications, 30, 1289–1314.

• Castelletti F, Mascaro A (2022). “BCDAG: An R package for Bayesian structural and Causal learning of Gaussian DAGs.” arXiv pre-print.

• Peluso S, Consonni G (2020). “Compatible priors for model selection of high-dimensional Gaussian DAGs.” Electronic Journal of Statistics, 14(2), 4110–4132.