# Perturbations and Causality in Gaussian Latent Variable Models

Armeen Taeb and Peter Bühlmann
Seminar for Statistics, ETH Zürich Correspondence email:
<EMAIL_ADDRESS>
###### Abstract

Causal inference is a challenging problem with observational data alone. The
task becomes easier when having access to data from perturbing the underlying
system, even when happening in a non-randomized way: this is the setting we
consider, encompassing also latent confounding variables. To identify causal
relations among a collections of covariates and a response variable, existing
procedures rely on at least one of the following assumptions: i) the response
variable remains unperturbed, ii) the latent variables remain unperturbed, and
iii) the latent effects are dense. In this paper, we examine a perturbation
model for interventional data, which can be viewed as a mixed-effects linear
structural causal model, over a collection of Gaussian variables that does not
satisfy any of these conditions. We propose a maximum-likelihood estimator –
dubbed _DirectLikelihood_ – that exploits system-wide invariances to uniquely
identify the population causal structure from unspecific perturbation data,
and our results carry over to linear structural causal models without
requiring Gaussianity. We illustrate the utility of our framework on synthetic
data as well as real data involving California reservoirs and protein
expressions.

##  1 Introduction

Identifying causal relations from observational data is challenging and one
can often only identify the corresponding _Markov equivalence class_ (MEC). At
the opposite pole are designed randomized experiments [28]: they are the gold
standard for causal inference but the feasibility to do the randomization is
hindered by cost or ethical reasons. It is possible though, under some
assumptions, to exploit non-specific and non-randomized interventions or
perturbations which frequently arise in many datasets: this is the topic of
the current paper.

In the context of observational data from structural causal models [24, 21],
one possibility is to find the MEC of directed acyclic graphs under the
faithfulness assumption [33] or the beta-min condition [32]. Some of the well-
known algorithms for structure learning of MECs with observational data
include the constraint based PC algorithm [30], score based greedy algorithm
GES [4] and hybrid methods that integrate constraint based and score based
methods such as ARGES [19]. In many applications though, we have available
both observational and unspecific interventional or perturbation data, where
the latter are coming from non-randomized experiments with unknown targets. In
genomics, for example, with the advance of gene editing technologies, high
throughput interventional gene expression data is being produced [9].
Interventional data can be viewed as _perturbations_ to components of the
system and can offer substantial gain in identifiability: [13] demonstrated
that combining interventional with observational data reduces ambiguity and
enhances identifiability to a smaller equivalence class than the MEC, known as
the I-MEC (Interventional MEC). A variety of methods have been proposed for
causal structure learning from observational and interventional data. This
includes the modified GES algorithm by [13] known as GIES, permutation-based
causal structure learning [34], penalized maximum-likelihood procedure in
Gaussian models [15], and methods based on a causal invariance framework [17,
23] building on a concept of stability [7, 8]. For a more comprehensive list,
see [10] and the references therein.

A common challenge for accurate structure learning is that there may be latent
variables for which it is expensive or impossible to obtain sample
observations. Such unobserved variables pose a significant difficulty as the
causal graphical model structure is not closed under marginalization;
therefore, the graphical structure corresponding to the marginal distribution
of the observed variables consists of potentially many confounding
dependencies that are induced due to the marginalization over the latent
variables. There are causal structural learning methods that account for the
presence of latent variables. In the observational setting, two prominent
examples are the Fast Causal Inference [30] and its variant RFCI [6] for DAG
learning and the two-stage deconfounding procedure [12] involving the sparse-
plus-low rank decomposition framework [3] as the first stage and the standard
DAG learning procedure in the second stage. As discussed earlier, the (R)FCI
algorithms or the two-stage deconfounding procedure will only enable to infer
a certain MEC but not the causal parameter and structure itself. In the joint
observational and interventional setting with unperturbed latent variables and
only shift interventions on the observed covariates, Causal Dantzig [25]
consistently estimates the causal relations of a response variable assuming
that the interventions do not directly affect the response variable. Such an
assumption is relaxed in the backShift procedure [27] while it still requires
that the latent variables remain unperturbed for identifying the causal
structure.

Guaranteed identifiability using these previous techniques for perturbation
data relies on at least one of the following assumptions: i) the response
variable remains unperturbed, ii) the latent variables remain unperturbed, and
iii) the latent effects are dense. In this paper, we propose a modeling
framework and an estimator that does not rely on any of these assumptions and
yet identifies the population DAG structure. Fig 1 demonstrates a toy example
of our setup among $4$ observed variables $X_{1},X_{2},X_{3},X_{4}$ and latent
variables $H$, and $\mathcal{A}$ represents external variables (to the
graphical structure among observed and latent variables) that provide
perturbations.

Figure 1: Toy example of 4 observed variables $X_{1},X_{2},X_{3},X_{4}$ and
latent variables $H$ where solid lines are connections among observed
variables and dotted lines are connections between observed and latent
variables; left: without perturbations, right: perturbations $\mathcal{A}$ on
all components indicated with red dotted lines.

We consider a Gaussian structural causal model (SCM) specifying the
perturbation model and the relationship between $p$ observed variables
$X\in\mathbb{R}^{p}$ and latent variables $H\in\mathbb{R}^{h}$. We consider
the setting with heterogeneous grouped data from different environments
$e\in{\cal E}$. Here $e$ denotes the index of an environment or a sub-
population and ${\cal E}$ is the space of different observed environments. As
we will formalize in Section 2.1, each group or environment $e$ corresponds to
some perturbations of the underlying SCM. The grouped data, across different
environments, is denoted by $(X^{e},H^{e})$ with $e\in{\cal E}$. The SCM is
parameterized by a connectivity matrix encoding the causal relationship among
the observed variables, a coefficient matrix encoding the latent variable
effects, and nuisance parameters involving the noise variances and
perturbation magnitudes among all of the variables. A key property of this
modeling framework is that the connectivity matrix and the latent variable
coefficient matrix remain _invariant_ across of all the perturbation
environments. With this insight, we propose a maximum-likelihood estimator –
dubbed _DirectLikelihood_ – to score a given DAG structure. _DirectLikelihood_
provides a flexible framework to incorporate additional knowledge including
do-interventions when their intervention-locations are known or additional
information on the perturbation structure (such as statistically identical
perturbations on all of the observed variables). Further, the framework can be
specialized to the setting considered by [27, 25] where the latent variables
are not perturbed across environments (i.e. $\mathcal{A}$ does not point to
$H$ in Fig 1), or to the setting where there is no latent confounding (i.e.
$H$ does not point to the covariates in Fig 1).

Besides the novel methodology, we provide conditions for which
_DirectLikelihood_ correctly identifies the population DAG structure. In
particular, we demonstrate that with at least two interventional environments,
where one of the environments consists of sufficiently large interventions on
each of the observed variables, and the latent effects satisfying a _latent
materiality_ assumption, _DirectLikelihood_ provides consistent estimates. The
_latent materiality_ assumption for an environment $e$ states that the latent
variables induce confounding dependencies among the observed variables;
formally, there exists at least one pair of variables $(X^{e}_{k},X^{e}_{l})$
such that $X^{e}_{k}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X^{e}_{l}~{}|~{}\\{X^{e}_{\backslash\\{k,l\\}},H^{e}\\}$
and $X^{e}_{k}\not\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X^{e}_{l}~{}|~{}X^{e}_{\backslash\\{k,l\\}}$,
where $X^{e}_{\backslash\\{k,l\\}}$ denotes the collections of variables
$X^{e}$ excluding $X^{e}_{k},X^{e}_{l}$. The _latent materiality_ assumption
is substantially weaker than the latent denseness assumption required in the
two-stage deconfounding procedure in [12] which insists that there are many
pairs of variables satisfying the condition above. Our theoretical results are
further specialized to the setting where the latent variables remain
unperturbed across all of the environments. When the latent variables are
unperturbed, _DirectLikelihood_ requires no assumption on the latent structure
for identifiability, whereas the two-stage deconfounding procedure still
requires latent denseness. We remark that the main focus of our analysis is on
identifiability guarantees, and we discuss in Section 7 future work on
understanding high-dimensional consistency properties of the
_DirectLikelihood_ procedure.

Further, we highlight a connection between distributional robustness and the
causal parameters in our perturbation model. Specifically, we prove that the
population causal parameters are minimizers of the worst-case risk over the
space of DAGs and distributional shifts from a certain perturbation class.
Here, the risk is measured by the Kullback-Libeler divergence between the
estimated and population Gaussian distributions. As with the DAG
identifiability, the relation between causality and distributional robustness
relies on the stringent assumption that the perturbations do not directly
affect the response variable or the latent variables [2, 26]. The results in
this paper provided a more complete picture on the connection between
perturbations, causality, and distributional robustness (see also Table 1).

As our final contribution, we propose an optimization procedure to solve
_DirectLikelihood_ in Section 5 and demonstrate the utility of our proposed
estimator with synthetic data and real data involving California reservoirs
and protein expression data in Section 6. The estimates provided by
_DirectLikelihood_ offer improvements over previous approaches in multiple
respects. First, the causal graphical structure that is obtained by
_DirectLikelihood_ is accurate even when there are interventions on the
response variable and the latent variables, or when the latent effects are not
dense across the observed variables. Previous methods, on the other hand, may
provide inaccurate estimates in such settings. Second, _DirectLikelihood_
produces models with few false positives and large number of true positives
(with respect to graphical structure) with moderate sample sizes, as compared
to competing methods like backShift that require much larger data. Finally, in
the analysis with real data, we demonstrate that accounting for latent effects
via the _DirectLikelihood_ procedure yields models that are more sensible
(with fewer spurious edges) than if latent variables are not taken into
account.

The outline of this paper is as follows. In Section 2, we describe the model
for observational and perturbation data and its representation as a mixed-
effects model, and then present the maximum-likelihood estimator
_DirectLikelihood_ to score a given DAG structure. In Section 3, we provide
theoretical guarantees for the optimally scoring DAGs (scored via
_DirectLikelihood_). In Section 4, the connection between the causal
parameters of the proposed perturbation model and distributional robustness is
explored. In Section 5, we present an optimization strategy for solving
_DirectLikelihood_ for a given DAG structure and how to use it to obtain the
best scoring DAGs. In Section 6, we demonstrate the utility of our approach
with real and synthetic data. We conclude with future research directions in
Section 7.

###  1.1 Related work

We have mentioned differences to backshift [27] and two-stage deconfounding
procedures and provide more comparisons throughout the paper.
_DirectLikelihood_ is similar in spirit to approaches based on invariance
principles [23, 25] as it exploits certain model parameters (such as the
connectivity matrix and latent variable effects) remaining unchanged across
perturbations. However, a key difference between _DirectLikelihood_ and these
other techniques – in addition to being able to incorporate perturbations on
the latent variables – is that _DirectLikelihood_ models the entire system of
observed variables as opposed to just the regression of the response variable
and the remaining observed variables. The virtue of this system-wide modeling
is that all of the variables can experience perturbations without sacrificing
consistency guarantees while the methods in [23, 25] assume that the
perturbations do not directly affect the response variable. This perspective
was also adopted in the backShift procedure [27], although _DirectLikelihood_
can allow for perturbations on the latent variables. For a summary of the
assumptions for _DirectLikelihood_ as compared to competing methods, see Table
1.

Method | Perturbed response variable | Unperturbed latent variables | Perturbed latent variables
---|---|---|---
IV, ICP, Causal Dantzig | x | ✓ | x
two-stage deconfounding | ✓ |  ✓ latent denseness |  ✓ latent denseness
backShift | ✓ | ✓ | x
_DirectLikelihood_ | ✓ | ✓ |  ✓ _latent materiality_

Table 1: Comparison of _DirectLikelihood_ with competing methods in the
following settings: response variable is perturbed, latent variables are
unperturbed, and the latent variables are perturbed. The methods are
Instrumental Variables IV [1], Invariant Causal Prediction ICP [23], two-stage
deconfounding [12] tailored for observational and interventional data,
backShift [27] and our proposal _DirectLikelihood_.

###  1.2 Notation

We denote the identity matrix by $\mathcal{I}$, with the size being clear from
context. The collection of $d\times d$ symmetric matrices are denoted by
$\mathbb{S}^{d}$ and positive-semidefinite matrices by $\mathbb{S}^{d}_{+}$
and the collection of strictly positive-definite matrices by
$\mathbb{S}^{d}_{++}$. The collection of non-negative vectors in
$\mathbb{R}^{d}$ is denoted by $\mathbb{R}_{+}^{d}$ and strictly positive
vectors by $\mathbb{R}_{++}^{d}$. Given a matrix $M\in\mathbb{R}^{d\times d}$
and a set $S\subseteq\\{1,2,\dots,d\\}$, we denote the restriction of $M$ to
rows and columns indexed by $S$ by $[M]_{S}$. We denote the number of nonzeros
in a matrix $M\in\mathbb{R}^{p\times p}$ by $\|M\|_{\ell_{0}}$. We apply a
similar notation to count the number of edges in a graph. We denote the the
index set of the parents of a random variable $X_{p}$ by $\text{PA}(p)$ and
the index sets for the descendants and ancestors by $\text{DES}(p)$ and
$\text{ANC}(p)$, respectively. Further, letting $\mathcal{D}$ be the DAG
underlying a collection of variables $(X,H)$, we denote the subgraph of
$\mathcal{D}$ restricted to the variables $X$ by $\mathcal{D}_{X}$ and
likewise for $\mathcal{D}_{H}$. Given a matrix $M\in\mathbb{R}^{d_{1}\times
d_{2}}$, we denote $\|M\|_{2}$ to be the largest singular value (spectral
norm). For two vectors $z_{1},z_{2}\in\mathbb{R}^{d}$, we denote $z_{1}\succeq
z_{2}$ to denote element-wise inequality. Finally, for random variables
$V_{1},V_{2}$ and random vectors $Z$, we use the notation
$\rho(V_{1},V_{2}|Z)$ to denote the partial correlation between $V_{1}$ and
$V_{2}$ given $Z$.

##  2 Modeling framework and maximum-likelihood estimator

In Section 2.1, we describe a data generation process associated with the
perturbation model in Fig 1. In Section 2.2, we propose _DirectLikelihood_ , a
maximum-likelihood estimator with respect to the marginal distribution of the
observed variables. _DirectLikelihood_ identifies estimates of the unknown
perturbation effects, the latent effects, and the causal relation among the
observed variables.

###  2.1 Modeling framework

We consider a directed acyclic graph $\mathcal{D}^{\star}$ whose $p+h$ nodes
correspond to jointly Gaussian and centered 111Without loss of generality, we
assume that the observed variables are centered. random variables
$(X,H)\subseteq\mathbb{R}^{p}\times\mathbb{R}^{h}$, where $X$ are observable
and $H$ are latent variables. As described in Section 1.1, our methodology is
also applicable in the setting where one may be primarily interested in the
causal effects of a response variable. As such, we distinguish $X_{p}$ as the
target or response variable. Owing to the joint Gaussianity of $(X,H)$, the
random pair $(X,H)$ satisfies the following (compactified) SCM:

$X=B^{\star}{X}+\Gamma^{\star}{H}+\epsilon.$ (1)

Here, $\epsilon=(\epsilon_{1},\epsilon_{2},\dots,\epsilon_{p})$ are
independent Gaussian random variables independent of $H$ where
$\epsilon\sim\mathcal{N}(0,\texttt{diag}(w^{1,\star}))$ for some
$w^{1,\star}\in\mathbb{R}^{p}_{++}$. The connectivity matrix
$B^{\star}\in\mathbb{R}^{p\times p}$ contains zeros on the diagonal and
encodes the causal relationship among the observables $X$, i.e.
$B^{\star}_{k,j}\neq 0$ if and only if
$j\in\text{PA}_{\mathcal{D}^{\star}_{X}}(k)$. The $p$-th row vector
$B^{\star}_{p,:}$ encodes the causal parents of the response variable and the
magnitude of their effects. The matrix $\Gamma^{\star}$ in (1) encodes the
effects of the latent variables on the observed variables where
$\Gamma^{\star}_{k,j}\neq 0$ if and only if
$j\in\text{PA}_{\mathcal{D}^{\star}_{H}}(k)$. For the sake of generality, we
do not immediately put any assumption on the number of latent variables $h$ or
the denseness of their effects.

The compact SCM (1) describes the generating process of $X$ in the
observational setting where there are no external perturbations on the system.
We next describe how the data generation process alters due to some type of
perturbations to the variables $(X,H)$. We consider perturbations that
directly shift the distributions of the random variables by some noise acting
additively to the system. Specifically, the perturbations $\mathcal{A}$
generate the random pair $(X^{e},H^{e})$ for each $e\in\mathcal{E}$ satisfying
the following SCM:

$\begin{gathered}X^{e}=B^{\star}{X^{e}}+\Gamma^{\star}{H}^{e}+\epsilon^{e}+\delta^{e}\\\
H^{e}\sim\mathcal{N}(0,\Psi^{\star,e}),\end{gathered}$ (2)

where for every $e\in\mathcal{E}$,
$\epsilon^{e}\stackrel{{\scriptstyle\text{dist}}}{{=}}\epsilon$,
$(H^{e},\delta^{e},\epsilon^{e})$ are jointly independent, and that the
collection $(X^{e},H^{e},\delta^{e},\epsilon^{e})$ is independent across $e$.
Further, $\delta^{e}\in\mathbb{R}^{p}$ is a Gaussian random vector
(independent across the coordinates) that represents the additive
perturbations, and $H^{e}\in\mathbb{R}^{h}$ is a Gaussian random vector that
represents the perturbed latent variables with covariance
$\Psi^{e}\in\mathbb{S}^{h}_{++}$. Notice that $\epsilon^{e},\delta^{e}$ are in
general not identifiable from the sum $\epsilon^{e}+\delta^{e}$ in (2); we
specify below in Section 2.1.1 an identifiable parametrization for the terms
$\epsilon^{e}+\delta^{e}$. The modeling framework (2) can also incorporate
information about do-interventions, as discussed in Section 2.1.3.

