Estimating Causal Effects with Hidden Confounding using Instrumental Variables and Environments

Recent works have proposed regression models which are invariant across data collection environments. These estimators often have a causal interpretation under conditions on the environments and type of invariance imposed. One recent example, the Causal Dantzig (CD), is consistent under hidden confounding and represents an alternative to classical instrumental variable estimators such as Two Stage Least Squares (TSLS). In this work we derive the CD as a generalized method of moments (GMM) estimator. The GMM representation leads to several practical results, including 1) creation of the Generalized Causal Dantzig (GCD) estimator which can be applied to problems with continuous environments where the CD cannot be fit 2) a Hybrid (GCD-TSLS combination) estimator which has properties superior to GCD or TSLS alone 3) straightforward asymptotic results for all methods using GMM theory. We compare the CD, GCD, TSLS, and Hybrid estimators in simulations and an application to a Flow Cytometry data set. The newly proposed GCD and Hybrid estimators have superior performance to existing methods in many settings.


Introduction
Causal inference is challenging because of confounding and reverse causality.One solution is to make strong assumptions about confounding (e.g.no hidden confounding) and the direction of causation (e.g.X is a cause of Y and not the reverse).Under these assumptions, causal parameters may be identifiable from observational data.
When these assumptions are not valid, instrumental variables (IV) are a classical method for identifying causal effects.The variable E is an instrument for the X → Y causal relation if 1) E is uncorrelated with the error term in the Y on X regression and 2) E is correlated with X (valid first stage).IV estimators such as Two Stage Least Squares (TSLS) remain consistent under hidden confounding and unknown direction of causality.IV methods date back to [28] and have more recently been generalized to high dimensional problems [14,9,2] and causal discovery applications where X ∈ R p is a vector and identifying the causes of each X i is of interest [6].
Recently, several works have proposed causal estimators based on the concept of data collection environment [24,20,11,16,8].[20] introduced the concept of data collection environment and developed a causal estimator, Invariant Causal Prediction (ICP).In this framework, each observation is collected in an environment.The environment may represent randomized experiments on some of the exposures of interest, shift, and/or noise interventions.Environments are typically discrete and often small in number (e.g. 2 or 3).
Estimators are constructed from environments based on the principle that parameters in a causal regression model Y on X should be invariant across environments while parameters in a merely associational model will vary.As a simple heuristic example, suppose we are interested in estimating the causal effect of X on Y .In truth Y is a cause of X and the true causal effect of X on Y is 0. Standard regression based estimators are inconsistent.More generally with purely observational data it will be impossible to determine the causal effect.However if we have data from two environments, e.g. an observational environment and an interventional environment where noise is added to X, then it is possible (under some conditions) to infer the causal effect as 0 by noting that the distribution of Y is identical in the observational and interventional environment.This would not be the case if X has a causal effect on Y .
The original environment estimator ICP has been generalized to problems with sequential data and non-linear models [21,11].[24] proposed the Causal Dantzig (CD) environment estimator to address two weaknesses of ICP: computational complexity and inconsistency under hidden confounding.While ICP requires fitting models on all subsets of the exposure variables (making the algorithm superexponential in the number of exposures), the CD estimator has computational burden similar to linear regression.Further, like instrumental variable estimators, the CD is consistent when hidden variables confound the X → Y causal relation.
In this work, we show that the Causal Dantzig can be represented as a generalized method of moments estimator (GMM).This immediately leads to a new estimator, termed the Generalized Causal Dantzig (GCD), which is equivalent to the CD in the two environment case but can be applied with continuous environments, a setting not handled by the original CD.The GMM representation facilitates straightforward asymptotic results based on GMM theory.GMM theory shows how to optimally weight the moment constraints in over-identified problems, which occurs whenever there are more than two data collection environments.Finally, the GMM representation of the CD facilitates construction of a Hybrid estimator which uses both the CD and IV moment constraints for estimating parameters.This Hybrid estimator is consistent in some settings where neither the CD nor TSLS are consistent.This work is organized as follows.In Section 2 we review the GMM representation of IV estimators and the environment invariance representation of the Causal Dantzig.In Section 3 we propose the Generalized Causal Dantzig (GCD), a GMM estimator, and construct a Hybrid estimator which use both IV and CD/GCD moment constraints.Section 4 derives asymptotic results for the GCD based on the GMM representation of the estimator.In Section 5, we assess consistency of the estimators in causal Structural Equation Models.Section 6 contains simulations which demonstrate some of the potential applications of the GCD and hybrid GCD-IV estimators.In Section 7 we apply IV, CD, GCD, and Hybrid estimators to Flow Cytometry data of [26].In several cases, we show that IV and Hybrid estimators identify more plausible causal relations than the CD alone.We conclude with a discussion in Section 8.All code and data for reproducing the computational aspects of this work is available.1

