# Discovery and inference of a causal network with hidden confounding††thanks:
L. Chen and C. Li contributed equally. Research is supported in part by NSF
grant DMS-1952539, NIH grants R01GM113250, R01GM126002, R01AG065636,
R01AG074858, R01AG069895, U01AG073079. The authors report there are no
competing interests to declare. The authors thank the editor, the associate
editor, and three anonymous referees for their helpful comments and
suggestions. C. Li would like to thank R. Oliver VandenBerg for the
suggestions on writing.

Li Chen111School of Statistics, University of Minnesota, Minneapolis, MN
55455. Chunlin Li222To whom correspondence should be addressed. Department of
Statistics, Iowa State University, Ames, IA 50011. Email<EMAIL_ADDRESS>Xiaotong Shen333School of Statistics, University of Minnesota, Minneapolis, MN
55455. Wei Pan444Division of Biostatistics, University of Minnesota,
Minneapolis, MN 55455.

###### Abstract

This article proposes a novel causal discovery and inference method called
GrIVET for a Gaussian directed acyclic graph with unmeasured confounders.
GrIVET consists of an order-based causal discovery method and a likelihood-
based inferential procedure. For causal discovery, we generalize the existing
peeling algorithm to estimate the ancestral relations and candidate
instruments in the presence of hidden confounders. Based on this, we propose a
new procedure for instrumental variable estimation of each direct effect by
separating it from any mediation effects. For inference, we develop a new
likelihood ratio test of multiple causal effects that is able to account for
the unmeasured confounders. Theoretically, we prove that the proposed method
has desirable guarantees, including robustness to invalid instruments and
uncertain interventions, estimation consistency, low-order polynomial time
complexity, and validity of asymptotic inference. Numerically, GrIVET performs
well and compares favorably against state-of-the-art competitors. Furthermore,
we demonstrate the utility and effectiveness of the proposed method through an
application inferring regulatory pathways from Alzheimer’s disease gene
expression data.

Keywords: Causal discovery, Gaussian directed acyclic graph, Invalid
instrumental variables, Uncertain interventions, Simultaneous inference, Gene
regulatory network.

##  1 Introduction

Understanding causal relations is part of the foundation of intelligence. A
directed acyclic graph (DAG) is often used to describe the causal relations
among multiple interacting units (Pearl,, 2009). Unlike classical causal
inference tasks where the DAG is determined a priori, causal discovery aims to
learn a graphical representation from data. It is useful for forming data-
driven conjectures about the underlying mechanism of a complex system,
including gene networks (Sachs et al.,, 2005), functional brain networks (Liu
et al.,, 2017), manufacturing pipelines (Kertel et al.,, 2022), and dynamical
systems (Li et al., 2020b, ). In such a situation, randomized experiments are
usually unethical or infeasible, and unmeasured confounders commonly arise in
practice. The presence of latent confounders can bias the causal effect
estimation and even distort causal directions, making causal discovery
challenging. To treat latent confounders, we use additive interventions as
instrumental variables (IVs), which are well-developed in conventional causal
inference (Angrist et al.,, 1996) yet are less explored in causal discovery of
a large-scale network. In this article, we focus on a Gaussian DAG model with
hidden confounders and develop methods that integrate the discovery and
inference of causal relations within the framework of uncertain additive
interventions (the targets of interventions are unknown).

Causal discovery has been extensively studied (Zheng et al.,, 2018; Aragam et
al.,, 2019; Gu et al.,, 2019; Lee and Li,, 2022; Zhao et al.,, 2022; Li et
al., 2023b, ); see Drton and Maathuis, (2017); Heinze-Deml et al., (2018);
Glymour et al., (2019); Vowels et al., (2021) for comprehensive reviews. For
observational data (without external interventions), some methods are able to
treat hidden confounding by either (a) producing less informative discoveries,
like a partial ancestral graph (Colombo et al.,, 2012) rather than a DAG, or
(b) employing a certain deconfounding strategy (Frot et al.,, 2019; Shah et
al.,, 2020) based on the pervasive confounding assumption. However, the former
may not reveal essential information, such as causal directions, while the
latter can be inconsistent in low-dimensional situations and may not
necessarily outperform the naive regression (Grimmer et al.,, 2020). Thus,
external interventions are useful to provide more information about causal
relations while relaxing the requirements on latent confounding.

As an example of external (additive) interventions, IVs have been well
developed in conventional causal inference to tackle unmeasured confounding;
see Lousdal, (2018) for a survey. In a classical bivariate setting where the
causal direction is known, an IV is required to influence the response
variable only through the cause variable, which is often fragile in practice
(Murray,, 2006). For instance, genetic variants like single nucleotide
polymorphisms (SNPs) are used as IVs in Mendelian randomization (MR) analysis
to discover putative causal genes of complex traits, where the IV conditions
are commonly violated due to the (horizontal) pleiotropy. Remedying these
invalid IVs has been the subject of recent work in causal inference (Kang et
al.,, 2016; Guo et al.,, 2018; Windmeijer et al.,, 2019; Burgess et al.,,
2020). The discussion of IV estimation in graphical modeling, however, remains
limited. The methods of Oates et al., (2016); Chen et al., (2018) estimate the
graph given valid IVs, while the work of Li et al., 2023a  propose the peeling
algorithm to construct the DAG in the case of uncertain interventions and
invalid IVs. None of these methods permit latent confounding. A recent work
(Xue and Pan,, 2020) discusses causal discovery of a bivariate mixed effect
graph where confounders and invalid IVs are allowed, but it remains unclear
how to extend it to a large-scale causal network.

Moreover, despite the progress in causal discovery, inference about the
discovered relations is often regarded as a separate task and has received
less attention in the literature. Notable exceptions include recent advances
in graphical modeling (Janková and van de Geer,, 2018; Li et al., 2020a, ; Shi
et al.,, 2023; Wang et al.,, 2023) and mediation analysis (Chakrabortty et
al.,, 2018; Shi and Li,, 2021; Li et al.,, 2022); however, these methods
cannot account for latent confounders. Indeed, due to unmeasured confounding,
the probability distribution of observed variables is no longer locally
Markovian with respect to the DAG (Pearl,, 2009), rendering these approaches
inappropriate. Consequently, there is a pressing need for new inference
methodologies.

This article contributes to the following aspects.

  * •

For modeling, we establish the identifiability conditions for a Gaussian DAG
with latent confounders utilizing additive interventions. To our knowledge,
this result is the first of its kind. Importantly, the conditions allow the
interventions to have unknown and multiple targets, which is suitable for
multivariate causal analysis (Murray,, 2006).

  * •

For methodology, we develop a novel method named the Graphical Instrumental
Variable Estimation and Testing (GrIVET), integrating order-based causal
discovery and likelihood-based inference. For causal discovery, we estimate
the ancestral relations and candidate IVs with a modified peeling algorithm to
treat unmeasured confounding. On this basis, we propose a sequential procedure
to estimate each direct effect using IVs, where a working response regression
is used to separate the direct effect from the mediation effects. Regarding
inference, we develop a new likelihood ratio test of multiple causal effects
to account for unmeasured confounders.

  * •

For theory, we show that GrIVET enjoys desired guarantees. In particular, it
consistently estimates the DAG structure and causal effects even when some
interventions do not meet the IV criteria. As for computation, only
$O((p+|\mathcal{E}^{+}|)\times\log(s)\times(q^{3}+nq^{2}))$ operations are
required almost surely, where $p$ and $q$ are the numbers of primary and
intervention variables, $s$ is sparsity, $|\mathcal{E}^{+}|$ is the size of
the ancestral relation set, and $n$ is the sample size. Moreover, under the
null hypothesis, we establish the convergence of the likelihood ratio
statistic to the null distribution in high-dimensional situations, ensuring
the validity of asymptotic inference.

  * •

The simulation studies and an application to the Alzheimer’s Disease
Neuroimaging Initiative dataset demonstrate the utility and effectiveness of
the proposed methods. The implementation of GrIVET is available at
https://github.com/chunlinli/grivet.

The rest of the article is structured as follows. Section 2 introduces a
linear structural equation model with hidden confounders and establishes its
identifiability. Section 3 presents a novel order-based method for causal
discovery and effect estimation. Section 4 develops a likelihood ratio test
for simultaneous inference of causal effects. Section 5 provides theoretical
justification of the proposed method. Section 6 performs simulation studies,
followed by an application to infer gene pathways with gene expression and SNP
data. Finally, Section 7 concludes the article. The Appendix contains
supporting lemmas, while the Supplementary Materials include illustrative
examples, technical proofs, and additional simulations.

##  2 Causal graphical model with confounders

###  2.1 Structural equations with confounders

We consider a structural equation model with $p$ primary variables
$\bm{Y}=(Y_{1},\ldots,Y_{p})^{\top}$ and $q$ intervention variables
$\bm{X}=(X_{1},\ldots,X_{q})^{\top}$,

$\begin{split}\bm{Y}=\mathbf{U}^{\top}\bm{Y}+\mathbf{W}^{\top}\bm{X}+\bm{\varepsilon},\quad\bm{\varepsilon}\sim
N(\bm{0},\bm{\Sigma}),\quad\operatorname{Cov}(\bm{\varepsilon},\bm{X})=\bm{0},\end{split}$
(1)

where $\mathbf{U}_{p\times p}$ is a matrix describing the causal influences
among $\bm{Y}$, $\mathbf{W}_{q\times p}$ is a matrix representing the
interventional effects of $\bm{X}$ on $\bm{Y}$, and $\bm{\varepsilon}$ is a
vector of possibly correlated errors. Specifically,

  * •

The parameter matrix $\mathbf{U}$, which is of primary interest, has a causal
interpretation in that $\mathrm{U}_{kj}\neq 0$ indicates that $Y_{k}$ is a
cause of $Y_{j}$, denoted by $Y_{k}\to Y_{j}$. Thus, $\mathbf{U}$ represents a
directed graph among primary variables. In what follows, we will focus on a
directed acyclic graph (DAG), where no directed cycle is permissible and
$\mathbf{U}$ is subject to the acyclicity constraint (Zheng et al.,, 2018;
Yuan et al.,, 2019).

  * •

The intervention variables $\bm{X}$ and errors $\bm{\varepsilon}$ are
uncorrelated by reparameterization. As a result, $\mathbf{W}$ is associational
instead of causal. Here, $\mathrm{W}_{lj}\neq 0$ indicates that $X_{l}$
intervenes on $Y_{j}$, denoted by $X_{l}\to Y_{j}$. As $\bm{X}$ represents
external interventions, no directed edge from a primary variable to an
intervention variable is allowed.

  * •

A non-diagonal $\bm{\Sigma}$ indicates the presence of unmeasured confounders.
For instance, $\bm{\varepsilon}=\bm{\Phi}^{\top}\bm{\eta}+\bm{e}$ can be (not
uniquely) written as a sum of correlated components
$\bm{\Phi}^{\top}\bm{\eta}$ and independent components $\bm{e}$ so that
$\bm{\Sigma}=\bm{\Phi}^{\top}\bm{\Phi}+\operatorname{Diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2})$,
where $\bm{\Phi}_{r\times p}$ is the matrix of confounding effects,
$\bm{\eta}\sim N(\bm{0},\mathbf{I}_{r\times r})$ represents $r$ independent
confounding sources, and $\bm{e}\sim
N(\bm{0},\operatorname{Diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2}))$
represents $p$ independent errors. Whenever $\Sigma_{jk}\neq 0$ for some
distinct $(j,k)$, we have $\Sigma_{jk}=\sum_{m=1}^{r}\Phi_{mj}\Phi_{mk}\neq
0$, implying that some confounding variable $\eta_{m}$ influences both $Y_{j}$
and $Y_{k}$.

As such, $(\mathbf{U},\mathbf{W})$ together represents a directed graph of $p$
primary variables and $q$ intervention variables, denoted as
$\mathcal{G}=(\bm{X},\bm{Y};\mathcal{E},\mathcal{I})$, where
$\mathcal{E}=\\{(k,j):\mathrm{U}_{kj}\neq 0\\}$ is the set of primary variable
edges and $\mathcal{I}=\\{(l,j):\mathrm{W}_{lj}\neq 0\\}$ is the set of
intervention edges. In $\mathcal{G}$, (a) if $Y_{k}\to Y_{j}$, then $Y_{k}$ is
a parent of $Y_{j}$, and $Y_{j}$ is a child of $Y_{k}$, (b) if
$Y_{k}\to\cdots\to Y_{j}$ (a directed path from $Y_{k}$ to $Y_{j}$), then
$Y_{k}$ is an ancestor of $Y_{j}$, and $Y_{j}$ is a descendant of $Y_{k}$, and
(c) if $Y_{k}\to\cdots\to Y_{m}\to\cdots\to Y_{j}$, then $Y_{m}$ is a mediator
of $Y_{k}$ and $Y_{j}$. In what follows, for a graph $\mathcal{G}$, denote the
parent set of $Y_{j}$ as $\textnormal{{pa}}_{\mathcal{G}}(j)=\\{k:Y_{k}\to
Y_{j}\\}$, the ancestor set of $Y_{j}$ as
$\textnormal{{an}}_{\mathcal{G}}(j)=\\{k:Y_{k}\to\cdots\to Y_{j}\\}$, and the
intervention set of $Y_{j}$ as $\textsc{in}_{\mathcal{G}}(j)=\\{l:X_{l}\to
Y_{j}\\}$. For $(k,j)$ such that $Y_{k}\to\cdots\to Y_{j}$, denote the
mediator set as $\textnormal{{me}}_{\mathcal{G}}(k,j)=\\{m:Y_{k}\to\cdots\to
Y_{m}\to\cdots\to Y_{j}\\}$.