The compactified SCM (2) characterizes the distribution among all of the
observed variables and encodes system-wide invariances. Specifically, (2)
insists that for every $k=1,2,\dots,p$, the regression coefficients when
regressing $X^{e}_{k}$ on the parent sets
$\\{X_{j}^{e}:j\in\text{PA}_{\mathcal{D}_{X}^{\star}}(k)\\}$ and
$\\{H_{l}^{e}:l\in\text{PA}_{\mathcal{D}_{H}^{\star}}(k)\\}$ remain invariant
for all environments $e\in\mathcal{E}$. This is a point of departure from
instrumental variable techniques or invariant causal prediction in two
significant ways: 1) such methods do not allow for the perturbations on the
latent variables or the response variable $X_{p}$ (i.e. they assume
$H^{e}\stackrel{{\scriptstyle\text{dist}}}{{=}}H$ and $\delta^{e}_{p}\equiv 0$
for all $e\in\mathcal{E}$) and 2) they only consider “local” invariances
arising from the distribution
$X^{e}_{p}~{}|~{}\\{(X_{j}^{e},H_{l}^{e}):j\in\text{PA}_{\mathcal{D}_{X}^{\star}}(p),l\in\text{PA}_{\mathcal{D}_{H}^{\star}}(p)\\}$.
The virtue of considering a joint model over all of the variables and
exploiting system-wide invariances is that we can propose a maximum-likelihood
estimator _DirectLikelihood_ which identifies the population DAG structure
even with perturbations on the response variable and the latent variables.

The SCM (2) is similar in spirit to previous modeling frameworks in the
literature. The authors [15] consider jointly observational and interventional
Gaussian data where the interventions are limited to do-interventions and
there are no latent variables. In the context of (2), this means that
$\delta^{e}\equiv 0$ and $\Gamma^{\star}\equiv 0$. As such, the framework
considered in this paper is a substantial generalization of [15]. Further, the
backShift [27] procedure considers the linear SCM (2) with the some
modifications: i) there are no do-interventions, ii) there are no
perturbations to the latent variables, i.e.
$H^{e}\stackrel{{\scriptstyle\text{dist}}}{{=}}H$ for all $e\in\mathcal{E}$,
and iii) $B^{\star}$ may be a cyclic directed graph. In addition, the
backShift algorithm relies on exploiting invariances of differences of
estimated covariance matrices across environments. Our _DirectLikelihood_
procedure is more in the ”culture of likelihood modeling and inference” and
has the advantage that it can cope well with having only a few observations
per group or environment. This likelihood perspective also fits much more into
the context of inference for mixed models as briefly discussed in Section
2.1.2.

####  2.1.1 Model specialization

Clearly, one cannot distinguish the parameters for $\epsilon^{e}$ and
$\delta^{e}$. We thus write, for all $e\in\mathcal{E}$:

$\epsilon^{e}+\delta^{e}\sim{\cal N}(0,\mathrm{diag}(w^{e,\star})),\
w^{e,\star}\in\mathbb{R}_{++}^{p}.$

Since we are mainly interested in the connectivity matrix $B^{\star}$, the
parameters $\Gamma^{\star}$, $w^{e,\star},\Psi^{e,\star}$ are nuisance
parameters and we may simplify the modeling framework by restricting the
parameter space for the covariances $\Psi^{e,*}$.. Our default proposal is to
model the latent variables as independent and identically distributed across
the environments. Specifically, we let
$\Psi^{e,\star}=\Psi^{\star}+\psi^{e,\star}\mathcal{I}$ where
$\psi^{e,\star}\in\mathbb{R}_{+}$. Further, without loss of generality,
$\Psi^{\star}$ can be taken to be the identity matrix by absorbing its effect
on $\Gamma^{\star}$ via the transformation
$\Gamma^{\star}\to\Gamma^{\star}{\Psi^{\star}}^{1/2}$ so that:

$\Psi^{e,\star}=(1+\psi^{e,\star})\mathcal{I},\
\psi^{e,\star}\in\mathbb{R}_{+}.$

Further, as an additional default setting, we assume that we have access to an
observational environment ($e=1$ without loss of generality) so that:

$w^{e,\star}\succeq w^{1,\star},\ \psi^{1,\star}=0.$

Here, the inequality $w^{e,\star}\succeq w^{1,\star}$ is element-wise.

In the setting where the latent variables are unperturbed across the
environments, one can take $\Psi^{e,\star}\equiv\mathcal{I}$ after the
transformation $\Gamma(\Psi^{e,\star})^{1/2}\to\Gamma^{\star}$. Fitting to a
model with equally distributed perturbations across the coordinates may be
attained by the reparametrization $w^{e,\star}=w^{1,\star}+\zeta^{e,\star}{\bf
1}$ for $\zeta^{e,\star}\in\mathbb{R}_{+}$. In general, other models for the
random terms $\epsilon^{e}+\delta^{e}$ and $H^{e}$ are possible. A connection
to random effects modeling is discussed next.

####  2.1.2 Interpretation as mixed-effects linear structural causal model

The framework in (2) bears some similarities to standard random effects mixed
models [16]. In particular, random effects mixed models are widely employed to
model grouped data, where some parameter components remain fixed and others
are random. In the context of our problem, the fixed parameters are the
matrices $B^{\star},\Gamma^{\star}$ and the random parameters are the shift
perturbations $\delta^{e}$.

For example, we can write for the response variable $Y=X_{p}$ and for
simplicity in the absence of latent variables: for each environment or group
$e$,

$Y^{e}=X^{e}\beta+Z^{e}b^{e}+\epsilon^{e},\ e=1,\ldots,m,$ (3)

where $Y^{e},\epsilon^{e}$ are $n^{e}\times 1$ vectors, $X^{e}$ is an
$n^{e}\times p$ design matrix, here $Z^{e}=\mathcal{I}_{n^{e}}$, $n^{e}$ is
the sample size within group or environment $e$, and the variables across $e$
are independent. The correspondence to (2) is as follows:
$\epsilon^{e}\sim{\cal N}(0,w^{1,\star}_{p}\mathcal{I}_{n^{e}})$,
$b^{e}\sim{\cal N}(0,v_{p}^{e,\star}\mathcal{I}_{n^{e}})$ (where
$v_{p}^{e,\star}=w_{p}^{e,\star}-w_{p}^{1,\star}$) and
$\beta={B^{\star}_{p,:}}^{T}$. There are three main differences to standard
mixed models. First, the distribution of $b^{e}\sim{\cal
N}(0,v_{p}^{e,\star}I_{n^{e}})$ changes with $e$ and the shrinkage effect
across groups is abandoned. Second, we take a multivariate view point for all
the variables $X_{j}^{e}\ (j=1,\ldots,p)$ in (2): they are all modelled with
random effects and can be individually written as in (3), but we allow for
dependence among all the $p$ variables. Finally, a difference between our
model in (1) and (2) and the standard mixed models is that the group specific
random effects, the random parameters $\delta^{e}$ in (2) or the random
parameters vector $b^{e}$ in (3), act in a _dynamic way_ on the system: the
effects of $\delta^{e}$ are _propagated_ through the structural equations; and
in practice, the order of propagation is usually unknown.

Thus, our model in (2) leads to a different way of describing group-specific
perturbations, calling also for a different likelihood calculation: in fact,
as we show, such dynamic perturbations allow to identify the causal structure.
The latter is not possible with standard mixed models but due to the
connection pointed out above, we refer to our formalization in (2) as ”mixed-
effects linear structural causal modeling”. We believe that the causal
inference literature has not much exploited this connection. We argue here
that our random effects approach is very useful and practical for modeling
perturbation data where the perturbations are believed to propagate further on
other variables in the system.

####  2.1.3 Incorporating do-interventions

The perturbation model (2) provides a flexible framework to incorporate
additional knowledge including do-interventions (eliminating the connections
between the perturbed variable and the corresponding parents) when their
intervention-locations are known. Specifically, in such setting, (2) can be
modified to:

$\begin{gathered}X^{e}=\mathcal{F}_{\texttt{do}(e)^{c}}(B^{\star}{X^{e}}+\Gamma^{\star}{H}^{e}+\epsilon^{e})+\delta^{e}\\\
H^{e}\sim\mathcal{N}(0,\Psi^{\star,e}),\end{gathered}$

where $\texttt{do}(e)\subseteq\\{1,\ldots,p\\}$ denotes do-locations in the
sub-graph of $\mathcal{D}^{\star}_{X}$ and
$\mathcal{F}_{S}\in\mathbb{R}^{p\times p}$ is a diagonal matrix with ones
corresponding to coordinates inside $S\subseteq\\{1,\ldots,p\\}$ and zeros
elsewhere. Accordingly, the _DirectLikelihood_ procedure described in Section
2.2 can be modified; see Section A in the supplementary material.

###  2.2 Scoring DAGs via _DirectLikelihood_

Let $\mathcal{D}$ be a given DAG structure among the observed variables (which
we can think of as the restriction $\mathcal{D}_{X}$ of a DAG among observed
and latent variables). In this section, we score this DAG via the maximum
likelihood procedure _DirectLikelihood_. We suppose that there are $m$
environments $|\mathcal{E}|=m$, and for every environment $e=1,2,\dots,m$, we
have samples of random pairs $(X^{e},H^{e})$:
$\\{X^{e}_{(i)}\\}_{i=1}^{n^{e}}$ for some positive integer $n^{e}$ which are
IID for each $e$ and independent across $e$. Thus, since the $X^{e}$’s are
independent for $e=1,2,\dots,m$ and the samples for each environment $e$ are
are IID, the maximum-likelihood estimator for the DAG structure $\mathcal{D}$
is given by:

$\displaystyle\operatornamewithlimits{arg\,min}_{\begin{subarray}{c}B\in\mathbb{R}^{p\times
p},\Gamma\in\mathbb{R}^{p\times\bar{h}}\\\
\\{\Psi^{e}\\}_{e=1}^{m}\subseteq\mathbb{S}^{\bar{h}}_{++},\\{w^{e}\\}_{e=1}^{m}\subseteq\mathbb{R}^{p}_{++}\end{subarray}}$
$\displaystyle\sum_{e=1}^{m}\hat{\pi}^{e}\sum_{i=1}^{n^{e}}-\log\texttt{prob}\left(X^{e}_{(i)}|B,\Gamma,\Psi^{e},w^{e}\right)$
(4) subject-to $\displaystyle~{}~{}~{}B\text{ compatible with }\mathcal{D}.$

Here, $\texttt{prob}\left(X^{e}_{(i)}|B,\Gamma,\Psi^{e},w^{e}\right)$
represents the Gaussian likelihood of $X^{e}_{(i)}$ given parameters
$B,\Gamma,\Psi^{e},w^{e}$; $\bar{h}\leq p$; the constraint $B\text{ compatible
with }\mathcal{D}$ ensures that the estimated $B$ has its support restricted
to the structure of $\mathcal{D}$, i.e. $B_{i,j}\neq 0$ if and only if $j\to
i$ in $\mathcal{D}$; $\hat{\pi}^{e}=\frac{n^{e}}{\sum_{e=1}^{m}n^{e}}$; and
$w^{e}$ is a surrogate for the variances of the sum $\delta^{e}+\epsilon$. The
maximum-likelihood estimator (4) can be rewritten as:

$\displaystyle(\hat{{B}},\hat{{\Gamma}},\\{(\hat{{\Psi}}^{e},\hat{{w}}^{e})\\}_{e=1}^{m})=\operatornamewithlimits{arg\,min}_{\begin{subarray}{c}B\in\mathbb{R}^{p\times
p},\Gamma\in\mathbb{R}^{p\times\bar{h}}\\\
\\{\Psi^{e}\\}_{e=1}^{m}\subseteq\mathbb{S}^{\bar{h}}_{++},\\{w^{e}\\}_{e=1}^{m}\subseteq\mathbb{R}^{p}_{++}\end{subarray}}$
$\displaystyle\sum_{e=1}^{m}\hat{\pi}^{e}\ell(B,\Gamma,\Psi^{e},w^{e};\hat{\Sigma}^{e}),$
(5) subject-to $\displaystyle~{}~{}~{}B\text{ compatible with }\mathcal{D}$

where

$\displaystyle\ell(\cdot)=\log\det\left(\texttt{diag}(w^{e})+\Gamma\Psi^{e}\Gamma^{T}\right)+\mathrm{trace}\left(\left[\texttt{diag}(w^{e})+\Gamma\Psi^{e}\Gamma^{T}\right]^{-1}(\mathcal{I}-B)\hat{\Sigma}^{e}(\mathcal{I}-B)^{T}\right),$

and $\hat{\Sigma}^{e}$ is the sample covariance matrix of the data
$\\{X^{e}_{(i)}\\}_{i=1}^{n^{e}}$. The input to the program (5) are the sample
covariance matrices $\hat{\Sigma}^{e}$ and the estimate $\bar{h}$ for the
number of latent variables. We note that the _DirectLikelihood_ estimator can
be specialized to different modeling options based on appropriate
reparametrization of the nuisance parameters $\Psi^{e},w^{e}$ in (5). For
example, in our default setting of IID latent variables with the environment
$e=1$ being observational (see Section 2.1.1), we add the following
constraints to (5):

$\displaystyle\Psi^{e}=(1+\psi^{e})\mathcal{I}\text{ with
}\psi^{e}\in\mathbb{R}_{+}\text{ for }e=1,\dots,m$ $\displaystyle w^{e}\succeq
w^{1}\text{ for }e=2,\dots,m;~{}\psi^{1}=0.$

Given estimates
$(\hat{{B}},\hat{{\Gamma}},\\{\hat{{\Psi}}^{e}\\}_{e=1}^{m},\\{\hat{{w}}^{e}\\}_{e=1}^{m})$,
our score for the DAG $\mathcal{D}$ is:

$\text{score}_{\lambda}(\mathcal{D})=\sum_{e=1}^{m}\hat{\pi}^{e}\ell(\hat{{B}},\hat{{\Gamma}},\hat{{\Psi}}^{e},\hat{{w}}^{e};\hat{\Sigma}^{e})+\lambda\|\text{moral}(\mathcal{D})\|_{\ell_{0}}.$
(6)

Here, $\text{moral}(\mathcal{D})$ denotes the moralization of $\mathcal{D}$
which forms an undirected graph of $\mathcal{D}$ by adding edges between nodes
that have commons children, $\lambda>0$ is a regularization parameter, and
$\lambda\|\text{moral}(\mathcal{D})\|_{\ell_{0}}$ is akin to the Bayesian
Information Criterion (BIC) score that prevents overfitting by incorporating
the denseness of the moral graph of $\mathcal{D}$ in the likelihood score. In
principle, a collection of DAGs can each be individually scored via (6) to
find the best fit to the data. We remark that regularization terms controlling
for complexity of estimated DAGs are commonly employed in structural causal
learning (see [10] and the references therein). A classically known fact is
that in a single environment setting, the moral graph of the DAGs in the
Markov equivalence class have the same cardinality [33]. In the context of
this paper with perturbations, incorporating the sparsity of the moral graph
plays a central role in our theoretical analysis for proving identifiability.

In comparison to the _DirectLikelihood_ procedure, backShift [27] fits the SCM
(2) (with some restrictions outlined in Section 2.1) by performing joint
diagonalization to the difference of sample covariance matrices.
_DirectLikelihood_ allows for much more modeling flexibility. First, in
contrast to backShift where the latent effects are subtracted by computing the
difference of covariances, _DirectLikelihood_ explicitly models these effects.
This feature of _DirectLikelihood_ enables the possibility of perturbations to
the latent variables and a manner to control the number of estimated latent
variables (as opposed to arbitrary number of latent variables with backShift).
We discuss in Section 3 that controlling the number of latent variables may
lead to identifiability using _DirectLikelihood_ with a single interventional
environment, whereas backShift is guaranteed to fail. Second,
_DirectLikelihood_ also models the perturbation magnitudes in each
environment, allowing for the flexibility of constraining the perturbation
magnitudes to improve estimation accuracy. Finally, _DirectLikelihood_ allows
to pool information over different environments $e$ for the parameter $B$ of
interest: this enable _DirectLikelihood_ to be used with only a few sample
points per environment.

###  2.3 Beyond Gaussianity

The _DirectLikelihood_ estimator (5) fits a Gaussian perturbation model (2) to
the data. However, the perturbation data of the observed variables may be non-
Gaussian but satisfy the linear SCM (2). In particular, the random variables
$H^{e},\delta^{e}$ may be non-Gaussian while still inducing a linear
relationship with the observed variables $X^{e}$. Nonetheless, since the
_DirectLikelihood_ estimator only operates on second moments, one may still
use the _DirectLikelihood_ procedure to find the best scoring DAGs and the
associated connectivity matrices without compromising identifiability
guarantees as shown in Section 3, still implying corresponding estimation
consistency. Further, we empirically explore the robustness of the
_DirectLikelihood_ procedure to non-Gaussianity as well as other model
misspecifications via numerical experiments in supplementary material Section
G.

##  3 Theoretical properties: identifiability via _DirectLikelihood_