Instrumental Variables and the Causal Dantzig
Let X ∈ R p be a set of endogenous exposures and Y ∈ R 1 be a response.The goal is to estimate the causal effect of X on Y .Consider a linear model of the form where E[δ Y ] = 0.Under a potential outcomes [1] or a structural equation modelling [19] framework, β 0j (j th element of β 0 ) can be given a causal interpretation as the average treatment effect (ATE) of X j on Y when shifting X j by 1 unit.In either of these frameworks, correlation between X and δ Y is induced by hidden confounders which exert a causal effect on X and Y .Straightforward regression of Y on X may result in inconsistent estimates of β 0 when the error term δ Y is correlated with X.
In this section we review two approaches to constructing consistent estimators with hidden confounding and reverse causality: the classical Two Stage Least Squares (TSLS) which uses Instrumental variables (IV) and the recently proposed Causal Dantzig which uses data collection environments.

Instrumental Variable Estimators
Instrumental variables techniques, dating back to [28], use instrumental variables (IVs) E ∈ R q to construct consistent estimates of β in the presence of hidden confounding.Suppose that 1) the instruments E are uncorrelated with the error term δ The latter condition implies that q ≥ p (i.e.there are at least as many instruments as exposures).We review the construction of IV estimators from a GMM perspective.See [15] for more background.
Let Z = (Y, X, E) and g IV (Z, β) = E(Y − X T β).Then the true causal parameter β 0 is the unique solution to This can be seen by noting The rank condition implies that the null-space of the matrix is 0 implying β 0 is the only solution.To construct a consistent estimator, the expectation is approximated with a sample.Define X ∈ R n×p , E ∈ R n×q , Y ∈ R n×1 to be a matrices of n i.i.d.observations.With i indexing observations, we have When q > p, the model is over-identified and there will typically be no β such that m IV ( β) = 0 in Equation (2.3).In this case, the standard GMM approach is to use estimator where W ≻ 0 is a positive definite weighting matrix.The TSLS IV estimator uses This weight matrix is chosen for asymptotic efficiency considerations which we discuss further in Section 4.2 (see also Section 1.3.4.2 of [15]).With W , Equation (2.4) has the form where X = E(E T E) −1 E T X.The TSLS estimator derives its name from the fact that it is computed by first regressing X on E (first stage) and then regressing Y on the predicted values from the first stage (second stage).Note that when p = q (just identified case) the unique β IV ( W ) does not depend on W and has the form

Causal Dantzig
The Causal Dantzig (CD) uses environments to estimate β 0 in Equation (2.1).Each observation belongs to one of a discrete set of environments.Let E = (1, 2, . ..) be a set of data collection environments with #E denoting the number of environments (at least 2).Let X e (Y e ) denote exposures (response) collected in environment e ∈ E. For any f, g ∈ E, the CD seeks a β which satisfies Solving for β CD one obtains assuming the inverse exists (see Equation 7of [24]).Equation (2.7) shows that the CD exploits how the covariance structure of X changes with the environment to construct a consistent estimate of β.This is in contrast to TSLS which exploits how the mean of X changes with the instrument E. The CD is a consistent and asymptotically normal estimator of β 0 under conditions specified in Sections 4 and 5.

New Estimators
We now construct the Generalized Causal Dantzig (GCD) and a Hybrid estimator.The GCD extends the CD to work in new settings (e.g.continuous environments).The Hybrid estimator imposes both TSLS and CD moment constraints to produce an estimator with potentially better asymptotic properties than either the CD or TSLS.Finally we discuss connections between the GCD and estimators proposed by Lewbel [13].

Generalized Causal Dantzig
Following notation used for the IV estimators in Section 2.1 where Z = (Y, X, E), define the Generalized Causal Dantzig (GCD) using the moment conditions Here EX T ∈ R q×p and vec(EX T ) vectorizes (column stacks) the matrix EX T [12].Specifically We See Section 9.2 for a proof.The result has conceptual and practical implications.On the conceptual side, comparing the IV moment constraints (Equation (2.2)) with the GCD moments constraints (Equation (3.1)) shows the environment and instrument (E) play similar roles in the estimators.TSLS requires each instrument/environment to be orthogonal to Y − X T β while CD requires each instrument/environment be orthogonal to each element of X(Y − X T β).On a practical side, the GCD provides a natural generalization of the CD to problems with continuous environments/instruments and GMM theory enables straightforward derivation of GCD asymptotics.This includes optimal weighting of invariance (equivalently moment) criteria in the overidentified case using a two-step estimator.We explore these ideas further in Sections 4 and 6.
We now construct the GCD estimator by imposing the constraints in Equation (3.1) on the sample.Note that vec(EX T ) ∈ R qp induces qp constraints on β.Define Define the GCD moment equations as and the resulting GCD estimator as where W ≻ 0 is some weighting matrix.When E is constructed from two environments, its dimension is 1 because q = #E − 1 = 1.In this case the qp = p constraints just identify β.Further β GCD ( W ) is invariant to different choices of W and identical to the original Causal Dantzig estimator β CD .The following theorem formalizes these results.
Theorem 3.2.Suppose W ≻ 0 and X T (E • X) ∈ R p×pq has column rank p. Then