###  2.2 Identifiability and instrumental variables

The causal parameter matrix $\mathbf{U}$ is generally non-identifiable555The
causal parameter $\mathbf{U}$ is said to be identifiable if for any
$(\mathbf{U},\mathbf{W},\bm{\Sigma})$ and
$(\mathbf{U}^{\prime},\mathbf{W}^{\prime},\bm{\Sigma}^{\prime})$, we have
$\mathbb{P}_{\mathbf{U},\mathbf{W},\bm{\Sigma}}=\mathbb{P}_{\mathbf{U}^{\prime},\mathbf{W}^{\prime},\bm{\Sigma}^{\prime}}$
implies $\mathbf{U}=\mathbf{U}^{\prime}$. Otherwise, it is said to be non-
identifiable. without further conditions on the Gaussian errors
$\bm{\varepsilon}$ or the interventions $\bm{X}$. Without invoking external
interventions ($\mathbf{W}\equiv\bm{0}$), $\mathbf{U}$ can be identified under
a certain error-scale assumption (Peters and Bühlmann,, 2014; Ghoshal and
Honorio,, 2018; Rajendran et al.,, 2021), which is sensitive to variable
scaling such as the common practice of standardizing variables (Reisach et
al.,, 2021). To overcome this limitation, interventions are introduced to
identify the causal parameters. With suitable interventions, $\mathbf{U}$ is
identifiable if no confounder is present in the model ($\bm{\Sigma}$ is
diagonal) (Oates et al.,, 2016; Chen et al.,, 2018; Li et al., 2023a, ). In
addition, it is worth mentioning that $\mathbf{U}$ can be estimated without
intervention if the errors $\bm{\varepsilon}$ are non-Gaussian (Shimizu et
al.,, 2006; Zhao et al.,, 2022); however, such methods are not applicable in
the case of unmeasured confounding.

This subsection establishes the identifiability of (1) in the presence of
unmeasured confounders using uncertain additive interventions (the targets of
interventions are unknown) as IVs. To proceed, we introduce the notion of IV
for our purpose.

######  Definition 1.

An intervention variable $X_{l}$ is said to be a valid IV of $Y_{k}$ in
$\mathcal{G}$ if (IV1) $X_{l}$ intervenes on $Y_{k}$, namely
$\mathrm{W}_{lk}\neq 0$, and (IV2) $X_{l}$ does not intervene on any other
primary variable $Y_{k^{\prime}}$, namely $\mathrm{W}_{lk^{\prime}}=0$ for
$k^{\prime}\neq k$. Otherwise, $X_{l}$ is called an invalid IV. Denote the
valid IV set of $Y_{k}$ as
$\textnormal{{iv}}_{\mathcal{G}}(k)=\\{l:X_{l}\rightarrow Y_{k},X_{l}\not\to
Y_{k^{\prime}},k^{\prime}\neq k\\}$.

######  Remark 1.

Consider a bivariate case where we are interested in the potential causal
effect $Y_{1}\to Y_{2}$. In causal inference literature (Angrist et al.,,
1996; Kang et al.,, 2016), a valid IV $X$ of $Y_{1}$ is required to satisfy
that (a) $X$ is related to the $Y_{1}$, referred to as relevance, (b) $X$ has
no directed edge to $Y_{2}$, called exclusion, and (c) $X$ is not related to
unmeasured confounders, called unconfoundedness. In (1), (IV1) is indeed the
relevance property, (IV2) generalizes the exclusion property for causal
discovery, and the requirement
$\operatorname{Cov}(\bm{\varepsilon},\bm{X})=\bm{0}$ corresponds to the
unconfoundedness.

To identify $\mathbf{U}$, two challenges emerge as the confounders arise.
First, determining causal directions in the graph becomes more challenging. In
(1), because of hidden confounding, the distribution
$\mathbb{P}(\bm{Y}\mid\bm{X})$ does not admit the causal Markov property
(Pearl,, 2009) according to $\mathcal{G}$, that is, $Y_{j}$ is not independent
of its non-descendants given
$(\bm{Y}_{\textnormal{{pa}}_{\mathcal{G}}(j)},\bm{X})$. As a result, the
existing methods based on this property can learn wrong causal directions due
to misspecification. To identify causal directions, we formalize the concept
of unmediated parents to highlight the causal relations that are critical in
identification.

######  Definition 2.

A primary variable $Y_{k}$ is an unmediated parent of $Y_{j}$ in $\mathcal{G}$
if $Y_{k}\to Y_{j}$ and there is no other directed path from $Y_{k}$ to
$Y_{j}$. In other words, $Y_{k}$ is an unmediated parent of $Y_{j}$ if no
mediator is between $Y_{k}$ and $Y_{j}$.

Another challenge comes from uncertain interventions and invalid IVs.
Assigning valid IVs for each primary variable can be difficult when the
targets of interventions are unknown. Thus, it may be effective to construct a
set of candidate IVs (including invalid IVs) for each primary variable, on
which we estimate the causal parameters $\mathbf{U}$. To this end, we define
$p$ candidate IV sets, one for each primary variable.

######  Definition 3.

An intervention variable $X_{l}$ is said to be a candidate IV of $Y_{k}$ in
$\mathcal{G}$ if (IV1’) $X_{l}$ intervenes on $Y_{k}$, and (IV2’) $X_{l}$ does
not intervene on any non-descendant of $Y_{k}$. Denote the candidate IV set of
$Y_{k}$ by $\textnormal{{ca}}_{\mathcal{G}}(k)=\\{l:X_{l}\rightarrow
Y_{k},X_{l}\rightarrow Y_{j}\text{ only if
}k\in\textnormal{{an}}_{\mathcal{G}}(j)\\}$.

The candidate IVs of $Y_{k}$ include all valid IVs of $Y_{k}$, but not vice
versa. A candidate IV of $Y_{k}$ may be invalid, as it could intervene on
descendants of $Y_{k}$.

######  Theorem 1 (Identifiability).

Suppose

  1. (A1)

$\operatorname{Cov}(\bm{X})$ is positive definite.

  2. (A2)

$\operatorname{Cov}(Y_{j},X_{l}\mid\bm{X}_{\\{1,\ldots,q\\}\setminus\\{l\\}})\neq
0$ whenever $X_{l}$ intervenes on an unmediated parent of $Y_{j}$.

  3. (A3)

_(Majority rule)_
$|\textnormal{{iv}}_{\mathcal{G}}(k)|>|\textnormal{{ca}}_{\mathcal{G}}(k)|/2$;
$k=1,\ldots,p$.

Then $(\mathbf{U},\mathbf{W},\bm{\Sigma})$ in (1) are identifiable in that if
$(\mathbf{U},\mathbf{W},\bm{\Sigma})$ and
$(\mathbf{U}^{\prime},\mathbf{W}^{\prime},\bm{\Sigma}^{\prime})$ encode the
same probability distribution, then
$(\mathbf{U},\mathbf{W},\bm{\Sigma})=(\mathbf{U}^{\prime},\mathbf{W}^{\prime},\bm{\Sigma}^{\prime})$.

To our knowledge, Theorem 1 is a new result for Gaussian DAG with hidden
confounding, establishing the identifiability of all parameters in (1). In
fact, if the causal parameter $\mathbf{U}$ is identifiable, then so are
parameters $\mathbf{W},\bm{\Sigma}$. Regarding the conditions, (A1) states
that $\operatorname{Cov}(\bm{X})$ has full rank, which is common in the IV
literature (Kang et al.,, 2016; Chen et al.,, 2018). Note that (A1) permits
discrete IV variables such as SNPs in data analysis. (A2) requires the
interventional effects through unmediated parents not to cancel out when an
invalid IV has multiple targets. (A3) requires valid IVs to dominate invalid
ones so that the causal effect can be identified in the presence of latent
confounders. Such a condition has been used in the causal inference literature
(Kang et al.,, 2016; Windmeijer et al.,, 2019). As shown in Supplementary
Materials Section 1, when (A3) fails, (1) can be non-identifiable. By
comparison, (A1)–(A2) together with (A4) are used for model identification in
the absence of unmeasured confounding (Li et al., 2023a, ).

  * (A4)

Each $Y_{k}$ is intervened by at least one valid IV.

Noting that (A4) is implied by (A3), treating hidden confounding demands
stronger conditions in view of Theorem 1.

##  3 Causal discovery

This section proposes a novel IV method to learn a DAG with unmeasured
confounders. First, we introduce the ancestral relation graph (ARG), which,
together with the candidate IV sets in Section 2.2, constitutes a basis for
the proposed method.

######  Definition 4 (Ancestral relation graph).

For a DAG $\mathcal{G}=(\bm{X},\bm{Y};\mathcal{E},\mathcal{I})$, its ancestral
relation graph is defined as
$\mathcal{G}^{+}=(\bm{X},\bm{Y};\mathcal{E}^{+},\mathcal{I}^{+})$, where

$\mathcal{E}^{+}=\Big{\\{}(k,j):k\in\textnormal{{an}}_{\mathcal{G}}(j)\Big{\\}},\qquad\mathcal{I}^{+}=\Big{\\{}(l,j):l\in\bigcup_{k\in\textnormal{{an}}_{\mathcal{G}}(j)\cup\\{j\\}}\textsc{in}_{\mathcal{G}}(k)\Big{\\}}.$

Here, $\mathcal{G}^{+}$ is a super-DAG of $\mathcal{G}$ in that
$\mathcal{E}^{+}\supseteq\mathcal{E}$ is the set of ancestral relations,
$\mathcal{I}^{+}\supseteq\mathcal{I}$ is a superset of interventional
relations, and $\mathcal{G}^{+}$ is acyclic. Note that $\mathcal{E}^{+}$
defines a partial order for the primary variables $\bm{Y}$ in that
$Y_{k}\prec_{\mathcal{G}}Y_{j}$ whenever $(k,j)\in\mathcal{E}^{+}$. Without
confounding, $\mathbf{U}$ can be consistently estimated via direct regressions
according to the known $\mathcal{G}^{+}$ (Shojaie and Michailidis,, 2010),
where $\mathcal{G}^{+}$ can be recovered by the peeling algorithm (Li et al.,
2023a, ). However, this approach no longer applies in the presence of hidden
confounders.

To address this obstacle, Sections 3.1–3.2 modify the peeling algorithm to
construct the ARG $\mathcal{G}^{+}$ and the candidate IV sets
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$, and then Sections
3.3–3.4 develop a method to estimate $\mathbf{U}$ assuming the ARG and
candidate IVs are known.

###  3.1 Identification of $\mathcal{G}^{+}$ and candidate IVs

In this subsection, we modify the peeling algorithm, originally designed for a
model without unmeasured confounders (Li et al., 2023a, ), to uncover
$\mathcal{G}^{+}$ and $\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq
p}$ in the presence of hidden confounders, of which the results can be
subsequently used as the inputs for identification of $\mathbf{U}$ in Section
3.3. The modified peeling algorithm essentially requires $p$ regressions to
identify the ARG and candidate IVs, which is suited for large-scale causal
discovery. Moreover, the produced ARG and candidate IV sets enjoy desirable
statistical properties; see Section 5.

Let us begin with an observation that (1) can be rewritten as

$\bm{Y}=\mathbf{V}^{\top}\bm{X}+(\mathbf{I}-\mathbf{U}^{\top})^{-1}\bm{\varepsilon},$
(2)

where $\mathbf{V}=\mathbf{W}(\mathbf{I}-\mathbf{U})^{-1}$ and
$\mathrm{V}_{lj}=\sum_{k=1}^{p}\mathrm{W}_{lk}(\mathrm{I}_{kj}+\mathrm{U}_{kj}+\cdots+(\bm{\mathbf{U}}^{p-1})_{kj})$.
Intuitively, $\mathrm{V}_{lj}\neq 0$ implies the dependence of $Y_{j}$ on
$X_{l}$ through a directed path $X_{l}\to Y_{k}\to\cdots\to Y_{j}$, and hence
that $X_{l}$ intervenes on $Y_{j}$ itself (when $k=j$) or its ancestor $Y_{k}$
(when $k\neq j$). In cases where $X_{l}$ intervenes exclusively on one primary
variable, the following proposition provides insights into the connection
between $\mathbf{V}$ and $\mathcal{G}^{+}$.

######  Proposition 1.

Suppose Assumptions (A1), (A2), and (A4) are satisfied. There exists at least
one intervention variable $X_{l}$ such that $\mathrm{V}_{lk}\neq 0$ and
$\mathrm{V}_{lk^{\prime}}=0$ for $k^{\prime}\neq k$ if and only if $Y_{k}$ is
a leaf node (has no descendant). Moreover, such $X_{l}$ is a valid IV of
$Y_{k}$ in $\mathcal{G}$.

Proposition 1 suggests that the leaves and their valid IVs in $\mathcal{G}$
can be identified by