We next investigate the theoretical properties of the _DirectLikelihood_
procedure. The main theorem in this section (Theorem 1) considers the general
setting with perturbed latent variables and establishes identifiability
properties under some population assumptions. Subsequently, Theorem 2 in
Section 3.1 analyze _DirectLikelihood_ under the specialization that the
latent variables are unperturbed. Throughout, the notation with ∗ indicates
the true underlying population objects which we aim to estimate from data.

Setup: We consider the perturbation model in (2) with a population
connectivity matrix $B^{\star}\in\mathbb{R}^{p\times p}$ and latent effects
coefficient matrix $\Gamma^{\star}\in\mathbb{R}^{p\times h}$. For every
environment $e$, the random vector $H^{e}$ has a covariance matrix
$\Psi^{e,\star}\in\mathbb{S}^{h}_{++}$ and the random vector
$\epsilon^{e}+\delta^{e}$ has a diagonal covariance matrix
$\texttt{diag}(w^{e,\star})$ for $w^{e,\star}\in\mathbb{R}^{p}_{+}$. In the
subsequent discussion, we allow for $H^{e},\delta^{e},$ and $\epsilon$ to be
non-Gaussian random vectors. As prescribed in Section 2.1.1 but not requiring
Gaussianity, we assume that the latent variables are independent and
identically distributed, i.e. $\Psi^{e,\star}=(1+\psi^{e,\star})\mathcal{I}$
with $\psi^{e,\star}\in\mathbb{R}_{+}$, and that for every environment
$e=1,2,\dots,m$, we have IID data
$\\{X^{e}_{(i)}\\}_{i=1}^{n^{e}}\subseteq\mathbb{R}^{p}$ where $e=1$ is an
observational environment (our theoretical results can be extended to the
settings with non-IID latent variables and perturbations in every
environment). To score a given DAG $\mathcal{D}$, we consider the modified
_DirectLikelihood_ estimator (5) in population:

$\displaystyle\min_{\begin{subarray}{c}B\in\mathbb{R}^{p\times
p},\Gamma\in\mathbb{R}^{p\times\bar{h}}\\\
\\{(\psi^{e},w^{e})\\}_{e=1}^{m}\subseteq\mathbb{R}_{+}\times\mathbb{R}^{p}_{++}\end{subarray}}$
$\displaystyle\sum_{e=1}^{m}\pi^{e,\star}\ell(B,\Gamma,(1+\psi^{e})\mathcal{I},w^{e};{\Sigma}^{e,\star}).$
(7) subject-to $\displaystyle~{}~{}~{}B\text{ compatible with
}\mathcal{D}~{}~{};\psi^{e}\leq C_{\psi}\text{ for }e=1,\ldots,m$
$\displaystyle~{}~{}~{}\psi^{1}=0~{};~{}w^{e}\succeq w^{1}\text{ for
}e=2,\dots,m.$

Comparing (7) to the _DirectLikelihood_ estimator (5), the reparametrization
$\Psi^{e}\to(1+\psi^{e})\mathcal{I}$ is to account for the latent variables
being IID and the constraints $\psi^{1}=0$ and $w^{e}\succeq w^{1}$ for
$e=2,\dots,m$ account for $e=1$ being an observational environment. Further,
the constraint $\|\psi\|_{\infty}\leq C_{\psi}$ bounds the strength of the
latent perturbations for some user-specified parameters $C_{\psi}\geq 0$.

We consider optimally scoring DAG(s) with their associated connectivity
matrices:

$\displaystyle{\mathcal{D}}_{\text{opt}}=\operatornamewithlimits{arg\,min}_{\text{DAG
}\mathcal{D}}\text{score}_{\lambda=0}(\mathcal{D})~{}~{};~{}~{}{B}_{\text{opt}}:\text{
associated connectivity matrix(ces)}.$ (8)

Here, $\text{score}_{\lambda=0}(\mathcal{D})$ is the achieved minimum in (7).
It is the analogue of (6) but using the population covariance matrix
$\Sigma^{e,*}$ and the population optimizers from (7). The sample
_DirectLikelihood_ procedure replaces ${\Sigma}^{e,\star}$ and $\pi^{e,\star}$
in (7) with the population covariance matrix $\hat{\Sigma}^{e}$ and the
mixture coefficients $\hat{\pi}^{e}$, respectively. Further, in the sample
setting, the regularization parameters $\lambda$ in the score evaluation (6)
must be tuned. Using the sample quantities as described above we denote by

$\displaystyle\hat{\mathcal{D}}_{\text{opt}},\ \ \hat{B}_{\text{opt}}$ (9)

the optimal scoring DAGs and connectivity matrices in the sample version. Our
objective is to demonstrate that under some assumptions,
$\mathcal{D}_{\text{opt}}=\mathcal{D}_{X}^{\star}$,
$B_{\text{opt}}=B^{\star}$, and in the limit of sample sizes for all
environments tending to infinity,
$\hat{\mathcal{D}}_{\text{opt}}\to\mathcal{D}_{X}^{\star}$ and
$\hat{B}_{\text{opt}}\to B^{\star}$ in probability for an appropriate choice
of $\lambda$.

Our consistency results are in the general setting where there are
perturbations on the latent variables and require an assumption on the latent
variable effects, dubbed _latent materiality_ , that is formalized below:

######  Definition 1 (_latent materiality_ for $e\in\mathcal{E}$).

The random variables $(X^{e},H^{e})$ satisfy _latent materiality_ if there
exists a pair $k,l$ such that:

$\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash\\{k,l\\}}^{e},H^{e})=0~{}\&~{}\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash\\{k,l\\}}^{e})\neq
0.$

In words, (1) states that the latent variables induce “some” confounding
dependencies among the observed variables in environment $e\in\mathcal{E}$. In
comparison, the latent denseness assumption needed for consistency of the two
stage deconfounding procedures [3, 12] require that the latent variables
induce “many” confounding dependencies. As such, _latent materiality_ is a
strictly (and often substantially) weaker condition than the denseness
assumption required for the success of the two stage deconfounding. We
investigate whether _DirectLikelihood_ is able to identify the population
connectivity matrix $B^{\star}$ under this weaker condition, and answer in the
affirmative under appropriate conditions on the strength and heterogeneity of
the interventions. We provide two sets of assumptions that lead to
identifiability. The first set requires two interventional environments that
have sufficiently large interventions on the observed variables, and the
second set requires two interventional environments with one interventional
environment consisting of much stronger interventions on the observed
variables than the other interventional environment. These assumptions are
described below where the observational environment is denoted by $e=1$ and
the two interventional environments are denoted by $e=2,3$:

Assumption 1 $\displaystyle-\text{the mixture effects are}\text{ \allowbreak
non-vanishing: }\pi^{e,\star}>0\text{ for }e=1,2,3$ (10) Assumption 2
$\displaystyle-\text{heterogeneity of perturbations for }e=2,3:~{}$
$\displaystyle~{}~{}~{}~{}\frac{w^{2,\star}_{k}-(1+\psi^{2,\star})w^{1,\star}_{k}}{w^{2,\star}_{l}-(1+\psi^{2,\star})w^{1,\star}_{l}}\neq\frac{w^{3,\star}_{k}-(1+\psi^{3,\star})w^{1,\star}_{k}}{w^{3,\star}_{l}-(1+\psi^{3,\star})w^{1,\star}_{l}}\text{
for all }k\neq{l}$ Assumption 3 $\displaystyle-\text{\emph{latent materiality
}}\text{ in Definition }\ref{defn:hidden_spurious}\text{ for environments
}e=2,3$ Assumption 4 $\displaystyle-\text{sufficiently large interventions on
variables for $e=2,3$: }$
$\displaystyle~{}\frac{\min_{k}(w^{e,\star}_{k})^{2}}{\max_{k}w^{e,\star}_{k}}>8{\kappa^{\star}}(1+2C_{\psi})^{2}(1+\|w^{1,\star}\|_{\infty})(1+\|\Gamma^{\star}\|_{2}^{2}+\|\Gamma^{\star}\|_{2}^{4})$
Assumption 2’ $\displaystyle-\text{heterogeneity of perturbations for }e=2,3:$
$\displaystyle~{}~{}~{}~{}\frac{w^{3,\star}_{k}-(1+\psi^{3,\star})w^{1,\star}_{k}}{w^{3,\star}_{l}-(1+\psi^{3,\star})w^{1,\star}_{l}}\neq\frac{w^{3,\star}_{k}-\frac{1+\psi^{3,\star}}{1+\psi^{2,\star}}w^{2,\star}_{k}}{w^{3,\star}_{l}-\frac{1+\psi^{3,\star}}{1+\psi^{2,\star}}w^{2,\star}_{l}}\text{
for all }k\neq{l}$ Assumption 3’ $\displaystyle-\text{\emph{latent materiality
}}\text{ in Definition }\ref{defn:hidden_spurious}\text{ for environments
}e=3$ Assumption 4’ $\displaystyle-\text{sufficiently large interventions on
variables in $S$ for $e=3$: }$
$\displaystyle~{}\frac{\min_{k}(w^{e,\star}_{k})^{2}}{\max_{k}w^{e,\star}_{k}}>8{\kappa^{\star}}(1+2C_{\psi})^{2}(1+\|w^{2,\star}\|_{\infty})(1+\|\Gamma^{\star}\|_{2}^{2}+\|\Gamma^{\star}\|_{2}^{4})$

where the quantity
$\kappa^{\star}\equiv\frac{1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2}}{1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2}}$
in Assumption 4 or 4’ of (10). Assumptions 1-4 or 1 $\&$ 2’-4’ in (10) impose
conditions on the population quantities associated with the environments
$e=1,2,3$. In particular, Assumption 1 in (10) require that the contribution
for each environment does not vanish in the large data limit; Assumptions 2
and 2’ in (10) ensure that the perturbations are heterogeneous. In principle,
the interventions on the observed variables in the environments $e=2,3$ may
come from identical distributions (i.e. $w^{2,\star}=w^{3,\star}$) or one of
them being even zero (i.e. $w^{2,\star}=w^{1,\star}$) with different latent
variable perturbations (i.e. $\psi^{2,\star}\neq\psi^{3,\star}$) without
compromising Assumption 2 or 2’ in (10). Additionally, one can show that if
the parameters $w^{3,\star},w^{2,\star},w^{1,\star}$ and
$\psi^{2,\star},\psi^{3,\star}$ are drawn from continuous distributions,
Assumption 2 and 2’ in (10) are satisfied almost surely. Assumptions 3 and 3’
in (10) insists that the _latent materiality_ in (1) is satisfied so that the
latent variables induce at least a single spurious dependency among the
observed variables. Finally, Assumptions 4 and 4’ in (10) require that the
perturbations on the observed variables are sufficiently large. This is akin
to strong instruments assumption in the instrumental variables literature [1].

Given Assumptions 1-4 or Assumptions 1 $\&$ 2’-4’ in (10), we first analyze
the theoretical properties of the population _DirectLikelihood_ procedure.

######  Theorem 1 (Identifiability in population: perturbed latent variables).

Suppose that the user-specified parameters $\bar{h}$ and $C_{\psi}$ in (7) are
chosen conservatively so that $\bar{h}\geq\text{dim}(H)$ and
$C_{\psi}\geq\psi^{e,\star}$ for all $e=1,2,\dots,m$. Under Assumptions 1-4 or
Assumptions 1 $\&$ 2’-4’ in (10), the following are satisfied for
$\mathcal{D}_{\text{opt}}$ in (8):

  1. 1.

$\mathcal{D}^{\star}_{X}\in{\mathcal{D}}_{\text{opt}}$ and any other optimum
$\mathcal{D}\in{\mathcal{D}}_{\text{opt}}$ satisfies:
$\text{moral}(\mathcal{D}^{\star}_{X})\subseteq\text{moral}(\mathcal{D})$.

  2. 2.

The optimum of $\arg\min_{D\in{D}_{\text{opt}}}\|\text{moral}(D)\|_{\ell_{0}}$
is unique and equal to $D^{\star}_{X}$. Further, the associated connectivity
matrix is equal to $B^{\star}$.

The proof is presented in the supplementary material Section B. The first
assertion in Theorem 1 states that the moral graph of any optimum
$\mathcal{D}\in\mathcal{D}_{\text{opt}}$ of the _DirectLikelihood_ procedure
is a superset of the moral graph of $\mathcal{D}^{\star}$, and the second
assertion states that the connectivity matrices yielding the sparsest moral
graphs among the optima are unique and equal to $B^{\star}$. These statements
do not guarantee recovering the other model parameters , viewed here as
nuisance part, including $\Gamma^{\star}$, and
$\\{(\psi^{e,\star},w^{e,\star})\\}_{e=1}^{m}$. However, under additional
assumptions namely: $\bar{h}=\text{dim}(H)$ and the incoherence of the
subspace $\text{col-space}(\Gamma^{\star})$, recovery of
$\Gamma^{\star}{\Gamma^{\star}}^{T}$ and
$\\{(\psi^{e,\star},w^{e,\star})\\}_{e=1}^{m}$ can be shown.

We note that Assumptions 1-4 or Assumptions 1 $\&$ 2’-4’ in (10) are
sufficient conditions for identifiability and are generally not necessary. As
an example, we show in supplementary material Section C that identifiability
cannot be achieved with only a single interventional environment if
$\bar{h}=p$ (e.g. most conservative choice for the number of latent
variables). However, we also demonstrate that if $\bar{h}<p$,
_DirectLikelihood_ will attain identifiability under certain configurations of
model parameters (i.e dense latent effects with sparse population DAG
$\mathcal{D}_{X}^{\star}$). Thus, Assumptions 1-4 or Assumptions 1 $\&$ 2’-4’
in (10) serve as protection for arbitrary population DAG structure and a class
of model parameters. We believe that relaxing these assumptions while
retaining identifiability guarantees is an interesting direction for future
research.

The virtue of incorporating the regularization term
$\lambda\|\mathcal{D}\|_{\ell_{0}}$ in (6) is that in the large data limit,
this penalty term encourages sparser moral graphs. Thus, in conjunction with
the results of Theorem 1, we demonstrate that in the large data limit, the set
$\hat{\mathcal{D}}_{\text{opt}}$ and $\hat{B}_{\text{opt}}$ asymptotically
converge to $\mathcal{D}^{\star}_{X}$ and $B^{\star}$, respectively. To appeal
to standard empirical process theory results, we constrain the parameter space
to be compact as described in the corollary below:

######  Corollary 1 (Asymptotic consistency for perturbed latent variables).

Consider the sample version of the _DirectLikelihood_ procedure in (7) with
the compactness constraints $\max\\{1/{\min_{k}w^{e}_{k}},\|B\|_{2}\\}\leq
C_{\text{comp}}$ for every every $e=1,2,\dots,m$ where
$C_{\text{comp}}>\max\\{{1}/{\min_{k}w^{e,\star}_{k}},\allowbreak\|B^{\star}\|_{2}\\}$
so that the true parameters are in the feasible region. Further, let
$\lambda\sim\mathcal{O}({\log(\sum_{e=1}^{m}n^{e})}\allowbreak/{\sum_{e=1}^{m}n^{e})}$
in (6). Under the conditions in Theorem 1, the following are satisfied for
$\hat{\mathcal{D}}_{\text{opt}}$ and $\hat{B}_{\text{opt}}$ in (9):
$\hat{\mathcal{D}}_{\text{opt}}\to\mathcal{D}^{\star}_{X}$ and
${B}_{\text{opt}}\to B^{\star}$ , in probability, as $n^{e}\to\infty$ for
$e=1,2,3$.

The proof of Corollary 1 is a straightforward consequence of Theorem 1 and
left out for brevity. The combined results of Theorem 1 and Corollary 1 state
that under perturbations that are sufficiently different across environments
and the _latent materiality_ condition, two interventional environments
suffice for consistent estimation.

_R_ emark 1: Assumptions 1-4 or Assumptions 1 $\&$ 2’-4’ in (10) needed for
identifiability suggest that perturbations on the latent variables can improve
identifiability. Specifically, the perturbations on the observed variables in
one interventional environment may be statistically identical to another
environment or even be completely equal to zero and still preserve
identifiability as long as the latent variables have been perturbed.

_R_ emark 2: As described in Section 2.1, the perturbation model (2) offers
flexibility with respect to many components of the model such as the structure
of the perturbations on the observed or latent variables. In particular, one
may fit to data the perturbation model (2) where the perturbation magnitudes
are equal in magnitude across the coordinates, e.g.
$\texttt{diag}(w^{e,\star})\propto\mathcal{I}$. We demonstrate in the
supplementary material Section D that _DirectLikelihood_ , under Assumptions
similar to (10), provides consistent estimators in this setting. Thus, in
principle one may have only two additional perturbation parameters per
environment: a scalar for the latent variables and a scalar for the observed
variables. As a point of contrast, in the setting where the perturbations
among the observed variables may vary, there are $p+1$ new variables for each
environment. The substantial reduction in the number of parameters can lead to
better statistical properties in practice.

###  3.1 Specializations: unperturbed latent variables

We next analyze the identifiability guarantees of the _DirectLikelihood_
procedure when the latent variables remain unperturbed across the
environments, i.e. the perturbation $\mathcal{A}$ does not point to $H$ in Fig
1. Specifically, we consider the setup described in the beginning of Section 3
with the modification that $\psi^{e,\star}=0$. Thus, we also modify the
_DirectLikelihood_ estimator (7) by setting $\psi\equiv 0$. We further
consider an arbitrary latent effects matrix $\Gamma^{\star}$, where the two-
stage deconfounding procedure will not perform well, since latent denseness
may not be satisfied. We demonstrate on the other hand, that the under
sufficient interventions, the connectivity matrix that attains the optimum
score via _DirectLikelihood_ in population is unique and equal to $B^{\star}$.