The unique minimizer in Equation
See Section 9.3 for a proof.With more than two environments, β GCD ( W ) will generally not be equivalent to the #E > 2 CD estimators proposed in [24].[24] proposed two approaches for adapting the CD to the case with more than two environments: 1) Merge environments to obtain two distinct environments.Fit the two environment CD estimator on the resulting data.2) Fit a minimax estimator which seeks to satisfy the CD invariance constraints in a one-versus-all environment approach (Equation 8 in [24]).No theory is given to guide merging of environments (if approach 1 is taken) or how to compute uncertainties on parameter estimates (if approach 2 is taken).The CD was previously applied to the Flow Cytometry Application studied in Section 6 of this work [17].In that application, a third strategy was used for adapting the CD to more than two environments which involves iteratively fitting the two environment CD to pairs of environments (see Section 6.2 of this work and [17] for a more detailed description).In contrast, the GCD with more than two environments uses the matrix W to weight moment conditions.GMM theory is then used to derive the optimal weight matrix (minimizes asymptotic variance) which can be estimated using standard two-step procedures.GMM theory provides (asymptotically) valid uncertainty estimates.This approach is used both in the simulations and the Flow Cytometry application of Section 6.We discuss optimal selection of W in Section 4.2.

Hybrid Estimator
The GMM representation of the GCD in Equation 3.1 and the GMM representation of the TSLS estimator in Equation 2.2 share a similar structure.The GCD enforces orthogonality between Y − X T β and vec(EX T ) while TSLS enforces orthogonality between Y − X T β and E. Using GMMs, it is straightforward to construct estimators which enforce both of these constraints.We define the Hybrid estimator moment conditions This estimator enforces orthogonality between Y − X T β and both vec(EX T ) and E, resulting in a total of q + qp constraints.Following standard GMM practice, the sample version of these orthogonality constraints is where W ∈ R (q+qp)×(q+qp) ≻ 0 is some weighting matrix.Since β IV , β GCD , and β H are all GMMs, asymptotic results for these estimators can be derived using GMM theory.We do this in the following section.

Relation to Lewbel [13]
Lewbel [13] constructed consistent estimators of causal effects with hidden confounding by exploiting heteroscedasticity in endogenous variables.The relationship between Lewbel's estimators and the CD was briefly discussed in the original CD paper [24].[24] claimed that the methods are different because Lewbel uses exogenous variables while the CD directly exploits endogenous covariance structure, "resulting in a different method."Complicating comparison of the methods is the fact that the original CD was constructed based on invariance and environments (Equation (2.6)) while Lewbel used GMM.
Our result that the CD can be represented as a GMM (Theorem (3.2)), facilitates comparison between the methods.In fact, for certain choices of Lewbel's variable Z, the CD and Lewbel are identical.See Supplementary Material Section 9.1 for details.Neither Lewbel's estimators nor the CD/GCD/Hybrid completely overlap in the cases they consider.Lewbel considers fully simultaneous systems and models additional exogenous variable which are not part of this work.Lewbel primarily considers univariate endogenous variable models, briefly developing an extension of the GMM to the p = 2 case (Section 3.3).Here we study the p > 1 case in depth.The Hybrid estimator is to our knowledge completely new and not considered in Lewbel.

Asymptotic Properties of Estimators
The GMM representation of the IV, GCD, and Hybrid estimators makes derivation of asymptotic properties straightforward.We first discuss consistency and then asymptotic normality.

Consistency
Assumptions 1 and Theorem 4.1 below may be specialized to the IV, GCD, or Hybrid estimators by substituting in the corresponding m, m, and β( W ). For example, for the GCD, let m = m GCD , m = m GCD , and β( W ) = β GCD ( W ).

Assumptions 1. Suppose (a) m(β) exists and is finite for all
See Theorem 1.1 and Section 1.3.4.1 of [15] for a proof and application to the i.i.d.case.Assumption 1 d) is the simplest to satisfy.For example, W = W = I satisfies the assumption.In general W will be chosen to be a consistent estimator of a W which minimizes the asymptotic variance of the estimator, as discussed in Section 4.2.Assumption 1 a) will generally hold whenever error terms have sufficient moments.While Assumption 1 c) is strong, it may be replaced with assumptions restricting β to a compact subset of R p and continuity of m(β) [10,18].Assumption 1 b) is closely related to the concepts of instrument/environment validity.We give conditions under which Causal Structural Equation Models (SEM) will satisfy Assumption 1 b) in Section 9.5.These results clarify the strengths and weaknesses of the IV, GCD, and Hybrid estimators.