$\begin{split}\textsc{leaf}(\mathcal{G})&=\\{k:\text{ for some $l$,
}\mathrm{V}_{lk}\neq 0\text{ and }\mathrm{V}_{lk^{\prime}}=0\text{ for all
}k^{\prime}\neq k\\}\\\
&=\\{k:k=\operatorname*{\arg\max}_{j}|\mathrm{V}_{lj}|\text{ for some
}l=\operatorname*{\arg\min}_{\|\mathbf{V}_{l,+}\|_{0}>0}\|\mathbf{V}_{l,+}\|_{0}\\},\\\
\textnormal{{iv}}_{\mathcal{G}}(k)&=\\{l:\mathrm{V}_{lk}\neq 0\text{ and
}\mathrm{V}_{lk^{\prime}}=0\text{ for all }k^{\prime}\neq k\\}\\\
&=\\{l:l=\operatorname*{\arg\min}_{\|\mathbf{V}_{l,+}\|_{0}>0}\|\mathbf{V}_{l,+}\|_{0}\text{
and }k=\operatorname*{\arg\max}_{j}|\mathrm{V}_{lj}|\\},\quad
k\in\textsc{leaf}(\mathcal{G}).\end{split}$ (3)

After the leaf nodes are learned, we can remove them to obtain a sub-DAG. If
$X_{l}$ is a valid IV of a non-leaf $Y_{k}$ in $\mathcal{G}$, its validity for
$Y_{k}$ is retained in the sub-DAG, implying (A4) continues to hold. Moreover,
Assumptions (A1)–(A2) are naturally upheld in the sub-DAG. Hence, the
requirements of Proposition 1 are satisfied in the sub-DAG, whose leaf
variables and their valid IVs can be learned in the same fashion. As a result,
we can successively identify and remove (i.e., peel) the leaf nodes from the
DAG and sub-DAGs. This yields a topological order of primary variables but
does not recover $\mathcal{G}^{+}$.

Next, we investigate how $\mathbf{V}$ can be further used to recover
$\mathcal{G}^{+}$ with $\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq
p}$. Subsequently, we use
$\mathcal{G}^{-}=(\bm{X}^{-},\bm{Y}^{-};\mathcal{E}^{-},\mathcal{I}^{-})$ to
denote a generic sub-DAG produced by peeling, where $\bm{Y}^{-}$ are the
primary variables in $\mathcal{G}^{-}$ and $\bm{Y}\setminus\bm{Y}^{-}$ are
peeled ones, $\bm{X}^{-}$ are intervention variables on $\bm{Y}^{-}$,
$\mathcal{E}^{-}$ is the set of causal relations among $\bm{Y}^{-}$, and
$\mathcal{I}^{-}$ is the set of interventional relations between $\bm{X}^{-}$
and $\bm{Y}^{-}$. Then each variable in $\bm{Y}^{-}$ is a non-descendant of
each in $\bm{Y}\setminus\bm{Y}^{-}$. Moreover,
$\textsc{leaf}(\mathcal{G}^{-})$ and
$\\{\textnormal{{iv}}_{\mathcal{G}^{-}}(k)\\}_{k\in\textsc{leaf}(\mathcal{G}^{-})}$
are identified by (3).

######  Proposition 2.

Suppose Assumptions (A1), (A2), and (A4) are satisfied. Let $Y_{k}$ be a leaf
node in $\mathcal{G}^{-}$ and $Y_{j}$ be in $\bm{Y}\setminus\bm{Y}^{-}$. Then
the following statements are true.

  1. (A)

If $\mathrm{V}_{lj}\neq 0$ for all
$l\in\textnormal{{iv}}_{\mathcal{G}^{-}}(k)$, we have
$(k,j)\in\mathcal{E}^{+}$.

  2. (B)

If $Y_{k}$ is an unmediated parent of $Y_{j}$, then $\mathrm{V}_{lj}\neq 0$
for all $l\in\textnormal{{iv}}_{\mathcal{G}^{-}}(k)$.

Proposition 2 outlines a method for identifying edges in $\mathcal{G}^{+}$
from the leaf variables of $\mathcal{G}^{-}$ to the peeled variables
$\bm{Y}\setminus\bm{Y}^{-}$ by

$\\{(k,j):Y_{k}\in\textsc{leaf}(\mathcal{G}^{-}),\
Y_{j}\in\bm{Y}\setminus\bm{Y}^{-}\text{ and }\mathrm{V}_{lj}\neq 0\text{ for
all }l\in\textnormal{{iv}}_{\mathcal{G}^{-}}(k)\\}.$ (4)

Specifically, (A) shows that any identified edge must be present in
$\mathcal{G}^{+}$, so no extra edges are identified. Meanwhile, (B) shows that
every directed edge from an unmediated parent must be correctly discovered.
Importantly, the collection of all such edges suffices to recover all
ancestral relationships, which guarantees that no edge in $\mathcal{E}^{+}$ is
overlooked. Upon the identification of $\mathcal{G}^{+}$, the candidate IV
sets can be learned by

$\textnormal{{ca}}_{\mathcal{G}}(k)=\\{l:(l,k)\in{\mathcal{I}}^{+}\text{ and
}(l,j)\in{\mathcal{I}}^{+},k\neq j\text{ only if
}(k,j)\in{\mathcal{E}}^{+}\\},\quad 1\leq k\leq p.$ (5)

Consequently, Propositions 1–2 enable the recovery of $\mathcal{G}^{+}$ and
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$.

###  3.2 Finite-sample estimation of $\mathcal{G}^{+}$ and candidate IVs

This subsection implements the modified peeling algorithm delineated in
Section 3.1 to estimate $\mathcal{G}^{+}$ and
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$. To proceed,
suppose data matrices $\mathbf{Y}_{p\times
n}=(\mathbf{Y}_{+,1},\ldots,\mathbf{Y}_{+,n})$ and $\mathbf{X}_{q\times
n}=(\mathbf{X}_{+,1},\ldots,\mathbf{X}_{+,n})$ are given, where
$(\mathbf{Y}_{+,i},\mathbf{X}_{+,i})_{i=1}^{n}$ are sampled from (1)
independently. We estimate $\mathbf{V}$ by
$\widehat{\mathbf{V}}=(\widehat{\mathbf{V}}_{+,1},\ldots,\widehat{\mathbf{V}}_{+,p})$
with sparse regressions

$\widehat{\mathbf{V}}_{+,j}=\operatorname*{\arg\min}_{\bm{\beta}}\
\sum_{i=1}^{n}(\mathrm{Y}_{j,i}-\bm{\beta}^{\top}\mathbf{X}_{+,i})^{2}\quad\text{s.t.}\quad\|\bm{\beta}\|_{0}\leq\kappa^{\prime}_{j}$
(6)

where $1\leq\kappa^{\prime}_{j}\leq q$ is tuned by BIC for $1\leq j\leq p$
Moreover, the truncated Lasso penalty (TLP) (Shen et al.,, 2012) is used as
the computational surrogate for $\|\cdot\|_{0}$, where TLP is defined as
$\operatorname{TLP}_{\tau}(\bm{\beta})=\sum_{j=1}^{r}\min(|\beta_{j}|/\tau,1)$
for $\bm{\beta}=(\beta_{1},\ldots,\beta_{r})$, and $\tau>0$ is a
hyperparameter in TLP; see Supplementary Materials Section 2 for details. The
modified peeling algorithm based on Section 3.1 is summarized in Algorithm
1.666In Algorithm 1 Step 7, the indices of $\mathbf{V}$ are kept so that
$\mathrm{V}_{lj}$ always represents the effect from $X_{l}$ to $Y_{j}$.

Input: Data $\mathbf{Y}_{p\times n}$ and $\mathbf{X}_{q\times n}$;

1

2Compute $\widehat{\mathbf{V}}$ via (6);

3

4Initialize $\mathbf{V}\leftarrow\widehat{\mathbf{V}}$,
$\widehat{\mathcal{E}}^{+}\leftarrow\emptyset$,
$\widehat{\mathcal{I}}^{+}\leftarrow\\{(l,k):\widehat{\mathrm{V}}_{lk}\neq
0\\}$;

5

6Initialize $\mathcal{G}^{-}$ by $\bm{Y}^{-}\leftarrow\bm{Y}$,
$\bm{X}^{-}\leftarrow\bm{X}$,
$\mathcal{E}^{-}\leftarrow\widehat{\mathcal{E}}^{+}$,
$\mathcal{I}^{-}\leftarrow\widehat{\mathcal{I}}^{+}$;

7

8while  _${\bm{Y}}^{-}$ is not empty_ do

9

10 Update $\textsc{leaf}(\mathcal{G}^{-})$ and
$\\{\textnormal{{iv}}_{\mathcal{G}^{-}}(k)\\}_{k\in\textsc{leaf}(\mathcal{G}^{-})}$
via (3);

11

12

13 Update $\widehat{\mathcal{E}}^{+}$ by adding (4);

14

15 Update $\mathcal{G}^{-}$ by removing $\textsc{leaf}(\mathcal{G}^{-})$ and
${\mathbf{V}}$ by keeping the columns in $\bm{Y}^{-}$;

16

17 end while

18

19Update $\widehat{\mathcal{E}}^{+}\leftarrow\\{(k,j):Y_{k}\to\cdots\to
Y_{j}\text{ in }\widehat{\mathcal{E}}^{+}\\}$;

20 Update
$\widehat{\mathcal{I}}^{+}\leftarrow\\{(l,j):(l,k)\in\widehat{\mathcal{I}}^{+}\text{
and }(k,j)\in\widehat{\mathcal{E}}^{+}\\}$;

21 Update $\widehat{\textnormal{{ca}}}_{\mathcal{G}}(k)$ by (5);

22

23return _$\widehat{\mathcal{E}}^{+}$ , $\widehat{\mathcal{I}}^{+}$, and
$\\{\widehat{\textnormal{{ca}}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$_;

Algorithm 1 Estimation of $\mathcal{G}^{+}$ and
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$

###  3.3 Identification of $\mathbf{U}$

In this subsection, we present a new method for identifying causal effects
$\mathbf{U}$, using the ARG $\mathcal{G}^{+}$ and candidate IV sets
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$ as inputs. Note
that $\\{\textnormal{{an}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$ and
$\\{\textnormal{{me}}_{\mathcal{G}}(k,j)\\}_{(k,j)\in\mathcal{E}^{+}}$ can be
derived from $\mathcal{G}^{+}$. Throughout this subsection, the subscript
$\mathcal{G}$ is dropped for brevity and $\bm{\alpha},\bm{\beta},\bm{\gamma}$
denote nuisance parameters in regression. Moreover, we assume that
$\bm{\varepsilon}$ and $\bm{X}$ are independent to simplify the derivation;
see Lemmas 1–2 in the Appendix for the case with $\bm{\varepsilon}$ and
$\bm{X}$ being uncorrelated.

Figure 1: Estimation of causal parameter $\mathrm{U}_{kj}$. (a) Display of the
relations among relevant variables. (b) Display of working response
regression.

#### The case with all IVs being valid.

We begin with a special case of (1) where all IVs are valid, that is,
$\textnormal{{ca}}(k)=\textnormal{{iv}}(k)$; $k=1,\ldots,p$.

To estimate $\mathbf{U}$, note that $\mathbf{U}$ is supported on
$\mathcal{E}^{+}$, namely $\mathbf{U}=(\mathbf{U}_{\mathcal{E}^{+}},\bm{0})$.
Here, we consider estimating $\mathrm{U}_{kj}$, as well as selecting nonzero
$\mathrm{U}_{kj}$ for graph recovery, for each $(k,j)\in\mathcal{E}^{+}$, as
described in Figure 1 (a).

To pinpoint the difficulties and motivate our approach, we make the following
observations. First, regression of $Y_{j}$ on $Y_{k}$ together with covariates
$(\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}},\bm{X})$ can bias the
estimation due to confounder $\eta$. Second, in hope of treating confounders
one might replace $Y_{k}$ with its surrogate
$\mathbb{E}(Y_{k}\mid\bm{Y}_{\textnormal{{an}}(k)},\bm{X})$ to regress $Y_{j}$
on $\mathbb{E}(Y_{k}\mid\bm{Y}_{\textnormal{{an}}(k)},\bm{X})$ with
$(\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}},\bm{X}_{\textnormal{{iv}}(k)^{c}})$
being covariates. However, this is also problematic. For explanation, note
that $\textnormal{{an}}(j)\setminus\\{k\\}$ can be partitioned into mediators
$\textnormal{{me}}(k,j)$ and non-mediators

$\textnormal{{nm}}(k,j)=\textnormal{{an}}(j)\setminus(\textnormal{{me}}(k,j)\cup\\{k\\}).$

In Figure 1 (a), $\bm{X}_{\textnormal{{iv}}(k)}$ can be associated with $\eta$
given
$\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}}=(\bm{Y}_{\textnormal{{me}}(k,j)},\bm{Y}_{\textnormal{{nm}}(k,j)})$,
violating the unconfoundedness of IVs (Remark 1) and causing an estimation
bias. This is because the mediators $\bm{Y}_{\textnormal{{me}}(k,j)}$ generate
additional associations after conditioning on them; see the Appendix for
technical discussion using the concept of d-separation (Pearl,, 2009).

Now, we propose a new method, which eliminates the impact of mediators
$\bm{Y}_{\textnormal{{me}}(k,j)}$ by introducing the working response
$\overline{Y}_{j}=Y_{j}-{\mathbf{U}}_{\textnormal{{me}}(k,j),j}^{\top}\bm{Y}_{\textnormal{{me}}(k,j)}$,
as depicted in Figure 1 (b). Of note, the definition of $\overline{Y}_{j}$
depends on $(k,j)$, which is dropped for simplicity. As in Angrist et al.,
(1996), we have

