learn_DAG {BCDAG}

MCMC scheme for Gaussian DAG posterior inference

Description

This function implements a Markov Chain Monte Carlo (MCMC) algorithm for structure learning of Gaussian DAGs and posterior inference of DAG model parameters

Usage

learn_DAG(
  S,
  burn,
  data,
  a,
  U,
  w,
  fast = FALSE,
  save.memory = FALSE,
  collapse = FALSE,
  verbose = TRUE
)

Arguments

S

integer final number of MCMC draws from the posterior of DAGs and parameters

burn

integer initial number of burn-in iterations, needed by the MCMC chain to reach its stationary distribution and not included in the final output

data

(n,q) data matrix

a

common shape hyperparameter of the compatible DAG-Wishart prior, a > q - 1

U

position hyperparameter of the compatible DAG-Wishart prior, a (q, q) s.p.d. matrix

w

edge inclusion probability hyperparameter of the DAG prior in [0,1]

fast

boolean, if TRUE an approximate proposal for the MCMC moves is implemented

save.memory

boolean, if TRUE MCMC draws are stored as strings, instead of arrays

collapse

boolean, if TRUE only structure learning of DAGs is performed

verbose

If TRUE, progress bars are displayed

Details

Consider a collection of random variables X_1, \dots, X_q whose distribution is zero-mean multivariate Gaussian with covariance matrix Markov w.r.t. a Directed Acyclic Graph (DAG). Assuming the underlying DAG is unknown (model uncertainty), a Bayesian method for posterior inference on the joint space of DAG structures and parameters can be implemented. The proposed method assigns a prior on each DAG structure through independent Bernoulli distributions, Ber(w), on the 0-1 elements of the DAG adjacency matrix. Conditionally on a given DAG, a prior on DAG parameters (D,L) (representing a Cholesky-type reparameterization of the covariance matrix) is assigned through a compatible DAG-Wishart prior; see also function rDAGWishart for more details.

Posterior inference on the joint space of DAGs and DAG parameters is carried out through a Partial Analytic Structure (PAS) algorithm. Two steps are iteratively performed for s = 1, 2, ... : (1) update of the DAG through a Metropolis Hastings (MH) scheme; (2) sampling from the posterior distribution of the (updated DAG) parameters. In step (1) the update of the (current) DAG is performed by drawing a new (direct successor) DAG from a suitable proposal distribution. The proposed DAG is obtained by applying a local move (insertion, deletion or edge reversal) to the current DAG and is accepted with probability given by the MH acceptance rate. The latter requires to evaluate the proposal distribution at both the current and proposed DAGs, which in turn involves the enumeration of all DAGs that can be obtained from local moves from respectively the current and proposed DAG. Because the ratio of the two proposals is approximately equal to one, and the approximation becomes as precise as q grows, a faster strategy implementing such an approximation is provided with fast = TRUE. The latter choice is especially recommended for moderate-to-large number of nodes q.

Output of the algorithm is a collection of S DAG structures (represented as (q,q) adjacency matrices) and DAG parameters (D,L) approximately drawn from the joint posterior. The various outputs are organized in (q,q,S) arrays; see also the example below. If the target is DAG learning only, a collapsed sampler implementing the only step (1) of the MCMC scheme can be obtained by setting collapse = TRUE. In this case, the algorithm outputs a collection of S DAG structures only. See also functions get_edgeprobs, get_MAPdag, get_MPMdag for posterior summaries of the MCMC output.

Print, summary and plot methods are available for this function. print provides information about the MCMC output and the values of the input prior hyperparameters. summary returns, besides the previous information, a (q,q) matrix collecting the marginal posterior probabilities of edge inclusion. plot returns the estimated Median Probability DAG Model (MPM), a (q,q) heat map with estimated marginal posterior probabilities of edge inclusion, and a barplot summarizing the distribution of the size of DAGs visited by the MCMC.

Value

An S3 object of class bcdag containing S draws from the posterior of DAGs and (if collapse = FALSE) of DAG parameters D and L. If save.memory = FALSE, these are stored in three arrays of dimension (q,q,S). Otherwise, they are stored as strings.

Author(s)

Federico Castelletti and Alessandro Mascaro

References

F. Castelletti and A. Mascaro (2021). Structural learning and estimation of joint causal effects among network-dependent variables. Statistical Methods and Applications, Advance publication.

F. Castelletti and A. Mascaro (2022). BCDAG: An R package for Bayesian structural and Causal learning of Gaussian DAGs. arXiv pre-print, url: https://arxiv.org/abs/2201.12003

F. Castelletti (2020). Bayesian model selection of Gaussian Directed Acyclic Graph structures. International Statistical Review 88 752-775.

Examples

# Randomly generate a DAG and the DAG-parameters
q = 8
w = 0.2
set.seed(123)
DAG = rDAG(q = q, w = w)
outDL = rDAGWishart(n = 1, DAG = DAG, a = q, U = diag(1, q))
L = outDL$L; D = outDL$D
Sigma = solve(t(L))%*%D%*%solve(L)
# Generate observations from a Gaussian DAG-model
n = 200
X = mvtnorm::rmvnorm(n = n, sigma = Sigma)

## Set S = 5000 and burn = 1000 for better results

# [1] Run the MCMC for posterior inference of DAGs and parameters (collapse = FALSE)
out_mcmc = learn_DAG(S = 50, burn = 10, a = q, U = diag(1,q)/n, data = X, w = 0.1,
                     fast = FALSE, save.memory = FALSE, collapse = FALSE)
# [2] Run the MCMC for posterior inference of DAGs only (collapse = TRUE)
out_mcmc_collapse = learn_DAG(S = 50, burn = 10, a = q, U = diag(1,q)/n, data = X, w = 0.1,
                              fast = FALSE, save.memory = FALSE, collapse = TRUE)
# [3] Run the MCMC for posterior inference of DAGs only with approximate proposal
# distribution (fast = TRUE)
# out_mcmc_collapse_fast = learn_DAG(S = 50, burn = 10, a = q, U = diag(1,q)/n, data = X, w = 0.1,
#                                    fast = FALSE, save.memory = FALSE, collapse = TRUE)
# Compute posterior probabilities of edge inclusion and Median Probability DAG Model
# from the MCMC outputs [2] and [3]
get_edgeprobs(out_mcmc_collapse)
# get_edgeprobs(out_mcmc_collapse_fast)
get_MPMdag(out_mcmc_collapse)
# get_MPMdag(out_mcmc_collapse_fast)

# Methods
print(out_mcmc)
summary(out_mcmc)
plot(out_mcmc)