Asymptotic Normality
Asymptotic normality of the GCD and Hybrid estimators can be derived from standard GMM theory.The asymptotic theory provides guidance on selection of the weight matrix W in the overidentified case (recall that overidentification will occur for the GCD whenever there are more than two environments or more generally when the dimension of E is greater than 1).The optimal weight matrix is chosen to minimize the asymptotic variance.For the GCD, this optimal weight matrix can be estimated using a GMM two-step procedure.
Assumptions 2 below are used for showing asymptotic normality of GMM estimators.These may be specialized to the IV, GCD, or Hybrid estimators by substituting in the corresponding g, m, m, and β.For example, for the GCD, and suppose β → P β 0 .Suppose there exists an M ∈ R k×p of full column rank such that M ( β) → P M.

Theorem 4.2. Under Assumptions 1 and 2,
with asymptotic variance See Section 1.3.4.1 of [15] for these results.The asymptotic variance is minimized by setting W = W ∝ V −1 which results in We now review weighting for the TSLS estimator.For the instrumental variables estimator is a consistent estimator of W and thus is an asymptotically optimal weighting for the IV estimator.The estimator β IV ( W ) can be computed using a set of two regressions, giving it the name Two Stage Least Squares.Since consistency of β IV only requires orthogonality of E and δ Y , there may be cases where β IV (W ) is a consistent, asymptotically normal estimate of β but W = E[EE T ] −1 is not asymptotically efficient.In these cases, an asymptotically efficient IV estimator can be constructed assuming one can construct a consistent estimate of W .One possibility is to use two-step procedures.First a pilot estimator β IV (W ) is computed using some initial weight matrix W . Possible choices include W = I (identity) or W = W (TSLS weight).Residuals are defined as Define Σ as a diagonal matrix with Σ ii = δ 2 Y,i and ) is asymptotically efficient.This two-step procedure is sometimes referred to as optimal GMM.See Section 6.4.2 of [4] for a discussion of two-step optimal GMM estimators and comparison with TSLS.
For the GCD, the asymptotically optimal estimators can be constructed from two-step procedures, following the same strategy as used for IV.The following theorem provides formal justification of this approach.Note that this is only necessary in the over-identified case since in the just identified case the estimator is invariant to different choices of W . Theorem 4.3.Suppose Assumptions 1 and 2 hold for the GCD estimator.Let W ≻ 0 be some initial weight matrix such that Let Σ be a diagonal matrix with Σ ii = δ 2 Y,i and define .
Then β GCD ( W GCD ) is asymptotically efficient with variance specified in Equation (4.1).
See Section 9.4 for a proof.For the GCD, the asymptotic variance can be estimated by first estimating M with

Causal Models and Consistency
Recall that Assumption 1 b) states We now consider Causal Structural Equation Models (SEM) which guarantee that Assumption 1 b) holds.In combination with regularity Assumptions 1 a), c), and d), these will ensure consistency of the estimators by Theorem 4.1.We also contrast these estimators and assumptions with Independent Instrumental Variable (IIV) methods.

Assumptions 3 (Causal SEM).
Let Z = (Y, X, E) be generated from a Causal SEM with independent exogenous variables ϵ E , ϵ H , ϵ X , ϵ Y and endogenous variables X ∈ R p and Y ∈ R 1 : where The corresponding Directed Acyclic Graph (DAG) is show in Figure 1 a).The dashed circle around H indicates that it is not measured and thus a hidden confounder.For GCD and Hybrid consistency, we need an additional assumption.
Assumptions 4 (Error Decomposition).We discuss this decomposition assumption further in Section 5.1.We have the following result.
X T has column rank p.
See Section 9.2 for a proof.Column rank conditions in Theorem 5.1 depend on observed random variables X and E and are empirically verifiable.The structure of f X dictates whether these column rank conditions hold and which estimator (IV, GCD, or Hybrid) will be most appropriate for a given problem.We discuss some specific cases now.For notational convenience let The IV estimator leverages how the instruments change the mean of X.If r(E) j = 0 for any j ∈ R p , then IV is not consistent.Further if q = dim(E) < dim(X) = p, then the IV is not consistent because the column rank of E[EX T ] is bounded by q < p.Thus with a single instrument, IV may only estimate the causal effect of a single exposure on the response.• GCD: Suppose r(E) = 0 so the IV estimator is not consistent.The GCD rank condition may be rewritten as The GCD leverages how E shifts higher moments of X (δ 1 X (E, ϵ X )).The GCD does not require q ≥ p because E[vec(EX T )X T ] ∈ R pq×p .With a single binary instrument/environment taking values s 1 or s 0 , the GCD will be consistent if • Hybrid: Consistency of the Hybrid estimator is weaker than for the GCD since The Hybrid estimator leverages how E shifts the mean (via the IV constraints) or higher moments (via the GCD constraints) of X.
Theorem 5.1 shows consistency of estimators under assumptions on the data generating model.In some settings, particular elements of the estimator may be consistent while others are not, e.g.β j ( W ) → P β 0j for some but not all j ∈ {1, . . ., p}.Such partial identifiability/consistency results have been derived for IV estimators under weaker conditions than considered here [25,20].Detailed consideration of these cases for the GCD and Hybrid estimators is beyond the scope of this work.