$\begin{split}&\operatorname{\mathbb{E}}\Big{(}\overline{Y}_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}\\\
\overset{\text{(i)}}{=}\
&\mathrm{U}_{kj}\operatorname{\mathbb{E}}\Big{(}Y_{k}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}+\sum_{k^{\prime}\in\textnormal{{nm}}(k,j)}\mathrm{U}_{k^{\prime}j}Y_{k^{\prime}}+\sum_{l\notin\textnormal{{iv}}(k)}\mathrm{W}_{lj}X_{l}+\operatorname{\mathbb{E}}\Big{(}\varepsilon_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}\\\
\overset{\text{(ii)}}{=}\
&\mathrm{U}_{kj}\widetilde{Y}_{k}+\bm{\gamma}^{\top}\bm{Z},\end{split}$ (7)

where
$\widetilde{Y}_{k}=\mathbb{E}(Y_{k}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$,
$\bm{Z}=(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{ca}}(k)^{c}})=(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{iv}}(k)^{c}})$,
equality (i) follows from (1), and equality (ii) holds because
$\mathbb{E}\left(\varepsilon_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\right)$
is a linear combination of
$(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{iv}}(k)^{c}})$ by Lemma
1 in Appendix. Observe that $\widetilde{Y}_{k}$ depends on
$\bm{X}_{\textnormal{{iv}}(k)}$ while $\bm{Z}$ does not. As a result, the
$\mathrm{U}_{kj}$ is identified through the working response regression.

This approach requires the knowledge of
$\mathbf{U}_{\textnormal{{me}}(k,j),j}$ prior to identifying
$\mathrm{U}_{kj}$. Given $\mathcal{G}^{+}$, we develop a sequential procedure
to learn $\bm{\mathbf{U}}$. First, we identify $\mathrm{U}_{kj}$ for each pair
$(k,j)$ such that the longest path in $\mathcal{G}^{+}$ between $k$ and $j$ is
equal to $d=1$. Then for $(k,j)$ such that the longest path in
$\mathcal{G}^{+}$ between $k$ and $j$ is $d=2$, the effects of mediators
$\mathbf{U}_{\textnormal{{me}}(k,j),j}$ are available. Thus, we can identify
$\mathrm{U}_{kj}$ in (7). Proceed similarly for $d=3,4,5,\ldots$ until all
pairs in $\mathcal{E}^{+}$ have been identified.

#### The case with invalid IVs.

In general, $\textnormal{{ca}}(k)\supseteq\textnormal{{iv}}(k)$ because of
invalid IVs, where $\textnormal{{ca}}(k)$ is known but $\textnormal{{iv}}(k)$
is unknown. Similar to Kang et al., (2016), we have

$\begin{split}&\operatorname{\mathbb{E}}\Big{(}\overline{Y}_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}\\\
{=}\
&\mathrm{U}_{kj}\operatorname{\mathbb{E}}\Big{(}Y_{k}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}+\sum_{k^{\prime}\in\textnormal{{nm}}(k,j)}\mathrm{U}_{k^{\prime}j}Y_{k^{\prime}}+\sum_{l\notin\textnormal{{iv}}(k)}\mathrm{W}_{lj}X_{l}+\operatorname{\mathbb{E}}\Big{(}\varepsilon_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}\Big{)}\\\
\overset{\text{(iii)}}{=}\
&\mathrm{U}_{kj}\widetilde{Y}_{k}+\bm{\gamma}^{\top}\bm{Z}+\sum_{l\in\textnormal{{ca}}(k)\setminus\textnormal{{iv}}(k)}\beta_{l}X_{l},\end{split}$
(8)

where
$\widetilde{Y}_{k}=\mathbb{E}(Y_{k}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$,
$\bm{Z}=(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{ca}}(k)^{c}})$,
equality (iii) holds by Lemma 1 in Appendix, and
$\beta_{l}=\mathrm{W}_{lj}\neq 0$ indicates $X_{l}$ is an invalid IV for
$Y_{k}$. However, since $\textnormal{{iv}}(k)$ has not been identified and
$\widetilde{Y}_{k}$ depends on $\bm{X}_{\textnormal{{ca}}(k)}$, the
representation of (iii) may not be unique. When the majority rule (A3) is
satisfied by the DAG, the term (iii) admits the unique expression as in (8),
providing the identification of $\mathrm{U}_{kj}$. This leads to a sparse
regression for an infinite sample

$\min_{\mathrm{U}_{kj},\bm{\beta},\bm{\gamma}}\
\mathbb{E}\left(\overline{Y}_{j}-\mathrm{U}_{kj}\widetilde{Y}_{k}-\bm{\gamma}^{\top}\bm{Z}-\bm{\beta}^{\top}\bm{X}_{\textnormal{{ca}}(k)}\right)^{2}\quad\text{s.t.}\quad\|\bm{\beta}\|_{0}\leq\kappa,$
(9)

where $0\leq\kappa<|\textnormal{{ca}}(k)|/2$ is an integer-valued
hyperparameter controlling the sparsity of $\bm{\beta}$.

###  3.4 Finite-sample estimation of U

Suppose $(\mathbf{Y}_{p\times n},\mathbf{X}_{q\times n})$ are given. To
estimate $\mathrm{U}_{kj}$, noting that $\widetilde{Y}_{k}$ is linear in
$(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$ by Lemma 1, we estimate
$\widetilde{\mathrm{Y}}_{k,i}$ by
$\widehat{\mathrm{Y}}_{k,i}=\widehat{\bm{\alpha}}_{1}^{\top}\bm{\mathbf{X}}_{+,i}+\widehat{\bm{\alpha}}_{2}^{\top}\bm{\mathbf{Y}}_{\textnormal{{nm}}(k),i}$,
where $(\widehat{\bm{\alpha}}_{1},\widehat{\bm{\alpha}}_{2})$ solves

$\begin{split}\min_{\bm{\mathbf{\alpha}}_{1},\bm{\mathbf{\alpha}}_{2}}\
\sum_{i=1}^{n}\Big{(}\mathrm{Y}_{k,i}-\bm{\alpha}_{1}^{\top}\bm{\mathbf{X}}_{+,i}+\bm{\alpha}_{2}^{\top}\bm{\mathbf{Y}}_{\textnormal{{nm}}(k),i}\Big{)}^{2}\quad\text{s.t.}\quad\|\bm{\mathbf{\alpha}}_{1}\|_{0}+\|\bm{\mathbf{\alpha}}_{2}\|_{0}\leq\nu_{1},\end{split}$
(10)

with $\nu_{1}$ being a tuning parameter. Let the final estimate
$\widehat{\mathrm{U}}_{kj}$ with
$(\widehat{\bm{\beta}},\widehat{\bm{\gamma}})$ be the solution to the working
response regression (provided that
$\widehat{\bm{\mathbf{U}}}_{\textnormal{{me}}(k,j),j}$ are available)

$\begin{split}\min_{\mathrm{U}_{kj},\bm{\beta},\bm{\gamma}}\
\sum_{i=1}^{n}\Big{(}\Big{(}\mathrm{Y}_{j,i}-\widehat{\bm{\mathbf{U}}}^{\top}_{\textnormal{{me}}(k,j),j}\bm{\mathbf{Y}}_{\textnormal{{me}}(k,j),i}\Big{)}-\mathrm{U}_{kj}\widehat{\mathrm{Y}}_{k,i}-\bm{\beta}^{\top}\bm{X}_{\textnormal{{ca}}(k),i}-\bm{\gamma}^{\top}\bm{Z}_{i}\Big{)}^{2}\\\
\text{s.t.}\quad\rho(\mathrm{U}_{kj})+\|\bm{\beta}\|_{0}\leq\kappa,\quad\|\bm{\gamma}\|_{0}\leq\nu_{2},\end{split}$
(11)

where $0\leq\kappa\leq|\textnormal{{ca}}(k)|/2$ and
$0\leq\nu_{2}\leq|\textnormal{{nm}}(k,j)|+|\textnormal{{ca}}(k)^{c}|$ are
tuning parameters. Depending on the purpose,
$\rho(\cdot)=\operatorname{I}(\cdot\neq 0)$ for graph recovery and
$\rho(\cdot)=0$ for effect estimation without selection. In (10)–(11),
$\nu_{1},\nu_{2}$ are added to treat possible high-dimensional situations and
the hyperparameters are tuned by BIC. Algorithm 2 summarizes the procedure.

Input: Data $\mathbf{Y}_{p\times n}$ and $\mathbf{X}_{q\times n}$, ARG
$\mathcal{G}^{+}$ and candidate IV sets
$\\{\textnormal{{ca}}_{\mathcal{G}}(k)\\}_{1\leq k\leq p}$;

1

2Initialize $\widehat{\mathbf{U}}\leftarrow\bm{0}$ and $d\leftarrow 1$;

3 while  _$d\leq$ the length of the longest directed path in
$\mathcal{G}^{+}$_ do

4 For $(k,j)\in\mathcal{E}^{+}$ so that the length of the longest directed
path from $Y_{k}$ to $Y_{j}$ is $d$, estimate $\widehat{\mathrm{U}}_{kj}$ with
(10)–(11);

5 Update $d\leftarrow d+1$;

6

7 end while

8return _$\widehat{\mathbf{U}}$_ ;

Algorithm 2 Estimation of $\mathbf{U}$

##  4 Likelihood inference

This section develops a likelihood ratio test for the presence of multiple
directed edges. Let $\mathcal{H}\subseteq\\{(k,j):k\neq j,\ 1\leq k,j\leq
p\\}$ be a hypothesized edge set for primary variables $\bm{Y}$, where
$(k,j)\in\mathcal{H}$ specifies a (hypothesized) directed edge
$Y_{k}\rightarrow Y_{j}$ in (1). Now consider simultaneous testing of directed
edges,

$H_{0}:\mathrm{U}_{kj}=0\text{ for all }(k,j)\in\mathcal{H}\quad\text{ versus
}\quad H_{a}:\mathrm{U}_{kj}\neq 0\text{ for some }(k,j)\in\mathcal{H}.$ (12)

The null hypothesis $H_{0}$ asserts that all hypothesized edges in
$\mathcal{H}$ are absent in the true graph $\mathcal{G}$. Rejecting $H_{0}$
indicates that at least one hypothesized edge in $\mathcal{H}$ presents in
$\mathcal{G}$.

#### The likelihood ratio.

Given $\mathcal{G}^{+}=(\bm{X},\bm{Y};\mathcal{E}^{+},\mathcal{I}^{+})$, let
$\bm{\theta}(\mathcal{G}^{+})=(\mathbf{U},\mathbf{W})$ encode the coefficient
parameters in $\mathcal{G}^{+}$, where
$\mathbf{U}=(\mathbf{U}_{\mathcal{E}^{+}},\bm{0})$ and
$\mathbf{W}=(\mathbf{W}_{\mathcal{I}^{+}},\bm{0})$. As such, the adjacency
matrix $\mathbf{U}$ automatically meets the acyclicity constraint. Given a
random sample ${(\mathbf{Y}_{+,i},\mathbf{X}_{+,i})}_{i=1}^{n}$, the log-
likelihood is written as (up to an additive constant)

$L(\bm{\theta}(\mathcal{G}^{+}),\bm{\Omega})=-\frac{1}{2}\sum_{i=1}^{n}\left\|\bm{\Omega}^{1/2}\left(\left(\mathbf{I}-\mathbf{U}^{\top}\right)\mathbf{Y}_{+,i}-\mathbf{W}^{\top}\mathbf{X}_{+,i}\right)\right\|_{2}^{2}+\frac{n}{2}\log\det(\bm{\Omega}),$
(13)

where $\bm{\Omega}=\bm{\Sigma}^{-1}$ is the inverse of $\bm{\Sigma}$ in (1).
Then the maximum likelihood estimation (MLE) of (1) can be written as

$\max_{(\mathcal{G}^{+},\bm{\Omega})}\max_{\bm{\theta}(\mathcal{G}^{+})}L(\bm{\theta}(\mathcal{G}^{+}),\bm{\Omega}).$
(14)

In view of (14), to obtain a likelihood ratio statistic for (12) we need to
compute the following quantities: (1) a consistent estimate
$\widehat{\mathcal{G}}^{+}$ of $\mathcal{G}^{+}$, (2) a consistent estimate
$\widehat{\bm{\Omega}}$ of $\bm{\Omega}$, and (3) two estimates,
$\widehat{\bm{\theta}}^{(0)}$ and $\widehat{\bm{\theta}}^{(1)}$, of
$\bm{\theta}(\mathcal{G}^{+})$ under $H_{0}$ and $H_{a}$, respectively. This
leads to the likelihood ratio defined as

$L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Omega}})-L(\widehat{\bm{\theta}}^{(0)},\widehat{\bm{\Omega}}),$
(15)

where $\mathcal{G}^{+}$ is estimated by Algorithm 1 and ${\bm{\Omega}}$ is
estimated from the residuals after fitting model (1) via Algorithm 2.

#### Inference subject to acyclicity.

In classical models, a likelihood ratio of form (15) has a nondegenerate and
tractable limiting distribution, typically a chi-squared distribution with
degrees of freedom $|\mathcal{H}|$. However, the likelihood ratio for (12) may
behave differently from classical ones since (15) may be degenerate or
intractable, as to be explained.

First, note that the maximum likelihood subject to a wrong ARG
$\widetilde{\mathcal{G}}^{+}\not\supseteq\mathcal{G}$ tends to be smaller than
that subject to the correct $\mathcal{G}^{+}$, that is,

$\max_{\widetilde{\mathcal{G}}^{+}\not\supseteq\mathcal{G}}\max_{\bm{\theta}(\widetilde{\mathcal{G}}^{+}),\bm{\Omega}}L(\bm{\theta}(\widetilde{\mathcal{G}}^{+}),\bm{\Omega})<\max_{\bm{\theta}({\mathcal{G}}^{+}),\bm{\Omega}}L(\bm{\theta}(\mathcal{G}^{+}),\bm{\Omega}),$

