Estimating Treatment Effects with Observed Confounders and Mediators
Shantanu Gupta Zachary C. Lipton David Childers
stat.ME cs.LG econ.EM stat.ML

Collared and non-collared manifold boundaries
    Mathieu Baillif[
Haute école du paysage, d'ingénierie et d'architecture (HEPIA),
département d'Informatique et systèmes de communications,
Genève – Suisse. mathieu.baillif@hesge.ch
]
    
==========================================================================================================================================================================================




Given a causal graph, the do-calculus
can express treatment effects as functionals 
of the observational joint distribution 
that can be estimated empirically.
Sometimes the do-calculus identifies multiple valid formulae,
prompting us to compare the statistical properties 
of the corresponding estimators.
For example, the backdoor formula applies  
when all confounders are observed and
the frontdoor formula applies 
when an observed mediator transmits the causal effect.
In this paper, we investigate the over-identified scenario
where both confounders and mediators are observed,
rendering both estimators valid.
Addressing the linear Gaussian causal model,
we demonstrate that either estimator 
can dominate the other by an unbounded constant factor.
Next, we derive an optimal estimator,
which leverages all observed variables,
and bound its finite-sample variance.
We
show that it
strictly outperforms the backdoor and frontdoor estimators
and that this improvement can be unbounded.
We also present a procedure for combining two datasets,
one with observed confounders and another with observed mediators.
Finally, we evaluate our methods on both 
simulated data and the IHDP and JTPA datasets.




§ INTRODUCTION


Causal effects are not, in general, identifiable from observational data alone.
The fundamental insight of causal inference is that
given structural assumptions on the data generating process,
causal effects may become expressible as functionals
of the joint distribution over observed variables. 
The do-calculus, introduced by <cit.>,
provides a set of three rules that can be used
to convert causal quantities into such functionals.
We are motivated by the observation that,
for some causal graphs, treatment effects may be over-identified.
Here, applications of the do-calculus produce distinct functionals, 
all of which, subject to positivity conditions,
yield consistent estimators of the same causal effect. 
Consider a causal graph (see Figure <ref>)
for which the treatment X, mediator M, confounder W, 
and outcome Y are all observable. 
Using the backdoor adjustment, we can express the 
average treatment effect of X on Y as a function of P(X,W,Y),
while the frontdoor adjustment expresses that 
same causal quantity via P(X,M,Y) <cit.>.
In our experiments, we work with a real-world dataset that contains both confounders and mediators.
Faced with the (fortunate) condition of overidentification,
our focus shifts from identification: is our effect estimable?, to optimality: which among multiple valid estimators dominates from a standpoint of statistical efficiency?

In this paper, we address this very graph,
focusing our analysis on the linear causal model <cit.>,
a central object of study in causal inference and econometrics
and also explore the semiparametric setting.
Over-identification can arise in many other causal graphs (e.g. multiple backdoor adjustment sets, multiple instrumental variables, etc.).
However, we focus on this graph because 
the frontdoor estimator is a canonical example 
of a novel identification result derived using graphical models. 
It is central in the causality literature <cit.> 
and is a natural first step in the study of over-identified causal models.
Deriving the finite sample variance of the backdoor and frontdoor estimators,
and precisely characterizing conditions under which each dominates,
we find that either may outperform the other to an arbitrary degree
depending on the underlying model parameters.
These expressions can provide guidance to practitioners
for assessing the suitability of each estimator.
For example, one byproduct of our analysis is to characterize
what properties make for the “ideal mediator”.
Moreover, in the data collection phase, if one has a choice
between collecting data on the mediator or the confounder,
these expressions, together with the practitioner's beliefs
about likely ranges for model parameters,
can be used to decide what data to collect.



Next, we propose techniques that leverage 
both observed confounders and mediators.
For the setting where we simultaneously observe
both the confounder and the mediator, 
we introduce an estimator 
that optimally combines all information.
We prove theoretically that this method achieves lower mean squared error (MSE) 
than both the backdoor and frontdoor estimators,
for all settings of the underlying model parameters.
Moreover, the extent to which this estimator can dominate the better 
of the backdoor and frontdoor estimators is unbounded.
Subsequently, we consider the partially-observed setting
in which two datasets are available,
one with observed confounders (but not mediators) {(X, W, Y)}_i=1^n, 
and another with observed mediators (but not confounders) {(X, M, Y)}_i=1^m. 
Interestingly, the likelihood is convex
given simultaneous observations
but non-convex under partially-observed data. 
We introduce
an estimator 
that is guaranteed to achieve higher likelihood
than either the backdoor or frontdoor estimators.
Finally,
we evaluate our methods on  
synthetic, semi-synthetic, and real datasets.
Our proposed estimators that combine confounders and mediators always exhibit lower MSE than the backdoor and frontdoor estimators
when our model assumptions are satisfied.

Our principal contributions are the following:


  * Derivation of the parameter regimes 
where either of the frontdoor and backdoor estimators dominate
vis-a-vis sample efficiency.

  * Demonstration
of strict (and unbounded) improvements
of the optimal (combined) estimator  
over both the frontdoor and backdoor estimators.

  * Adaptation of a semi-parametric estimator to our graph,
showing the benefits of our approach in non-linear settings.

  * Analysis for the partially observed case,
where mediators and confounders are observed separately (but never simultaneously).




§ RELATED WORK


The backdoor adjustment formalizes the practice 
of controlling for known confounders
and is widely applied in statistics and econometrics <cit.>.
The frontdoor adjustment,
which leverages observed mediators to identify causal effects
even amid unobserved confounding, 
has seen increasing application in real-world datasets
<cit.>. 

In the most similar work to ours, 
<cit.> compare the frontdoor and backdoor adjustments,
computing bias (but not variance) formulas for each 
and performing sensitivity analysis.
Exploring a real-world job training dataset,
they demonstrate that the frontdoor estimator 
outperforms its backdoor counterpart (in terms of bias).
The finite sample variance of the frontdoor estimator for the linear Gaussian case was previously derived by <cit.>.
<cit.> compare the frontdoor and backdoor estimators based on their asymptotic variances and also show that the combined estimator's variance cannot be higher than the other two.
<cit.> derive the finite-sample variances of two possible adjustments in a three-variable binary causal graph
and show that the optimal estimator depends on the model parameters.
<cit.> introduce a graphical criterion for comparing the asymptotic variances of adjustment sets for the backdoor criterion in linear causal models.
<cit.> extend this work,
showing that the same graphical criterion is valid 
for non-parametric causal models.
They also present a semi-parametric efficient estimator 
that exploits the conditional independencies in a causal graph.

Researchers have also worked 
to generalize the frontdoor criterion.
<cit.> introduce 
the conditional frontdoor criterion,
allowing for both treatment-mediator confounders
and mediator-outcome confounders. 
<cit.> propose a method 
for including observed confounders along with a mediator 
with discrete treatments.

The study of overidentified models 
dates at least back to  <cit.>. 
<cit.> formalized the result 
that in the presence of overidentification, 
multiple estimators can be combined to improve efficiency. 
This was extended to the non-parametric setting by <cit.>.
A related line of work considers methods 
for combining multiple datasets for causal inference. 
<cit.> study the problem 
of handling biases while combining heterogeneous datasets, 
while <cit.> present Bayesian methods 
for combining datasets with different covariates and some common covariates.



§ PRELIMINARIES


In this work, we work within the structural causal model (SCM) framework due to  <cit.>,
formalizing causal relationships via directed acyclic graphs (DAGs).
Each X → Y edge in this DAG indicates that the variable X 
is (potentially) a direct cause of variable Y.
All measured variables are deterministic functions of their parents and a set of jointly independent per-variable noise terms. 

Linear Gaussian SCM    In linear Gaussian SCMs,
each variable is assumed to be a linear function of its parents. 
The noise terms are assumed to be additive and Gaussian. 
In this paper,
the finite sample results are derived for
the linear Gaussian SCM 
for the overidentified confounder-mediator graph 
(Figure <ref>),
where the structural equations can be written as

    w_i    = u^w_i, 
    
        x_i    = d w_i + u^x_i, 
    
        m_i    = c x_i + u^m_i, 
    
        y_i    = a m_i + b w_i + u^y_i,
    
        u^w_i    ∼𝒩(0, σ^2_u_w) 
    
        u^x_i    ∼𝒩(0, σ^2_u_x) 
    
        u^m_i    ∼𝒩(0, σ^2_u_m) 
    
        u^y_i    ∼𝒩(0, σ^2_u_y).

Here, w_i, x_i, m_i, and y_i are realized values 
of the random variables W, X, M, Y, respectively,
and u^w_i, u^x_i, u^m_i, u^y_i are realized values 
of the corresponding noise terms.
The zero mean assumption in Eq. <ref> simplifies analysis, but is not necessary for the results presented in this paper.



 §.§ The Backdoor and Frontdoor Adjustments


The effect of a treatment X is expressible 
in terms of the post-intervention distributions of 
the outcome Y for different values of the treatment X=x.
An intervention do(X=x) in a causal graph
can be expressed via the mutilated graph that results 
from deleting all incoming arrows to X, 
setting X's value to X=x for all instances, 
while keeping the SCM otherwise identical. 
This distribution is denoted as P(Y|do(X=x)).

The backdoor and frontdoor adjustments <cit.> 
express treatment effects
as functionals of the observational distribution. 
Consider our running example of the causal model in Figure <ref>. 
We denote X as the treatment, Y as the outcome,
W as a confounder, and M as a mediator. 
Our goal is to estimate the causal quantity P(Y|do(X=x)).

Backdoor Adjustment     
When all confounders of both X and Y are observed—in 
our example, W—then 
the causal effect of X on Y, i.e., P(Y|do(X=x)) can be written as

    P(Y|   do(X=x)) 
       = ∑_w P(Y|X=x,W=w)P(W=w).

Frontdoor Adjustment     
This technique applies even when the confounder W is unobserved.
Here we require access to a mediator M 
that (i) is observed;
(ii) transmits the entire causal effect from X to Y;
and (iii) is not influenced by the confounder W given X. 
The effect of X on Y is computed in two stages. 
We first find the effect of X on M, then the effect of M on Y as:

    P(M=m|do(X=x)) = P(M=m|X=x) 
       
            P(Y|   do(M=m)) 
       = ∑_x P(Y|M=m,X=x)P(X=x).