Decomposition Condition and Do Operators
Assumption 4 will hold for many models including an additive hidden confounder model where δ 2 X (H, ϵ X ) = γ(H) for some function γ.The assumption will not hold when the instrument/environment E represents a do operator.Consider the structural model for X for some f X (H, ϵ X ).When E = −1, the data is observational and generated from the DAG in Figure 1 b).When E = 1, the data is interventional and generated from the DAG in Figure 1 c).For this model the error term is and cannot be decomposed according to Assumption 4. Inconsistency of the CD under do interventions was discussed in [24].Since the Hybrid estimator also requires the error decomposition assumption, it will also be inconsistent for instruments/environments which are used to represent do operators.

Independent Instrumental Variables
Independent instrumental variable (IIV) models assume that δ Y ⊥ ⊥ E. This can be leveraged to construct consistent estimates of β 0 even when r(E) = 0, much like the GCD [27,7,22].Note that δ Y ⊥ ⊥ E is equivalent to for all bounded, continuous functions η and ϕ.Thus estimates of β 0 may be constructed by finding β which satisfy empirical versions of Equation (5.1) (up to sampling variability) for a large set of η and ϕ.See [22] for specific implementations of this approach.Assumptions 3 do not require δ Y ⊥ ⊥ E because δ Y (H, X, ϵ X ) is a function of X which is a function of E. One common setting where Assumptions 3 holds but δ Y ̸⊥ ⊥ E is heteroskedastic response models such as where s : R p → R + and E[ϵ Y ] = 0.Here IIV models should not be used.

Simulations
We demonstrate some of the applications of the GCD and Hybrid estimators in simulations.

Continuous Environments
We fit the GCD to a model with continuous environments/instruments.
where exogenous variables ϵ X and ϵ Y are standard normal.The true causal parameter is β 0 = 1.We simulate n = 100 samples.IV is not consistent for this model because E shifts the variance of X (larger E implies larger variance for X), but not the mean of X.With univariate E and X the GCD is The CD is not directly applicable here because the environment is continuous.
One could discretize E with some function e : R → {0, 1} and then apply the CD using environment e(E).We consider this approach with e(E i ) = N = 1000 times and plot sampling distributions of the Generalized Causal Dantzig (GCD), the Causal Dantzig (CD), and Ordinary Least Squares (OLS) in Figure 2. Note that E is centered (mean shifted to 0 in the sample) before the GCD is fit.OLS is inconsistent due to hidden confounding.The GCD and CD sampling distributions are both centered at the true causal effect of 1.The GCD empirical sampling distribution is more concentrated around the causal effect.The GCD is also easier to fit because it does not require selection of the function e to binarize the continuous variable E into discrete environments.

Overidentified Models
The CD/GCD constraints overidentify β 0 whenever there is more than one instrument/environment.We consider data generating model: All exogenous variables (H, ϵ 1 , ϵ 2 , ϵ 3 , ϵ Y ) are standard normal.The true parameter value is β 0 = (0, 1, 0) because only X 2 has a causal effect on Y .The hidden confounder H will cause OLS to be inconsistent.The sample size is n = 200.We simulate N = 500 runs.We consider four estimators: GCD (using both E 1 and E 2 ), GCDE1 (uses only environment E 1 ), GCDE2 (uses only environment E 2 ), and OLS (ordinary least squares).Note that GCDE1 is equivalent to the CD using E 1 to indicate the environment of the observation.GCDE2 is not equivalent to a CD estimator because E 2 is a continuous environment.For the GCD, the two-step estimator is used since the two environments overidentify β 0 .For initial GCD weight matrix we use Table 1 shows the empirical coverage probabilities of 95% confidence intervals for each parameter and the median CI width.The rows of the table correspond to different estimators.OLS is inconsistent because of hidden confounding.This results in the coverage probabilities being well below nominal levels (0 in the case of β 2 ).The three GCD methods (GCD, GCDE1, and GCDE2) all have empirical coverage probabilities near or above 95%.However GCDE1 and GCDE2 obtain this coverage by producing extremely wide confidence intervals.For example, the median CI width for GCDE2 for β 2 is 10 times that for GCD.Similarly, the median width of the β 3 CI for GCDE1 is about 10 times the median width for GCD.By only using the information in one of the environment variables, GCDE1 and GCDE2 produce highly uncertain estimators with very wide confidence intervals.

Hybrid Estimator
Consider a model in which the instrument E shifts the mean and variance of X: IV is inconsistent when R = 0.It will be a poor estimator when R is near 0 because the instrument is weak.In similar fashion GCD is inconsistent when α v = 0.It will be a poor estimator when α v is near 0. The Hybrid estimator is consistent whenever IV or GCD is consistent because it can leverage changes in the mean or variance induced by the instrument/environment E.Here we consider two simulation settings: In Model 1, R = 5 and α v = 1 so that E has a strong effect on the mean of X and only a weak effect on the variance.In Model 2, R = 1 and α v = 5 so that E has a weak effect on the mean of X and a strong effect on the variance.We fit IV, GCD, and the Hybrid estimator on these two models.IV and GCD do not require specification of a weight matrix because the number of parameters equals the number of constraints.For the Hybrid estimator, we perform a two step procedure to estimate the optimal weight matrix, using for an initial weighting.We let β 0 = 1 and simulate n = 200 samples.Empirical sampling distributions with N = 1000 simulations are shown in Figure 3.The Hybrid estimator performs well for both models while IV and the GCD each only perform well for one of the models.