######  Theorem 2 (Identifiability in population: unperturbed latent
variables).

Let $\bar{h}\geq\text{dim}(H)$ in the _DirectLikelihood_ estimator (7).
Letting $S\subseteq\\{1,2,\dots,p\\}$ encode the location of perturbations,
suppose that Assumption 2 in (10) is modified to
$\frac{w^{2,\star}_{k}-w^{1,\star}_{k}}{w^{2,\star}_{l}-w^{1,\star}_{l}}\neq\frac{w^{3,\star}_{k}-w^{1,\star}_{k}}{w^{3,\star}_{l}-w^{1,\star}_{l}}$
for all $k,l\in S,k\neq l$ and Assumption 4 in (10) is modified to
$w^{e,\star}_{k}>w^{1,\star}_{k}$ for $e=2,3$ and all $k\in S$. Then, under
Assumptions 1 in (10) and modified Assumptions 2 and 4, we have for
$\mathcal{D}_{\text{opt}}$ and $B_{\text{opt}}$ as in (8):

  1. (a)

${\mathcal{D}}_{\text{opt}}=\mathcal{D}^{\star}_{X}$ and
${B}_{\text{opt}}=B^{\star}$ if ${S}=\\{1,\ldots,p\\}$.

  2. (b)

Any optimum $B\in{B}_{\text{opt}}$ satisfies $B_{p,:}=B^{\star}_{p,:}$ for the
sets $\texttt{ANC}(p)\subseteq S$ or $\texttt{DES}(p)\cup p\subseteq S$.

  3. (c)

$\bar{B}=\arg\min_{B\in{B}_{\text{opt}}}\|B\|_{\ell_{0}}$ satisfies
$\bar{B}_{p,:}=B^{\star}_{p,:}$ if $\texttt{PA}(p)\cup p\subseteq S$ and
$\mathcal{D}^{\star}_{X}$ is faithful to the distribution of $X|H$.

The proof of Theorem 2 is similar in nature to that of Theorem 1 and can be
found in the supplementary material Section E. Further, analogous to Corollary
1, one can readily show the large limit convergence of the population
_DirectLikelihood_ to the sample _DirectLikelihood_ , although we omit this
for brevity.

Remark 2: The conditions needed for identifiability of the unperturbed latent
variable setting (Theorem 2) differ from the perturbed setting (Theorem 1) in
multiple ways. First, there are no conditions on the strength of perturbations
on the observed variables. Further, the latent coefficient matrix
$\Gamma^{\star}$ may be arbitrary without needing conditions like _latent
materiality_. Finally, the setting with unperturbed latent variables requires
two interventional environments where all observed variables are perturbed,
whereas the setting with perturbed latent variables only requires a single
environment with perturbations on all the observed variables and another
environment where the latent variables are perturbed, highlighting that
perturbations on the latent variables is useful for identifiability.

Remark 3: Theorem 2(a) is similar in nature to the backShift procedure [27].
Nonetheless, _DirectLikelihood_ provides additional flexibility such as
controlling the number of latent variables, incorporating do-interventions,
and structure in the strength of shift interventions that lead to more
desirable statistical properties. As an example, a necessary condition for
identifiability using the backShift procedure is that there are at least two
interventional environments. We demonstrate in supplementary material Section
C that this is also a necessary condition with _DirectLikelihood_ if
$\bar{h}=p$. However, under $\bar{h}<p$, _DirectLikelihood_ may attain
identifiability with only a single interventional environment. As another
example, a single interventional environment consisting of the same magnitude
perturbation across the coordinates is sufficient for consistency via
_DirectLikelihood_(see Section D of the supplementary material for the
theoretical statement).

##  4 Connections to distributional robustness

Recent works have demonstrated an intrinsic connection between distributional
robustness and causal inference. Specifically, in the setting where the
response variable is not directly perturbed and there is no latent
confounding, the causal parameter $B^{\star}_{p,:}$ linking the covariates
$X_{{\backslash}p}$ to the response variable $X_{p}$ in the SCM (2) satisfies
the following max-risk optimization problem:

$B^{\star}_{p,:}=\min_{\begin{subarray}{c}\beta\in\mathbb{R}^{p}\\\
\beta_{p}=0\end{subarray}}\max_{\begin{subarray}{c}\mathcal{P}_{e}\in\mathcal{P}\\\
X^{e}\sim\mathcal{P}_{e}\end{subarray}}\|X_{p}^{e}-X^{e}\beta\|_{2}^{2},$ (11)

for a certain perturbation distribution class $\mathcal{P}$ consisting of
distributions $\mathcal{P}_{e}$ indexed by environments $e$ [2]. In
particular, the causal coefficients $B^{\star}_{p,:}$ are solutions to a
robust optimization problem subject to distributional changes to the system
which do not act directly on $X_{p}$. Given access to exogenous variables or
different environments, [26] allow for non-perturbed latent variables and
possibly direct action of change on the target of interest, and prove a
relation between the causal parameters and a particular robust optimization
program.

In this section, we demonstrate that the joint causal parameters $B^{\star}$
minimize a certain worst-case risk in the setting where there may be
perturbations to all the variables including the latent variables, further
strengthening the connection between causal inference and distributional
robustness. We consider the following perturbation distribution class
parameterized by the quantities $C_{\zeta},C_{\psi}\geq 0$:

$\displaystyle\mathcal{P}_{C_{\zeta},C_{\psi}}$ $\displaystyle=$
$\displaystyle\Bigg{\\{}\text{distribution }\mathcal{P}_{e}\text{ over random
pairs }(X^{e},H^{e})\text{ satisfying default SCM }\eqref{eq:SCM_env}\text{
and}$ $\displaystyle w^{e,\star}=w^{1,\star}+\zeta^{e,\star}{\bf{1}}\text{
with
}\zeta^{e,\star}\in[0,C_{\zeta}],\psi^{e,\star}\in[0,C_{\psi}]\Bigg{\\}},$

where the default SCM is the setting with IID latent variables, i.e.
$\Psi^{e,\star}=(1+\psi^{e,\star})\mathcal{I}$. Recall that the sum
$\epsilon^{e}+\delta^{e}$ in the SCM (2) is distributed as follows:
$\epsilon^{e}+\delta^{e}\sim\mathcal{N}(0,\texttt{diag}(w^{e,\star}))$. The
constraints on $w^{e,\star}$ ensure that the perturbations on the observed
variables are IID with magnitude less than a pre-specified level $C_{\zeta}$;
finally, the constraints on $\psi_{e}^{\star}$ ensure that the perturbations
on the latent variables have magnitude less than a pre-specified level of
$C_{\psi}$. We note that the distributions inside $\mathcal{P}$ are specified
by parameters that are invariant, namely the population connectivity matrix
$B^{\star}$, the latent effects matrix $\Gamma^{\star}$, and noise variable
$\epsilon$ with variance of its coordinates encoded in $w^{1,\star}$. We
consider the following worst-case optimization program that identifies
parameters $B,\Gamma,v$ that are robust to perturbations from the class
$\mathcal{P}_{C_{\zeta},C_{\psi}}$:

$(B_{\text{robust}},\Gamma_{\text{robust}},w^{1}_{\text{robust}})=\operatornamewithlimits{arg\,min}_{\begin{subarray}{c}B\text{
is a DAG}\\\
\Gamma\in\mathbb{R}^{p\times\bar{h}},w^{1}\in\mathbb{R}^{p}_{++}\end{subarray}}\max_{\begin{subarray}{c}\mathcal{P}_{e}\in\mathcal{P}_{C_{\zeta},C_{\psi}}\\\
(X^{e},H^{e})\sim\mathcal{P}_{e}\end{subarray}}~{}\texttt{KL}(\Sigma^{e,\star},{\Sigma}_{B,\Gamma,w^{1}}(\bar{\zeta}^{e},\bar{\psi}^{e})),$
(12)

where $\Sigma_{B,\Gamma,w^{1}}(\cdot,\cdot)$ is an estimated covariance model
with definition shown below and KL is the Gaussian Kullback-Leibler divergence
between the estimated and population covariance models. Here,
$\bar{\zeta}^{e},\bar{\psi}^{e}\in\mathbb{R}_{+}$ are estimates for the
nuisance perturbation parameters $\zeta^{e,\star},\psi^{e,\star}$ that vary
across the perturbation distributions. For a given $B,\Gamma,w^{1}$, the
quantities $(\bar{\zeta}^{e},\bar{\psi}^{e})$ are obtained by finding the best
fit to data: $(\bar{\zeta}^{e},\bar{\psi}^{e})=\arg\min_{0\leq\zeta^{e}\leq
C_{\zeta},0\leq\psi^{e}\leq
C_{\psi}}\allowbreak\texttt{KL}(\Sigma^{e,\star},{\Sigma}_{B,\Gamma,w^{1}}(\zeta^{e},\psi^{e}))$
where
${\Sigma}_{B,\Gamma,d}(\zeta^{e},{\psi}^{e})=(\mathcal{I}-B)^{-1}\allowbreak(\texttt{diag}(w^{1})+\zeta^{e}\mathcal{I})+(1+\psi^{e})\Gamma\Gamma^{T})(\mathcal{I}-B)^{-T}$
is the covariance specified by the model parameters.

In comparison to (11), the risk in (12) is measured jointly over the entire
collection of observed variables (via the covariance matrix). As observed
previously, this system-wide perspective is crucial for allowing perturbations
on all of the variables. The following theorem connects the max-risk solutions
$B_{\text{robust}}$ to the causal parameter $B^{\star}$.

######  Theorem 3.

Suppose that the estimated number of latent variables $\bar{h}$ in (12) is
chosen conservatively, i.e. $\bar{h}\geq\text{dim}(H)$. Let the maximum
perturbation size on the observed variables in the perturbation class satisfy
$C_{\zeta}\geq{\kappa^{\star}}(1+2C_{\psi})^{2}(1+\|w^{1,\star}\|_{\infty})(1+\|\Gamma^{\star}\|_{2}^{2}+\|\Gamma^{\star}\|_{2}^{4})$
where
$\kappa^{\star}\equiv\frac{1+\max_{i}\|B^{\star}_{:,i}\|_{2}}{1+\min_{i}\|B^{\star}_{:,i}\|_{2}}$.
Suppose there exists a perturbation distribution
$\mathcal{P}_{e}\in\mathcal{P}_{C_{\zeta},C_{\psi}}$ with parameters
$\zeta^{e,\star}=C_{\zeta}$, $\psi^{e,\star}\neq 0$ such that the random pairs
$(X^{e},H^{e})$ drawn from this distribution satisfy the _latent materiality_
assumption in Definition 1. Then:

  1. 1.

Any optimal connectivity matrix $B\in B_{\text{robust}}$ satisfies
$\text{moral}(B^{\star})\subseteq\text{moral}(B)$

  2. 2.

The optimum of
$\arg\min_{B\in{B}_{\text{robust}}}\|\text{moral}(B)\|_{\ell_{0}}$ is unique
and equal to $B^{\star}$.

Remark 6: The proof of Theorem 3 is presented in the supplementary material
Section F. This theorem result states that the causal parameter $B^{\star}$ is
a minimizer of the max-risk optimization problem over the perturbation class
$\mathcal{P}_{C_{\zeta},C_{\psi}}$ (and produces the sparsest moral graph
among the optimum), establishing a fundamental relation between causality and
distributional robustness. Further, under similar assumptions as required in
Theorem 2, analogous connections can be established for the setting with
unperturbed latent variables.

##  5 Computing the _DirectLikelihood_ estimates

Solving the _DirectLikelihood_ estimator (5) for a DAG $\mathcal{D}$ is a
challenging task, as the problem is non-convex over the decision variables
$B,\Gamma,\\{\Psi^{e}\\}_{e=1}^{m},\\{w^{e}\\}_{e=1}^{m}$. Further, searching
over the space of DAGs is super-exponential in the number of variables. These
computational challenges are common in causal structure learning problems and
are made worse with the presence of multiple environments and latent
confounding. In this section, we propose some heuristics for computing
_DirectLikelihood_ based on perturbation data to find optimal scoring DAGs; we
discuss open questions regarding computations involving the _DirectLikelihood_
in Section 7. The outline of this section is as follows: in Section 5.1, we
describe a method to compute _DirectLikelihood_ for a given DAG structure,
that is, when the support of $B$ is pre-specified. Building on this, in
Section 5.2, we describe some computational heuristics for structure search
over different DAGs.

###  5.1 Scoring a DAG

As announced above, we first assume that a DAG hence also the support of $B$
are pre-specified. The goal, for a given DAG, is to estimate the unknown
parameters. As prescribed in Section 2.1, we employ the _DirectLikelihood_
procedure in the default setting (see Section 2.1.1) with IID latent variables
and jointly observational and interventional data. While the optimization
program (5) is jointly non-convex, solving for the connectivity matrix $B$
with the nuisance parameters $\psi,\Gamma,$ and $\\{w^{e}\\}_{e=1}^{m}$ fixed
is a convex program. Since we are mainly interested in an accurate estimate
for the connectivity matrix, we propose the following alternating minimization
strategy: starting with an initialization of all of the model parameters, we
fix $B$ and perform gradient updates to find updated estimates for the
nuisance parameters, and then update $B$ – by solving a convex program to
optimality with the remaining parameters fixed. We find that the alternating
method described above is relatively robust to the initialization scheme, but
we nonetheless propose the following concrete strategy:

$\displaystyle 1){}$ $\displaystyle B_{(0)}\text{ via linear regression with
observational data}$ (13) $\displaystyle 2){}$
$\displaystyle{w}^{1}_{(0)}=\texttt{diag}\left\\{\mathcal{I}-B_{(0)})\hat{\Sigma}^{1}(\mathcal{I}-B_{(0)})^{T}\right\\}$
$\displaystyle 3){}$ $\displaystyle\Gamma_{(0)}=UD^{1/2}\text{ where
}UDU^{T}\text{ is SVD of
}(\mathcal{I}-B_{(0)})\hat{\Sigma}^{1}(\mathcal{I}-B_{(0)})^{T}$
$\displaystyle 4){}$ $\displaystyle\text{ initialize
}w^{e}_{(0)}={w}^{1}_{(0)}+\zeta^{e}{\bf 1}\text{ and solve
}\zeta^{e},{\psi^{e}}\text{ by 2-dimensional gridding},$

where the first step follows since the DAG structure is known, the fourth step
is based on the observation that for a fixed $B,\Gamma,{w}^{1}_{0}$, the
optimization problems for $\zeta^{e},$ and $\psi^{e}$ decouples across the
environments $e=2,3,\dots,m$. The entire procedure, involving the
initialization step and the parameter updates is presented in Algorithm 1.

Algorithm 1 Scoring $\mathcal{D}$ via _DirectLikelihood_

1: Input: $\hat{\Sigma}^{e}$ for $e=1,2,\dots,m$; regularization $\lambda\geq
0$; number of latent variables $\bar{h}$

2: Initialize parameters: via relation (13)

3: Alternating minimization:

  * (a)

Fixing $(\Gamma_{(t)},\\{(\psi_{(t)}^{e},w^{e}_{(t)})\\}_{e=1}^{m})$, update
$B_{(t+1)}$ by solving the convex optimization program (5). Fixing
$B_{(t+1)}$, perform gradient updates until convergence to find
$(\Gamma_{(t+1)},\\{(\psi^{e}_{(t+1)},w^{e}_{(t+1)})\\}_{e=1}^{m})$

  * (b)

Perform alternating iterates for positive integers $t$ until convergence at
iteration $T$

4: Compute $\text{score}_{\lambda}(\mathcal{D})$: plug-in the estimates
$(B_{(T)},\Gamma_{(T)},\\{(\psi^{e}_{(T)},w^{e}_{(T)})\\}_{e=1}^{m})$ into (6)

5: Output: $\text{ score}_{\lambda}(\mathcal{D})$ and the connectivity matrix
$B_{(T)}$

Step 3 of Algorithm 1 involves two convergence criteria: the convergence of
the gradient steps for the parameters
$(\Gamma_{(t)},\\{(\psi^{e}_{(t)},w^{e}_{(t)})\\}_{e=1}^{m})$ as well as the
convergence of the alternating procedure. For the first criterion, we
terminate the gradient descent when the relative change in the likelihood
score is below $\epsilon_{1}$. For the second criterion, we terminate the
alternating minimization at step $T$ when
$\|B_{(T)}-B_{(T-1)}\|_{\infty}\leq\epsilon_{2}$. In the numerical experiments
in Section 6, we set $\epsilon_{1}=10^{-6}$ and $\epsilon_{2}=10^{-2}$.
Finally, for all our experiments, we select the regularization parameter
$\lambda$ via holdout-validation.

###  5.2 Identifying candidate DAGs

We have discussed how to score a given DAG using the _DirectLikelihood_
estimator. Searching over all DAGs is typically not possible unless the number
of observed variables $p$ is small. In fact, performing a combinatorial search
is known to be very challenging and in some sense NP-hard [5]. One could rely
on greedy strategies [4]; we discuss below a strategy which exploits a
reasonable set of candidate DAGs. In some applications with domain expertise,
a set of plausible DAGs may be considered as candidate DAGs to be scored by
the _DirectLikelihood_. Without this knowledge however, this candidate set
must be obtained from data. In this section, we propose a heuristic to
identify a collection of candidate DAGs to be scored via Algorithm 1. Our
approach is to identify these DAGs by assuming no latent confounding. In
general, fitting a DAG without taking into account the effect of latent
variables yields a denser graph (compared to the population or Markov
equivalent DAGs) since marginalization of latent variables induces confounding
dependencies. As such, scoring such DAGs using Algorithm 1 may yield
connectivity matrices that are more dense than the population connectivity
matrix, although the magnitude of the spurious edges will be small. In our
numerical experiments in Section 6, we find the optimally scoring DAG(s)
(using _DirectLikelihood_) among the candidate DAGs. Then, for each optimal
DAG, we perform backward deletion by removing each of its edges (in the
reverse order of their edge strength) and computing the likelihood score of
the resulting DAGs (using _DirectLikelihood_). We then choose the DAG(s) that
obtain the smallest likelihood score along the entire path.