We can then write the causal effect of X on Y as

    P(   Y|do(X=x)) 
       = ∑_m P(M=m|do(X=x)) P(Y|do(M=m)).




§ VARIANCE OF BACKDOOR & FRONTDOOR ESTIMATORS


In this section, we analyze 
the backdoor and frontdoor estimators and characterize the regimes where each dominates. 
We work with the linear SCM described in Eq. <ref>. 
Throughout, our goal is to estimate 
the causal effect of X on Y.
In terms of the underlying parameters of the linear SCM, 
the quantity that we wish to estimate is ac. 
Absent measurement error, both estimators are unbiased (see proof in Appendix <ref>)
and thus we focus our comparison on their respective variances.

Variance of the Backdoor Estimator   
The backdoor estimator requires only that we observe {X, Y, W}
(but not necessarily the mediator M).
Say we observe the samples {x_i, y_i, w_i}^n_i=1.
We can estimate the causal effect ac 
by taking the coefficient on X in an OLS regression of Y on {X, W}. 
This controls for the confounder W 
and corresponds naturally to the adjustment described 
in Eq. <ref>.

The finite sample and asymptotic variances of the backdoor estimator are (see proof in Appendix <ref>)

    (ac)_backdoor = a^2 σ^2_u_m + σ^2_u_y/(n-3)σ^2_u_x, 
       lim_n →∞(√(n)(ac - ac))_backdoor =
            a^2 σ^2_u_m + σ^2_u_y/σ^2_u_x.


Variance of the Frontdoor Estimator   
The frontdoor estimator is used when {X, Y, M} samples
are observed.
Say we observe the samples {x_i, y_i, m_i}^n_i=1.
First, we estimate c by taking the coefficient on X in an OLS regression of M on X. Let the estimate be c. This corresponds to the adjustment in Eq. <ref>. Then, we estimate a by taking the coefficient on M in an OLS regression of Y on {M, X}. Let the estimate be a_f. This corresponds to the adjustment in Eq. <ref>.

The finite sample variances of c and a_f are (see proof in Appendix <ref>)

    (c) = σ^2_u_m/(n-2)(d^2 σ^2_u_w + σ^2_u_x), 
       (a_f) = b^2 σ^2_u_wσ^2_u_x + σ^2_u_y (d^2 σ^2_u_w + σ^2_u_x)/(n-3) (d^2 σ^2_u_w + σ^2_u_x) σ^2_u_m.

Using the facts that Cov(a_f, c) = 0 and (a_f^2, c^2) = (a_f)(c), the finite sample variance of the frontdoor estimator is (see proof in Appendix <ref>)

    (a_fc)    = c^2 (a_f) + a^2 (c) + 2 (a_f) (c).

And the asymptotic variance, which does not require Gaussianity, is (see proof in Appendix <ref>)

    lim_n →∞   (√(n)(a_fc - ac)) =  c^2 (b^2 σ^2_u_wσ^2_u_x + σ^2_u_y D)/D σ^2_u_m + a^2 σ^2_u_m/D, 
       where   D = d^2 σ^2_u_w + σ^2_u_x.


The Ideal Frontdoor Mediator   
A natural question then arises:
what properties of a mediator 
make the frontdoor estimator most precise?
We can see that (a_f c) is non-monotonic in the mediator noise σ_u_m.
Eq. <ref> provides us with guidance.
(a_f c) is a convex function of σ^2_u_m.
The ideal mediator
will have noise variance σ^2*_u_m which minimizes Eq. <ref>. That is,

    σ^2*_u_m   = _σ^2_u_m[ (a_fc) ] 
    σ^2*_u_m   = |c|√(b^2 σ^2_u_wσ^2_u_x + σ^2_u_y D )/|a|√(n-2/n-3),

where D = d^2 σ^2_u_w + σ^2_u_x.

Comparison of Backdoor and Frontdoor Estimators   
The relative performance of the backdoor and frontdoor estimators depend on the underlying SCM's parameters. 
Using Eqs. <ref> and <ref>, the ratio of the backdoor to frontdoor variance is

    R_Var = Var(ac)_backdoor/Var(a_fc)
       = (n-2)σ^2_u_m D^2 (a^2 σ^2_u_m + σ^2_u_y)/σ^2_u_x((n-3)a^2σ^4_u_mD + (2 σ^2_u_m + c^2 (n-2) D ) E),

where D = (d^2 σ^2_u_w + σ^2_u_x) and E = (b^2 σ^2_u_wσ^2_u_x + σ^2_u_y D).
The backdoor estimator dominates 
when R_Var < 1
and vice versa when R_Var > 1.
Note that there exist parameters
that cause any value of R_Var > 0. 
In particular, as σ^2_u_x→ 0, R_Var→∞ and as σ^2_u_x→∞, R_Var→ 0,
regardless of the sample size n. 
Thus, either estimator can dominate the other by any arbitrary constant factor.



§ COMBINING MEDIATORS & CONFOUNDERS


Having characterized the performance of each estimator separately, we now consider optimal strategies for estimating treatment effects in the overidentified regime, where we observe both the confounder and the mediator simultaneously.
Say we observe n samples {x_i, y_i, w_i, m_i}_i=1^n. 
We show that the maximum likelihood estimator (MLE) 
is strictly better than the backdoor and frontdoor estimators. The MLE will be optimal since our model satisfies 
the necessary regularity conditions for MLE optimality (by virtue of being linear and Gaussian).
The combined estimator is unbiased (see Appendix <ref>) and thus we focus on the variance.

Let the vector 𝐬_𝐢 = [x_i, y_i, w_i, m_i] denote the i^th sample.
Since the data is multivariate Gaussian, 
the log-likelihood of the data is ℒℒ = -n/2[log(Σ) + (ΣΣ^-1)],
where Σ = Cov([X, Y, W, M])
and Σ = 1/n∑_i=1^n𝐬_𝐢𝐬_𝐢^⊤.
The MLE for a Gaussian graphical model is Σ^MLE = Σ <cit.>. Let the MLE estimates for parameters c and a be c and a_c, respectively. Then

    c   = Σ_1,4/Σ_1,1,    a_c = Σ_1,4Σ_3,3 - Σ_1,3Σ_3,4/Σ_3,3Σ_4,4 - Σ^2_3,4.

The MLE estimate for c in Eq. <ref> 
is the same as for the frontdoor—the
coefficient of X in an OLS regression of M on X. 
The MLE estimate for a in Eq. <ref> 
is the coefficient of M in an OLS regression of Y on {M, W}. 
The finite sample variance of a_c is (see proof in Appendix <ref>)

    (a_c)    = σ^2_u_y/(n-3)(c^2 σ^2_u_x + σ^2_u_m).

The variance of c is the same 
as the frontdoor case as in Eq. <ref>.
Let r_1 = √(n-3/n-5), r_2 = √(3(n-2)/n-4), and L = ( c^2 σ^2_u_y/c^2 σ^2_u_x + σ^2_u_m + a^2 σ^2_u_m/d^2 σ^2_u_w + σ^2_u_x). We can bound the finite sample variance of the combined estimator as

    L/n   ≤(a_c c) ≤( c^2 (a_c) + a^2 (c) +  r_1 . 
       . ( 2 |c| (a_c) √((c)) + r_2 (a_c) (c) ) ).

The lower bound is derived using the Cramer-Rao theorem (since the estimator is unbiased) and for the upper bound, we use the Cauchy-Schwarz inequality. The complete proof is in Appendix <ref>.
And the asymptotic variance, which does not require Gaussianity, is (see proof in Appendix <ref>)

    lim_n →∞(√(n)(a_cc - ac))    = L.


The Ideal Mediator   
Just as with the frontdoor estimator, 
we can ask what makes for an ideal mediator 
in this case. 
Eq. <ref> shows that lim_n →∞(√(n)a_cc)
is a convex function of σ^2_u_m.
The ideal mediator will have noise variance σ^2*_u_m which minimizes the variance in Eq. <ref>. 
We use the asymptotic variance here since we only have finite-sample bounds on the variance of the combined estimator.
This means that

    σ^2*_u_m = _σ^2_u_m[ lim_n →∞(√(n)a_cc) ] 
       σ^2*_u_m = max{0, |c| σ_u_y√(d^2 σ^2_u_w + σ^2_u_x)/|a| - c^2 σ^2_u_x}.




 §.§ Comparison with Backdoor and Frontdoor Estimators


We can compare Eqs. <ref> and <ref> too see that, asymptotically, the combined estimator has lower variance than the backdoor estimator for all values of model parameters. That is, as n →∞, (√(n)a_c c) ≤(√(n)ac)_backdoor.
Similarly, we can compare Eqs. <ref> and <ref> to see that, asymptotically, the combined estimator is always better than the frontdoor estimator for all values of model parameters. That is, as n →∞, (√(n)a_c c) ≤(√(n)a_f c).

In the finite sample case, using Eqs. <ref> and <ref>, we can see
that for all model parameters,
for a large enough n,
the combined estimator will dominate the backdoor.
That is, ∃ N, s.t., ∀ n > N, (a_c c) ≤(ac)_backdoor,
where the dependence of N on the model parameters is stated in Appendix <ref>. We can make a similar argument for the dominance of the combined estimator over the frontdoor estimator. Using Eqs <ref> and <ref>, it can be shown that ∃ N, s.t., ∀ n > N, (a_c c) ≤(a_f c),
where the dependence of N on the model parameters is stated in Appendix <ref>.

Next, we show that the combined estimator can dominate the better of the backdoor and frontdoor estimators by an arbitrary amount. That is, we show that the quantity R = min{(ac)_backdoor, (a_fc) }/(a_cc)
is unbounded. Consider the case when (ac)_backdoor = (a_fc).
This condition holds for certain settings of the model parameters
(see Appendix <ref> for an example). 
Here,

    R    = (ac)_backdoor/(a_cc)
       ≥(n-2) D E (a^2 σ^2_u_m + σ^2_u_y)/σ^2_u_x( F + σ^2_u_y( σ^2_u_m + √(3)σ^2_u_m H  ) ) ,

