We consider the problem of jointly estimating multiple related directed
acyclic graph (DAG) models based on high-dimensional data from each graph. This
problem is motivated by the task of learning gene regulatory networks based on
gene expression data from different tissues, developmental stages or disease
states. We prove that under certain regularity conditions, the proposed
β0β-penalized maximum likelihood estimator converges in Frobenius norm to
the adjacency matrices consistent with the data-generating distributions and
has the correct sparsity. In particular, we show that this joint estimation
procedure leads to a faster convergence rate than estimating each DAG model
separately. As a corollary, we also obtain high-dimensional consistency results
for causal inference from a mix of observational and interventional data. For
practical purposes, we propose \emph{jointGES} consisting of Greedy Equivalence
Search (GES) to estimate the union of all DAG models followed by variable
selection using lasso to obtain the different DAGs, and we analyze its
consistency guarantees. The proposed method is illustrated through an analysis
of simulated data as well as epithelial ovarian cancer gene expression data