as $n\to\infty$ under some regularity conditions for consistency. Thus, we
assume $\widehat{\mathcal{G}}^{+}=\mathcal{G}^{+}$ in this paragraph. Then
$\widehat{\bm{\theta}}^{(0)}$ is the MLE subject to ${\mathcal{G}}^{+}$ and
$\mathbf{U}_{\mathcal{H}}=\bm{0}$, which is equal to the MLE subject to the
graph
$\mathcal{G}^{+}_{0}=(\bm{X},\bm{Y};\mathcal{E}^{+}\setminus\mathcal{H},\mathcal{I}^{+})$.
Meanwhile, to test whether any edge in $\mathcal{H}$ exists,
$\widehat{\bm{\theta}}^{(1)}$ is the MLE subject to an augmented graph
$\mathcal{G}^{+}_{1}=(\bm{X},\bm{Y};\mathcal{E}^{+}\cup\mathcal{H},\mathcal{I}^{+})$
with hypothesized edges being added, namely,
$\widehat{\mathbf{U}}^{(1)}=(\widehat{\mathbf{U}}^{(1)}_{\mathcal{E}^{+}\cup\mathcal{H}},\bm{0})$
and
$\widehat{\mathbf{W}}^{(1)}=(\widehat{\mathbf{W}}^{(1)}_{\mathcal{I}^{+}},\bm{0})$.
Of note, since $\mathcal{H}$ is pre-specified by the user,
$\mathcal{G}^{+}_{1}$ is not necessarily acyclic, and thus, not all edges in
$\mathcal{H}$ could present in $\widehat{\mathbf{U}}^{(1)}$. Furthermore, if a
hypothesized edge $(k,j)$ is present in $\widehat{\mathbf{U}}^{(1)}$, then
$\\{(k,j)\\}\cup\mathcal{E}^{+}$ must have no directed cycle and (15) is
strictly positive (nondegenerate). However, even if (15) does not degenerate
to zero, its limiting distribution can be complicated when there exist
multiple ways of augmenting $\mathcal{G}^{+}$ with the edges in $\mathcal{H}$
while maintaining the resulting graph as a DAG. Therefore, a regularity
condition for $\mathcal{H}$ is necessary to rule out intractable situations.

On the ground of the foregoing discussion, we introduce the concepts of
nondegeneracy and regularity to characterize the behavior of (15) as in Li et
al., 2023a .

######  Definition 5 (Nondegeneracy and regularity with respect to
$\mathcal{G}^{+}$).

  1. (A)

An edge $(k,j)\in\mathcal{H}$ is said to be nondegenerate with respect to an
ancestral graph
$\mathcal{G}^{+}=(\bm{Y},\bm{X};\mathcal{E}^{+},\mathcal{I}^{+})$ if
$\\{(k,j)\\}\cup\mathcal{E}^{+}$ contains no directed cycle. Otherwise,
$(k,j)$ is said to be degenerate. Let $\mathcal{D}\subseteq\mathcal{H}$ be the
set of all nondegenerate edges with respect to $\mathcal{G}^{+}$. A null
hypothesis $H_{0}$ is said to be nondegenerate with respect to
$\mathcal{G}^{+}$ if $\mathcal{D}\neq\emptyset$. Otherwise, $H_{0}$ is said to
be degenerate.

  2. (B)

A null hypothesis $H_{0}$ is said to be regular with respect to
$\mathcal{G}^{+}$ if $\mathcal{D}\cup\mathcal{E}^{+}$ contains no directed
cycle. Otherwise, $H_{0}$ is called irregular.

Suppose $H_{0}$ is nondegenerate and regular. Then
$\widehat{\bm{\theta}}^{(0)}$ is the MLE subject to the graph
$\mathcal{G}^{+}_{0}=(\bm{X},\bm{Y};\mathcal{E}^{+}\setminus\mathcal{D},\mathcal{I}^{+})$
and $\widehat{\bm{\theta}}^{(1)}$ is the MLE subject to the graph
$\mathcal{G}^{+}_{1}=(\bm{X},\bm{Y};\mathcal{E}^{+}\cup\mathcal{D},\mathcal{I}^{+})$.

Now, we investigate the limiting distribution of (15) and derive an asymptotic
test based on it. To this end, define the statistic

$T(\mathcal{D})=\begin{cases}2\left(L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Omega}})-L(\widehat{\bm{\theta}}^{(0)},\widehat{\bm{\Omega}})\right)&\text{if
}|\mathcal{D}|\text{ is fixed},\\\
\left(2\left(L(\widehat{\bm{\theta}}^{(1)},\widehat{\bm{\Omega}})-L(\widehat{\bm{\theta}}^{(0)},\widehat{\bm{\Omega}})\right)-|\mathcal{D}|\right)/\sqrt{2|\mathcal{D}|}&\text{if
}|\mathcal{D}|\to\infty.\end{cases}$ (16)

######  Theorem 2 (Limiting distribution).

Assume the null hypothesis $H_{0}$ is nondegenerate and regular. Suppose
$\mathbb{P}(\widehat{\mathcal{G}}^{+}={\mathcal{G}}^{+})\to 1$ as
$n\to\infty$. Then we have $\mathbb{P}(\widehat{\mathcal{D}}=\mathcal{D})\to
1$. In addition, if
$\|\widehat{\bm{\Omega}}-\bm{\Omega}\|_{2}^{2}=O_{\mathbb{P}}({|S|\log(p\vee
n)/n})$ where $S=\\{(k,j):\Omega_{kj}\neq 0\\}$, then under $H_{0}$,

$\begin{split}T(\widehat{\mathcal{D}})\overset{d}{\longrightarrow}\begin{cases}\chi^{2}_{|\mathcal{D}|},&\text{if
}|\mathcal{D}|\text{ is fixed and }{|S|\log(p\vee n)}/{n}\rightarrow 0,\\\
N(0,1),&\text{if }|\mathcal{D}|\to\infty\text{ and
}{|{\mathcal{D}}|}|S|\log(p\vee n)/{n}\rightarrow 0.\end{cases}\end{split}$

On the basis of Theorem 2, we conduct inference by substituting
$|\mathcal{D}|$ by its estimate $|\widehat{\mathcal{D}}|$ and proceed with the
empirical rule: (1) use the chi-squared test when
$|\widehat{\mathcal{D}}|<50$, and (2) use the normal test when
$|\widehat{\mathcal{D}}|\geq 50$.

Theorem 2 requires a good estimator $\widehat{\bm{\Omega}}$ of
$\bm{\Omega}=\bm{\Sigma}^{-1}$ to account for the confounding effects, where
$\bm{\Sigma}=\operatorname{Cov}(\bm{\varepsilon})$. To estimate $\bm{\Omega}$,
let
$\widehat{\bm{\varepsilon}}_{+,i}=(\mathbf{I}-\widehat{\mathbf{U}})^{\top}\mathbf{Y}_{+,i}-\widehat{\mathbf{W}}^{\top}\mathbf{X}_{+,i}$;
$i=1,\ldots,n$ be the estimated residuals after fitting (1) with Algorithm 2.
Here we use the neighborhood selection method (Meinshausen and Bühlmann,,
2006) with an additional refitting to obtain a positive definite estimate
$\widehat{\bm{\Omega}}$. In Supplementary Materials, we include the
computational details and show that this estimator satisfies
$\|\widehat{\bm{\Omega}}-\bm{\Omega}\|_{F}^{2}=O_{\mathbb{P}}(|S|\log(p\vee
n)/n)$ so that Theorem 2 applies.

######  Remark 2.

In Theorem 2, we focus on nondegenerate and regular hypotheses. For a
degenerate case, we define the p-value as one. For an irregular case where
$\mathcal{D}\cup\mathcal{E}^{+}$ contains a directed cycle, we decompose
$H_{0}$ into sub-hypotheses $H_{0}^{(1)},\ldots,H_{0}^{(r)}$, each of which is
regular. Then testing $H_{0}$ is reduced to multiple testing for
$H_{0}^{(1)},\ldots,H_{0}^{(r)}$.

Finally, we discuss two aspects of likelihood estimation and inference in the
presence of unmeasured confounding. First, when $\bm{\Sigma}$ is non-diagonal,
the likelihood in (13) cannot be factorized according to $\mathcal{G}$ (or
$\mathcal{G}^{+}$). This implies that, unlike the case without latent
confounders (Shojaie and Michailidis,, 2010), the parameters of each equation
in (1) cannot be estimated separately given $\mathcal{G}^{+}$. Indeed, the
likelihood estimation of $(\mathbf{U},\mathbf{W})$ in (1) requires a
preliminary estimate of $\bm{\Omega}$ to account for correlations arising from
hidden confounding. Furthermore, compared to Li et al., 2023a , the likelihood
ratio (15) is no longer a sum of likelihood ratios of equations associated
with nondegenerate hypothesized edges, rendering inference more challenging in
both computation and theory when hidden confounders are present.
Computationally, the likelihood ratio (15) requires maximization of the full
likelihood, which is costly for a large-scale graph. Theoretically, estimating
$\bm{\Omega}$ and $(\mathbf{U},\mathbf{W})$ in high-dimensional situations may
suffer from the curse of dimensionality.

Second, to mitigate the challenges in inference, we may conduct inference with
respect to a sub-DAG to achieve dimensionality reduction. Specifically, let
$\mathcal{D}$ be the nondegenerate edges of $H_{0}$. Given ARG
$\mathcal{G}^{+}$, we perform likelihood inference using a sub-DAG (of ARG)
${\mathcal{G}}^{+}_{\text{sub}}=({\bm{X}}_{\text{sub}},{\bm{Y}}_{\text{sub}};{\mathcal{E}}_{\text{sub}}^{+},{\mathcal{I}}_{\text{sub}}^{+})$,
where all edges specified in $\mathcal{D}$ are among primary variables
${\bm{Y}}_{\text{sub}}$, and ${\bm{Y}}_{\text{sub}}$ are non-descendants of
$\bm{Y}\setminus{\bm{Y}}_{\text{sub}}$ in the graph
$(\bm{X},\bm{Y};\mathcal{E}^{+}\cup\mathcal{D},\mathcal{I}^{+})$,
${\bm{X}}_{\text{sub}}$ is the set of intervention variables of
${\bm{Y}}_{\text{sub}}$, ${\mathcal{E}}_{\text{sub}}^{+}$ is the set of
ancestral relations among ${\bm{Y}}_{\text{sub}}$, and
${\mathcal{I}}_{\text{sub}}^{+}$ is the set of interventional relations
between ${\bm{X}}_{\text{sub}}$ and ${\bm{Y}}_{\text{sub}}$ in ARG
$\mathcal{G}^{+}$. Then the test statistic (16) is computed within the sub-DAG
${\mathcal{G}}_{\text{sub}}^{+}$, which reduces computation. Furthermore,
Theorem 2 holds true when the estimator of the smaller precision matrix
${\bm{\Omega}}_{\text{sub}}$ enjoys the desired convergence rate
$O_{\mathbb{P}}(\sqrt{|{S}_{\text{sub}}|\log({p}_{\text{sub}}\vee n)/n})$ in
operator norm, where the subscript ${}_{\text{sub}}$ denotes the quantities
corresponding to the structural equations of ${\bm{Y}}_{\text{sub}}$.

##  5 Theory

In this section, we develop a theory to quantify the finite sample performance
as well as the complexities of Algorithms 1–2 when TLP is used for
computation.

To proceed, we introduce some technical conditions for casual discovery
consistency. For $(k,j)\in\mathcal{E}^{+}$, let
$\widetilde{\bm{\Sigma}}^{(k,j)}$ be the covariance matrix of
$(\operatorname{\mathbb{E}}(Y_{k}\mid\bm{Y}_{\textnormal{{nm}}_{\mathcal{G}}(k,j)},\bm{X}),\bm{Y}_{\textnormal{{nm}}_{\mathcal{G}}(k,j)},\bm{X})$.
Moreover, let
$s=\max_{(k,j)\in\mathcal{E}^{+}}(\kappa+\nu_{2},\nu_{1})\vee\max_{1\leq k\leq
p}\|\mathbf{V}_{+,k}\|_{0}$ be the maximum sparsity-level in the estimation
procedure, where $\nu_{1},\nu_{2},\kappa$ depends on $(k,j)$ which is dropped
for conciseness. Assume there exist constants $c_{0},c_{1},c_{2},c_{3}>0$ such
that

  1. (C1)

$\min_{(k,j)\in\mathcal{E}^{+}}\min_{B:|B|\leq
2s}\min_{\bm{\mathbf{v}}:\|\bm{\mathbf{v}}\|_{2}=1,\|\bm{\mathbf{v}}_{B^{c}}\|_{1}\leq
3\|\bm{\mathbf{v}}_{B}\|_{1}+c_{0}s\sqrt{\log(p)/n}}\langle\bm{\mathbf{v}},\widetilde{\bm{\Sigma}}^{(k,j)}\bm{\mathbf{v}}\rangle\geq
c_{1}$.

  2. (C2)

$\min_{\mathrm{V}_{kj}\neq 0}|\mathrm{V}_{kj}|\geq c_{2}\sqrt{{\log(q\vee
n)}/{n}}$.

  3. (C3)

$\min_{\mathrm{U}_{kj}\neq 0}|\mathrm{U}_{kj}|\geq c_{3}\sqrt{\log(p\vee
n)/n}$.

  4. (C4)