In Section 1, we outlined procedures to identify DAGs without latent
confounding, with the constraint based PC algorithm, score based GES, and the
hybrid method ARGES being among the most popular for structure learning in the
observational setting. In principal, when domain expertise is not available,
many of these methods can be used to find candidate DAGs. For simplicity, in
our synthetic illustrations in Section 6, we perform GES on pooled
environmental data. The GES procedure greedily adds or deletes edges in the
space of Markov equivalent DAGs based on $\ell_{0}$ regularized likelihood
score and is asymptotically consistent [4]. We select the regularization
parameter to be twice the analogue from the BIC score (as was suggested in
[20]). Algorithm 2 presents the entire procedure of finding candidate DAGs,
scoring them, and selecting the final output.

Algorithm 2 Optimizing _DirectLikelihood_

1: Input: $\hat{\Sigma}^{e}$ for $e\in 1,2,\dots,m$; regularization parameter
$\lambda>0$; number of latent variables $\bar{h}$

2: Find candidate DAGs:
$\tilde{\mathcal{D}}_{\text{cand}}=\mathcal{D}_{1},\mathcal{D}_{2},\dots,\mathcal{D}_{q}$
using domain expertise, GES with pooled data, or some structure learning
algorithm

3: Score each DAG: For each $\mathcal{D}_{i}$, compute
$\text{score}_{\lambda}(\mathcal{D}_{i})$ via Algorithm 1 and obtain
$\tilde{\mathcal{D}}_{\text{opt}}=\arg\max_{\mathcal{D}\in\tilde{\mathcal{D}}_{\text{cand}}}\text{score}_{\lambda}(\mathcal{D})$
and associated connectivity matrices $\tilde{B}_{\text{opt}}$

4: Backward deletion: set
${\mathcal{D}}_{\text{cand}}=\tilde{\mathcal{D}}_{\text{opt}}$ and for each
${\mathcal{D}}\in\tilde{\mathcal{D}}_{\text{opt}}$, perform for
$i=1,2,\dots,\\#\text{edges}({\mathcal{D}})$

  1. 1.

Let ${\mathcal{D}}_{i}$ be the DAG after deleting smallest $i$ edges in
magnitude in ${\mathcal{D}}$

  2. 2.

Compute $\text{score}_{\lambda}({\mathcal{D}}_{i})$ via Algorithm 1

  3. 3.

Add ${\mathcal{D}}_{i}$ to ${\mathcal{D}}_{\text{cand}}$

5: Output: Compute
$\hat{\mathcal{D}}_{\text{opt}}=\arg\max_{\mathcal{D}\in\mathcal{D}_{\text{cand}}}\text{score}_{\lambda}(\mathcal{D})$
and the associated $\hat{B}_{\text{opt}}$.

We remark that the $\arg\max$ in steps 3 and 5 of Algorithm 2 may not be
unique in the infinite sample limit, due to potential non-identifiability. In
practice, the optimization is done to find all optimal DAGs within a relative
tolerance value from the minimum (set at $10^{-3}$ in our experiments), and
outputs also its several associated parameter estimates.

##  6 Experiments

In this section, we illustrate the utility of _DirectLikelihood_ with
simulated and real data. In Section 6.1.1, we study the accuracy of the
_DirectLikelihood_ procedure in estimating the population causal graph
underlying the observed variables. In Section 6.1.2, we provide comparisons of
_DirectLikelihood_ with other methods, including Invariant Causal Predictions,
Causal Dantzig, backShift, and the two-stage deconfounding procedure [12].
Finally in Section 6.2, we evaluate the utility of _DirectLikelihood_ for
learning causal networks on two real datasets, one involving California
reservoir volumes [31] and the other involving protein mass spectroscopy [29].
In supplementary material Section G, we examine _DirectLikelihood_ under model
miss-specifications, namely: non-Gaussian variables in a linear structure
equation model, dependent latent variables, and non-linear SCMs.

Algorithm 2 requires as input the regularization parameter $\lambda$ and the
number of latent variables $\bar{h}$. We select the regularization parameter
$\lambda$ via holdout-validation. Specifically, we partition in the data in
each setting into a training set and a validation set, where the validation
set compromises of some portion of the data in the observational environment.
Unless specified otherwise, the validation set in all numerical experiments is
taken to be $20\%$ of the samples in the observational data. Given estimates
$(\hat{B},\hat{\Gamma},\hat{w}^{1})$ after supplying training data into
_DirectLikelihood_ , we then compute the validation performance as the
negative log-likelihood
$\ell(\hat{B},\hat{\Gamma},\mathcal{I},\hat{w}^{1},\Sigma_{\texttt{valid}}^{1})$,
where $\Sigma_{\texttt{valid}}^{1}$ is the sample covariance of the validation
data. As smaller negative log-likelihood is indicative of better fit to data,
we select $\lambda$ that minimizes the negative log-likelihood on validation
data. We observe that our procedure is generally robust to the choice
$\bar{h}$ and furthermore, _DirectLikelihood_ procedure is consistent as long
as $\bar{h}\geq\text{dim}(H)$ (see Section 3). Thus, we select $\bar{h}$ to be
moderately large (relative to the ambient dimension) so that it is an
overestimate of the true number of latent variables, although holdout-
validation can also be performed to select $\bar{h}$.

###  6.1 Synthetic experiments

####  6.1.1 DAG structural recovery

_Setup_ : we consider a collection of $p=10$ observed variables influenced by
$h\in\\{1,2\\}$ latent variables. To generate the connectivity matrix
$B^{\star}\in\mathbb{R}^{p\times p}$, we sample from an Erdös Rényi graph with
edge probabilities $0.1$ until we find a DAG structure, and form $B^{\star}$
by setting edge strengths equal to $-0.7$. The resulting DAG and connectivity
matrix consists of $10$ nonzero entries. The entries of the latent coefficient
matrix $\Gamma^{\star}\in\mathbb{R}^{p\times h}$ are generated IID distributed
uniformly from the interval $[0,\sqrt{0.3/\sqrt{h}}]$ and the entries below
$0.5\sqrt{0.3/\sqrt{h}}$ are set to zero. The noise term $\epsilon$ is
distributed according to $\epsilon\sim\mathcal{N}(0,0.5\mathcal{I}_{p})$.
Unless otherwise specified, the latent variables $H$ are generated as
$H\sim\mathcal{N}(0,\mathcal{I}_{h})$. These parameters specify the
distribution of the observed and latent variables when there are no
perturbations and we denote this environment by $e=1$. In addition to this
observational environment, we suppose there are $m-1$ interventional
environments. The number of samples generated in the observational environment
is set to $n^{1}=300$ and $n^{e}=5t$ for positive integer $t$ is the sample
size for environment $e$. The values for $t$, the number of environments, and
the magnitude of perturbations on the observed and latent variables is
specified later.

For each environment $e$, we set
$\delta^{e}_{k}\sim\mathcal{N}(0,\zeta+\text{Unif}(0,1))$ for $k=1,2,\dots,p$
and certain values of $\zeta$, and
$H^{e}\sim(1+\psi^{e,\star})\mathcal{N}(0,\mathcal{I}_{h})$ with
$\psi^{e,\star}\sim\frac{1}{2}(1+\text{Unif}(0,1))$. We generate data from
$m=7$ environments, one observational environment with no perturbations and
six interventional environments, and consider the following five settings
$(a)$ $h=1,\zeta=5$, $(b)$ $h=1,\zeta=2$, $(c)$ $h=2,\zeta=5$, and $(d)$
$h=1,\zeta=5$ and the last five environments have two observed variables that
are chosen randomly to receive do-interventions with values set identically
equal to $5$. The perturbation data for each setting is supplied to the
_DirectLikelihood_ procedure to score each DAG in a collection of candidate
DAGs. We set $\bar{h}=h+1$ in the _DirectLikelihood_ estimator (5) and
constrain the latent variable perturbation $\psi^{e}\leq\psi_{\text{max}}=2$
for interventional environments $e=2,3,\dots,7$. We then evaluate the accuracy
of the _DirectLikelihood_ procedure (Algorithm 2) for DAG structural recovery
in each of the settings $(a-d)$ averaged across $10$ independent trials. The
accuracy of DAG recovery is computed with respect to false positives (edges
produced in the estimated DAG that are missing or in the reverse direction in
the population DAG) and true positives (edges in the estimated DAG present in
the correct direction in the population DAG). The set of candidate DAGs to
score via Algorithm 2 is obtained by performing the GES algorithm on pooled
data. Since _DirectLikelihood_ always finds a single graph as the optimum in
these numerical experiments, we compute for comparison the average size of the
observational Markov equivalence class obtained after the pooled GES step in
setting $(a)$: $9$ DAGs for $t=64$, $8.8$ for $t=16$, $9.3$ for $t=4$, $6$ for
$t=2$, $6.4$ for $t=1$.

(a) $h=1,\zeta=5$ (b) $h=2,\zeta=5$ (c) $h=1,\zeta=2$ (d) $h=1,\zeta=5$, do-
interventions

Figure 2: Structure estimation accuracy of Algorithm 2 (best scoring DAG)
using candidate DAGs obtained by the GES algorithm on pooled data for
different problem settings. Total number of true discoveries equals to $10$.
The curve for each $t$ corresponds to $5t$ samples for each interventional
environment, with $t\in\\{1,2,4,16,64\\}$. For each curve, the accuracy of the
estimated DAG in comparison to the population DAG is calculated by ordering
the edges according to their strengths and sequentially counting (from
strongest edge to weakest edge) an edge to be a false discovery if it is
missing or in a reverse direction in the population DAG, and otherwise count
as a true discovery. Each curve is averaged across 10 independent trials.

####  6.1.2 Comparison to previous methods

_DirectLikelihood_ _vs Invariant Causal Predictions, causal Dantzig, and
backShift_ : We compare the performance of _DirectLikelihood_ to Invariant
Causal Predictions [23], causal Dantzig [25], and backShift [27] for finding
the causal parents of a response variable. Consistency guarantees for these
previous methods require at least one of these assumptions i) there are no
latent effects, ii) the latent variables remain unperturbed across
environments, iii) the response variable remains unperturbed across
environments. Specifically, while Invariant Causal Predictions (ICP) does not
impose any conditions on the specific relationship among the covariates,
assumptions i) and iii) are needed for consistently estimating the causal
parents. Causal Dantzig allows for latent effects, although it requires
assumptions ii) and iii) for consistency. Finally, guarantees using the
backShift procedure require assumption ii). Via numerical simulations, we
illustrate the impact of using these previous approaches when assumptions
i-iii) are not satisfied. We leave out the comparison to Instrumental
Variables [1] as the number of instruments (environments) must be larger than
the number of covariates. One in principle could apply Anchor regression [26],
although this method does not obtain causal parameters.

We consider a causal structure among $p=10$ variables and $h=1$ latent
variable with $X_{p}$ denoting the response variable and $X_{1:p-1}$ denoting
the covariates. We modify the parents and children of the DAG in Section 6.1.1
(so that the response variable in the population DAG has more parents and
children): $X_{3},X_{4}$ are parents of the response variable and
$X_{7},X_{8},X_{9}$ are children of the response variable. We set all the edge
weights in the DAG to be $-0.7$. The entries of the latent coefficient matrix
$\Gamma^{\star}\in\mathbb{R}^{p\times 1}$ are generated IID distributed
uniformly from the interval $[0,\sqrt{0.3}]$. The noise term $\epsilon$ is
distributed according to $\epsilon\sim\mathcal{N}(0,0.5\mathcal{I}_{p})$. We
generate an observational environment $e=1$ and four interventional
environments $e=2,3,4,5$ and consider the following four settings:

  1. Setting 1.

no perturbations on the response variable and the latent variables:
$\delta_{k}^{e}\sim\mathcal{N}(0,5+\text{Unif}(0,1))$ for all
$k=1,2,\dots,p-1$ and $\psi^{e}=0$ for all $e$

  2. Setting 2.

no perturbations on the latent variables and perturbations on the response
variable: $\delta_{k}^{e}\sim\mathcal{N}(0,5+\text{Unif}(0,1))$ for all
$k=1,2,\dots,p$ and $\psi^{e}=0$ for all $e$

  3. Setting 3.

no perturbations on the response variable and perturbations on the latent
variables : $\delta_{k}^{e}\sim\mathcal{N}(0,5+\text{Unif}(0,1))$ for all
$k=1,2,\dots,p-1$ and $\psi^{e}\sim 1+\text{Unif}(0,1)$ for all $e$

  4. Setting 4.

perturbations on the response and latent variables:
$\delta_{k}^{e}\sim\mathcal{N}(0,5+\text{Unif}(0,1))$ for all $k=1,2,\dots,p$
and $\psi^{e}\sim 1+\text{Unif}(0,1)$ for all $e$

We obtain $1000$ IID observational data and interventional data and supply
this data to the algorithms for each of the procedures. For the ICP and causal
Dantzig methods, we set the significance threshold at $0.01$ and for the
backShift procedure, we perform stability selection (with stability parameter
$0.70$) as is prescribed in [27]. We produce the set of candidate DAGs for the
_DirectLikelihood_ procedure using pooled GES, and set $\bar{h}=2$. Table 2
compares the false positives and true positives associated with identifying
the causal parents of the response variable (across $10$ independent trials)
of _DirectLikelihood_ and competing methods. The population has two causal
parents for the response variable, so that the total size of true positives is
at most 2.

A few remarks are in order. First, in all of the settings, ICP returns the
empty set due to the latent effects. Further, backShift performs poorly in all
settings, even when there are no perturbations on the latent variables (the
setting where [27] prove identifiability guarantees). We do observe however
that if we increase the strength and dynamic range of the perturbations,
backShift is able to accurately estimate the causal parents when the
perturbations do not affect the latent variables. Namely, we consider Setting
2 with $\delta_{k}^{e}\sim\mathcal{N}(0,5+5*\text{Unif}(0,1))$, set the
stability threshold at $0.51$ and find that $\text{TP}=0.9;\text{FP}=0$ when
$n^{e}=1000$ and $\text{TP}=1.5;\text{FP}=0$ when $n^{e}=10000$ for
$e=1,2,\dots,5$. Next, as supported by theoretical guarantees, causal Dantzig
estimates the causal parameters accurately in Setting 1 when there are no
perturbations on the response or the latent variables. However, in settings
where the latent variables or the response variables are perturbed, causal
Dantzig yields many false positives and often incorrectly identifies the
causal children of the response variable as the estimated causal parents. The
_DirectLikelihood_ procedure on the other hand, does not yield many false
positives and has comparable power performance. We note that the power
performance of _DirectLikelihood_ in Setting 4 is negatively affected by the
performance of pooled GES to select candidate DAGs. Specifically, the largest
number of true positives among the candidate DAG (without scoring) is on
average equal to $1.2$ and thus _DirectLikelihood_ is performing as well as
possible given the candidate DAGs that are supplied as input. In Section 7, we
discuss future directions for more rigorous techniques to obtain and score
candidate DAGs.

Method | Setting 1 | Setting 2 | Setting 3 | Setting 4
---|---|---|---|---
_DirectLikelihood_ | TP = 2, FP = 0.3 | TP = 2, FP = 0.1 | TP = 2, FP = 0.8 | TP = 1.2, FP = 0.5
causal Dantzig | TP = 2, FP = 0 | TP = 2, FP = 5 | TP = 2, FP = 3 | TP = 2, FP = 5.1
backShift | TP = 0, FP = 1.6 | TP = 0, FP = 0 | TP = 0, FP = 1.4 | TP = 0, FP = 0
ICP | empty set | empty set | empty set | empty set

Table 2: Comparison of _DirectLikelihood_ with other methods for identifying
the causal parents of the response variable. Maximum possible number of true
discoveries is equal $2$. There are $1000$ samples in the observational and
each of the four interventional environments.

_DirectLikelihood_ _vs two-stage deconfounding_ : The two-stage deconfounding
procedure first employs a sparse+low-rank decomposition on data from each
environment to deconfound the latent effects and then employs the
_DirectLikelihood_ procedure with $\Gamma\equiv 0$ (i.e. as latent effects are
in principle removed) in the second stage. As described in Section 1, the
accuracy of the first step heavily relies on the denseness of the latent
effects. We generate the following synthetic example to compare the
performance of these algorithms. We set $h=3$ and consider the synthetic setup
described earlier in Section 6.1.1 with the following modifications: the first
two columns of $\Gamma^{\star}\in\mathbb{R}^{p\times 3}$ consist of standard
basis elements with the coordinate corresponding to $X_{6}$ and $X_{5}$
nonzero, and a third column with entries sampled IID from the uniform
distribution with entries less than $0.5$ set to zero. We generate an
observational environment $e=1$ and four interventional environments
$e=2,3,4,5$ where $\delta_{k}^{e}\sim\mathcal{N}(0,\zeta+\text{Unif}(0,1))$
for $k=1,2,\dots,p$ with $\zeta=2$. We generate the latent perturbation
coefficient $\psi^{e}\sim\text{Uniform}(0,0.5)$. We obtain $n^{1}=1000$ IID
observational data and $5t$ IID interventional data for each interventional
environment with $t\in\\{60,200\\}$. The number of latent variables included
in the model must be selected by the user in the _DirectLikelihood_ and two-
stage deconfounding procedures ($\bar{h}$ in _DirectLikelihood_ and two
regularization parameters in the first step of the two-stage deconfounding
procedure). Since we are interested in comparing identifiability properties of
these procedures, we choose $\bar{h}=3$ in _DirectLikelihood_. Further, we
chose the regularization parameters in the deconfounding step of the two-stage
deconfounding procedure by choosing the best predictive model on a validation
set with number of latent variables less than or equal to $h=3$.