where D = d^2 σ^2_u_w + σ^2_u_x, E = c^2 σ^2_u_x + σ^2_u_m, r_1 = √(n-3/n-5), r_2 = √(n-2/n-4), F = (n-3) a^2 σ^2_u_m E, G = r_1 σ_u_m/√((n-2) D), H = r_1 r_2 + |c| (n-2) D ( |c| + G ) and, in Eq. <ref>, we used Eq. <ref>. We can see that as σ_u_x→ 0, R →∞ and thus R is unbounded. This shows that, even in finite samples, combining confounders and mediators can lead to an arbitrarily better estimator than the better of the backdoor and frontdoor estimators.



 §.§ Semi-Parametric Estimators
 

<cit.> derive the efficient influence function
and semi-parametric efficiency bound for a generalized model with discrete treatment and non-linear relationships between the variables. While they allow for confounding of the treatment-mediator link and the mediator-outcome link, the graph in Figure <ref> has additional restrictions.
As per <cit.>, this graph is locally overidentified.
This suggests that it is possible to improve the estimator by <cit.> (which we refer to as IF-Fulcher).
In our model, there are two additional conditional independences compared to the graph studied in <cit.>: Y  X | (M, W), and M  W | X.
We incorporate these conditional independences in IF-Fulcher
by using [Y|M,W,X]=[Y|M,W], and f(M|X,W)=f(M|X) to create an estimator we refer to as
IF-Restricted:

    Ψ = 1/n∑_i=1^n (Y_i-𝔼[Y|M_i,W_i])f(M|x^*)/f(M|X_i) + 
       1{X_i=x^*}/P(X_i=x^*|W_i)×{𝔼[Y|M_i,W_i]- . 
       . ∑_m𝔼[Y|m,W_i]f(m|X_i)} + ∑_m𝔼[Y|m,W_i]f(m|x^*),

where, if f, P, and  are consistent estimators, then Ψp→𝔼[Y|do(X=x^*)].
By double robustness of the given estimator, if f, P, and  are correctly specified, then IF-Restricted has identical asymptotic distribution as IF-Fulcher.
But using the additional restrictions improves estimation of nuisance functions.
Thus we expect the proposed semi-parametric estimator to perform better in finite samples.
<cit.>, in contemporaneous work, analyzed the same graph and showed that, in addition, the efficient influence function is also changed when imposing these conditional independences (see Example 10 in their paper) (we refer to the estimator for this influence function as IF-Rotnitzky).
For our experiments with binary treatments, we use linear regression for f, 𝔼 and logistic regression for P.
Another way to adapt IF-Fulcher is for the case when we do not observe the confounders (as in the frontdoor adjustment). In this case, we can set W_i = ∅ and apply Ψ. We call this special case IF-Frontdoor.



§ COMBINING REVEALED-CONFOUNDER AND REVEALED-MEDIATOR DATASETS


We now consider a situation in which
the practitioner has access to two datasets.
In the first one, the confounders are observed but the mediators are unobserved. 
In the second one, the mediators are observed but the confounders are unobserved. 
This situation might arise if data is collected by two groups,
the first selecting variables to measure to apply the backdoor adjustment
and the second selecting variables to apply the frontdoor adjustment.
Given the two datasets,
we wish to optimally 
leverage all available data to estimate the effect of X on Y.

A naive approach would be to apply the backdoor and frontdoor estimator
to the first and second dataset, respectively, 
and take a weighted average of the two estimates. 
However, in this case, the variance will be between that of 
the frontdoor and backdoor estimator.
We analyze the MLE,
showing
that this estimator 
has lower asymptotic variance 
than both the backdoor and frontdoor estimators.

Combined Log-Likelihood under Partial Observability   
Say we have P samples of {x_i, y_i, w_i}^P_i=1. 
Let each such sample be denoted by the vector 𝐩_i = [x_i, y_i, w_i].
Moreover, say we have Q samples of {x_i, y_i, m_i}^Q_i=1. 
Let each such sample be denoted using the vector 𝐪_j = [x_j, y_j, m_j]. 
Let the observed data be represented as D. That is, D = {𝐩_1, 𝐩_2, , 𝐩_P, 𝐪_1, 𝐪_2, , 𝐪_Q }.
Let N = P + Q and let k = P/N.
Since the data is multivariate Gaussian, the conditional log-likelihood given k can be written as

    ℒℒ(D|k)    = -N/2[k  ( logΣ_p + (Σ_p Σ^-1_p)) + 
        (1-k) ( logΣ_q + (Σ_q Σ^-1_q)) ],

where Σ_p = Cov([X, Y, W]), Σ_q = Cov([X, Y, M]), 
Σ_p = ∑_i=1^P 𝐩_i 𝐩_i^⊤/P and Σ_q = ∑_i=1^Q 𝐪_i 𝐪_i^⊤/Q.



Cramer-Rao Lower Bound   
To compute the variance of the estimate of
e = ac, 
we compute the Cramer-Rao variance lower bound. 
We first compute the Fisher information matrix 𝐈 
as 𝐈 = - [ ∇^2_θℒℒ], where θ represents the eight model parameters.
Let e be the MLE. Since regularity holds for our model (due to linearity and Gaussianity), the MLE is asymptotically normal. Using the Cramer-Rao theorem, for constant k, as N →∞, we have √(N)(e - e) d→𝒩(0, V_e), where V_e is a function of 𝐈^-1.
The closed form expression for V_e is given in Appendix <ref>.

For any fixed k ∈ (0, 1),
(V_e - (√(P)ac)_backdoor) < 0 and (V_e - (√(Q)a_f c)) < 0, where
 is asymptotic variance.
This shows that the combined estimator always has lower 
asymptotic variance than that of the backdoor and frontdoor estimators on the individual datasets. Moreover, we also find cases where the combined estimator outperforms both the backdoor and frontdoor estimators even when the total number of samples are the same.
That is, there exist model parameters such that (V_e - (√(N)ac)_backdoor) < 0 and (V_e - (√(N)a_f c)) < 0 for some k ∈ (0, 1).
This means that is these cases, it is better to collect a mix of confounders and mediators rather than only collecting mediators or confounders.
Despite having access to the same number of samples, a mix of confounders and mediators can lead to lower variance.
This happens when the variances of the backdoor and frontdoor estimators are close to each other.
In Figure <ref>, 
we present two examples of causal graphs 
where having a mix of confounders and mediators leads to the lowest asymptotic variance
(see Appendix <ref> 
for parameter values).

The Maximum Likelihood Estimator   
Computing an analytical solution for the model parameters 
that maximizes the log-likelihood turns out to be intractable. 
As a result, we update our estimated parameters 
to maximize the likelihood numerically.
The likelihood in Eq. <ref> is non-convex. 
So we intialize the parameters using the two datasets (see Appendix <ref> for details) and run the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm <cit.> to maximize the likelihood.
In our experiments, 
the non-convexity of the likelihood 
never proved a practical problem.
When we find the global minimum, 
this estimator is optimal and dominates 
both the backdoor and frontdoor estimators.



§ EXPERIMENTS


Synthetic Data    
To show that the empirical variance of the various estimators 
is close to the theoretical variance
(Table <ref>), we randomly initialize parameters and for each instance, we compute the Mean and Standard Deviation of Absolute Percentage Error
of theoretical variance as a predictor of empirical variance
(see Appendix <ref>). Next, we compare the estimators under different settings of the model parameters. Unless stated otherwise, the model parameter values we use for experiments are a = 10, b = 4, c = 5, d = 5, σ^2_u_w = 1, σ^2_u_x = 1, σ^2_u_m = 1, σ^2_u_y = 1.





The quantity of interest is the causal effect ac = 50.
For Figure <ref>, we set {σ^2_u_x = 0.05, σ^2_u_m = 0.05},
which makes the backdoor estimator better 
as predicted by Eq. <ref>. 
For Figure <ref>,
we set {σ^2_u_w = 2, σ^2_u_x = 0.01, σ^2_u_m = 0.1} which makes the frontdoor estimator better as predicted by Eq. <ref>.
The plots in Figure <ref> and <ref> corroborate these predictions at different sample sizes. Furthermore, the optimal combined estimator always outperforms 
both the backdoor and frontdoor estimators.

Next, we evaluate the procedure for combining datasets described in Section <ref>,
generating two datasets with equal numbers of samples. 
In the first, only {X, Y, W} are observed. 
In the second, only {X, Y, M} are observed. 
We set {σ^2_u_x = 0.05, σ^2_u_m = 0.05},
which makes the backdoor estimator better
(Figure <ref>),
and then set {σ^2_u_w = 2, σ^2_u_x = 0.01, σ^2_u_m = 0.1},
which makes the frondoor estimator better
(Figure <ref>).
The plots show that the combined estimator 
has lower MSE than either for various sample sizes
(Figure <ref>),
supporting our theoretical claims.

IHDP Dataset   
<cit.> constructed a dataset from 
the Infant Health and Development Program (IHDP). 
This semi-synthetic dataset has been used 
for benchmarking in causal inference <cit.>.
The dataset is
is based on a randomized experiment to measure
the effect of home visits from a specialist 
on future cognitive test scores of children.
The treatment is binary
and
the covariates contain both continuous and categorical variables
representing measurements on the child and their mother.
We use samples from the NPCI package <cit.>.
We converted the randomized data into an observational study
by removing a biased subset of the treated group.
This set contains 747 samples with 25 covariates.

We use the covariates and the treatment assignment 
from the real study. 
We use a procedure similar to <cit.> 
to simulate the mediator and the outcome. 
The mediator M takes the form M ∼𝒩(cX, σ^2_u_m), where X is the treatment. 
The response Y takes the form Y ∼𝒩(aM + w^⊤𝐛, 1) where w is the vector of standardized 
(zero mean and unit variance) covariates 
and values in the vector 𝐛 are randomly sampled
(0, 1, 2, 3, 4) with probabilities (0.5, 0.2, 0.15, 0.1, 0.05). The ground truth causal effect is c× a.