$\max_{1\leq k\leq
p}\\{|\textnormal{{an}}_{\mathcal{G}}(k)|,|\textsc{in}_{\mathcal{G}}(k)|,\|\mathbf{U}_{+,k}\|_{1}\\}=O(1)$,
and
$\max_{(k,j)\in\mathcal{E}^{+}}(\operatorname{Diag}(\widetilde{\bm{\Sigma}}^{(k,j)}))=O(1)$.

Condition (C1) is a restricted eigenvalue condition, which is common in high-
dimensional estimation (Bickel et al.,, 2009) and can be viewed as a stronger
version of (A1) in Theorem 1. (C2) and (C3) impose restrictions on the minimal
signal strengths of $\mathbf{V}$ and $\mathbf{U}$ so that the ARG
$\mathcal{G}^{+}$ and DAG $\mathcal{G}$ can be consistently recovered,
respectively. They are similar to the beta-min condition (Meinshausen and
Bühlmann,, 2006) and the degree of separation condition (Shen et al.,, 2012)
in the variable selection literature.

######  Theorem 3.

Suppose Assumptions (A1)–(A3) in Theorem 1 are satisfied and assume $\bm{X}$
is sub-Gaussian with mean zero and parameter $\varsigma^{2}$.

  1. (A)

(Parameter estimation) Suppose (C1), (C2), (C4) are met with sufficiently
large $c_{0},c_{1},c_{2}$. Suppose the tuning parameters are suitably chosen
such that

    1. (i)

In Algorithm 1, $0.01c_{2}\sqrt{\log(q\vee n)/n}\leq\tau^{\prime}\leq
0.4\min_{\mathrm{V}_{kj}\neq 0}|\mathrm{V}_{kj}|$,
$\kappa^{\prime}_{j}=\|\bm{\mathbf{V}}_{+,j}\|_{0}$ for $1\leq j\leq p$.

    2. (ii)

In Algorithm 2, $0.5c_{3}\sqrt{{\log(p\vee n)}/{n}}\leq\tau$,
$\nu_{1}=\lceil\operatorname{TLP}_{\tau}((\bm{\mathbf{\alpha}}_{1},\bm{\mathbf{\alpha}}_{2}))\rceil$,
$\nu_{2}=\lceil\operatorname{TLP}_{\tau}(\bm{\gamma})\rceil$, and
$\kappa=\lceil\operatorname{TLP}_{\tau}(\bm{\beta})\rceil$ for any
$(k,j)\in\mathcal{E}^{+}$.

Then there exists constant $C_{1}>0$ such that when $n$ is sufficiently large

$|\widehat{\mathrm{U}}_{kj}-{\mathrm{U}_{kj}}|\leq C_{1}\sqrt{{\log(p\vee
n)}/{n}},$

almost surely under $\mathbb{P}_{(\mathbf{U},\mathbf{W},\bm{\Sigma})}$.
Moreover, Algorithms 1 and 2 respectively terminate in
$O(p\times\log(s)\times(q^{3}+nq^{2}))$ and
$O(|\mathcal{E}^{+}|\times\log(s)\times(q^{3}+nq^{2}))$ operations almost
surely.

  2. (B)

(Graph recovery) Additionally, if (C3) is satisfied with $c_{3}>C_{1}>\tau$,
then when $n$ is sufficiently large we have
$\widehat{\mathcal{G}}=\mathcal{G}$ almost surely.

By Theorem 3, the proposed method achieves causal discovery consistency in
terms of consistent parameter estimation and structure recovery. Moreover,
Algorithms 1–2 enjoy low-order polynomial time complexity almost surely
provided that the data are randomly sampled from (1).

##  6 Numerical examples

###  6.1 Simulations

This subsection investigates via simulations the operating characteristics of
GrIVET, including the qualities of structure learning, parameter estimation,
and statistical inference.

To generate an observation $(\bm{Y},\bm{X})$, we first introduce hidden
variables $\bm{\eta}\sim N(\bm{0},\mathbf{I}_{r\times r})$ as unmeasured
confounders. Then, we sample $\bm{X}$ from $N(\bm{0},\mathbf{I}_{q\times q})$
for continuous interventions or from $\\{-1,1\\}^{q}$ with equal probability
for discrete interventions. Given $\bm{X}$ and $\bm{\eta}$, we generate
$\bm{Y}$ according to

$\bm{Y}=\mathbf{U}^{\top}\bm{Y}+\mathbf{W}^{\top}\bm{X}+\bm{\Phi}^{\top}\bm{\eta}+\bm{e},\quad\bm{e}\sim
N\left(\mathbf{0},\operatorname{Diag}\left(\sigma_{1}^{2},\ldots,\sigma_{p}^{2}\right)\right).$
(17)

We conduct simulations with the following settings.

  * •

Hub graph. Let $p=101$, $q=252$, and $r=10$. For $\mathbf{U}$,
$(\mathrm{U}_{1,j})_{2\leq j\leq p}$ are independently sampled from
$\\{-1,1\\}$ with equal probability, while the rest are set to $0$. This
generates a sparse graph with the dense neighborhood of the first node. Let
$\mathbf{W}_{q\times p}=(\mathbf{I}_{p\times p},\mathbf{I}_{p\times
p},\mathbf{F}^{\top})^{\top}$ where the entries
$(\mathrm{F}_{j,2j},\mathrm{F}_{j,2j+1})_{1\leq j\leq q-2p}$ are set to $1$,
while other entries of $\mathbf{F}$ are zero. Then $X_{j},X_{2j}$ are IVs of
$Y_{j}$ for $j=1,\ldots,p$ and $X_{2p+1},\ldots,X_{q}$ are invalid IVs with
two intervention targets. For the confounders, $\mathrm{\Phi}_{1,1}$ and
$(\mathrm{\Phi}_{jk})^{1\leq j\leq r}_{10j-8\leq k\leq 10j+1}$ are sampled
uniformly from $(-0.4,-0.6)\cup(0.4,0.6)$, while other entries of
$\mathbf{\Phi}$ are zero. We generate $(\sigma_{1},\ldots,\sigma_{p})$
uniformly from $(0.4,0.6)$.

  * •

Random graph. Let $p=100$, $q=250$, and $r=10$. For $\mathbf{U}$, the upper
off-diagonals $(\mathrm{U}_{kj})_{k<j}$ are sampled independently from
$\\{0,1\\}$ according to $\text{Bernoulli}({1}/{10p})$ while other entries are
zero. Set $\mathbf{W}_{q\times p}=(\mathbf{I}_{p\times p},\mathbf{I}_{p\times
p},\mathbf{F}^{\top})^{\top}$ where
$(\mathrm{F}_{j,2j-1},\mathrm{F}_{j,2j})_{1\leq j\leq 1-2p}$ are set to $1$,
while other entries of $\mathbf{F}$ are zero. Then $X_{j},X_{2j}$ are IVs of
$Y_{j}$ for $j=1,\ldots,p$ and $X_{2p+1},\ldots,X_{q}$ are invalid IVs with
two intervention targets. For the confounders, $(\mathrm{\Phi}_{jk})^{1\leq
j\leq r}_{10j-9\leq k\leq 10j}$ are sampled uniformly from
$(-0.4,-0.6)\cup(0.4,0.6)$, while other entries of $\mathbf{\Phi}$ are zero.
We generate $(\sigma_{1},\ldots,\sigma_{p})$ uniformly from $(0.4,0.6)$.

#### Structure learning.

After obtaining ancestral relations from Algorithm 1, we implement Algorithm 2
to confirm parental relations but with constraints also imposed on the
parameter of interest. Four graph metrics are used for evaluation: the false
discovery rate (FDR), the true positive rate (TPR), the Jaccard index (JI),
and the structural Hamming distance (SHD). The results in Table 1 demonstrate
the strong performance of GrIVET in structure learning. Note that a high TPR
indicates GrIVET’s capability to detect the true existing edges, while the FDR
remains low, signifying the high specificity of GrIVET. In Supplementary
Materials Section 3.3, we further compare GrIVET with RFCI (Colombo et al.,,
2012) and LRpS-GES (Frot et al.,, 2019) in terms of structural learning
accuracy. GrIVET compares favorably against the competitors.

Table 1: False discovery rate (FDR), true positive rate (TPR), structural Hamming distance (SHD), and Jaccard index (JI) of GrIVET for causal discovery over 1000 simulation replications. To compute the metrics, let TP, RE, FP, and FN be the numbers of identified edges with correct directions, those with wrong directions, estimated edges not in the skeleton of the true graph, and missing edges compared to the true skeleton. Then $\mathrm{FDR}=(\mathrm{RE}+\mathrm{FP})/(\mathrm{TP}+\mathrm{RE}+\mathrm{FP})$, $\mathrm{TPR}=\mathrm{TP}/(\mathrm{TP}+\mathrm{FN})$, $\mathrm{SHD}=\mathrm{FP}+\mathrm{FN}+\mathrm{RE}$, and $\mathrm{JI}=\mathrm{TP}/(\mathrm{TP}+\mathrm{SHD})$. Graph | Intervention | $n$ |  FDR($\%$) |  TPR($\%$) | SHD |  JI($\%$)
---|---|---|---|---|---|---
Hub | Continuous | 500 | 0.000 | 100.000 | 0.000 | 100.000
|  | 400 | 0.000 | 99.998 | 0.002 | 99.998
|  | 300 | 0.000 | 99.998 | 0.002 | 99.998
| Discrete | 500 | 0.000 | 99.999 | 0.001 | 99.999
|  | 400 | 0.000 | 99.998 | 0.002 | 99.998
|  | 300 | 0.000 | 99.999 | 0.001 | 99.999
Random | Continuous | 500 | 0.011 | 98.600 | 0.001 | 98.589
|  | 400 | 0.000 | 98.600 | 0.000 | 98.600
|  | 300 | 0.018 | 98.590 | 0.003 | 98.575
| Discrete | 500 | 0.000 | 98.600 | 0.000 | 98.600
|  | 400 | 0.024 | 98.600 | 0.002 | 98.576
|  | 300 | 0.000 | 98.600 | 0.000 | 98.600

#### Parameter estimation.

Table 2: Parameter estimation: the average of largest absolute difference (Max AD), the average absolute differences (Mean AD), and the average squared differences (Mean SqD) between the estimated parameters and the true parameters for two competing methods over 1000 simulation replications. Graph | Intervention | $n$ | GrIVET |  |  Direct regression (Li et al., 2023a, )
---|---|---|---|---|---
|  |  | (Max AD, Mean AD, Mean SqD) |  | (Max AD, Mean AD, Mean SqD)
Hub | Continuous | 500 | (0.06107, 0.01808, 0.00052) |  | (0.12817, 0.02448, 0.00142)
|  | 400 | (0.06863, 0.02037, 0.00066) |  | (0.13196, 0.02637, 0.00156)
|  | 300 | (0.07922, 0.02347, 0.00087) |  | (0.13395, 0.02873, 0.00170)
| Discrete | 500 | (0.06119, 0.01803, 0.00051) |  | (0.12770, 0.02434, 0.00141)
|  | 400 | (0.06932, 0.02030, 0.00065) |  | (0.13041, 0.02621, 0.00153)
|  | 300 | (0.08046, 0.02355, 0.00088) |  | (0.13334, 0.02867, 0.00169)
Random | Continuous | 500 | (0.02836, 0.01445, 0.00034) |  | (0.04254, 0.01791, 0.00076)
|  | 400 | (0.03245, 0.01660, 0.00045) |  | (0.04390, 0.01899, 0.00079)
|  | 300 | (0.03760, 0.01939, 0.00060) |  | (0.04709, 0.02150, 0.00091)
| Discrete | 500 | (0.02910, 0.01505, 0.00037) |  | (0.04287, 0.01808, 0.00075)
|  | 400 | (0.03272, 0.01686, 0.00046) |  | (0.04432, 0.01962, 0.00081)
|  | 300 | (0.03619, 0.01879, 0.00057) |  | (0.04756, 0.02146, 0.00094)

We compare the proposed IV estimation method in Section 3.3 with the
regression method without any adjustment for confounding (Li et al., 2023a, ).
To evaluate the quality of estimation, we consider three metrics, the average
maximum absolute deviation, the mean absolute deviation, and the mean square
deviation between true coefficients and estimates over 1000 runs. As
demonstrated in Table 2, GrIVET enhances parameter estimation by accounting
for latent confounding. As anticipated, GrIVET’s estimation improves with
increasing sample size $n$, while the naive regression method (Li et al.,
2023a, ) remains inconsistent. Furthermore, GrIVET’s advantages become more
pronounced when stronger confounding effects are present, as evidenced by
additional simulations in the Supplementary Materials.

#### Inference.

We now evaluate the empirical performance of the proposed tests in terms of
size and power. For the empirical size, we calculate the percentage of times
$H_{0}$ is rejected out of 1000 simulations when $H_{0}$ is true. For the
power, we consider three alternative hypotheses $H_{a}$, where all the edges
in $H_{0}$ exist. The empirical power of a test is the percentage of times
$H_{0}$ is rejected out of 1000 simulations when $H_{a}$ is true. The
adjacency matrix $\mathbf{U}$ is modified according to the null and
alternative hypotheses.

  * •