Both the _DirectLikelihood_ procedure and the second stage of the two-stage
deconfounding score a set of candidate DAGs. Noticing that the sparseness of
the latent effects induces a spurious edges between the pairs
$(X_{5},X_{10}),(X_{8},X_{10}),\allowbreak(X_{5},X_{3})$, we generate $8$
candidate DAGs adding $5$ edges at random to all $8$ DAGs in the Markov
equivalence class of the population DAG $\mathcal{D}^{\star}_{X}$ and a final
candidate DAG that adds the directed edges $X_{5}\to X_{10}$, $X_{8}\to
X_{10}$, $X_{5}\to X_{3}$ to the population DAG $\mathcal{D}^{\star}_{X}$.
Thus, the total number of candidate DAGs is equal to $9$. Table 3 compares the
structural recovery (across $10$ independent trials) of _DirectLikelihood_ and
the two-stage deconfounding procedure. We observe that since the denseness
assumption is violated, the two-stage deconfounding produces a DAG with false
positives, even when the number of samples in each environment is large (i.e.
$1000$ samples). Furthermore, in the low-sample regime, the two-stage
procedure yields fewer true discoveries than _DirectLikelihood_. It is worth
noting that in addition to the superior performance of _DirectLikelihood_ as
compared to two-stage deconfounding, the _DirectLikelihood_ solution is faster
to compute since it does not involve tuning two regularization parameters.

Method |  $300$ samples/ interven. environment |  $1000$ samples/ interven. environment
---|---|---
_DirectLikelihood_ | TP = 10, FP = 0 | TP = 10, FP = 0
two-stage deconfounding | TP = 9.2, FP = 1.9 | TP = 10, FP = 2.9

Table 3: Comparison of _DirectLikelihood_ with two-stage deconfounding and
backShift procedures. Maximum possible number of true discoveries is equal
$10$. There are $1000$ samples in the observational environment and
$\\{300,1000\\}$ samples in each of the four interventional environments.

###  6.2 Experimental results on real datasets

####  6.2.1 California reservoirs

The California reservoir network consists of approximately $1530$ reservoirs
that act as buffer against severe drought conditions and are a major source of
water for agricultural use, hydropower generation, and industrial use. Water
managers of these reservoirs have to assess likelihood of system-wide failure
and effectiveness of potential policies. Due to similarities in hydrological
attributes (e.g. altitude, drainage area, spatial location), the reservoir
network is highly interconnected. Thus, effective reservoir management
requires understanding of reservoir interdependencies. [31] used historical
data of volumes of the largest 55 reservoirs to obtain an undirected graphical
model of the California reservoir network. The previous analysis does not
provide causal implications, namely how change in management of one reservoir
affects the entire system. As such, we seek a causal network among the
reservoirs.

We consider the 10 largest reservoirs (with respect to capacity) in
California, where daily volume data (downloaded from
https://github.com/armeentaeb/WRR-Reservoir) are available during the period
of study (January 2003–December 2015). Following the preprocessing steps in
[31], we average the data from daily down to 156 monthly observations. A
seasonal adjustment step is performed to remove predictable seasonal patterns.
The resulting data was demonstrated in [31] to be well-approximated by a
multivariate Gaussian distribution. The data naturally be categorized to an
observational environment and three interventional environments. Specifically,
the observational environment is during a normal period (2003-2006, 2010-2012)
with no drought conditions, one interventional environment corresponding to an
abnormally dry period (2007, 2013) with small changes to management, moderate
drought period (2008-2009) with significant changes to management, and a
severe drought period (2014-2015) with extreme changes to management. We take
as training data the periods 2003-2006, 2010 as well as all the interventional
data, and take the validation data to be the period 2011-2012 from
observational data.

We include two latent variables in the model as was discovered in [31] (e.g
$\bar{h}=2$) and supply the observational and interventional data to the
_DirectLikelihood_ procedure. After holdout-validation, we identify a causal
graph with $9$ edges as shown in Figure 3(a). The connections are between
pairs of reservoirs with at least of these commonalities: a) similar
hydrological attributes (e.g. hydrological zone and elevation), b) coordinated
management by a district or a state-wide project, and c) similar usages (e.g.
hydropower generation). For example, the reservoirs New Melones (NML), Don
Pedro (NP), New Exchequer (EXC) and Pine Flat (PNF) are all in the San Joaquin
district. Further, Shasta (SHA), Trinity (CLE), Oroville (ORO) and Folsom
(FOL) are in the network of Central Valley and State Water projects and their
reservoir operations are coordinated. Finally, Almanor (ALM) and Trinity (CLE)
are primarily used for hydroelectric power generation. Specifically, the
Pacific Gas $\&$ Electric Company owns Almanor and has historically negotiated
with Trinity Public Utilities District that use water from Trinity to generate
electricity.

For comparison, we obtain _DirectLikelihood_ estimates when no latent
variables are included in the model. The model we obtain after holdout-
validation contains $14$ edges as shown in Figure 3(b). Unlike the structure
with latent variables, the model without latent variables contains many
spurious edges: namely connections between pairs of reservoirs that are
geographically far apart (e.g. Oroville - Pine Flat and Trinity - New
Melones). The same phenomena was also noted in [31] in the context of
undirected graphical models.

(a) with latent variables (b) without latent variables

Figure 3: Causal graphical structure among the volumes of 10 largest
reservoirs in California (with respect to capacity) using _DirectLikelihood_
procedure: a) incorporating latent variables, and b) without latent variables.

####  6.2.2 Protein expressions

We next analyze the protein mass spectroscopy dataset [29]. This dataset
(downloaded from
http://www.sciencemag.org/content/suppl/2005/04/21/308.5721.523.DC1/Sachs.SOM.Datasets.zip)
contains a large number of measurements of the abundance of $11$
phosphoproteins and phospholipids recorded under different experimental
conditions in primary human immune system cells. The different experimental
conditions are characterized by associated reagents that inhibit or activate
signaling nodes, corresponding to interventions at different points of the
protein-signaling network. Following the previous works [18, 17], we take 8
environments consisting of an observational environment and 7 interventional
environments. Knowledge about some of the “ground truth” is available, which
helps verification of results.

To identify a set of candidate DAGs to score using our _DirectLikelihood_
procedure, we consider all DAGs that are Markov equivalent to the ground truth
DAG reported in [29] (total of $176$ DAGs). We include two latent variables
(e.g. set $\bar{h}=2$) in the _DirectLikelihood_ estimator, and after holdout-
validation, we select a causal graphical structure with six total edges. We
compare our findings to the direct causal relations reported in the literature
[11, 17, 18, 29] in Table 4.

Edge |  [29](a) |  [29](b) | [18] | [11] | [17]
---|---|---|---|---|---
PKC $\to$ p38 | ✓ | ✓ | ✓ | ✓ | ✓
Akt $\to$ Erk |  |  | ✓ |  | ✓
Mek $\to$ Raf |  |  | ✓ | ✓ | ✓
PKC $\to$ JNK | ✓ | ✓ | ✓ | ✓ |
PIP2 $\to$ PIP3 |  |  | ✓ |  |
PLCg $\to$ PIP2 | ✓ | ✓ |  | ✓ | ✓

Table 4: Comparing the findings of _DirectLikelihood_(ordered by edge
strength) with different causal discovery methods. Here, we are only including
edges found by _DirectLikelihood_ and note that additional edges have been
identified by the other methods. The consensus network according to [29] is
denoted by “[29](a)” and their reconstructed network by “[29](b)”.

We next compare our findings when we account for latent effects to the setting
where we do not account for latent effects, namely by setting $\Gamma\equiv 0$
in the _DirectLikelihood_ procedure. The causal graphical model we obtain
without accounting for latent effects also consists of six edges, but three of
those are different than the model that incorporates latent variables. These
edges (in the order of strength) are Akt $\to$ PKA, PIP3 $\to$ PIP2, PLCg
$\to$ PIP3. The edge Akt $\to$ PKA was never reported in previous work, and
the edge PLCg $\to$ PIP3 was not reported in methods that accounted for latent
effects [17, 18]. Thus, these two edges appear to be spurious dependencies due
to common latent variables. The edge PIP3 $\to$ PIP2 in the causal structure
without latent variables is also reported in [11, 17, 29] while the reverse
direction PIP2 $\to$ PIP3 is discovered in our causal structure with latent
effects (see Table 4) and was also reported in [18].

##  7 Discussion

In this paper, we proposed a framework to model unspecific perturbation data
among a collection of Gaussian observed and latent variables. It can be
represented as a certain mixed-effects linear structural causal model where
the interventions are modeled as random effects which propagate dynamically
through the structural equations. This framework allows for perturbations on
all components of the system, including a response variable of interest or the
latent variables. We demonstrated the utility of _DirectLikelihood_ for
identifying causal relationships on both synthetic and real datasets.

There are several interesting directions for further investigation that arise
from our work. In Section 5, we proposed heuristics for searching over the
space of DAGs and for solving the _DirectLikelihood_ estimator (5) to score
DAGs. While the empirical results in Section 6 support the utility of our
heuristics, there is much room for more rigorous optimization techniques (e.g.
provably consistent greedy methods). Next, the theoretical results in Section
3 are based on analysis in the large data limit. However, our empirical
results in Section 6 suggest that the _DirectLikelihood_ procedure may provide
accurate estimates with moderate data size. As such, an exciting research
direction is to develop high-dimensional consistency guarantees for
_DirectLikelihood_. Further, in the setting where the perturbations are
limited, there may be multiple DAGs that are equally representative of the
data, or equivalently, multiple DAGs that yield the same exact likelihood
score in the population case (known as the interventional Markov equivalence
class). The characterization of these equivalence classes will be central to
developing greedy algorithms as well as constructing active learning
strategies for maximally informative interventions [13, 14]. Finally, the
perturbation model (2) assumes a linear relationship between the observed and
latent variables. It would be of practical interest to explore extensions of
our framework to non-linear settings, or alternatively, characterize the
extent to which linear models capture the causal effects.

## Acknowledgements

A. Taeb and P. Bühlmann both acknowledge scientific interaction and exchange
at “ETH Foundations of Data Science”. The authors would like to thank Mona
Azadkia, Yuansi Chen, Juan Gamella, and Marloes Maathius for helpful
discussions and feedback on the manuscript. The dataset and the code to
produce the results of this paper can be found at
https://github.com/armeentaeb/perturbations-and-causality.

## References

  * [1] J. Angrist, G. Imbens, and D. Rubin, Identification of causal effects using instrumental variables, Journal of the American Statistical Association, 91 (1996), pp. 444–455.
  * [2] P. Bühlmann, Invariance, causality and robustness, Statistical Science, 35 (2020), pp. 404–426.
  * [3] V. Chandrasekaran, P. Parillo, and A. Willsky, Latent variable graphical model selection via convex optimization, Annals of Statistics, 40 (2012), pp. 1935–1967.
  * [4] D. Chickering, Optimal structure identification with greedy search, Journal of Machine Learning Research, 3 (2002), pp. 507–554.
  * [5] D. Chickering, C. Meek, and D. Heckerman, Large-sample learning of bayesian networks is np-hard, Journal of Machine Learning Research, 5 (2004), pp. 1287–1330.
  * [6] D. Colombo, M. Maathuis, M. Kalisch, and T. Richardson, Learning high-dimensional directed acyclic graphs with latent and selection variables, Annals of Statistics, 40 (2012), pp. 294–321.
  * [7] A. Dawid, Decision-theoretic foundations for statistical causality, arXiv 2004.12493, (2020).
  * [8] A. Dawid and V. Didelez, Identifying the consequences of dynamic treatment strategies: a decision-theoretic overview, Statistical Surveys, 4 (2010), pp. 184–231.
  * [9] A. Dixit, O. Parnas, and B. Li, Perturb-seq: dissecting molecular circuits with scalable single-cell rna profiling of pooled genetic screens, Cell, 167 (2016), pp. 1853–1866.
  * [10] M. Drton and M. Maathius, Structure learning in graphical modeling, Annual Review of Statistics and Its Application, 4 (2017), pp. 365–393.
  * [11] D. Eaton and K. Murphy, Exact bayesian structure learning from uncertain interventions, In Proceedings of the 10th International Conference on Artificial Intelligence and Statistics (AISTATS), (2007).
  * [12] B. Frot, P. Nandy, and M. Maathius, Robust causal structure learning with hidden variables, Journal of Royal Statistical Society, Series B, 81 (2019), pp. 459–487.
  * [13] A. Hauser and P. Bühlmann, Characterization and greedy learning of interventional markov equivalence classes of directed acyclic graphs, Journal of Machine Learning Research, 13 (2012), pp. 2409–2464.
  * [14] A. Hauser and P. Bühlmann, Two optimal strategies for active learning of causal models from interventional data, International Journal of Approximate Reasoning, 4 (2014), pp. 926–939.
  * [15] A. Hauser and P. Bühlmann, Jointly interventional and observational data: estimation of interventional markov equivalence classes of directed acyclic graphs, Journal of Royal Statistical Society, Series B, 77 (2015), pp. 291–318.
  * [16] A. McLean, L. Sanders, and W. Walter, A unified approach to mixed linear models, Journal of American Statistical Association, 45 (1991), pp. 54–64.
  * [17] N. Meinshausen, A. Hauser, J. Mooij, J. Peters, P. Versteeg, and P. Bühlmann, Methods for causal inference from gene perturbation experiments and validation, Proceeding of National Academy of Sciences, 113 (2016), pp. 7361–7368.
  * [18] J. Mooij and T. Heskes, Cyclic causal discovery from continuous equilibrium data, In Proceedings of the 29th Annual Conference on Uncertainty in Artificial Intelligence, (2013), pp. 431–439.
  * [19] P. Nandy, A. Hauser, and M. Maathius, High-dimensional consistency in score-based and hybrid structure learning, Annals of Statistics, 46 (2018), pp. 3151–3183.
  * [20] C. Nowzohour and P. Bühlmann, Score-based causal learning in additive noise models, Statistics, 50 (2016), pp. 471–485.
  * [21] J. Pearl, Causality: Models, reasoning, and inference, Cambridge University Press, 2nd edition, 2009.
  * [22] J. Peters and P. Bühlmann, Identifiability of Gaussian structural equation models with equal error variances, Biometrika, 101 (2014), pp. 219–228.
  * [23] J. Peters, P. Bühlmann, and N. Meinshausen, Causal inference using invariant prediction: identification and confidence intervals, Journal of the Royal Statistical Society, Series B, 78 (2016), pp. 947–1012.
  * [24] J. Robins, M. Hernan, and B. Brumback, Marginal structural models and causal inference in epidemiology, Epidemiology, 11 (2000), pp. 550–560.
  * [25] D. Rothenhäusler, P. Bühlmann, and N. Meinshausen, Causal dantzig: fast inference in linear structural equation models with hidden variables under additive interventions, Annals of Statistics, 47 (2019), pp. 1688–1722.
  * [26] D. Rothenhäusler, N. Meinshausen, P. Bühlmann, and J. Peters, Anchor regression: heterogeneous data meets causality, arXiv:1801.06229, (2020).
  * [27] D. Rothhenhäusler, C. Heinze, J. Peters, and N. Meinshausen, backshift: Learning causal cyclic graphs from unknown shift interventions, In Advances in Neural Information Processing Systems, (2016).
  * [28] D. Rubin, Causal inference using potential outcomes, Journal of the American Statistical Association, 100 (2015), pp. 322–331.
  * [29] K. Sachs, O. Perez, D. Lauffenburger, and G. Nolan, Causal protein-signaling networks derived from multiparameter single-cell data, Science, 308 (2005), pp. 523–529.
  * [30] P. Spirites, C. Glymour, and R. Scheines, Causation, Prediction, and Search, Cambridge: MITPress, 2000.
  * [31] A. Taeb, J. Reager, M. Turmon, and V. Chandrasekaran, A statistical graphical model of the california reservoir system, Water Resources Research, 53 (2017), pp. 9721–9739.
  * [32] S. van de Geer and P. Bühlmann, $\ell_{0}$-penalized maximum likelihood for sparse directed acyclic graphs, Annals of Statistics, 41 (2013), pp. 536–567.
  * [33] T. Verma and J. Pearl, Equivalence and synthesis of causal models, In Proceedings of the 6th Conference on Uncertainty in Artificial Intelligence (UAI), (1991), pp. 255–270.
  * [34] Y. Wang, L. Solus, K. Yang, and C. Uhler, Permutation based causal inference algorithms with interventions, In Advances in Neural Information Processing Systems, (2017), pp. 5822–5831.

## Supplementary Material

The proof of the theoretical results in the supplementary material are based
on the following population quantities that we summarize. Let
$B^{\star}\in\mathbb{R}^{p\times p}$ be the population connectivity matrix,
$\Gamma^{\star}\in\mathbb{R}^{p\times h}$ be the matrix encoding the effects
of latent variables on the observed variables,
$w^{1,\star}\in\mathbb{R}^{p}_{++}$ encode the variance of the coordinates of
$\epsilon$, and $w^{e,\star}_{k}=w^{1,\star}_{k}+\text{var}(\delta^{e}_{k})$
with $w^{e,\star}\in\mathbb{R}^{p}_{++}$. Let
$\\{\psi^{e,\star}\\}_{e=1}^{m}\subseteq\mathbb{R}_{+}$ be the perturbation on
the latent variables. Let
$\kappa^{\star}=\frac{1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2}}{1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2}}$.
Finally, for a matrix $M\in\mathbb{R}^{d\times d}$, we denote $\|M\|_{2}$ as
the spectral norm (largest singular value) of $M$.