We evaluate our estimators and the four IF estimators: IF-Fulcher (IF-Fulc), IF-Restricted (IF-Restr), IF-Frontdoor (IF-FD), and IF-Rotnitzky (IF-Rotz) (Section <ref>).
We test the estimators on two settings of the model parameters (Table <ref>, the Complete dataset setting). 
The MSE values are computed across 1000 instantiations of the dataset created 
by simulating the mediators and outcomes. 
We first evaluate the estimators on the complete dataset of 747 samples. We see that for Setting 1 (S1): a=10, c=5, σ_u_m = 1, the backdoor estimator dominates the frontdoor estimator whereas
for Setting 2 (S2): a=10,c=1,σ_u_m = 2,
the frontdoor estimator is better.
In both cases, the combined estimator (Section <ref>) outperforms both estimators. Furthermore, we see that IF-Restricted outperforms IF-Frontdoor, showing the value of leveraging the covariates. Moreover, IF-Restricted also outperforms IF-Fulcher, suggesting
that incorporating the additional model restrictions improves performance. Next, we randomly split the data into two sets, 
one with the confounder observed 
and the other with the mediator observed,
finding that the estimator that combines the datasets (Section <ref>)
outperforms the frontdoor and backdoor estimator
(Table <ref>, the Partial dataset setting).
We compute the MSE over 1000 realizations of the dataset.





National JTPA Study   
The National Job Training Partnership Act (JTPA) Study evaluates the effect of a job training program on future earnings.
We use the dataset from <cit.>.
The binary treatment X represents if a participant signed up for the program. The outcome Y represents future earnings. The collected covariates (like race, study location, age) are the confounders W. 
The covariates contain both categorical and continuous variables.
There was non-compliance among the treated units. The binary mediator M represents compliance, that is, whether the participant make use of JTPA services after signing up.
The study contained a randomized component which allowed us to compute the ground truth treatment effect, which was 862.74. <cit.> showed that the backdoor estimator has high bias, suggesting that
there was unmeasured confounding, so we omit the backdoor estimator in our results.
They also justify the assumptions required for the frontdoor estimator and show that it works well for this study.

A comparison of the frontdoor estimator, the combined estimator (Section <ref>), IF-Restricted, IF-Frontdoor and IF-Fulcher (Section <ref>) is shown in Table <ref> (the “C” data setting). For IF-Restricted, we only use the f(M|X,W)=f(M|X) restriction and do not use the [Y|M,W,X]=[Y|M,W] restriction since it is not valid. We compute the variance and MSE using 1000 bootstrap iterations. The combined estimator has lower variance and MSE than the frontdoor estimator. 
IF-Restricted outperforms IF-Frontdoor, 
reinforcing the utility of combined estimators. 
Furthermore, IF-Restricted outperforms IF-Fulcher, showing that using model restrictions is valuable.
Next, we evaluate our procedure
for the partially-observed setting (Section <ref>). 
We compute variance and MSE across 1000 bootstrap iterations. At each iteration, we randomly split our dataset into two datasets of equal size, one with revealed confounders, one with revealed mediators.
The combined estimator does not outperform the frontdoor estimator (Table <ref>, the “P” data setting).
This is expected since the backdoor adjustment works poorly and 
the revealed-confounder data is unlikely to help.
Despite this, the combined estimator does not suffer too badly
and has low bias despite the required assumptions for one of the identification strategies not holding.



§ DISCUSSION


In this paper, we studied over-identified graphs
with confounders and mediators,
showing that the two identification strategies 
can lead to estimators with arbitrarily different variances. 
We show that having access to both confounders and mediators
(either simultaneously or in separate datasets)
can give (unbounded) performance gains.
We also show that our results qualitatively apply to general non-linear settings.



  
Future Work
We see several promising lines for future work, including
(i) extensions to more general graphs;
(ii) online data collection subject to some cost structure over the observations;
and (iii) leveraging overidentification to mitigate errors due to measurement and confounding.
Our experiments show the applicability of our methods in the frontdoor-backdoor graph,
with combined estimators yielding gains in both linear and non-linear settings. 
We expect these insights to extend to other over-identified settings (e.g. graphs with multiple instrumental variables, multiple confounders, etc.)
and we hope next to extend the results to more general over-identified causal graphs.
Additionally, we plan to analyze the online data collection setting.
Here, subject to budget constraints, a practitioner must choose 
which variables to observe at each time step.
This direction seems especially important in medical applications (where each test may be costly) 
and survey studies (with a cap on the number of questions).
Our current results suggest that the optimal strategy 
must depend on the model parameters. 
At each step, the revealed data will improve 
our estimates of the model parameters,
in turn impacting what we collect in the future.

One potential limitation of the method is that situations where there exist multiple valid identification formulas may be uncommon in practice, when finding a single source of identification can already be difficult. However, we believe that in reality, many identification approaches are often available,  
but members of the community find flaws in each of the proposed estimators. 
In these cases, with multiple imperfect estimators of the same causal effect,
we believe that overidentification might be leveraged 
to create robust combined estimators.







§ BRIEF REVIEW OF OLS REGRESSION


Since we use OLS regression for our results, 
we briefly review OLS estimators.
We consider the following setup:

    𝐲   = 𝐗β + 𝐞,

where 𝐲 and 𝐞 are n × 1 vectors, 
𝐗 is an n × d matrix of observations, and 
β is the d × 1 coefficient vector that we want to estimate. 
If 𝐞𝐗 
and 𝐞∼𝒩(0, σ^2_e 𝐈_𝐧), where 𝐈_𝐧 is the n × n identity matrix, then the OLS estimate of β is

    β   = (𝐗^⊤𝐗)^-1𝐗^⊤𝐲
       = β + (𝐗^⊤𝐗)^-1𝐗^⊤𝐞,

with 𝔼[β] = β and Var(β) = σ^2_e 𝔼[(𝐗^⊤𝐗)^-1].
If each row X_i of 𝐗 is sampled from X_i i.i.d.∼𝒩(0, Σ), then the distribution of (𝐗^⊤𝐗)^-1 is an Inverse-Wishart distribution. 
Then the variance of β is

    Var(β) = σ^2_e Σ^-1/n - d - 1.




§ COVARIANCE OF A_HAT AND C_HAT




 §.§ Frontdoor estimator
 

We prove that Cov(a_f, c) = 0 for the frontdoor estimator. The expressions for a_f and c are

    c   [t]
           = ∑ x_i m_i/∑ x^2_i
       = c + ∑ x_i u^m_i/∑ x^2_i
    a_f    [t]
           = ∑ x^2_i ∑ m_i y_i - ∑ x_i m_i ∑ x_i y_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2
       = a + ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2,

where e_i = -b/du^x_i + u^y_i. Using the fact the (u^x, x) is bivariate normally distributed, we get

    [e|x]    = b σ^2_u_x/d (d^2 σ^2_u_w + σ^2_u_x) x 
       = F x,