Hub graph, fixed $\mathcal{H}$. For the size, consider
$\mathcal{H}=\\{(2,7)\\}$, $\mathcal{H}=\\{(2,7),(7,12),(12,17)\\}$, and
$\mathcal{H}=\\{(2,7),(7,12),(12,17),(17,22),(22,27)\\}$. For the power,
consider $\mathcal{H}=\\{(1,2)\\}$, $\mathcal{H}=\\{(1,2),(1,12),(1,22)\\}$,
and $\mathcal{H}=\\{(1,2),(1,12),(1,22),(1,32),(1,42)\\}$.

  * •

Random graph, fixed $\mathcal{H}$. We consider $\mathcal{H}=\\{(1,6)\\}$,
$\mathcal{H}=\\{(1,6),(6,11),(11,16)\\}$, and
$\mathcal{H}=\\{(1,6),(6,11),(11,16),(16,21),(21,26)\\}$ for both size and
power.

  * •

Random graph, random $\mathcal{H}$. We also consider testing 50 randomly
selected edges individually. Here, a random graph is generated so that 20 of
these selected edges are present in the true DAG (i.e., $H_{a}$ is valid). As
a result, for every selected edge, $H_{0}$ holds in roughly $600$ repetitions
and $H_{a}$ holds in roughly $400$ repetitions.

As shown in Table 3 for fixed $\mathcal{H}$, empirical sizes are close to the
nominal $\alpha=0.05$ under $H_{0}$, and the proposed test enjoys desirable
power under $H_{a}$. Figure 2 presents similar results for testing random
$\mathcal{H}$. The Supplementary Materials display that the sampling
distribution of the test statistic is close to the derived asymptotic
distribution in Theorem 2. Additional simulation details and results are also
available in Supplementary Materials.

Table 3: Empirical size for GrIVET at nominal level $\alpha=0.05$, respectively for $|\mathcal{D}|=1$, $|\mathcal{D}|=3$ and $|\mathcal{D}|=5$, over 1000 simulation replications. Graph | Intervention | $n$ |  Size ($|\mathcal{D}|=1,3,5$) |  Power ($|\mathcal{D}|=1,3,5$)
---|---|---|---|---
Hub | Continuous | 500 | (0.028,0.026,0.029) | (1.000,1.000,1.000)
|  | 400 | (0.043,0.038,0.035) | (1.000,1.000,1.000)
|  | 300 | (0.037,0.030,0.034) | (1.000,1.000,1.000)
| Discrete | 500 | (0.036,0.040,0.027) | (1.000,1.000,1.000)
|  | 400 | (0.051,0.040,0.040) | (1.000,1.000,1.000)
|  | 300 | (0.052,0.041,0.035) | (1.000,1.000,1.000)
Random | Continuous | 500 | (0.038,0.037,0.026) | (1.000,1.000,1.000)
|  | 400 | (0.033,0.031,0.028) | (1.000,1.000,1.000)
|  | 300 | (0.033,0.025,0.030) | (1.000,1.000,1.000)
| Discrete | 500 | (0.040,0.029,0.027) | (1.000,1.000,1.000)
|  | 400 | (0.042,0.034,0.040) | (1.000,1.000,1.000)
|  | 300 | (0.029,0.033,0.034) | (1.000,1.000,1.000)

(a) $H_{0}$ holds

(b) $H_{a}$ holds

Figure 2: The boxplots of the empirical rejection probabilities for testing
randomly selected edges. The nominal level is $\alpha=0.05$.

###  6.2 ADNI data analysis

In this subsection, GrIVET is applied to analyze the Alzheimer’s Disease
Neuroimaging Initiative (ADNI) dataset (available at
https://adni.loni.usc.edu). The goal is to infer gene pathways related to
Alzheimer’s Disease (AD) in order to elucidate the gene-gene interactions in
AD/cognitive impairment patients and healthy individuals, respectively.

#### Dataset.

The dataset comprises gene expression levels adjusted for five covariates:
gender, handedness, education level, age, and intracranial volume. For data
analysis, we select genes with at least one SNP at a marginal significance
level below $10^{-14}$, resulting in $p=21$ genes as primary variables. For
these genes, we further extract their marginally most correlated two SNPs,
yielding $q=42$ SNPs as unspecified intervention variables for subsequent data
analysis. All gene expression levels are normalized.

The dataset initially categorizes individuals into four groups: Alzheimer’s
Disease (AD), Early Mild Cognitive Impairment (EMCI), Late Mild Cognitive
Impairment (LMCI), and Cognitive Normal (CN). For our analysis, we treat 247
CN individuals as controls and the remaining 462 individuals as cases (AD-
MCI). We then use the gene expressions and the SNPs to infer gene pathways for
the 462 AD-MCI and 247 CN control cases, respectively.

(a) AD-MCI

(b) CN

Figure 3: Display of the genes associated with proposed tests. (a) and (b):
Solid/dashed arrows indicate significant/insignificant edges at $\alpha=0.05$
after adjustment for multiplicity by the Bonferroni-Holm correction.

#### Hypotheses.

We focus on statistical inferences related to genes APP and CASP3 (Julia and
Goate,, 2017; Su et al.,, 2001). As in Figure 3, for each edge $(k,j)$, we
consider testing $H_{0}:\mathrm{U}_{kj}=0\text{ versus
}H_{a}:\mathrm{U}_{kj}\neq 0$.

#### Results.

Figure 3 displays the p-values and significant results under the level
$\alpha=0.05$ after the Holm-Bonferroni adjustment for $2\times 7=14$ tests.
The tests exhibit strong evidence for the presence of
$\\{\mathrm{LRP1}\rightarrow\mathrm{CASP3},\
\mathrm{APP}\rightarrow\mathrm{APOE}\\}$ in the AD-MCI group, but no evidence
in the CN group. Meanwhile, this result suggests the presence of connections
$\\{\mathrm{CAPN1}\rightarrow\mathrm{CASP3},\
\mathrm{ATP5F1}\rightarrow\text{CASP3}\\}$ in the CN group but not so in the
AD-MCI group. In both groups, we identify directed connection
$\mathrm{APP}\rightarrow\text{APBB1}$. Figure 4 shows the residual correlation
matrices for both groups, suggesting the existence of unmeasured confounding.
The Supplementary Materials include normal Q-Q plots of residuals,
demonstrating that the normality assumption is approximately satisfied for
both groups.

(a) AD-MCI

(b) CN

Figure 4: Display of residual correlation matrices for AD-MCI and CN groups.

Some of our discoveries agree with the existing findings. Specifically, our
result indicates the presence of connection APP $\rightarrow$ APOE for the AD-
MCI group, but not for the CN group, which seems consistent with the knowledge
that APP and APOE are functionally linked in brain cholesterol metabolism (Liu
et al.,, 2017) and the contributions of APOE to the pathophysiology of AD
(Bu,, 2009). The connection LRP1 $\rightarrow$ CASP3 also differs in AD-MCI
and CN groups, which may serve to support the conclusion that activated CASP3
may be a factor in functional decline and may have an important role in
neuronal cell death and plaque formation in AD brain (Su et al.,, 2001) given
the finding that both APOE and its receptor LRP1 are present in amyloid
plaques (Poirier,, 1996). Moreover, the connection CAPN1 $\rightarrow$ CDK5R1
discovered in both groups can be found in the AlzNet database (interaction ID
24614).

##  7 Discussion

This article proposes a novel instrumental variable procedure that integrates
causal discovery and inference for a Gaussian directed acyclic graph with
hidden confounders. One future research direction is to develop methodologies
for analyzing discrete/mixed-type (primary variable) data. Additionally, the
present work uses individual-level data from a single study for causal
discovery and inference. In many real applications, due to privacy concerns
and ownership restrictions, the data are only available in the form of summary
statistics (e.g., GWAS summary data) or in other privatized forms. Extending
GrIVET to leverage these data is an important topic. Furthermore,
multisource/decentralized data are ubiquitous, raising new challenges in
communication, privacy, and handling of corrupted data. It would be promising
to employ modern machine learning techniques, such as federated learning
(Xiong et al.,, 2021; Gao et al.,, 2021), to address these challenges and
fully unleash the potential of large-scale causal discovery and inference.

Finally, we discuss two limitations of the present work.

  * •

GrIVET necessitates the availability of valid IVs for each primary variable
due to the hardness of causal identification in the presence of hidden
confounding. In genetic research, there is an ample supply of genetic variants
(e.g., SNPs) serving as IVs. Nonetheless, obtaining valid IVs can be
challenging in certain applications. It is thus crucial to investigate the
potential for causal discovery even when faced with an insufficient number of
IVs.

  * •

For inference, Theorem 2 requires that
$\mathbb{P}(\widehat{\mathcal{G}}^{+}=\mathcal{G}^{+})\to 1$, which is
guaranteed by Condition (C2) in Theorem 3. Fulfilling this requirement can be
challenging; in such cases, one might turn to the post-selection inference
framework (Berk et al.,, 2013) by concentrating on the parameters within the
selected model. However, the test results should be meticulously interpreted,
as these parameters cease to be causal or structural (Berk et al.,, 2013)
unless $\mathbb{P}(\widehat{\mathcal{G}}^{+}=\mathcal{G}^{+})\to 1$. In
essence, (C2) enables the causal meaning of the tested parameters to be
carried over to finite-sample inference. Exploring ways to lift the signal
strength condition while preserving the causal interpretation for statistical
inference after DAG structure learning (Wang et al.,, 2023) is an important
research topic.

##  Appendix A Appendix

#### Definition of d-separation (Pearl,, 2009).

Consider a DAG $\mathcal{G}$ with node variables
$(Z_{1},\ldots,Z_{d})^{\top}$. Nodes $Z_{k}$ and $Z_{j}$ are adjacent if
$Z_{k}\to Z_{j}$ or $Z_{k}\leftarrow Z_{j}$. An undirected path between
$Z_{k}$ and $Z_{j}$ in $\mathcal{G}$ is a sequence of distinct nodes
$(Z_{k},\ldots,Z_{j})$ such that all pairs of successive nodes in the sequence
are adjacent. A non-endpoint node $Z_{m}$ on an undirected path
$(Z_{k},\ldots,Z_{m-1},Z_{m},Z_{m+1},\ldots,Z_{j})$ is called a collider if
$Z_{m-1}\to Z_{m}\leftarrow Z_{m+1}$. Otherwise, it is called a non-collider.
Let $A\subseteq\\{1,\ldots,d\\}$, where $A$ does not contain $k$ and $j$. Then
$\bm{Z}_{A}$ is said to block an undirected path $(Z_{k},\ldots,Z_{j})$ if at
least one of the following holds: (1) the undirected path contains a non-
collider that is in $\bm{Z}_{A}$, or (2) the undirected path contains a
collider that is not in $\bm{Z}_{A}$ and has no descendant in $\bm{Z}_{A}$. A
node $Z_{k}$ is d-separated from $Z_{j}$ given $\bm{Z}_{A}$ if $\bm{Z}_{A}$
block every undirected path between $Z_{k}$ and $Z_{j}$; $k\neq j$.

#### Additional discussion of Figure 1 (a).

Let $(k,j)\in\mathcal{E}^{+}$ and suppose all IVs are valid. We explain why
$\bm{X}_{\textnormal{{ca}}(k)}$ may not be valid IVs after conditioning on
$\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}}$, as mentioned in Section 3.3.
Let $l\in\textnormal{{ca}}(k)$ and $m\in\textnormal{{me}}(k,j)$ such that
$Y_{k}$ is an unmediated parent of $Y_{m}$. Note that in Figure 1 (a) of the
main text, whenever $\eta\to Y_{m}$, then
$\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}}$ does not d-separate
$\bm{X}_{\textnormal{{ca}}(k)}$ and $\eta$, since $Y_{m}$ is a collider in the
undirected path $(X_{l},Y_{k},Y_{m},\eta,Y_{j})$. As a result,
$\bm{X}_{\textnormal{{ca}}(k)}$ and $\eta$ can be associated conditioned on
$\bm{Y}_{\textnormal{{an}}(j)\setminus\\{k\\}}$.

#### Additional discussion on identification of $\mathbf{U}$.

We have the following result.

######  Lemma 1.

In (1), assume $\bm{X}$ and $\bm{\varepsilon}$ are independent.

  1. (A)

$\mathbb{E}(Y_{k}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$ is a linear
combination of $(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$.

  2. (B)

$\mathbb{E}(\varepsilon_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})$ is a
linear combination of
$(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{ca}}(k)^{c}})$.

###### Proof.

Here, (A) follows directly from (1). For (B), we have

$\mathbb{E}(\varepsilon_{j}\mid\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X})=\mathbb{E}(\varepsilon_{j}\mid\bm{\varepsilon}_{\textnormal{{nm}}(k,j)},\bm{X})=\mathbb{E}(\varepsilon_{j}\mid\bm{\varepsilon}_{\textnormal{{nm}}(k,j)})=\bm{\pi}^{\top}\bm{\varepsilon}_{\textnormal{{nm}}(k,j)},$

where the last equality is due to the normality of $\bm{\varepsilon}$.
Finally, in (1), we immediately have
$\bm{\varepsilon}_{\textnormal{{nm}}(k,j)}$ is linear in
$(\bm{Y}_{\textnormal{{nm}}(k,j)},\bm{X}_{\textnormal{{ca}}(k)^{c}})$. ∎

Now, we show that $\operatorname{Cov}(\bm{\varepsilon},\bm{X})=\bm{0}$ is
sufficient to derive the identification results in Section 3.3. Given random
variables $\zeta$ and $\bm{\xi}$, let $\mathbb{L}(\zeta\mid\bm{\xi})$ be the
best linear approximation of $\zeta$ using $\bm{\xi}$, namely
$\mathbb{L}(\zeta\mid\bm{\xi})=\widetilde{\bm{\omega}}^{\top}\bm{\xi}$ where