Application to Flow Cytometry Data
[26] measured the abundance of 11 biochemical agents in thousands of cells using flow cytometry.These data were collected under several conditions (or environments) in which external reagents were added to the system.Each reagent has the effect of stimulating or inhibiting particular agents in the system.Five conditions used in this work (1 observational and 4 interventional) are described in Table 2.In the observational condition, only a general perturbation (CD3+CD28) is applied.For the observational condition, the expression of the 11 biochemical agents was measured in 853 cells.In Condition 3, the reagent Psitectorignin, an inhibitor of PIP2 was added to the system in addition to the general perturbation.The effect of this perturbation should reduce the abundance of PIP2 (one of the 11 measured agents) as well as alter the abundance of any agents which PIP2 itself effects.[17] fit the CD to this data to infer a causal signalling network (graph) among the 11 agents.Here we compare the performance of the Causal Dantzig with IV, GCD, and Hybrid estimators.Before fitting any models, we hyperbolic arcsine transform the data.This technique is used to approximately normalize the flow cytometry data to better satisfy modelling assumptions and reduce the influence of outliers [23].

Univariate Analysis
We first consider the problem of determining whether a particular agent X has a total causal effect on another agent Y .In general, simple regression of Y on X will not consistently estimate the total causal effect of X on Y because hidden confounding and reverse causality will induce an association between X and Y even when X has no causal effect on Y .To address this problem, we consider (X, Y ) data from two environments: an observational environment involving only a general system perturbation and an interventional environment which includes the general perturbation plus a reagent designed to perturb X. Presence or absence of the additional reagent is modeled using an instrument (environment) variable E. We focus on two causes (different X variables), PIP2 and MEK, in order to demonstrate similarities and differences in IV and CD modeling results.

PIP2
Using the Causal Dantzig, [17] (Figure 3) did not find that PIP2 is a direct cause of changes in any of the other 10 biochemical agents.This implies that PIP2 should not have a total effect on any of the agents in the system.To investigate this, we consider abundance measures from cells in two conditions: observational (the general perturbation only) and Condition 3 (Psitectorignin plus general perturbation).The condition is treated as a binary instrumental variable / environment.This is justified by the fact that Psitectorignin is meant to directly target PIP2 and any effects on other agents should occur by way of PIP2.Since the CD and GCD are identical in this case, we compare the CD, IV, and Hybrid estimators.decreasing its mean.Plcg and PIP3 abundances are also strongly influenced by the intervention.This suggests that PIP2 is a cause (either directly or possibly indirectly through other agents) of both Plcg and PIP3.Note that the intervention primarily effects the mean, rather than the variance, of PIP2.Thus the CD may struggle to identify an effect because it is sensitive to variance, not mean, shifts.In contrast, IV is better suited to settings where the instrument/environment effects the exposure mean.This is a possible explanation for why the CD did not identify PIP2 as a cause of changes of other agents in [17].Table 3 contains parameter estimates, confidence intervals, and p-values for the CD, IV and Hybrid estimators fit to the data in Figure 4.As expected, the CD does not find a significant causal effect while IV does.The Hybrid estimator, which can leverage changes in mean or variance, produces estimates very similar to the IV.

MEK
We now consider estimating the total causal effect of MEK on RAF, using data from the observational condition and condition 4. Since the condition 4 reagent targets MEK, this serves as a good instrument.Using the Causal Dantzig, [17] found that MEK has a direct effect on RAF.
Figure 5 shows the scatter plot of MEK versus RAF.The mean of MEK in condition 4 is higher than in the observational condition.Further the variance of MEK has increased in condition 4 relative to the observational condition.Thus both CD and IV are likely suitable to estimating a causal effect in this situation.Table 4 shows IV, CD and Hybrid parameter estimates, confidence intervals, and p-values fit using the data in Figure 5.All three estimators identify a causal effect.