where F =  b σ^2_u_x/d (d^2 σ^2_u_w + σ^2_u_x). The covariance then is

    Cov(a_f, c)    = 𝔼[(a_f - a)(c - c)] 
       = 𝔼[𝔼[(a_f - a)(c - c)|x, m]] 
       =𝔼[ 𝔼[( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ m_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ( ∑ x_i u^m_i/∑ x^2_i) |x,m ] ] 
       = 𝔼[ (∑ x^2_i ∑ m_i 𝔼[e_i|x] - ∑ m_i x_i ∑ x_i 𝔼[e_i|x]/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) (∑ u^m_i x_i/∑ x^2_i) ] 
       = 𝔼[ F (∑ x^2_i ∑ m_i x_i - ∑ m_i x_i ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) (∑ u^m_i x_i/∑ x^2_i) ] 
       = 0,

where in Eq. <ref> we used the expression from Eq. <ref>. Also, in Eq. <ref>, we took u^m_i out of the conditional expectation because u^m_i is given x_i and m_i (because u^m_i = m_i - c x_i).



 §.§ Combined estimator
 

We prove that Cov(a_c, c) = 0 for the combined estimator from Section <ref>. The expressions for a_c and c are

    c   [t]
           = ∑ x_i m_i/∑ x^2_i
       = c + ∑ x_i u^m_i/∑ x^2_i
    a_c    [t]
           = ∑ w^2_i ∑ m_i y_i - ∑ w_i m_i ∑ w_i y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2
       = a + ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2.

The covariance is

    Cov(a_c, c)    = 𝔼[(a_c-a)(c-c)] 
       = 𝔼[𝔼[(a_c-a)(c-c)|x, m, w]] 
       = 𝔼[ (∑ w^2_i ∑ m_i 𝔼[u^y_i] - ∑ m_i w_i ∑ w_i 𝔼[u^y_i]/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2) ( ∑ u^m_i x_i/∑ x^2_i) ] 
       = 0 ,

where in <ref> we used the fact that [u^y_i] = 0.



§ UNBIASEDNESS OF THE ESTIMATORS





 §.§ Backdoor estimator



Recall that for the backdoor estimator, we take the coefficient of X in an OLS regression of Y on { X, W }. The outcome y_i can be written as

    y_i = ac x_i + b w_i + a u^m_i + u^y_i.

The error term a u^m_i + u^y_i is independent of (x_i, w_i). In this case, the OLS estimator is unbiased. Therefore, [ac_backdoor] = ac.



 §.§ Frontdoor estimator


For the frontdoor estimator, we first compute c by taking the coefficient of X in an OLS regression of M on X. The mediator m_i can be written as

    m_i = c x_i + u^m_i.

The error term u^m_i is independent of x_i. In this case, the OLS estimator is unbiased and hence, [c] = c.

We then compute a_f by taking the coefficient of M in an OLS regression of Y on { M, X }. The outcome y_i can be written as

    y_i    = a m_i + b/d x_i - b/d u^x_i + u^y.

In this case, the error term - b/d u^x_i + u^y is correlated with x_i. The expression for a is given in Eq. <ref>. The expectation [a_f] is

    [a_f]    = a + [ ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2] 
       = a + [ [ ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2| x,m ] ] 
       = a + [ ∑ x^2_i ∑ m_i [e_i|x] - ∑ x_i m_i ∑ x_i [e_i|x]/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2] 
       = a + [ ∑ x^2_i ∑ m_i (F x_i) - ∑ x_i m_i ∑ x_i (F x_i)/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2] 
       = a,

where, in Eq. <ref>, the expression for E[e_i|x] is taken from Eq. <ref>. Using the fact that (a_f, c) = 0 (see proof in Appendix <ref>), we can see that the frontdoor estimator is unbiased as

    [a_fc]    = [a_f][c] + (a_f, c) 
       = ac.




 §.§ Combined estimator



In the combined estimator, the expression for c is the same as for the frontdoor estimator. Therefore, as shown in Appendix <ref>, [c] = c. We compute a by taking the coefficient of M in an OLS regression of Y on { M, W }. The outcome y_i can be written as

    y_i = a m_i + b w_i + u^y_i.

The error term u^y_i is independent of (m_i, w_i). In this case, the OLS estimator is unbiased. Therefore, [a_c] = a. Using the fact that (a_c, c) = 0 (see proof Appendix <ref>), we can see that the combined estimator is unbiased as

    [a_cc]    = [a_c][c] + (a_c, c) 
       = ac.




§ VARIANCE RESULTS FOR THE FRONTDOOR, BACKDOOR, AND COMBINED ESTIMATORS




 §.§ Backdoor estimator


The outcome y_i can be written as 

    y_i = ac x_i + b w_i + a u^m_i + u^y_i.

We estimate the causal effect ac 
by taking the coefficient on X in an OLS regression of Y on {X, W}.
Let Σ = ([X, W]). Using Eq. <ref>, 
the finite sample variance of the backdoor estimator is

    (ac)_backdoor   = (a u^m + u^y) (Σ^-1)_1,1/n - 3
       = a^2 σ^2_u_m + σ^2_u_y/(n-3)σ^2_u_x.


OLS estimators are asymptotically normally for arbitrary error distributions (and hence, Gaussianity is not needed). Therefore, the asymptotic variance of the backdoor estimator is

    lim_n →∞(√(n)(ac - ac))_backdoor   = (a u^m + u^y) (Σ^-1)_1,1 = a^2 σ^2_u_m + σ^2_u_y/σ^2_u_x.




 §.§ Frontdoor estimator





  §.§.§ Variance of c_hat


The regression of M on X can be written as m_i = c x_i + u^m_i. Let Σ_c = (X). Using Eq. <ref>, (c) is

    (c)    = (u^m) (Σ^-1_c)/n-2
            = σ^2_u_m/(n-2)(d^2 σ^2_u_w + σ^2_u_x).




  §.§.§ Variance of a_hat frontdoor


The regression of Y on {M, X} can be written as y_i = a m_i + b/d x_i + e_i, where e_i = -b/d u^x_i + u^y_i. 
In this case, the error e_i is not independent of the regressor x_i. 
Using the fact that (u^x, x) has a bivariate normal distribution, Var(e|x) is

    Var(e|x)    = b^2 σ^2_u_wσ^2_u_x + σ^2_u_y (d^2 σ^2_u_w + σ^2_u_x)/(d^2 σ^2_u_w + σ^2_u_x)
       := V_e .

Note that V_e is a constant and does not depend on x.
From Eqs. <ref> and <ref>, we know that

    c   = c + ∑ x_i u^m_i/∑ x^2_i
    a_f    = a + ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2.

where e_i = -b/du^x_i + u^y_i.
Let

    A    = ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2
    
        C    = ∑ x_i u^m_i/∑ x^2_i.


First, we derive the expression for (a_f) as follows,

    (a_f)    = (a + A) 
       = (A) 
       = ([A|x, m]) + [(A|x,m)] 
       = ([∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2|x, m]) + [(A|x,m)] 
       = ([∑ x^2_i ∑ m_i [e_i|x] - ∑ x_i m_i ∑ x_i [e_i|x]/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2]) + [(A|x,m)] 
       = ([∑ x^2_i ∑ m_i (F x_i) - ∑ x_i m_i ∑ x_i (F x_i)/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2]) + [(A|x,m)] 
       = [(A|x,m)] 
       = [(∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2|x,m)] 
       = [ 1/( ∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2 )^2( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i |x,m)] 
       = [ 1/( ∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2 )^2(e_i|x_i) ∑ x^2_i ( ∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2 ) ] 
       = V_e [∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2] 
       = V_e [D],

where, in Eq. <ref>, we used the result from Eq. <ref>, and D = ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2. Using the fact that D has the distribution of a marginal from an inverse Wishart-distributed matrix, that is, if the matrix M ∼ℐ𝒲(([M, X])^-1, n), then D = M_1,1, in Eq. <ref>, we get

    (a_f)    = V_e [D] 
       = V_e ([M, X])^-1_1,1/n - 2 - 1
       =  b^2 σ^2_u_wσ^2_u_x + σ^2_u_y (d^2 σ^2_u_w + σ^2_u_x)/(n-3) (d^2 σ^2_u_w + σ^2_u_x) σ^2_u_m,

where the expression for V_e is taken from Eq. <ref>.



  §.§.§ Covariance of a_hat square frontdoor and c_hat square



We prove that (a_f^2, c^2) = (a_f)(c). This covariance can be written as

    (a_f^2, c^2)    = [(a_f^2 - [a_f^2])(c^2 - [c^2])] 
       = [(a_f^2 - (a_f) - ^2[a_f])(c^2 - (c) - ^2[c])] 
       = [a_f^2c^2] - (a_f)(c) - a^2 (c) - c^2 (a_f) - a^2 c^2 .

We can write [a_f^2c^2] as

    [a_f^2c^2]    = [(a + A)^2 (c + C)^2] 
       = [a^2 c^2 + c^2 A^2 + a^2 C^2 + A^2 C^2 + 2 aA C^2 + 2cC A^2] 
       = a^2 c^2 + c^2 (a) + a^2 (c) + [A^2 C^2] + [2 aA C^2] + [2cC A^2] .

Substituting the result from Eq. <ref> in Eq. <ref>, we get

    (a_f^2, c^2)    = [A^2 C^2] + [2 aA C^2] + [2cC A^2] .

Now we expand each term in Eq. <ref> separately. [2 aA C^2] is

    [2 aA C^2]    = 2 a [ (∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 2 a [ [ (∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) | x, m] ] 
       = 2 a [ (∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ x^2_i ∑ m_i [e_i|x] - ∑ x_i m_i ∑ x_i [e_i|x]/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 2 a [ (∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ x^2_i ∑ m_i (F x_i) - ∑ x_i m_i ∑ x_i (F x_i)/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 0,

where, in Eq. <ref>, the expression for [e|x] is taken from Eq. <ref>. 

Next, we simplify [2cC A^2] as

    [   2cC A^2] 
       = 2c [ (∑ x_i u^m_i/∑ x^2_i) ( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2)^2 ] 
       = 2c [ [ (∑ x_i u^m_i/∑ x^2_i) ( ∑ x^2_i ∑ m_i e_i - ∑ x_i m_i ∑ x_i e_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2)^2 | x,m ] ] 
       = 2c [ [ (∑ x_i u^m_i/∑ x^2_i) (∑ x^2_i)^2 (∑ m_i e_i)^2 + (∑ x_i m_i)^2 (∑ x_i e_i)^2 - 2 ∑ x^2_i ∑ m_i e_i ∑ x_i m_i ∑ x_i e_i/(∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2 )^2| x,m ] ] 
       = 2c [ (∑ x_i u^m_i/∑ x^2_i) ( (e|x) (∑ x^2_i)/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2 + F^2 (2 (∑ x^2_i)^2 (∑ x_i m_i)^2 - 2 (∑ x^2_i)^2 (∑ x_i m_i)^2)/(∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2)^2) ] 
       = 2c [ (∑ x_i u^m_i/∑ x^2_i) ( (e|x) (∑ x^2_i)/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 2c  V_e [ (∑ x_i u^m_i/∑ x^2_i) ( ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 2c  V_e [ (c - c ) ( ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] 
       = 2c  V_e ( [ c( ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] - c [ ( ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2) ] ) 
       = 2c  V_e ( [ c D ] - c [ D ] ),

where D = ∑ x^2_i/∑ x^2_i ∑ m^2_i - (∑ x_i m_i)^2.
Using the fact that c and D are independent of each other (see proof at the end of this section), we get

    [c D]    = [c][D] 
       = c [D].

Substituting the result from Eq. <ref> in Eq. <ref>, we get

    [2cC A^2]    = 2c V_e ( c [D ] - c [ D ] ) 
       = 0.

We proceed similarly to Eq. <ref> to write [A^2 C^2] as

    [A^2 C^2]    = V_e [C^2 D].

Then we further simplify [A^2 C^2] as

    [A^2 C^2]    = V_e [C^2 D] 
       = V_e [(c - c)^2 D] 
       = V_e ( [c^2] [D] - c^2 [D]  ) 
       = V_e ( (c) [D] + ^2[c] [D] - c^2 [D]  ) 
       = V_e (c) [D] 
       = V_e (c) ([M, X])^-1_1,1/n - 2- 1
       = V_e/(n-3)σ^2_u_m(c) 
       = (a_f) (c),

where, in Eq. <ref>, we used the fact that if the matrix M ∼ℐ𝒲(([M, X])^-1, n), then D = M_1,1 (that is, D has the distribution of a marginal from an inverse Wishart-distributed matrix), and in Eq. <ref>, the expression for V_e is taken from Eq. <ref>.

Substituting the results from Eqs. <ref>, <ref>, and <ref> in Eq. <ref>, we get

    (a_f^2, c^2)    = (a_f) (c).




  
Proof that c and D are independent. Let Σ be the following sample covariance matrix:

    Σ = 1/n[   ∑ m^2_i ∑ m_i x_i; ∑ m_i x_i   ∑ x^2_i ].

The distribution of Σ is a Wishart distribution. That is, Σ∼𝒲(([M, X]), n). Then (Σ_1,1 - Σ_1,2Σ_2,2^-1Σ_2,1) and (Σ_2,1, Σ_2,2) are independent <cit.>. We can see that

    Σ_1,1 - Σ_1,2Σ_2,2^-1Σ_2,1   = ∑ m^2_i ∑ x^2_i - (∑ x_i m_i)^2/∑ x^2_i
       = 1/D.

Therefore, we get

    1/D   (∑ x^2_i, ∑ x_i m_i) 
    ∴   D    (∑ x^2_i, ∑ x_i m_i) 
    ∴   D    ∑ x_i m_i/∑ x^2_i
    ∴   D    c.




  §.§.§ Finite Sample Variance of a_hat frontdoor * c_hat


The variance of the product of two random variables can be written as

    (a_fc)    = (a_f^2, c^2) + ((a_f) + ^2[a_f])((c) + ^2[c]) - ((a_f, c) + [a_f][c])^2 
       = (a_f^2, c^2) + ((a_f) + a^2)((c) + c^2) - ((a_f, c) + ac)^2,

where in Eq. <ref> we used the facts that [a_f] = a, and [c] = c (see Appendix <ref>). Using the facts that (a_f^2, c^2) = (a_f) (c) (from Eq. <ref>) and (a_f, c) = 0 (from Appendix <ref>), we get

    (a_fc)    = a^2 (c) + c^2 (a_f) + 2 (c) (a_f).




  §.§.§ Asymptotic Variance of a_hat frontdoor * c_hat


Using asymptotic normality of OLS estimators, which does not require Gaussianity, we have

    √(n)( [ a_f;   c ] -
        [ a; c ])    d→𝒩(0, 
        lim_n →∞[      _∞(a_f) (√(n)a_f, c); (√(n)a_f, c)        _∞(c) ]) 
    ∴  √(n)( [ a_f;   c ] -
        [ a; c ])    d→𝒩(0, 
        [ _∞(a_f)       0;       0   _∞(c) ]),

where _∞(a_f) and _∞(c) are the asymptotic variances of a_f and c, respectively. The expressions for asymptotic variances are

    _∞(a_f)    = V_e ([M, X])^-1_1,1 =  b^2 σ^2_u_wσ^2_u_x + σ^2_u_y (d^2 σ^2_u_w + σ^2_u_x)/(d^2 σ^2_u_w + σ^2_u_x) σ^2_u_m
    _∞(c)    = (u^m) (Σ^-1_c)
            = σ^2_u_m/d^2 σ^2_u_w + σ^2_u_x.

In order to compute the asymptotic variance of a_f c, we use the Delta method:

    √(n)(a_f c - ac)    d→𝒩(0,
        [ c a ][ _∞(a_f)       0;       0   _∞(c) ][ c; a ]) 
    ∴  √(n)(a_f c - ac)    d→𝒩(0,
         c^2 _∞(a_f) + a^2 _∞(c) ).




 §.§ Combined estimator





  §.§.§ Finite sample variance of a_hat combined
 

We can write the regression of Y on {M, W} 
as y_i = a m_i + b w_i + u^y_i. 
Let Σ_a_c = Cov([M, W]). 
Using Eq. <ref>, 
we get

    (a_c)    = (u^y_i)(Σ^-1_a_c)_1,1/n-3 = σ^2_u_y/(n-3)(c^2 σ^2_u_x + σ^2_u_m).




  §.§.§ Bounding the finite sample variance


We first compute the lower bound of the combined estimator. Since the estimator is unbiased (see Appendix <ref>), we can apply the Cramer-Rao theorem to lower bound the finite sample variance.

Let the vector 𝐬_𝐢 = [x_i, y_i, w_i, m_i] denote the i^th sample.
Since the data is multivariate Gaussian, 
the log-likelihood of the data is

    ℒℒ = -n/2[log(Σ) + (ΣΣ^-1)],

where Σ = Cov([X, Y, W, M])
and Σ = 1/n∑_i=1^n𝐬_𝐢𝐬_𝐢^⊤.
Let e = ac and e = a_c c. Since we want to lower bound the variance of e, we reparameterize the log-likelihood by replacing c with e/a to simplify calculations. Next, we compute the Fisher Information Matrix for the eight model parameters:

    𝐈   = -E [      ∂^2 ℒℒ/∂ e^2    ∂^2 ℒℒ/∂ e ∂ a    ∂^2 ℒℒ/∂ e ∂ b                   ∂^2 ℒℒ/∂ e ∂σ_u_y;    ∂^2 ℒℒ/∂ a ∂ e      ∂^2 ℒℒ/∂ a^2    ∂^2 ℒℒ/∂ a ∂ b                   ∂^2 ℒℒ/∂ a ∂σ_u_y;                 ⋮                 ⋮                 ⋮                 ⋱                 ⋮;  ∂^2 ℒℒ/∂σ_u_y∂ e  ∂^2 ℒℒ/∂σ_u_y∂ a  ∂^2 ℒℒ/∂σ_u_y∂ b                    ∂^2 ℒℒ/∂^2 σ_u_y ].

Therefore, using the Cramer-Rao theorem, we have

    (e)    = (a_c c) 
       ≥ (𝐈^-1)_1,1
       = 1/n( c^2 σ^2_u_y/c^2 σ^2_u_x + σ^2_u_m + a^2 σ^2_u_m/d^2 σ^2_u_w + σ^2_u_x).


Next, we compute a finite sample upper bound for (a_c^2, c^2). We derive this in a similar manner as the frontdoor estimator in Appendix <ref>.
From Eqs. <ref> and <ref>, we know that

    a_c    = a + ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2
    c   = c + ∑ x_i u^m_i/∑ x^2_i.

Let

    A    = ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2
    
        C    = ∑ x_i u^m_i/∑ x^2_i.

Then, similarly to Eq. <ref>, we get

    (a_c^2, c^2)    = [A^2 C^2] + [2 aA C^2] + [2cC A^2].

Now we simplify each term in Eq. <ref> separately. [2 aA C^2] can be simplified as

    [2 aA C^2]    = 2 a [ ( ∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2) ] 
       = [ [ ( ∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2) |x,m,w ] ] 
       = [ ( ∑ x_i u^m_i/∑ x^2_i)^2 ( ∑ w^2_i ∑ m_i [u^y_i] - ∑ w_i m_i ∑ w_i [u^y_i]/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2) ] 
       = 0,

where, in Eq. <ref>, we used the fact that [u^y] = 0.

Next, we simplify [2cC A^2] as

    [2cC A^2]    = 2c [ ( ∑ x_i u^m_i/∑ x^2_i) ( ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2)^2 ] 
       = 2c [ [ ( ∑ x_i u^m_i/∑ x^2_i) ( ∑ w^2_i ∑ m_i u^y_i - ∑ w_i m_i ∑ w_i u^y_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2)^2 | x,m,w ] ] 
       = 2c [ ( ∑ x_i u^m_i/∑ x^2_i) (u^y) ( ∑ w^2_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2) ] 
       = 2c σ^2_u_y[ C D ],

where D = ∑ w^2_i/∑ w^2_i ∑ m^2_i - (∑ w_i m_i)^2. We can upper bound the expression in Eq. <ref> as

    [2cC A^2]    = 2c σ^2_u_y[ C D ] 
       ≤ 2 |c| σ^2_u_y[ C D ] 
       ≤ 2 |c| σ^2_u_y√([C^2] [D^2])
       = 2 |c| σ^2_u_y√((c) ((D) + [D]^2))
       = 2 |c| σ^2_u_y√((c) (2 [(([M, W]))^-1_1,1]^2/(n-2-1)^2 (n-2-3) + ( (([M, W]))^-1_1,1/n-2-1)^2 ))
       = 2 |c| σ^2_u_y√((c))√(2/(n-3)^2 (n-5) (c^2 σ^2_u_x + σ^2_u_m)^2  + 1/(n-3)^2 (c^2 σ^2_u_x + σ^2_u_m)^2)
       = 2 |c| σ^2_u_y/(n-3) (c^2 σ^2_u_x + σ^2_u_m)√((c))√(n-3/n-5)
       = 2 |c| (a) √((c))√(n-3/n-5),

where, in Eq. <ref>, we used the Cauchy–Schwarz inequality, and in Eq. <ref>, we used the fact that if the matrix M ∼ℐ𝒲(([M, W])^-1, n), then D = M_1,1 (that is, D has the distribution of a marginal from an inverse Wishart-distributed matrix).

Similarly to Eq. <ref>, we simplify [A^2 C^2] as

    [A^2 C^2]    = σ^2_u_y[C^2 D].

The expression in Eq. <ref> can be upper bounded using the Cauchy-Schwarz inequality as

    [A^2 C^2]    = σ^2_u_y[C^2 D] 
       ≤σ^2_u_y√([C^4] [D^2])
       = σ^2_u_y√([C^4] ((D) + ^2[D]))
       = σ^2_u_y√([C^4] (2 [(([M, W]))^-1_1,1]^2/(n-2-1)^2 (n-2-3) + ( (([M, W]))^-1_1,1/n-2-1)^2 ))
       = σ^2_u_y/(n-3) (c^2 σ^2_u_x + σ^2_u_m)√(n-3/n-5)√([C^4])
       = (a_c) √(n-3/n-5)√([C^4]).

We can simplify E[C^4] as follows,

    E[C^4]    = [ ( ∑ x_i u^m_i/∑ x^2_i)^4 ] 
       = [ [ ( ∑ x_i u^m_i/∑ x^2_i)^4 | x ] ] 
       = [ 1/( ∑ x^2_i )^4[ ( ∑ x_i u^m_i )^4 | x ] ] 
       = [ 1/( ∑ x^2_i )^4{( ( ∑ x_i u^m_i )^2 | x ) + [ ( ∑ x_i u^m_i )^2 | x ]^2 }] 
       = [ 1/( ∑ x^2_i )^4{( ( ∑ x_i u^m_i )^2 | x ) + σ^4_u_m( ∑ x^2_i )^2 }] 
       = [ 1/( ∑ x^2_i )^4{( σ^2_u_m∑ x^2_i ( ∑ x_i u^m_i )^2/σ^2_u_m∑ x^2_i| x ) + σ^4_u_m( ∑ x^2_i )^2 }] 
       = [ 1/( ∑ x^2_i )^4{σ^4_u_m(∑ x^2_i )^2 ( ( ∑ x_i u^m_i )^2/σ^2_u_m∑ x^2_i| x ) + σ^4_u_m( ∑ x^2_i )^2 }] 
       = [ 1/( ∑ x^2_i )^4{σ^4_u_m(∑ x^2_i )^2 2 + σ^4_u_m( ∑ x^2_i )^2 }] 
       = [ 1/( ∑ x^2_i )^4{3 σ^4_u_m( ∑ x^2_i )^2 }] 
       = 3 σ^4_u_m[ 1/( ∑ x^2_i )^2] 
       = 3 σ^4_u_m[ ( 1/∑ x^2_i) + [ 1/∑ x^2_i]^2 ] 
       = 3 σ^4_u_m[ 2/(n-2)^2 (n-4) (d^2 σ^2_u_w + σ^2_u_x)^2 + 1/(n-2)^2 (d^2 σ^2_u_w + σ^2_u_x)^2] 
       = 3 σ^4_u_m/(n-2)^2 (d^2 σ^2_u_w + σ^2_u_x)^2[ n-2/n-4] 
       = 3 ((c^2))^2 [ n-2/n-4],

where, in Eq. <ref>, we used the fact that ( ∑ x_i u^m_i )^2/σ^2_u_m∑ x^2_i| x has a Chi-squared distribution, that is, ( ∑ x_i u^m_i )^2/σ^2_u_m∑ x^2_i| x ∼χ^2(1), and in Eq. <ref>, we used the fact that 1/∑ x^2_i has a scaled inverse Chi-squared distribution, that is, 1/∑ x^2_i∼Scale-inv-χ^2(n, (d^2 σ^2_u_w + σ^2_u_x)^2/n).

Substituting the result from Eq. <ref> in Eq. <ref>, we get

    [A^2 C^2]    ≤(a_c) (c) √(3 (n-3)(n-2)/(n-5)(n-4)).


Substituting the results from Eqs. <ref>, <ref>, and <ref> in Eq. <ref>, we get

    (a_c^2, c^2) ≤√(n-3/n-5)( 2 |c| (a_c) √((c)) + √(3)√(n-2/n-4)(a_c) (c) ).

The variance of the product of two random variables can be written as

    (a_cc)    = (a_c^2, c^2) + ((a_c) + ^2[a_c])((c) + ^2[c]) - ((a_c, c) + [a_c][c])^2 
       = (a_c^2, c^2) + ((a_c) + a^2)((c) + c^2) - ((a_c, c) + ac)^2,

where we used the facts that [a_c] = a, and [c] = c (see Appendix <ref>). Using the fact that (a_c, c) = 0 (see Appendix <ref>) and the upper bound for (a_c^2, c^2), we get

    (   a_c c) ≤
        c^2 (a_c) + a^2 (c) + √(n-3/n-5)( 2 |c| (a_c) √((c)) + √(3)√(n-2/n-4)(a_c) (c) ).




  §.§.§ Asymptotic variance


Using asymptotic normality of OLS estimators, which does not require Gaussianity, we have

    √(n)( [ a_c;   c ] -
        [ a; c ])    d→𝒩(0, 
        lim_n →∞[      _∞(a_c) (√(n)a_c, c); (√(n)a_c, c)        _∞(c) ]) 
    ∴  √(n)( [ a_c;   c ] -
        [ a; c ])    d→𝒩(0, 
        [ _∞(a_c)       0;       0   _∞(c) ]),

where _∞(a_c) and _∞(c) are the asymptotic variances of a_c and c, respectively. The expressions for the asymptotic variances are

    _∞(a_c)    = (u^y_i)(Σ^-1_a_c)_1,1 = σ^2_u_y/c^2 σ^2_u_x + σ^2_u_m
    _∞(c)    = (u^m) (Σ^-1_c)
            = σ^2_u_m/d^2 σ^2_u_w + σ^2_u_x.

In order to compute the asymptotic variance of a_c c, we use the Delta method:

    √(n)(a_c c - ac)    d→𝒩(0,
         [ c a ][ _∞(a_c)       0;       0   _∞(c) ][ c; a ]) 
    √(n)(a_c c - ac)    d→𝒩(0,
         c^2 _∞(a_c) + a^2 _∞(c) ).




§ COMPARISON OF COMBINED ESTIMATOR WITH BACKDOOR AND FRONTDOOR ESTIMATORS


In this section, we provide more details on the comparison of the combined estimator presented in <ref> to the backdoor and frontdoor estimators.



 §.§ Comparison with the backdoor estimator
 

In Section <ref>, we made the claim that

    ∃ N, s.t., ∀ n > N, (a_c c) ≤(ac)_backdoor.

In this case, by comparing Eqs. <ref> and <ref>, we have

    N = 2 ( σ^4_u_x F + d^2 σ^2_u_w (σ^2_u_mD + σ^2_u_xF) + c^2 σ^6_u_x√( c^2 σ^6_u_x D^2 (F + 2√(3)σ^2_u_m)))/σ^2_u_m D^2,

where D = d^2 σ^2_u_w + σ^2_u_x, E = c^2 σ^2_u_x + σ^2_u_m, and F = E + (1 + 2√(3))σ^2_u_m. Thus, for a large enough n, the combined estimator has lower variance than the backdoor estimator for all model parameter values.



 §.§ Comparison with the frontdoor estimator
 

In Section <ref>, we made the claim that

    ∃ N, s.t., ∀ n > N, (a_c c) ≤(a_f c).

In this case, by comparing Eqs. <ref> and <ref>, we have

    N = 2 ( σ^6_u_m + 2√(3)c^2 σ^4_u_mσ^2_u_x - c^4 σ^2_u_mσ^4_u_x + ( D + σ^4_u_m√(σ^4_u_m + 4√(3) c^2 σ^2_u_mσ^2_u_x - 2 c^4 σ^4_u_x)) )/D^2,

where D = c^6 σ^6_u_x. Thus, for a large enough n, the combined estimator has lower variance than the frontdoor estimator for all model parameter values.



 §.§ Combined estimator dominates the better of backdoor and frontdoor



In this section, we provide more details for the claim in Section <ref> that the combined estimator can dominate the better of the backdoor and frontdoor estimators by an arbitrary amount. We show that the quantity

    R = min{(ac)_backdoor, (a_fc) }/(a_cc)

is unbounded.

We do this by considering the case when (ac)_backdoor = (a_fc). Note that

    (ac)_backdoor = (a_fc) 
        b = √(- D (-a^2 σ^4_u_m ( (n-2)d^2 σ^2_u_w + σ^2_u_x ) + (-(n-2)d^2 σ^2_u_w (σ^2_u_m - c^2 σ^2_u_x) + σ^2_u_x E ) σ^2_u_y ) /σ^2_u_wσ^4_u_x (2 σ^2_u_m + (n-2)c^2 D )),

where D = d^2 σ^2_u_w + σ^2_u_x, and E = (n-2)c^2 σ^2_u_x - (n-4)σ^2_u_m. Hence, if the parameter b is set to the value given in Eq. <ref>, the backdoor and frontdoor estimators will have equal variance. We have to ensure that the value of b is real. b will be a real number if

    |c|    ≤σ_u_m/σ_u_x√(1 - 2 σ^2_u_x/ (n-2) D ),  and
    
        n    > 2.

For the value of b in Eq. <ref>, the quantity R becomes

    R    = (ac)_backdoor/(a_cc)
       ≥(n-2) D E (a^2 σ^2_u_m + σ^2_u_y)/σ^2_u_x( (n-3) a^2 σ^2_u_m E + σ^2_u_y( σ^2_u_m + √(3)σ^2_u_m( r_1 r_2 + |c| (n-2) D ( |c| + r_1 σ_u_m/√((n-2) D)) )  ) ) ,

where D = d^2 σ^2_u_w + σ^2_u_x, E = c^2 σ^2_u_x + σ^2_u_m, r_1 = √(n-3/n-5) and r_2 = √(n-2/n-4).

R does not depend on the parameter b. It is possible to set the other model parameters in a way that allows R to take any positive value. In particular, it can be seen that as σ_u_x→ 0, R →∞, which shows that R is unbounded.



§ COMBINING PARTIALLY OBSERVED DATASETS




 §.§ Cramer-Rao Lower Bound



We are interested in estimating the value of the product ac. 
Let e = ac. We reparameterize the likelihood in Eq. <ref> 
by replacing c with e / a. 
This simplifies the calculations and improves numerical stability. 
Now, we have the following eight unknown model parameters: {e,a,b,d,σ^2_u_w,σ^2_u_x,σ^2_u_m,σ^2_u_y}.

In order to compute the variance of the estimate of parameter e = ac, 
we compute the Cramer-Rao variance lower bound. 
We first compute the Fisher information matrix (FIM) 𝐈 
for the eight model parameters:

    𝐈   = -E [      ∂^2 ℒℒ/∂ e^2    ∂^2 ℒℒ/∂ e ∂ a    ∂^2 ℒℒ/∂ e ∂ b                   ∂^2 ℒℒ/∂ e ∂σ_u_y;    ∂^2 ℒℒ/∂ a ∂ e      ∂^2 ℒℒ/∂ a^2    ∂^2 ℒℒ/∂ a ∂ b                   ∂^2 ℒℒ/∂ a ∂σ_u_y;                 ⋮                 ⋮                 ⋮                 ⋱                 ⋮;  ∂^2 ℒℒ/∂σ_u_y∂ e  ∂^2 ℒℒ/∂σ_u_y∂ a  ∂^2 ℒℒ/∂σ_u_y∂ b                    ∂^2 ℒℒ/∂^2 σ_u_y ]

Let e be the MLE. Since standard regularity conditions hold for our model (due to linearity and Gaussianity), the MLE is asymptotically normal. We can use the Cramer-Rao theorem to get the asymptotic variance of e. That is, for constant k, as N →∞, we have

    √(N)(e - e)    d→𝒩(0, V_e),    and
    
        V_e    = (𝐈^-1)_1,1.


Below, we present the closed form expression for V_e.
Let V_e = X/Y. Then

    X =    (a^2 σ^2_u_m+σ^2_u_y) (-a^8 d^2 (k-1) σ^2_u_w (σ^2_u_m)^5 (d^2 σ^2_u_w+σ^2_u_x)^2 +a^6 (σ^2_u_m)^3 (d^2 σ^2_u_w+σ^2_u_x) 
        (b^2 σ^2_u_wσ^2_u_x (c^2 d^4 (σ^2_u_w)^2+d^2 σ^2_u_w (c^2 (k+1) σ^2_u_x+(-2 k^2+2 k+1) σ^2_u_m)+σ^2_u_x (c^2 k σ^2_u_x+ 
       σ^2_u_m)) + σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (c^2 d^4 (σ^2_u_w)^2+d^2 σ^2_u_w (c^2 (k+1) σ^2_u_x+(3-2 k) σ^2_u_m)+ 
        k σ^2_u_x (c^2 σ^2_u_x+σ^2_u_m)))-4 a^5 b c d (k-1) k σ^2_u_wσ^2_u_x (σ^2_u_m)^3 (d^2 σ^2_u_w+σ^2_u_x) 
        (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x))+a^4 (σ^2_u_m)^2 (b^4 (σ^2_u_w)^2 (σ^2_u_x)^2 (c^2 d^4 (σ^2_u_w)^2+ 
        d^2 σ^2_u_w (2 c^2 (-k^2+k+1) σ^2_u_x+(k+1) σ^2_u_m)+σ^2_u_x (c^2 (-2 k^2+2 k+1) σ^2_u_x+2 σ^2_u_m))- 
        b^2 σ^2_u_wσ^2_u_xσ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (2 c^2 d^4 (k-2) (σ^2_u_w)^2+d^2 σ^2_u_w (c^2 (4 k^2-3 k-5) σ^2_u_x+2 
        (k^2-2 k-1) σ^2_u_m)+σ^2_u_x (c^2 (4 k^2-5 k-1) σ^2_u_x-2 (k+1) σ^2_u_m))-(σ^2_u_y)^2 (d^2 σ^2_u_w+σ^2_u_x)^2 
        (c^2 d^4 (2 k-3) (σ^2_u_w)^2+d^2 σ^2_u_w (c^2 (2 k^2-k-3) σ^2_u_x+(k-3) σ^2_u_m)+k σ^2_u_x (c^2 (2 k-3) σ^2_u_x- 
        2 σ^2_u_m)))-4 a^3 b c d (k-1) k σ^2_u_wσ^2_u_x (σ^2_u_m)^2 σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x))+
        a^2 σ^2_u_m (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x)) (b^4 (σ^2_u_w)^2 (σ^2_u_x)^2 (c^2 (d^2 σ^2_u_w-(k-2) σ^2_u_x)+σ^2_u_m)+ 
        b^2 σ^2_u_wσ^2_u_xσ^2_u_y (2 c^2 d^4 (σ^2_u_w)^2+2 d^2 σ^2_u_w (c^2 (-k^2) σ^2_u_x+k (c^2 σ^2_u_x+σ^2_u_m)+2 c^2 σ^2_u_x)+σ^2_u_x (2 c^2 
        (-k^2+k+1) σ^2_u_x+(k+1) σ^2_u_m))+(σ^2_u_y)^2 (d^2 σ^2_u_w+σ^2_u_x) (d^2 σ^2_u_w+k σ^2_u_x) (σ^2_u_m-c^2 (2 k-3) 
        (d^2 σ^2_u_w+σ^2_u_x)))+c^2 (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x))^2 (b^4 (σ^2_u_w)^2 (σ^2_u_x)^2+b^2 σ^2_u_wσ^2_u_xσ^2_u_y
        (2 d^2 k σ^2_u_w+k σ^2_u_x+σ^2_u_x)+(σ^2_u_y)^2 (d^2 σ^2_u_w+σ^2_u_x) (d^2 σ^2_u_w+k σ^2_u_x))),

and

    Y =     (a^2 σ^2_u_m (d^2 σ^2_u_w+σ^2_u_x)+b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x)) (a^6 (-d^2) (k-1) 
       σ^2_u_w (σ^2_u_m)^4 (d^2 σ^2_u_w+σ^2_u_x)^2+a^4 (σ^2_u_m)^2 (d^2 σ^2_u_w+σ^2_u_x) (b^2 σ^2_u_wσ^2_u_x (d^2 k σ^2_u_w
        (c^2 σ^2_u_x-2 k σ^2_u_m+2 σ^2_u_m)+σ^2_u_x (c^2 k σ^2_u_x+σ^2_u_m))+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (d^2 σ^2_u_w
        (c^2 k σ^2_u_x-3 k σ^2_u_m+3 σ^2_u_m)+k σ^2_u_x (c^2 σ^2_u_x+σ^2_u_m)))-4 a^3 b c d (k-1) k σ^2_u_wσ^2_u_x (σ^2_u_m)^2 
        (d^2 σ^2_u_w+σ^2_u_x) (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x))+a^2 σ^2_u_m (b^4 (σ^2_u_w)^2 (σ^2_u_x)^2 (σ^2_u_x
        (σ^2_u_m-2 c^2 (k-1) k σ^2_u_x)-d^2 (k-1) σ^2_u_w (2 c^2 k σ^2_u_x+σ^2_u_m))+2 b^2 σ^2_u_wσ^2_u_xσ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) 
        (σ^2_u_x (σ^2_u_m-2 c^2 (k-1) k σ^2_u_x)-2 d^2 (k-1) k σ^2_u_w (c^2 σ^2_u_x+σ^2_u_m))+(σ^2_u_y)^2 (d^2 σ^2_u_w+σ^2_u_x)^2 
        (-d^2 (k-1) σ^2_u_w (2 c^2 k σ^2_u_x+3 σ^2_u_m)-k σ^2_u_x (2 c^2 (k-1) σ^2_u_x+(k-2) σ^2_u_m)))-4 a b c d (k-1) 
        k σ^2_u_wσ^2_u_xσ^2_u_mσ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (b^2 σ^2_u_wσ^2_u_x+σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x))+b^6 c^2 k (σ^2_u_w)^3 (σ^2_u_x)^4+b^4
        (σ^2_u_w)^2 (σ^2_u_x)^2 σ^2_u_y (d^2 σ^2_u_w+σ^2_u_x) (3 c^2 k σ^2_u_x-k σ^2_u_m+σ^2_u_m)+b^2 σ^2_u_wσ^2_u_x (σ^2_u_y)^2 (d^2 σ^2_u_w+σ^2_u_x) 
        (d^2 k σ^2_u_w  (3 c^2 σ^2_u_x-2 (k-1) σ^2_u_m)+σ^2_u_x (3 c^2 k σ^2_u_x-k^2 σ^2_u_m+σ^2_u_m))+(σ^2_u_y)^3 (d^2 σ^2_u_w+σ^2_u_x)^2 
        (d^2 σ^2_u_w (c^2 k σ^2_u_x-k σ^2_u_m+σ^2_u_m)+k σ^2_u_x (c^2 σ^2_u_x-k σ^2_u_m+σ^2_u_m))),

where k = P/N.



 §.§ Comparison with frontdoor and backdoor estimators





In this section, we show some examples of regimes where the combining partially observed datasets results in lower variance than applying either of the backdoor or frontdoor estimator even when the total number of samples are the same. In other words, there exist settings of model parameters such that, for some k ∈ (0, 1), we have

    V_e ≤(√(N)ac)_backdoor,  and  
        V_e ≤(√(N)a_f c).


Figure <ref> shows three examples where the optimal value of k is between 0 and 1. We plot the variance as predicted by the expression for V_e versus the value of k. The plots show that in some cases, it is better to collect a mix of confounders and mediators rather than only mediators or only confounders. The expression for V_e in the previous section allows us to verify that. This happens when the variance of the frontdoor and backdoor estimators do not differ by too much. 

In Figure <ref>, the model parameters are {a=10, b=3.7, c=5, d=5, σ^2_u_w=1, σ^2_u_x=1, σ^2_u_m=0.64, σ^2_u_y=1 }. In this case, the variance of the frontdoor estimator is lower than the backdoor estimator. Despite this, it is not optimal to only collect mediators. The optimal value of k is 0.303, that is, 30% of the collected samples should be confounders and the rest should be mediators to achieve lowest variance.

In Figure <ref>, the model parameters are {a=10, b=3.955, c=5, d=5, σ^2_u_w=1, σ^2_u_x=1, σ^2_u_m=0.64, σ^2_u_y=1 }. In this case, the variance of the frontdoor estimator is almost equal to that of the backdoor estimator. The optimal ratio k is 0.505, that is, we should collect the same of amount of confounders as mediators.

In Figure <ref>, the model parameters are {a=10, b=4.3, c=5, d=5, σ^2_u_w=1, σ^2_u_x=1, σ^2_u_m=0.64, σ^2_u_y=1 }. In this case, the variance of the frontdoor estimator is greater than the backdoor estimator. The optimal ratio k is 0.735, that is, we should collect the more confounders than mediators.



 §.§ Parameter initialization for finding the MLE



The likelihood in Eq. <ref> is non-convex. 
As a result, we cannot start with arbitrary initial values 
for model parameters because we might encounter a local minimum. 
To avoid this, we use the two datasets to initialize our parameter estimates. 
Each of the eight parameters can be identified 
using only data from one of the datasets.
For example, d can be initialized using 
the revealed-confounder dataset 
(via OLS regression of X on W). 
The parameter e is can be identified using either dataset, 
so we pick the value with lower bootstrapped variance.

After initializing the eight model parameters, 
we run the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm <cit.>
to find model parameters that minimize the negative log-likelihood.



§ MORE DETAILS ON EXPERIMENTS


Here we provide more details for how results in Table <ref> are generated.
We initialize the model parameters 
by sampling 200 times from the following distributions:

    a,b,c,d    ∼Unif[-10, 10]
    σ^2_u_w,σ^2_u_x,σ^2_u_m,σ^2_u_y   ∼Unif[0.01, 2].

For each initialization, we compute the Mean Absolute Percentage Error (MAPE) of the theoretical variance as a predictor of empirical variance:

    MAPE = | _theoretical - _empirical|/_empirical * 100%

We report the mean and standard deviation of the MAPE across 1000 realizations of datasets
sampled from Eq. <ref>.  We find that the theoretical variance is close 
to the empirical variance even for small sample sizes
(Table <ref>).