$\widetilde{\bm{\omega}}=\operatorname*{\arg\min}_{\bm{\omega}}\
\operatorname{\mathbb{E}}(\zeta-\bm{\omega}^{\top}\bm{\xi})^{2}.$

For random variables $\zeta$, $\zeta^{\prime}$, and $\bm{\xi}$, we have that
(a)
$\mathbb{L}(\zeta+\zeta^{\prime}\mid\bm{\xi})=\mathbb{L}(\zeta\mid\bm{\xi})+\mathbb{L}(\zeta^{\prime}\mid\bm{\xi})$,
(b) $\mathbb{L}(c\zeta\mid\bm{\xi})=c\mathbb{L}(\zeta\mid\bm{\xi})$ for
$c\in\mathbb{R}$, (c) $\mathbb{L}(\zeta\mid\bm{\xi})=0$ if
$\operatorname{Cov}(\zeta,\bm{\xi})=\bm{0}$, (d)
$\mathbb{L}(\zeta\mid\bm{\xi})=\zeta$ if
$\zeta\in\operatorname{Span}(\bm{\xi})$, and (e)
$\mathbb{L}(\zeta\mid\bm{\xi})=\mathbb{L}(\zeta\mid\mathbf{A}\bm{\xi})$ for
invertible $\mathbf{A}$. Thus, $\mathbb{L}(\cdot\mid\star)$ mimics
$\mathbb{E}(\cdot\mid\star)$, and Lemma 2 holds. The proof is similar to that
of Lemma 1.

######  Lemma 2.

In (1), Lemma 1 holds with $\mathbb{E}(\cdot\mid\star)$ being replaced by
$\mathbb{L}(\cdot\mid\star)$.

As a result, if $\bm{X}$ and $\bm{\varepsilon}$ are uncorrelated as in (1),
the derivation in Section 3.3 holds with $\mathbb{E}(\cdot\mid\star)$ being
replaced by $\mathbb{L}(\cdot\mid\star)$.

## References

  * Angrist et al., (1996) Angrist, J. D., Imbens, G. W., and Rubin, D. B. (1996).  Identification of causal effects using instrumental variables.  Journal of the American Statistical Association, 91(434):444–455.
  * Aragam et al., (2019) Aragam, B., Amini, A. A., and Zhou, Q. (2019).  Globally optimal score-based learning of directed acyclic graphs in high-dimensions.  In Proceedings of the 33rd International Conference on Neural Information Processing Systems, pages 4450–4462.
  * Berk et al., (2013) Berk, R., Brown, L., Buja, A., Zhang, K., and Zhao, L. (2013).  Valid post-selection inference.  The Annals of Statistics, 41(2):802–837.
  * Bickel et al., (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of Lasso and Dantzig selector.  The Annals of Statistics, 37(4):1705–1732.
  * Bu, (2009) Bu, G. (2009).  Apolipoprotein E and its receptors in Alzheimer’s disease: pathways, pathogenesis and therapy.  Nature Reviews Neuroscience, 10(5):333–344.
  * Burgess et al., (2020) Burgess, S., Foley, C. N., Allara, E., Staley, J. R., and Howson, J. M. (2020).  A robust and efficient method for Mendelian randomization with hundreds of genetic variants.  Nature Communications, 11(1):1–11.
  * Chakrabortty et al., (2018) Chakrabortty, A., Nandy, P., and Li, H. (2018).  Inference for individual mediation effects and interventional effects in sparse high-dimensional causal graphical models.  arXiv preprint arXiv:1809.10652.
  * Chen et al., (2018) Chen, C., Ren, M., Zhang, M., and Zhang, D. (2018).  A two-stage penalized least squares method for constructing large systems of structural equations.  Journal of Machine Learning Research, 19(1):40–73.
  * Colombo et al., (2012) Colombo, D., Maathuis, M. H., Kalisch, M., and Richardson, T. S. (2012).  Learning high-dimensional directed acyclic graphs with latent and selection variables.  The Annals of Statistics, 40(1):294–321.
  * Drton and Maathuis, (2017) Drton, M. and Maathuis, M. H. (2017).  Structure learning in graphical modeling.  Annual Review of Statistics and Its Application, 4:365–393.
  * Frot et al., (2019) Frot, B., Nandy, P., and Maathuis, M. H. (2019).  Robust causal structure learning with some hidden variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(3):459–487.
  * Gao et al., (2021) Gao, E., Chen, J., Shen, L., Liu, T., Gong, M., and Bondell, H. (2021).  FedDAG: Federated DAG structure learning.  Transactions on Machine Learning Research.
  * Ghoshal and Honorio, (2018) Ghoshal, A. and Honorio, J. (2018).  Learning linear structural equation models in polynomial time and sample complexity.  In International Conference on Artificial Intelligence and Statistics, pages 1466–1475. PMLR.
  * Glymour et al., (2019) Glymour, C., Zhang, K., and Spirtes, P. (2019).  Review of causal discovery methods based on graphical models.  Frontiers in Genetics, 10.
  * Grimmer et al., (2020) Grimmer, J., Knox, D., and Stewart, B. M. (2020).  Naïve regression requires weaker assumptions than factor models to adjust for multiple cause confounding.  arXiv preprint arXiv:2007.12702.
  * Gu et al., (2019) Gu, J., Fu, F., and Zhou, Q. (2019).  Penalized estimation of directed acyclic graphs from discrete data.  Statistics and Computing, 29(1):161–176.
  * Guo et al., (2018) Guo, Z., Kang, H., Tony Cai, T., and Small, D. S. (2018).  Confidence intervals for causal effects with invalid instruments by using two-stage hard thresholding with voting.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(4):793–815.
  * Heinze-Deml et al., (2018) Heinze-Deml, C., Maathuis, M. H., and Meinshausen, N. (2018).  Causal structure learning.  Annual Review of Statistics and Its Application, 5:371–391.
  * Janková and van de Geer, (2018) Janková, J. and van de Geer, S. (2018).  Inference in high-dimensional graphical models.  In Handbook of Graphical Models, pages 325–350. CRC Press.
  * Julia and Goate, (2017) Julia, T. and Goate, A. M. (2017).  Genetics of $\beta$-amyloid precursor protein in Alzheimer’s disease.  Cold Spring Harbor Perspectives in Medicine, 7(6).
  * Kang et al., (2016) Kang, H., Zhang, A., Cai, T. T., and Small, D. S. (2016).  Instrumental variables estimation with some invalid instruments and its application to Mendelian randomization.  Journal of the American Statistical Association, 111(513):132–144.
  * Kertel et al., (2022) Kertel, M., Harmeling, S., and Pauly, M. (2022).  Learning causal graphs in manufacturing domains using structural equation models.  arXiv preprint arXiv:2210.14573.
  * Lee and Li, (2022) Lee, K.-Y. and Li, L. (2022).  Functional structural equation model.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 84(2):600–629.
  * (24) Li, C., Shen, X., and Pan, W. (2020a).  Likelihood ratio tests for a large directed acyclic graph.  Journal of the American Statistical Association, 115(531):1304–1319.
  * (25) Li, C., Shen, X., and Pan, W. (2023a).  Inference for a large directed acyclic graph with unspecified interventions.  Journal of Machine Learning Research, 24(73):1–48.
  * (26) Li, C., Shen, X., and Pan, W. (2023b).  Nonlinear causal discovery with confounders.  Journal of the American Statistical Association, pages 1–32.
  * Li et al., (2022) Li, L., Shi, C., Guo, T., and Jagust, W. J. (2022).  Sequential pathway inference for multimodal neuroimaging analysis.  Stat, 11(1):e433.
  * (28) Li, Y., Torralba, A., Anandkumar, A., Fox, D., and Garg, A. (2020b).  Causal discovery in physical systems from videos.  In Proceedings of the 34th International Conference on Neural Information Processing Systems, pages 9180–9192.
  * Liu et al., (2017) Liu, Z., Zhang, M., Xu, G., Huo, C., Tan, Q., Li, Z., and Yuan, Q. (2017).  Effective connectivity analysis of the brain network in drivers during actual driving using near-infrared spectroscopy.  Frontiers in Behavioral Neuroscience, 11:211.
  * Lousdal, (2018) Lousdal, M. L. (2018).  An introduction to instrumental variable assumptions, validation and estimation.  Emerging Themes in Epidemiology, 15(1):1–7.
  * Meinshausen and Bühlmann, (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the Lasso.  The Annals of Statistics, 34(3):1436–1462.
  * Murray, (2006) Murray, M. (2006).  Avoiding invalid instruments and coping with weak instruments.  Journal of Economic Perspectives, 20(4):111–132.
  * Oates et al., (2016) Oates, C. J., Smith, J. Q., and Mukherjee, S. (2016).  Estimating causal structure using conditional DAG models.  Journal of Machine Learning Research, 17(1):1880–1903.
  * Pearl, (2009) Pearl, J. (2009).  Causality.  Cambridge University Press.
  * Peters and Bühlmann, (2014) Peters, J. and Bühlmann, P. (2014).  Identifiability of Gaussian structural equation models with equal error variances.  Biometrika, 101(1):219–228.
  * Poirier, (1996) Poirier, J. (1996).  Apolipoprotein E in the brain and its role in Alzheimer’s disease.  Journal of Psychiatry and Neuroscience, 21(2):128–134.
  * Rajendran et al., (2021) Rajendran, G., Kivva, B., Gao, M., and Aragam, B. (2021).  Structure learning in polynomial time: Greedy algorithms, Bregman information, and exponential families.  In Advances in Neural Information Processing Systems, volume 34, pages 18660–18672.
  * Reisach et al., (2021) Reisach, A., Seiler, C., and Weichwald, S. (2021).  Beware of the simulated DAG! Causal discovery benchmarks may be easy to game.  Advances in Neural Information Processing Systems, 34:27772–27784.
  * Sachs et al., (2005) Sachs, K., Perez, O., Pe’er, D., Lauffenburger, D. A., and Nolan, G. P. (2005).  Causal protein-signaling networks derived from multiparameter single-cell data.  Science, 308(5721):523–529.
  * Shah et al., (2020) Shah, R. D., Frot, B., Thanei, G.-A., and Meinshausen, N. (2020).  Right singular vector projection graphs: fast high dimensional covariance matrix estimation under latent confounding.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(2):361–389.
  * Shen et al., (2012) Shen, X., Pan, W., and Zhu, Y. (2012).  Likelihood-based selection and sharp parameter estimation.  Journal of the American Statistical Association, 107(497):223–232.
  * Shi and Li, (2021) Shi, C. and Li, L. (2021).  Testing mediation effects using logic of Boolean matrices.  Journal of the American Statistical Association, pages 1–14.
  * Shi et al., (2023) Shi, C., Zhou, Y., and Li, L. (2023).  Testing directed acyclic graph via structural, supervised and generative adversarial learning.  Journal of the American Statistical Association, pages 1–24.
  * Shimizu et al., (2006) Shimizu, S., Hoyer, P. O., Hyvärinen, A., and Kerminen, A. (2006).  A linear non-Gaussian acyclic model for causal discovery.  Journal of Machine Learning Research, 7:2003–2030.
  * Shojaie and Michailidis, (2010) Shojaie, A. and Michailidis, G. (2010).  Penalized likelihood methods for estimation of sparse high-dimensional directed acyclic graphs.  Biometrika, 97(3):519–538.
  * Su et al., (2001) Su, J. H., Zhao, M., Anderson, A. J., Srinivasan, A., and Cotman, C. W. (2001).  Activated caspase-3 expression in Alzheimer’s and aged control brain: correlation with Alzheimer pathology.  Brain Research, 898(2):350–357.
  * Vowels et al., (2021) Vowels, M. J., Camgoz, N. C., and Bowden, R. (2021).  D’ya like DAGs? A survey on structure learning and causal discovery.  ACM Computing Surveys (CSUR).
  * Wang et al., (2023) Wang, Y. S., Kolar, M., and Drton, M. (2023).  Confidence sets for causal orderings.  arXiv preprint arXiv:2305.14506.
  * Windmeijer et al., (2019) Windmeijer, F., Farbmacher, H., Davies, N., and Davey Smith, G. (2019).  On the use of the lasso for instrumental variables estimation with some invalid instruments.  Journal of the American Statistical Association, 114(527):1339–1350.
  * Xiong et al., (2021) Xiong, R., Koenecke, A., Powell, M., Shen, Z., Vogelstein, J. T., and Athey, S. (2021).  Federated causal inference in heterogeneous observational data.  arXiv preprint arXiv:2107.11732.
  * Xue and Pan, (2020) Xue, H. and Pan, W. (2020).  Inferring causal direction between two traits in the presence of horizontal pleiotropy with GWAS summary data.  PLoS Genetics, 16(11):e1009105.
  * Yuan et al., (2019) Yuan, Y., Shen, X., Pan, W., and Wang, Z. (2019).  Constrained likelihood for reconstructing a directed acyclic Gaussian graph.  Biometrika, 106(1):109–125.
  * Zhao et al., (2022) Zhao, R., He, X., and Wang, J. (2022).  Learning linear non-Gaussian directed acyclic graph with diverging number of nodes.  Journal of Machine Learning Research, 23(269):1–34.
  * Zheng et al., (2018) Zheng, X., Aragam, B., Ravikumar, P., and Xing, E. P. (2018).  DAGs with NO TEARS: continuous optimization for structure learning.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 9492–9503.