Multivariate Analysis
For the multivariate analysis, we fit models using the 5 conditions specified in Table 2.Each condition is treated as a data generating environment.Thus there are 5 environments.An instrument/environment variable E ∈ R 4 is created following the procedure outlined in Equation (3.2).We iteratively treat each of the 11 agents as a response and regress it on the other 10 agent abundance measurements.Due to the limited number of conditions relative to reagents (4 instruments and 10 exposures), the IV can not be directly used in this situation.However, we would like the estimator to be sensitive to shifts in mean induced by the interventions.Thus we fit a Hybrid estimator specified in Equation (3.5).
We compare the performance of the Hybrid estimator with the CD and GCD.Note that, if an agent is used as the target or response variable, then intervention on that agent is not allowed in either CD, GCD or Hybrid.This is equivalent to an instrument having a direct effect on the response and would violate Assumption 1b.So when one of Akt, PIP2, Mek or PKC is used as the response, the total number of environments is 4 (observational and 3 conditions), while for all the other cases, the number of environments is 5 (observational and 4 conditions).For constructing the weight matrices, we use two-step estimators proposed Theorem 4.3 with initial weight matrices as described in Section 6.2 for the GCD and Section 6.3 for the Hybrid.
When the number of environments is greater than two, the GCD and the CD are not equivalent.[24] proposed two methods for handling the case with greater than 2 environments.We use the hiddenICP version of the CD from the R package InvariantCausalPrediction. With K > 2 environments, hiddenICP iteratively fits the CD with one environment versus all the other environments (K fits).This produces K point estimates and K confidence intervals for each parameter.The parameter estimates across the K fits are averaged to create a single point estimate.The confidence intervals lower limit is the smallest of the lower limits of the individual CD fit confidence intervals.Likewise the confidence interval upper limit is the largest of the upper limits of the individual confidence intervals.Thus the intervals are conservative.Figure 6 compares 95% confidence intervals for the CD, GCD, and Hybrid estimators with a) Plcg and b) PIP3 as the response.CD confidence intervals are very wide relative to GCD and Hybrid.The GCD and the Hybrid estimators perform similarly.In our univariate analysis, we found that PIP2 had a total effect on both Plcg and PIP3.In the multivariate analysis here, GCD and Hybrid identify PIP2 as direct causes of changes in both Plcg and PIP3.The CD fails to identify this effect, due to suboptimal merging of environments and an excessively conservative strategy for constructing confidence intervals.
We define a causal effect as strong if the entire 95% confidence interval is outside the range of (−.2, .2). Figure 7 shows the strong causal relations found by the Hybrid estimator in a graph.Compared with the results from [17] (Figure 3), our hybrid estimator finds more strong causal relations: 24 strong relations for hybrid versus 13 for CD.Some relations are found by both methods, e.g. the causal effects and the reverse causal effects between Raf-Mek, P38-Jnk, and Erk-Akt.The Hybrid model identifies the entire Raf->Mek->Erk path which was seen as a major validation of the original Bayesian network applied to the system [26] while the original application of the CD failed to identify the Mek->Erk edge.Further, the CD did not discover any causal effects from PIP2 to Plcg and PIP2 to PIP3, which are illustrated to exist from Figure 4 and discovered by our hybrid estimator.In addition, the hybrid estimator also found similar causal relations with other methods (ICP from [20] and the Bayesian network from [26]).For example, ICP and the hybrid both discovered the causal relations from PIP2 to Plcg and PKA to Erk, while the Bayesian network and the hybrid both discovered the relations from Plcg to PKC and PIP2 to PKC.In this work, we proposed two new methods, the Generalized Causal Dantzig (GCD) and Hybrid, for estimating causal effects in the presence of hidden confounding and reverse causality.The GCD generalizes the Causal Dantzig estimator of [24] to problems with continuous environments.Further we developed theory (based on GMM estimation) for the GCD in the over-identified case.This was not present in the original CD paper.The Hybrid estimator enforces both GCD environment and instrumental variable (IV) moment constraints, further illustrating connections between the concepts of environment and instrument which have been previously noted [11,20].We demonstrated the utility of these estimators in simulations and an application to Flow Cytometry data.

Discussion
In this work, we did not discuss high-dimensional estimation.[24] proposed using the Dantzig selector [5] to regularize the Causal Dantzig in high dimensional problems.This required extensive theoretical development.Since both the GCD and Hybrid are are GMM estimators, existing methods for fitting high dimensional GMM estimators (penalty terms, tuning parameter selection algorithms, and theory) may be applicable for the GCD and Hybrid estimators [3].This represents a direction for future research.The performance of these high-dimensional estimators may be compared with existing high-dimensional instrumental variable estimators [2,9,14].

CD and Lewbel [13]
We demonstrate equivalence between Lewbel's GMM (Equation 12 of [13]) and the CD/GCD for the case of univariate X and some restrictions on Lewbel's model (Equations 1 and 2 of [13]).With univariate X, the GMM defining the GCD (Equation (3.1)) becomes Translate from Lewbel notation to our notation according to: Lewbel's X is a vector of observed exogenous regressors which we do not consider (thus 0).Assume Lewbel's γ 2 = 0 (triangular system) and µ = E[Z] = 0 (environment/instrument E is mean 0).Under these restrictions and notational changes, Q 1 , Q 2 and Q 3 of Lewbel Section 3.1 are 0 and Thus Equation (9.1) equals Equation 12of Lewbel.Hence the estimators are equivalent.