##  Appendix A Incorporating do-interventions

Recall from Section 2.1.3 (main paper) that the structural causal model with
do-interventions is modified to be:

$\begin{gathered}X^{e}=\mathcal{F}_{\texttt{do}(e)^{c}}(B^{\star}{X^{e}}+\Gamma^{\star}{H}^{e}+\epsilon^{e})+\delta^{e}\\\
H^{e}\sim\mathcal{N}(0,\Psi^{\star,e}),\end{gathered}$

Given data cross environments $e=1,2,\dots,m$, we can optimize the parameters
of the model via the negative log-likelihood (4) (main paper). It is
straightforward to see that the negative log-likelihood
$\log\texttt{prob}(\cdot)$ decouples across the parameters
$(B,\Gamma,\\{\Psi^{e}\\}_{e=1}^{m},\\{{w^{e}_{\texttt{do}(e)^{c}}}\\}_{e=1}^{m})$
and $\\{{w^{e}_{\texttt{do}(e)}}\\}_{e=1}^{m}$. In other words, the structure
of the DAG $\mathcal{D}$ only plays a role in the term involving the
parameters
$(B,\Gamma,\\{\Psi^{e}\\}_{e=1}^{m},\\{{w^{e}_{\texttt{do}(e)^{c}}}\\}_{e=1}^{m})$,
and we thus focus on that component of the likelihood:

$\displaystyle(\hat{{B}},\hat{{\Gamma}},\\{(\hat{{\Psi}}^{e},\hat{{w}}^{e})\\}_{e=1}^{m})=\operatornamewithlimits{arg\,min}_{\begin{subarray}{c}B\in\mathbb{R}^{p\times
p},\Gamma\in\mathbb{R}^{p\times\bar{h}}\\\
\\{\Psi^{e}\\}_{e=1}^{m}\subseteq\mathbb{S}^{\bar{h}}_{++},\\{w^{e}\\}_{e=1}^{m}\subseteq\mathbb{R}^{p}_{++}\end{subarray}}$
$\displaystyle\sum_{e=1}^{m}\hat{\pi}^{e}\ell(B,\Gamma,\Psi^{e},w^{e};\hat{\Sigma}^{e},\texttt{do}(e)),$
(14) subject-to $\displaystyle~{}~{}~{}B\text{ compatible with }\mathcal{D}$

where

$\displaystyle\ell(\cdot)=\log\det\left(\left[\texttt{diag}(w^{e})+\Gamma\Psi^{e}\Gamma^{T}\right]_{\texttt{do}(e)^{c}}\right)$
$\displaystyle+\mathrm{trace}\left(\left[\texttt{diag}(w^{e})+\Gamma\Psi^{e}\Gamma^{T}\right]_{\texttt{do}(e)^{c}}^{-1}\left[(\mathcal{I}-\mathcal{F}_{\texttt{do}(e)^{c}}B)\hat{\Sigma}^{e}(\mathcal{I}-\mathcal{F}_{\texttt{do}(e)^{c}}B)^{T}\right]_{\texttt{do}(e)^{c}}\right),$

Here, we assume that the location of the do-interventions are known so that
the input to the program (14) are the sample covariance matrices
$\hat{\Sigma}^{e}$, the do-intervention locations $\texttt{do}(e)$, and the
estimate $\bar{h}$ for the number of latent variables.

##  Appendix B Proof of Theorem 1 (main paper)

Recall that the _DirectLikelihood_ estimator (5) (main paper) scores candidate
DAGs and the best scoring DAGs are chosen as output (there is no penalty term
in the score function as $\lambda=0$ in the population setting). As stated in
the theoretical results in Section 3 (main paper), we assume that all possible
DAGs among the observed variables may be scored. Thus, we consider the
_DirectLikelihood_ estimator (5) (main paper) specialized to IID latent
variables and $e=1$ denoting observational environment with the additional
decision variable over the space of DAGs to find optimal DAGs with associated
parameter estimates:

$\displaystyle(\hat{B},\hat{\Gamma},\hat{\psi},\\{\hat{{w}}^{e}\\}_{e=1}^{m})=\operatornamewithlimits{arg\,min}_{\begin{subarray}{c}B,\Gamma,\psi\in\mathbb{R}^{m}_{+}\\\
\\{w^{e}\\}_{e=1}^{m}\subseteq\mathbb{R}^{p}_{++}\\\ \text{DAG
}\mathcal{D}\end{subarray}}$
$\displaystyle\sum_{e=1}^{m}\pi^{e,\star}\ell(B,\Gamma,(1+\psi^{e})\mathcal{I},w^{e};\Sigma^{e,\star})$
(15) subject-to. $\displaystyle B\text{ compatible with
}\mathcal{D}~{};~{}\|\psi\|_{\infty}\leq C_{\psi}$
$\displaystyle\psi^{1}=0~{};~{}w^{e}\succeq w^{1}\text{ for }e=2,\dots,m$

Here the decision variable $\psi$ encodes the latent perturbations and
consists of coordinates $\psi=(\psi^{1},\psi^{2},\ldots,\psi^{m})$. As stated
in Theorem 1 (main paper), we assume that the number of latent variables in
the model is a conservative estimate of the true number of latent variables,
i.e. $\bar{h}\geq\text{dim}(H)$. The proof strategy for proving Theorem 1
(main paper) is based on appealing to the following three lemmas:

######  Lemma 1.

Optimal solutions of (15) satisfy the following equivalence:

$\displaystyle(B,\Gamma,\psi,\\{{w}^{e}\\}_{e=1}^{m})\text{ optimum of
}\eqref{eqn:estim_cons_exp_pert}$ $\displaystyle\iff B\text{ compatible with a
DAG},\Gamma\in\mathbb{R}^{p\times\bar{h}},\\{w^{e}\\}_{e=1}^{m}\subseteq\mathbb{R}^{p}_{++},\psi\in\mathbb{R}^{m}_{+}$
$\displaystyle\|\psi\|_{\infty}<C_{\psi},\psi^{1}=0,w^{e}\succeq w^{1}\text{
for }e=2,\dots,m~{}\&~{}\text{ for every }e=1,2,\dots,m$
$\displaystyle\Sigma^{e,\star}=(\mathcal{I}-{B})^{-1}(\texttt{diag}(w^{e})+(\psi^{e}+1)\Gamma\Gamma^{T})(\mathcal{I}-{B})^{-T}.$

######  Lemma 2.

Let $(\tilde{B},\tilde{\Gamma},\tilde{\psi},\\{\tilde{w}^{e}\\}_{e=1}^{m})$ be
an optimal solution of (15). The following statements hold:

  1. 1.

Suppose $\tilde{\psi}^{e}\neq\psi^{e,\star}$ for some $e\in\\{2,3\\}$. Under
Assumptions 1-4 in (10) or Assumptions 1 $\&$ 2’-4’ in (10) (main paper),
$\text{moral}(B^{\star})\subset\text{moral}(\tilde{B})$.

  2. 2.

Suppose $\tilde{\psi}^{e}=\psi^{e,\star}$ for $e=2,3$. Under Assumptions 1-4
in (10) or Assumptions 1 $\&$ 2’-4’ in (10) (main paper),
$\text{moral}(\tilde{B})=\text{moral}(B^{\star})$.

######  Lemma 3.

Let $(\tilde{B},\tilde{\Gamma},\tilde{\psi},\\{\tilde{w}^{e}\\}_{e=1}^{m})$ be
an optimal solution of (15). Suppose $\tilde{\psi}^{e}=\psi^{e,\star}$ for
$e=2,3$. Then, $\tilde{B}=B^{\star}$.

Combining Lemma 2 and 3 will conclude the proof of Theorem 1 (main paper), due
to the fact that
$(B^{\star},\Gamma^{\star},\psi^{\star},\\{w^{e,\star}\\}_{e=1}^{m})$ are
optimal solutions of (15). We now prove each lemma.

###  B.1 Useful notations

We introduce some notations that will be repeatedly used. Specifically, we
define for $e\in\mathcal{E}$:

$\displaystyle{\kappa}_{\text{cond}}^{e}$
$\displaystyle\equiv~{}~{}\min_{k,l}|(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star})|_{k,l}$
$\displaystyle~{}~{}\text{s.t.}~{}~{}\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash{\\{k,l\\}}}^{e},H^{e})\neq
0$ $\displaystyle{\kappa}_{\text{latent}}^{e}$
$\displaystyle\equiv\max_{k,l}~{}~{}|(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}({\Gamma^{\star}}^{T}\texttt{diag}(w^{e,\star})^{-1}\Gamma^{\star}+\frac{1}{1+\psi^{e,\star}}\mathcal{I})^{-1}$
$\displaystyle~{}~{}~{}~{}~{}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star})|_{k,l}$
$\displaystyle~{}~{}\text{s.t.}~{}~{}\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash{\\{k,l\\}}}^{e},H^{e})=0$

The intuition behind the quantities $\kappa_{\text{cond}}^{e}$ and
$\kappa_{\text{latent}}^{e}$ is based on the decomposition of
$(\Sigma^{e,\star})^{-1}$. Specifically, from the Woodbury inversion lemma:

$\displaystyle(\Sigma^{e,\star})^{-1}$ $\displaystyle=$
$\displaystyle(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star})$
$\displaystyle-$
$\displaystyle(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}({\Gamma^{\star}}^{T}\texttt{diag}(w^{e,\star})^{-1}\Gamma^{\star}+\frac{1}{1+\psi^{e,\star}}\mathcal{I})^{-1}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star}).$

Standard multivariate analysis states that for any pair of indices $(k,l)$
with $k\neq l$, $\left[(\Sigma^{e,\star})^{-1}\right]_{k,l}\neq 0$ if and only
if $\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash\\{k,l\\}}^{e})\neq 0$. Similarly,
since the precision matrix
$\text{cov}(X^{e}|H^{e})^{-1}=(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star})$,
we have that
$\left[(\mathcal{I}-B^{\star})^{T}\texttt{diag}(w^{e,\star})^{-1}(\mathcal{I}-B^{\star})\right]_{k,l}\neq
0$ if and only if
$\rho(X_{k}^{e},X_{l}^{e}|X_{\backslash\\{k,l\\}}^{e},H^{e})\neq 0$. Thus, by
definition, $\kappa_{\text{cond}}^{e}>0$ and by the _latent materiality_ in
Definition 1 (main paper), $\kappa_{\text{latent}}^{e}>0$. Then,

$\displaystyle\kappa_{\text{cond}}^{e}\geq\frac{(1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2})}{2\max_{k}w^{e,\star}_{k}}.$
(16)

Similarly, we have due to $\min_{k}w^{e,\star}_{k}\geq
8\|\Gamma^{\star}\|_{2}^{2}(1+C_{\psi})$ and $\min_{k}w^{e,\star}_{k}\geq
8\|w^{1,\star}\|_{\infty}$:

$\displaystyle\kappa_{\text{latent}}^{e}\geq\frac{8^{3}(1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2})(1+\psi^{e,\star})}{9^{3}(\max_{k}w^{e,\star}_{k})^{2}}.$
(17)

###  B.2 Proof of Lemma 1

###### Proof.

Let $\mathcal{M}(B,\Gamma,\psi^{e},w^{e})$ denote a model associated with each
equation in the SCM (2) (main paper). For notational convenience, we use the
short-hand notation $\mathcal{M}^{e}$ for this model. We let
$\Sigma(\mathcal{M}^{e})=(\mathcal{I}-\mathcal{F}_{\texttt{do}(e)^{c}}{B})^{-1}(\texttt{diag}(w^{e})+(1+\psi^{e})[\Gamma\Gamma^{T}]_{\texttt{do}(e)^{c}})(\mathcal{I}-\mathcal{F}_{\texttt{do}(e)^{c}}{B})^{-T}$
be the associated covariance model parameterized by the parameters
$(B,\Gamma,\psi^{e},w^{e})$. The optimal solution of the population
_DirectLikelihood_ can be equivalently reformulated as:

$\displaystyle\operatornamewithlimits{arg\,min}_{\\{\mathcal{M}^{e}\\}_{e=1}^{m}}$
$\displaystyle\sum_{e=1}^{m}\pi^{e,\star}\texttt{KL}(\Sigma^{\star,e},\Sigma(\mathcal{M}^{e})).$
(18)

Notice that for the decision variables
$\mathcal{M}^{e,\star}=(B^{\star},\Gamma^{\star},\psi^{e,\star},w^{e,\star})$
for each $e=1,2,\dots,m$, (18) achieves zero loss. Hence, any other optimal
solution of (18) must yield zero loss, or equivalently,
$\Sigma(\mathcal{M}^{e})=\Sigma^{e,\star}$ for any optimal collection
$\\{\mathcal{M}^{e}\\}_{e=1}^{m}$. ∎

###  B.3 Proof of Lemma 2

###### Proof.

We first provide the proof of Lemma 2 under Assumptions 1-4 in (10) (main
paper). Lemma 1 implies that for every $e=2,3$,

$\displaystyle\Sigma^{e,\star}-(1+\tilde{\psi}^{e})\Sigma^{1,\star}$ (19)
$\displaystyle=(\mathcal{I}-B^{\star})^{-1}\left(\texttt{diag}\left(w^{e,\star}-(1+\tilde{\psi}^{e})w^{1,\star}\right)+(\psi^{e,\star}-\tilde{\psi}^{e})\Gamma^{\star}{\Gamma^{\star}}^{T}\right)(\mathcal{I}-B^{\star})^{-T}$
$\displaystyle=(\mathcal{I}-\tilde{B})^{-1}\texttt{diag}\left(\tilde{w}^{e}-(1+\tilde{\psi}^{e})\tilde{w}^{1}\right)(\mathcal{I}-\tilde{B})^{-T}$

Since $\min_{k}w^{e,\star}_{k}\geq
C_{\psi}(\|w^{1,\star}\|_{\infty}+C_{\psi}\|\Gamma^{\star}\|_{2}^{2})$ from
Assumption 4 in (10) (main paper), we conclude that the matrix
$\Sigma^{e,\star}-(1+\tilde{\psi}^{e})\Sigma^{1,\star}$ is invertible for
$e=2,3$.

To establish the first component of Lemma 2, consider $e\in\\{2,3\\}$ for
which $\psi^{e,\star}\neq\tilde{\psi}^{e}$. After an inversion of (19) for
this environment, we obtain:

$\displaystyle(\mathcal{I}-B^{\star})^{T}({M}^{e}+{L}^{e})(\mathcal{I}-B^{\star})$
(20)
$\displaystyle=(\mathcal{I}-\tilde{B})^{T}\texttt{diag}\left(\tilde{w}^{e}-(1+\tilde{\psi}^{e})\tilde{w}^{1}\right)^{-1}(\mathcal{I}-\tilde{B})$

where,

$\displaystyle{M}^{e}=\texttt{diag}\left(w^{e,\star}-(1+\tilde{\psi}^{e})w^{1,\star}\right)^{-1}~{}~{};~{}~{}{L}^{e}={M}^{e}\Gamma^{\star}\left[{\Gamma^{\star}}^{T}{M}^{e}\Gamma^{\star}+\frac{1}{\Delta\psi^{e}}\mathcal{I}\right]^{-1}{\Gamma^{\star}}^{T}{M}^{e},$

Here, we have introduced a short-hand notation:
$\Delta\psi^{e}=(\psi^{e,\star}-\tilde{\psi}^{e})$. Notice that the nonzero
entries of
$(\mathcal{I}-\tilde{B})^{T}\allowbreak\texttt{diag}\allowbreak\left(\tilde{w}^{e}-(1+\tilde{\psi}^{e})\tilde{w}^{1}\right)^{-1}(\mathcal{I}-\tilde{B})$
encode the moral graph induced by $\tilde{B}$. Our strategy is to use
Assumptions 1-4 in (10) (main paper) to show
$(\mathcal{I}-B^{\star})^{T}(M^{e}+L^{e})(\mathcal{I}-B^{\star})$ has non-
zeros in the entries corresponding to the moral graph of $B^{\star}$ and at
least one nonzero outside of the moral graph. To conclude this, we consider
the following intermediate terms close to $M^{e}$ and $L^{e}$:

$\displaystyle\bar{{M}}^{e}=\texttt{diag}(w^{e,\star})^{-1}~{};~{}\bar{L}^{e}=\bar{M}^{e}\Gamma^{\star}\left[{\Gamma^{\star}}^{T}\bar{M}^{e}\Gamma^{\star}+\frac{1}{1+\psi^{e,\star}}\mathcal{I}\right]^{-1}{\Gamma^{\star}}^{T}\bar{M}^{e}$

Notice that:

$\|{\bar{M}}^{e}-M^{e}\|_{2}\leq\frac{5(1+C_{\psi})\|w^{1,\star}\|_{\infty}}{4(\min_{k}w^{e,\star}_{k})^{2}}~{}~{}~{};~{}~{}~{}\|M^{e}\|_{2}\leq\frac{5}{4(\min_{k}w^{e,\star}_{k})}$
(21)