Proof of Theorem 3.1
Recall q = #E − 1.We show that for any k ∈ {1, . . ., p} This implies the result because Noting that Equation (3.2) implies (s e1 P (E e = s e1 ))/(s e0 P (E e = s e0 )) = −1, we have So sufficient to show The r.h.s.implies that b e = b f for all e, f and c = b e for all e.Since the rows of A sum to 1, this implies the l.h.s.set condition is satisfied.Thus the r.h.s. is a subset of the l.h.s.Now we show set equivalence by arguing that if β is not in the r.h.s., it is not in the l.h.s.Note that all elements of A are non-negative and the rows of A sum to 1.If β is not in the r.h.s.set then either: A me b e + A m,q+1 c.Thus the l.h.s.equality constraints are not satisfied.

Proof of Theorem
where C does not depend on β.This is a quadratic function in β.The matrix M is positive definite under the assumptions that W ≻ 0 and X T (E • X) has column rank p.Thus the unique minimizer is given by 2. Define P = X T (E • X) W 1/2 ∈ R p×pr .When q = 1, P is invertible, so we have

3.
Recall that E i = 0 which implies that n 1 s 11 + n 2 s 10 = 0. Then we have

Proof of Theorem 4.3
Sufficient to show that W GCD → P W .
Let γ i be the i th row of E • X.Then
The quantity B 1 → P W −1 by the LLN.Sufficient to show B 2 , B 3 → P 0. Note that Similarly ∞ → P 0.

Proof of Theorem 5.1
Note that for all x we have First note Next we have • For the Hybrid estimator

Sufficient to show E E vec(EX T )
δ Y = 0.Each element of this vector was shown to be mean 0 in the proofs for m IV and m GCD .

Fig 1 .
Fig 1. a) Causal DAG model with instrument/environment E b) Causal DAG for observational data.c) Causal DAG for interventional data (do(X = x 1 ))

Fig 3 .
Fig 3. a) Model 1: IV and Hybrid dominate the GCD when the mean shift is strong but noise shift is weak.b) Model 2: GCD and Hybrid dominate IV when the mean shift is weak but the noise shift is strong.

Fig 5 .
Fig 5. Scatter plot that shows the relation between Mek and Raf.

Fig 6 .
Fig 6. a) Hybrid, GCD and CD estimation results when Plcg is the response.b) Hybrid, GCD and CD estimation results when PIP3 is the response.

Fig 7 .
Fig 7. Causal network found by the Hybrid estimator.The red circles represent the agents on which interventions were applied.

1 .A
max b > c.Select m ∈ argmax e∈{1,...,q} b e .Then b m = q e=1 A me b m + A m,q+1 b m > q e=1 me b e + A m,q+1 c.Thus the l.h.s.equality constraints are not satisfied.2. min b < c.Select m ∈ argmin e∈{1,...,q} b e .Then b m = q e=1 A me b m + A m,q+1 b m < q e=1 ϵ Y )]where the last equality holds because there are no backdoor paths from X to (H, ϵ Y ) or from X to (X, ϵ Y ).Thus there is constantc such that c = E[δ 2 Y (x, ϵ Y )] for all x and −c = E[δ 1 Y (H, ϵ Y )].Next note that Z 1 ⊥ ⊥ ϵ Y |Z 2for all sets of variables Z 1 and Z 2 which do not contain Y because all paths from Z 1 to ϵ Y must contain collider Y which is not in Z 2 (by assumption) and has no descendants ([19] d-Separation Definition 1.2.3).Thus we haveE ⊥ ⊥ ϵ Y |X E, H ⊥ ⊥ ϵ Y |X, ϵ X ϵ X ⊥ ⊥ ϵ Y |X • For the IV estimator m IV (β) = E[E(Y − X T β)] = E[EX T ](β 0 − β) + E[Eδ Y (H, X, ϵ Y )]Sufficient to show that E[Eδ Y (H, X, ϵ Y )] = 0. We have

= 0 •
For the GCD estimatorm GCD (β) = E[vec(EX T )(Y − X T β)] = E[vec(EX T )X T ](β 0 − β) + E[vec(EX T )δ Y (H, X, ϵ X )]Sufficient to show that E[vec(EX T )δ Y (H, X, ϵ X )] = 0. Let i ∈ {1, . . ., q} and j ∈ {1, . . ., p}.Elements of the vector are of the form .6)Under conditions specified in Section 4, Equation (2.6) will have a unique solution which equals the causal estimand β 0 .The CD estimator is constructed by enforcing sample based versions of constraints in Equation (2.6).Since X e ∈ R p , when #E = 2 (two data collection environments), the invariances specified in Equation (2.6) produce p sample constraints on β ∈ R p .Suppose there are n e observations from environment e.Let X e ∈ R ne×p (Y e ∈ R ne ) represent the design matrix (response vector) for environment e.Then the sample version of constraints in Equation (2.6) (with f = 2 and g = 1) is

Table 1
Coverage and Median Width of Confidence Intervals for Different Estimators

Table 2
Description of the conditions we used in our data application

Table 3
Estimation results using PIP2 as the input variable.The left shows the results when Plcg is the response, while the right being PIP3 as response.

Table 4
Estimation result for Mek→ Raf.