where the inequalities follow by noting that
$5(1+C_{\psi})\|w^{1,\star}\|_{\infty}\leq\min_{k}w^{e,\star}_{k}$ from
Assumption 4 in (10) (main paper). Now let $(k,l)$ be any pair of indices
connected in the moral graph of $B^{\star}$. Then:

$\displaystyle|(\mathcal{I}-B^{\star})^{T}M^{e}(\mathcal{I}-B^{\star})|_{k,l}$
$\displaystyle\geq$
$\displaystyle|(\mathcal{I}-B^{\star})^{T}\bar{M}^{e}(\mathcal{I}-B^{\star})|_{k,l}$
(22) $\displaystyle-$
$\displaystyle(1+\max_{i}\|B^{\star}_{:,i}\|_{2})^{2}\|\bar{M}^{e}-M^{e}\|_{2}$
$\displaystyle\geq$
$\displaystyle{{\kappa}^{e}_{\text{cond}}}-\frac{5(1+C_{\psi})\|w^{1,\star}\|_{\infty}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})}{4(\min_{k}w_{k}^{e,\star})^{2}}$
$\displaystyle\geq$
$\displaystyle\frac{(1+\min_{i}\|B^{\star}_{:,i}\|_{2})^{2}}{2(\max_{k}w^{e,\star}_{k})}-\frac{5(1+C_{\psi})\|w^{1,\star}\|_{\infty}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})}{4(\min_{k}w_{k}^{e,\star})^{2}}$
$\displaystyle\geq$
$\displaystyle\frac{(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})}{4\max_{k}w^{e,\star}_{k}}.$

Here, the second to last inequality follows from the relation (16) and the
last inequality follows from
$\frac{\min_{k}w^{e,\star}_{k}}{\text{cond}(\texttt{diag}(w^{e,\star}))}\geq
5{\kappa^{\star}}(1+C_{\psi})\|w^{1,\star}\|_{\infty}$. Next, we control
$\|(\mathcal{I}-B^{\star})^{T}L^{e}(\mathcal{I}-B^{\star})\|_{\infty}$. Using
the inequality
$\left[{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}+\frac{1}{\Delta\psi^{e}}\mathcal{I}\right]^{-1}\preceq({\Delta\psi^{e}})\mathcal{I}$,
we have that:

$\displaystyle\|(\mathcal{I}-B^{\star})^{T}L^{e}(\mathcal{I}-B^{\star})\|_{\infty}\leq\frac{25C_{\psi}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})\|\Gamma^{\star}\|_{2}^{2}}{16(\min_{k}w^{e,\star}_{k})^{2}}$
(23)

Since
$\frac{\min_{k}w^{e,\star}_{k}}{\text{cond}(\texttt{diag}(w^{e,\star}))}\geq
7C_{\psi}\|\Gamma^{\star}\|_{2}^{2}$, comparing (23) and (22), we conclude
that for any indices $(k,l)$ connected in the moral graph of $B^{\star}$

$\displaystyle|(\mathcal{I}-B^{\star})^{T}M^{e}(\mathcal{I}-B^{\star})|_{k,l}>\|(\mathcal{I}-B^{\star})^{T}L^{e}(\mathcal{I}-B^{\star})\|_{\infty}.$

To finish the proof of the first assertion of Lemma 2, we have to show that
for indices $(k,l)$ attaining the optimum $\kappa_{\text{latent}}$,
$|(\mathcal{I}-B^{\star})^{T}L^{e}(\mathcal{I}-B^{\star})|_{k,l}>0$ or
equivalently,
$|(\mathcal{I}-B^{\star})^{T}(\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e})(\mathcal{I}-B^{\star})|_{k,l}>0$.
Notice that:

$\displaystyle\left|(\mathcal{I}-B^{\star})^{T}\left(\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}\right)(\mathcal{I}-B^{\star})\right|_{k,l}$
$\displaystyle\geq$
$\displaystyle|(\mathcal{I}-B^{\star})^{T}\bar{L}^{e}(\mathcal{I}-B^{\star})|_{k,l}$
$\displaystyle-$
$\displaystyle|(\mathcal{I}-B^{\star})^{T}\left(\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\bar{L}^{e}\right)(\mathcal{I}-B^{\star})|_{k,l}$
$\displaystyle\geq$
$\displaystyle\kappa_{\text{latent}}-\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\bar{L}^{e}\right\|_{2}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})$
$\displaystyle\geq$
$\displaystyle\frac{8^{3}(1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2})(1+\psi^{e,\star})}{9^{3}(\max_{k}w^{e,\star}_{k})^{2}}$
$\displaystyle-$
$\displaystyle\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\bar{L}^{e}\right\|_{2}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})$

Where the last inequality follows from the relation (17). Thus, it suffices to
show that:

$\displaystyle\frac{8^{3}(1+\min_{i}\|B^{\star}_{:,i}\|_{2}^{2})(1+\psi^{e,\star})}{9^{3}(\max_{k}w^{e,\star}_{k})^{2}}-\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}{L}^{e}-\bar{L}^{e}\right\|_{2}(1+\max_{i}\|B^{\star}_{:,i}\|_{2}^{2})>0.$
(24)

To that end, we control the term
$\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\bar{L}^{e}\right\|_{2}$.

$\displaystyle\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\bar{L}^{e}\right\|_{2}$
(25)
$\displaystyle\leq{2}\underbrace{\left\|(M^{e}-\bar{M}^{e})\Gamma^{\star}\left[\frac{\Delta\psi^{e}}{1+\psi^{e,\star}}{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}+\frac{1}{1+\psi^{e,\star}}\mathcal{I}\right]^{-1}{\Gamma^{\star}}^{T}{M}^{e}\right\|_{2}}_{\text{Term
1}}$
$\displaystyle+\underbrace{\left\|(M^{e}-\bar{M}^{e})\Gamma^{\star}\left[\frac{\Delta\psi^{e}}{1+\psi^{e,\star}}{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}+\frac{1}{1+\psi^{e,\star}}\mathcal{I}\right]^{-1}{\Gamma^{\star}}^{T}(M^{e}-\bar{M}^{e})\right\|_{2}}_{\text{Term
2}}$
$\displaystyle+\underbrace{\Bigg{\|}\bar{M}^{e}\Gamma^{\star}\Bigg{\\{}\left[\frac{\Delta\psi^{e}{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}+\mathcal{I}}{1+\psi^{e,\star}}\right]^{-1}-\left[{\Gamma^{\star}}^{T}\bar{M}^{e}\Gamma^{\star}+\frac{\mathcal{I}}{1+\psi^{e,\star}}\mathcal{I}\right]^{-1}\Bigg{\\}}{\Gamma^{\star}}^{T}\bar{M}^{e}\Bigg{\|}_{2}}_{\text{Term
3}}.$

We bound each of the individual terms in (25). Using the inequalities
$\left[{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}+\frac{1}{\Delta\psi^{e}}\mathcal{I}\right]^{-1}\preceq({\Delta\psi^{e}})\mathcal{I}$
and $\|\tilde{M}^{e}\|_{2}\leq\frac{1}{\min_{k}w_{k}^{e,\star}}$ and the
relation (21), Term 1 and Term 2 can be bounded as follows:

Term 1
$\displaystyle\leq\frac{10(1+C_{\psi})^{2}\|w^{1,\star}\|_{\infty}\|\Gamma^{\star}\|_{2}^{2}}{4(\min_{k}w^{e,\star}_{k})^{3}}$
(26) Term 2
$\displaystyle\leq\frac{25\|\Gamma^{\star}\|_{2}^{4}(1+C_{\psi})^{3}\|w^{1,\star}\|_{\infty}}{16(\min_{k}w^{e,\star}_{k})^{4}}.$

To bound Term 3, we use Taylor series expansion yielding

$\displaystyle(A+E)^{-1}-A^{-1}=A^{-1}\sum_{k=1}^{\infty}(EA^{-1})^{k}.$

Further, if $\|E\|_{2}\|A^{-1}\|_{2}<1$, we can bound the spectral norm of the
difference $(A+E)^{-1}-A^{-1}$ as follows:

$\displaystyle\|(A+E)^{-1}-A^{-1}\|_{2}\leq\|A^{-1}\|_{2}\sum_{k=1}^{\infty}\|E\|_{2}^{k}\|A^{-1}\|_{2}^{k}=\frac{\|A^{-1}\|_{2}^{2}\|E\|_{2}}{1-\|E\|_{2}\|A^{-1}\|_{2}}$
(27)

In the context of Term 3,
$E=\frac{\Delta\psi^{e}}{1+\psi^{e,\star}}{\Gamma^{\star}}^{T}M^{e}\Gamma^{\star}-{\Gamma^{\star}}^{T}\bar{M}^{e}\Gamma^{\star}$
and $A=(1+\psi^{e,\star})\mathcal{I}$. One can check that
$\|E\|_{2}\leq\frac{(\frac{5}{4}C_{\psi}+1)\|\Gamma^{\star}\|_{2}^{2}}{(\min_{k}w^{e,\star}_{k})^{2}}$.
Thus, employing the relation $(\min_{k}w_{k}^{e,\star})^{2}\geq
5(\frac{5}{4}C_{\psi}+1)\|\Gamma^{\star}\|_{2}^{2}$, we have that:

$\displaystyle\text{Term
3}\leq\frac{5\|\Gamma^{\star}\|_{2}^{4}(\frac{5}{4}C_{\psi}+1)}{4\min_{k}(w^{e,\star}_{k})^{3}}$
(28)

Combining the bounds in (26) and (28) with (25), we find that:

$\displaystyle\left\|\frac{1+\psi^{e,\star}}{\Delta\psi^{e}}L^{e}-\tilde{L}^{e}\right\|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{10(1+C_{\psi})^{2}\|w^{1,\star}\|_{\infty}\|\Gamma^{\star}\|_{2}^{2}}{4(\min_{k}w^{e,\star}_{k})^{3}}+\frac{25\|\Gamma^{\star}\|_{2}^{4}(1+C_{\psi})^{3}\|w^{1,\star}\|_{\infty}}{16(\min_{k}w^{e,\star}_{k})^{4}}$
$\displaystyle+$
$\displaystyle\frac{5\|\Gamma^{\star}\|_{2}^{4}(\frac{5}{4}C_{\psi}+1)}{4\min_{k}(w^{e,\star}_{k})^{3}}$
$\displaystyle\leq$
$\displaystyle\left(\frac{10}{4}+\frac{25}{16}+\frac{5}{4}\right)\frac{(1+2C_{\psi})^{2}\max\\{\|\Gamma^{\star}\|_{2}^{2},\|\Gamma^{\star}\|_{2}^{4}\\}\max\\{1,\|w^{1,\star}\|_{\infty}\\}}{(\min_{k}w^{e,\star}_{k})^{3}},$

where the second inequality follows from
$\min_{k}w^{e,\star}_{k}\geq\|w^{1,\star}\|_{\infty}(1+C_{\psi})$. Thus, since
$\frac{\min_{k}w^{e,\star}_{k}}{\text{cond}(\texttt{diag}(w^{e,\star}))}\geq
8{\kappa^{\star}}(1+C_{\psi})^{2}\max\\{\|\Gamma^{\star}\|_{2}^{2},\|\Gamma^{\star}\|_{2}^{4}\\}\max\\{1,\|w^{1,\star}\|_{\infty}\\}$,
the sufficient condition in (24) is satisfied. This concludes that if
$\tilde{\psi}^{e}\neq\psi^{e,\star}$ for $e\in\\{2,3\\}$,
$(\mathcal{I}-B^{\star})^{T}(M^{e}+L^{e})(\mathcal{I}-B^{\star})$ will have a
non-zero outside of the moral graph of $B^{\star}$ and thus according to (20),
$\text{moral}(B^{\star})\subset\text{moral}(\tilde{B})$. We have established
the first component of Lemma 2. The second component (where
$\tilde{\psi}^{e}=\psi^{e,\star}$ for $e=2,3$) follows from (19).

We next provide a proof of Lemma 2 under Assumptions 1 $\&$ 2’-4’ in (10)
(main paper). Lemma 1 implies the following relations:

$\displaystyle(\mathcal{I}-B^{\star})^{-1}\left(\texttt{diag}\left(w^{3,\star}-(1+\tilde{\psi}^{3})w^{1,\star}\right)+(\psi^{3,\star}-\tilde{\psi}^{3})\Gamma^{\star}{\Gamma^{\star}}^{T}\right)(\mathcal{I}-B^{\star})^{-T}$
(29)
$\displaystyle=(\mathcal{I}-\tilde{B})^{-1}\texttt{diag}\left(\tilde{w}^{3}-(1+\tilde{\psi}^{3})\tilde{w}^{1}\right)(\mathcal{I}-\tilde{B})^{-T}$
$\displaystyle(\mathcal{I}-B^{\star})^{-1}\Bigg{(}\texttt{diag}\left(w^{3,\star}-\frac{1+\tilde{\psi}^{3}}{1+\tilde{\psi}^{2}}w^{2,\star}\right)$
$\displaystyle+\left(1+\psi^{3,\star}-\frac{(1+\psi^{2,\star})(1+\tilde{\psi}^{3})}{1+\tilde{\psi}^{2}}\right)\Gamma^{\star}{\Gamma^{\star}}^{T}\Bigg{)}(\mathcal{I}-B^{\star})^{-T}$
$\displaystyle=(\mathcal{I}-\tilde{B})^{-1}\texttt{diag}\left(\tilde{w}^{3}-\frac{(1+\tilde{\psi}^{3})}{1+\tilde{\psi}^{2}}\tilde{w}^{2}\right)(\mathcal{I}-\tilde{B})^{-T}$

Using the relation (29) and a similar analysis as with the proof under
Assumptions 1-4 in (10) (main paper), one can arrive at the conclusion of
Lemma 2 with Assumptions 1 $\&$ 2’-4’ in (10) (main paper). ∎

###  B.4 Proof of Lemma 3

###### Proof.

The proof technique of this lemma is similar in spirit to the proof of Theorem
1 in [27]. We consider the setup with Assumptions 1-4 in (10) (main paper) and
for brevity, leave out the proof with Assumptions 1 $\&$ 2’-4’ in (10) (main
paper). We have from (19) that for $e=2,3$:

$\displaystyle(\mathcal{I}-B^{\star})^{-1}\texttt{diag}\left(w^{e,\star}-{(1+{\psi}^{e,\star})}w^{1,\star}\right)(\mathcal{I}-B^{\star})^{-T}$
(30)
$\displaystyle=(\mathcal{I}-\tilde{B})^{-1}\texttt{diag}\left(\tilde{w}^{e}-{(1+{\psi}^{e,\star})}\tilde{w}^{1}\right)(\mathcal{I}-\tilde{B})^{-T}$

From relation (30), we have:

$\displaystyle(\mathcal{I}-\tilde{B})(\mathcal{I}-B^{\star})^{-1}\texttt{diag}\left(w^{2,\star}-{(1+{\psi}^{e,\star})}w^{1,\star}\right)(\mathcal{I}-B^{\star})^{-T}(\mathcal{I}-\tilde{B})^{T}$
$\displaystyle=$
$\displaystyle\texttt{diag}\left(\tilde{w}^{2}-{(1+{\psi}^{e,\star})}\tilde{w}^{1}\right)$
$\displaystyle(\mathcal{I}-\tilde{B})(\mathcal{I}-B^{\star})^{-1}\texttt{diag}\left(w^{3,\star}-{(1+{\psi}^{e,\star})}w^{1,\star}\right)(\mathcal{I}-B^{\star})^{-T}(\mathcal{I}-\tilde{B})^{T}$
$\displaystyle=$
$\displaystyle\texttt{diag}\left(\tilde{w}^{3}-{(1+{\psi}^{e,\star})}\tilde{w}^{1}\right)$

Let
$\phi^{e}_{k}:=\left[(\mathcal{I}-\tilde{B})(\mathcal{I}-B^{\star})^{-1}\texttt{diag}(w^{e,\star}-w^{1,\star}{(1+{\psi}^{e,\star})}\right]_{k,}$
for any $k=1,2,\dots,p$. Let
$\xi_{k}:=\left[(\mathcal{I}-\tilde{B})(\mathcal{I}-B^{\star})^{-1}\right]_{k,}$.
Then

$\displaystyle\phi_{k}^{e}\perp\xi_{l}~{}~{}\text{for any }k\neq l.$ (31)

Notice that for any $k=1,2,\dots,p$, $\\{\xi_{l}\\}_{l\neq k}$ are linearly
independent. The condition above means that $\phi_{k}^{2}$ and $\phi_{k}^{3}$
(where neither would be exactly a zero vector because of Assumption 4 in (10)
(main paper) ensuring that
$w^{2,\star}-w^{1,\star}{(1+{\psi}^{e,\star})},w^{3,\star}-w^{1,\star}{(1+{\psi}^{e,\star})}\neq
0$) live inside the one-dimensional null-space of the matrix formed by
concatenating the vectors $\\{\xi_{l}\\}_{l\neq k}$. In particular, for every
$k$, we have that for some constant $c\neq 0$:
$\xi_{k}\texttt{diag}(w^{2,\star}-w^{1,\star}{(1+{\psi}^{e,\star})}))=c\xi_{k}\texttt{diag}(w^{3,\star}-w^{1,\star}{(1+{\psi}^{e,\star})}))$.
It is straightforward to check that Assumption 2 in (10) (main paper)