11affiliationtext: Department of Computer Science
ETH Zürich
Switzerland 22affiliationtext: Max Planck Institute for Intelligent Systems,
Tübingen, Germany 33affiliationtext: Department of Engineering, University of
Cambridge, United Kingdom {kkladny, jvk, bs<EMAIL_ADDRESS>
# Causal Effect Estimation from Observational and Interventional Data
Through Matrix Weighted Linear Estimators

Klaus-Rudolf Kladny  Julius von Kügelgen  Bernhard Schölkopf  Michael
Muehlebach

###### Abstract

We study causal effect estimation from a mixture of observational and
interventional data in a confounded linear regression model with multivariate
treatments. We show that the statistical efficiency in terms of expected
squared error can be improved by combining estimators arising from both the
observational and interventional setting. To this end, we derive methods based
on matrix weighted linear estimators and prove that our methods are
asymptotically unbiased in the infinite sample limit. This is an important
improvement compared to the pooled estimator using the union of interventional
and observational data, for which the bias only vanishes if the ratio of
observational to interventional data tends to zero. Studies on synthetic data
confirm our theoretical findings. In settings where confounding is substantial
and the ratio of observational to interventional data is large, our estimators
outperform a Stein-type estimator and various other baselines.

##  1 Introduction

Estimating the causal effect of a treatment variable on an outcome of interest
is a fundamental scientific problem that is central to disciplines such as
econometrics, epidemiology, and social science (Angrist and Pischke, 2009;
Morgan and Winship, 2014; Imbens and Rubin, 2015; Hernán and Robins, 2020). A
fundamental obstacle to this task is the possibility of hidden confounding:
unobserved variables that influence both the treatment and the outcome may
introduce additional associations between them (Reichenbach, 1956). As a
result, estimators purely based on observational (passively collected) data
can be biased and typically do not recover the true causal effect.

This contrasts experimental studies such as randomized controlled trials
(RCTs; Neyman, 1923; Fisher, 1936), where the treatment assignment mechanism
is modified through an external intervention, thus breaking potential
influences of confounders on the treatment. For this reason, RCTs have become
the gold standard for causal effect estimation. However, obtaining such
interventional data is difficult in practice because the necessary experiments
are often infeasible, unethical, or very costly to perform.

(a) observational

(b) interventional

Figure 1: Overview. We estimate the causal effect $\bm{\alpha}$ of a
continuous multi-variate treatment $\mathbf{X}$ on a scalar outcome $Y$ in a
linear Gaussian structural equation model with unobserved confounders
$\mathbf{Z}$. Given a large sample from the observational setting in (a) and a
small sample from the interventional setting in (b), we establish an optimal,
sample-size dependent matrix weighting scheme for combining the biased, low-
variance observational estimator with the unbiased, high-variance
interventional estimator.

In contrast, observational data is usually cheap and abundant, motivating the
study of causal inference from observational data (Rubin, 1974; Pearl, 2009).
In fact, in certain situations causal effects can be identified and estimated
from purely observational data, even under hidden confounding, e.g., in the
presence of natural experiments (instrumental variables; Angrist et al., 1996)
or observed mediators (front-door adjustment; Pearl, 1995). However, this does
not apply to the general case in which a treatment $\mathbf{X}$ and an outcome
$Y$ are confounded by an unobserved variable $\mathbf{Z}$ as shown in Fig.
1(a).

In the present work, we study treatment effect estimation in this general
setting under the assumption that we have access to both observational and
interventional data. The latter can be viewed as sampled from the setting
shown in Fig. 1(b), where the arrow $\mathbf{Z}\to\mathbf{X}$ has been removed
as a result of the intervention on $\mathbf{X}$ (graph surgery; Spirtes et
al., 2000), and is thus unbiased for our task. Due to small sample size,
however, the estimator based only on interventional data may exhibit high
variance. Our main idea is therefore to use the (potentially large amounts of)
observational data for variance reduction—at the cost of introducing some
bias. This is achieved by forming a combined estimator, which is superior to
the purely interventional one in terms of mean squared error.

We make the key assumption that both the treatment $\mathbf{X}\to Y$ and
confounding $\mathbf{Z}\to\\{\mathbf{X},Y\\}$ effects are linear, but allow
for treatment $\mathbf{X}$ and unobserved confounder $\mathbf{Z}$ to be
continuous and multi-variate. We then consider a class of estimators of the
causal effect parameter vector that combine the unbiased, but high-variance
interventional estimator and the biased, but low-variance observational
estimator through weight matrices—akin to a multi-variate convex combination.
We study the statistical properties of these estimators, establish theoretical
optimality results, and investigate their empirical behavior through
simulations.

In summary, we highlight the following contributions:

  * •

We introduce a new framework of weighing linear estimators using matrices and
show that several existing approaches fall into this category (§ 4).

  * •

We prove that, unlike pooling observational and interventional data (Prop.
4.1), our matrix weighting approaches achieve vanishing mean squared error in
the interventional sample limit (Props. 4.3 and 4.4) if the ratio between
observational and interventional data is non-vanishing.

  * •

We discuss two practical approaches for variance reduction in estimating
optimal weight matrices (§ 4.4; Prop. 4.5), and demonstrate through
simulations that our estimators outperform baselines and existing methods in
situations where confounding is substantial (§ 5).

##  2 Related Work

Causal reasoning, i.e., inferring a causal query such as a causal effect, can
be split up into the tasks of (i) identification and (ii) estimation. Step (i)
operates at the population level and seeks to answer whether a causal question
can—at least in principle—be answered given infinite data. If the answer is
positive and a valid estimand is provided, step (ii) then aims to construct a
statistically efficient estimator.

A causal query is identified from a set of assumptions if it can be expressed
in terms of the available distributions (e.g., a mixture of different
observational and interventional distributions). To this end, Pearl’s do-
calculus (1995; 2009) provides an axiomatic set of rules for manipulating
causal expressions based on graphical criteria. The identification task has
been studied extensively (Tian and Pearl, 2002; Pearl and Bareinboim, 2014;
Bareinboim and Pearl, 2016) and has by now been solved for many settings of
interest: In these cases, the do-calculus—and its extensions (Correa and
Bareinboim, 2020)—are sound and complete in that they provide a valid estimand
if and only if one exists (Huang and Valtorta, 2006; Shpitser and Pearl, 2006;
Bareinboim and Pearl, 2012; Lee et al., 2020).

In our setting from Fig. 1, the causal effect $\bm{\alpha}$ is not
identifiable from observational data, but is trivially identified by
intervening on $\mathbf{X}$. Yet, this leaves open the question of how to
estimate $\bm{\alpha}$ from finite data in the best possible way. In contrast
to the plethora of works on identification, there is much less prior
literature about statistical efficiency of causal parameter estimation,
particularly for confounded settings.

Figure 2: Relation between shrinkage and estimator weighting in 2D. Classical
shrinkage methods like ridge regression (yellow) shrink the interventional
estimator $\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ toward the origin. Scalar
estimator weighting (red) instead shrinks toward the observational estimator
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$. Dashed circles show the covariances
of $\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ and
$\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$, here assumed isotropic.
$\bm{\Delta}$ (green) is the confounding-induced bias of
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$.

A common source of inspiration for both prior work and our approach is that of
shrinkage estimation. In light of the bias-variance decomposition of the mean
squared error (e.g., Hastie et al., 2009, p. 24), shrinkage can yield a
strictly better (“dominating”) estimator by reducing variance, at the cost of
introducing some bias. These ideas were first introduced in frequentist
statistics by Stein (1956); James and Stein (1961) who showed that the maximum
likelihood estimate of a multivariate mean is dominated by shrinking towards a
fixed point such as the origin. Similar ideas are also at the heart of
empirical Bayes analysis (Robbins, 1964; Efron and Morris, 1973; Efron, 2012).
For estimating a parameter vector $\bm{\alpha}$ in a linear model, as is the
focus of the present work, a classical shrinkage method is ridge regression
(Hoerl, 1970).

Instead of shrinking towards the origin, an intuitive idea for causal effect
estimation is to shrink towards the observational estimator. The hope is that
the latter constitutes a better attractor if the confounding bias is not too
large—despite a slight increase in variance compared to shrinking toward a
constant. We refer to this approach as scalar estimator weighting. Fig. 2
shows a visual comparison to classical shrinkage estimation. The most closely
related work on estimator weighting is that of Green and Strawderman (1991);
Green et al. (2005) and Rosenman et al. (2020). The former two consider
general biased and unbiased estimators. The latter propose weighting schemes
for estimating vectors of multiple binary treatment effects. These works are
strongly inspired by James-Stein shrinkage estimators and minimize a
generalized version of Stein’s unbiased risk estimate (Wasserman, 2006, p.
150). Rosenman et al. (2020) show optimality among scalar weights with respect
to minimizing the true risk as the dimensionality of the estimated treatment
effects goes to infinity. However, these theoretical results rely on knowledge
of the true covariance matrix of the interventional estimator (which is
typically unknown in practice), and the behavior of their estimators in the
infinite sample limit is not analyzed.

Other work that focuses on combining observational and interventional data to
estimate causal effects of binary treatments includes, e.g., Kallus et al.
(2018); Cheng and Cai (2021); Ilse et al. (2021); Rosenman et al. (2022); Hatt
et al. (2022), see Colnet et al. (2020) for a comprehensive survey.

Yang and Ding (2020) also study combining estimators of binary treatment
effects. However, in their framework an estimator with less bias in addition
to a second error-prone estimator is computed from a second observational
“validation set”, in which all confounders are measured. Our framework, in
contrast, does not require measurements of the confounders.

In the present work, we consider a general linear regression setting with
continuous (rather than binary) multi-variate treatments. To combine
observational and interventional data, we introduce a new class of matrix
(rather than scalar) weighted estimators, of which ridge regression and data
pooling are special cases. Instead of employing Stein’s unbiased risk
estimate, we develop and analyze estimates for the theoretically optimal
weight matrix, without making assumptions about the covariance structure of
estimators.

Most approaches to causal estimation, including the present work, assume that
the causal structure among variables is known and takes the general form of
the directed acyclic graph in Fig. 1. For prior work on leveraging
observational and interventional data for causal discovery, or structure
learning, see, e.g., Wang et al. (2017).

##  3 Setting & Preliminaries

#### Notation.

Upper case $Y$ denotes a scalar random variable, lower-case $y$ a scalar, bold
lower-case $\mathbf{x}$ a vector, and bold upper-case $\mathbf{X}$ either a
matrix or random vector. The spectral norm of a matrix $\mathbf{X}$ is denoted
by $\left\lVert\mathbf{X}\right\rVert_{2}$.

#### Causal Model.

To formalize our problem setting, we adopt the structural causal model
framework of Pearl (2009). Specifically, we assume that the causal
relationships between the $d$-dimensional confounder $\mathbf{Z}$, the
$p$-dimensional treatment $\mathbf{X}$, and the scalar outcome $Y$ are
captured by the following linear Gaussian structural equation model (SEM):

$\displaystyle\mathbf{Z}$
$\displaystyle\;\leftarrow\;\mathbf{N}_{\mathbf{Z}},$
$\displaystyle\mathbf{N}_{\mathbf{Z}}\sim\mathcal{N}(\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}},\bm{\Sigma}_{\mathbf{N}_{\mathbf{Z}}})$
(1) $\displaystyle\mathbf{X}$
$\displaystyle\;\leftarrow\;\mathbf{B}\mathbf{Z}+\mathbf{N}_{\mathbf{X}},$
$\displaystyle\mathbf{N}_{\mathbf{X}}\sim\mathcal{N}(\bm{\mu}_{\mathbf{N}_{\mathbf{X}}},\bm{\Sigma}_{\mathbf{N}_{\mathbf{X}}})$
(2) $\displaystyle Y$
$\displaystyle\;\leftarrow\;\mathbf{Z}^{\top}\bm{\gamma}+\mathbf{X}^{\top}\bm{\alpha}+N_{Y},$
$\displaystyle N_{Y}\sim\mathcal{N}(\mu_{N_{Y}},\sigma^{2}_{N_{Y}})$ (3)

with $\mathbf{B}\in\mathbb{R}^{p\times d}$, $\bm{\gamma}\in\mathbb{R}^{d}$,
$\bm{\alpha}\in\mathbb{R}^{p}$, and
$(\mathbf{N}_{\mathbf{Z}},\mathbf{N}_{\mathbf{X}},N_{Y})$ mutually independent
exogenous noise variables. The SEM in 1–(3) induces an observational
distribution over $(\mathbf{Z},\mathbf{X},Y)$ which is referred to as
$\mathbb{P}_{\text{obs}}$, see Fig. 1(a).

To model the interventional setting, we consider a soft intervention
(Eberhardt and Scheines, 2007), which randomizes the treatment $\mathbf{X}$ by
replacing the assignment in 2 with

$\mathbf{X}\;\leftarrow\;\widetilde{\mathbf{N}}_{\mathbf{X}},\qquad\qquad\widetilde{\mathbf{N}}_{\mathbf{X}}\sim\mathbb{P}_{\widetilde{\mathbf{N}}_{\mathbf{X}}},\qquad$
(4)

where $\widetilde{\mathbf{N}}_{\mathbf{X}}$ is mutually independent of
$\mathbf{N}_{\mathbf{Z}}$ and $N_{Y}$. We note that
$\widetilde{\mathbf{N}}_{\mathbf{X}}$ may be non-Gaussian. The modified
interventional SEM consisting of 1, 4 and 3 induces a different,
interventional distribution over $(\mathbf{Z},\mathbf{X},Y)$, which we refer
to as $\mathbb{P}_{\text{int}}$, see Fig. 1(b).

For ease of notation and for the remainder of this work, we assume without
loss of generality that all noise variables are zero-mean. Details on how to
extend our method to non zero-mean noise variables are provided in App. D.

#### Data.

We assume access to two separate datasets of observations of $(\mathbf{X},Y)$
of size $n$ and $m$, each sampled independently from the observational and
interventional distributions (i.i.d.), respectively:

$\displaystyle(\mathbf{x}_{i},y_{i})\;$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\;\mathbb{P}_{\text{obs}},\quad
i=1,...,n,$ $\displaystyle(\mathbf{x}_{i},y_{i})\;$
$\displaystyle\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\;\mathbb{P}_{\text{int}},\quad
i=n+1,...,n+m,$

where $\mathbb{P}_{\text{obs}}$ and $\mathbb{P}_{\text{int}}$ denote the
distributions of $(\mathbf{X},Y)$ in the observational and interventional
settings, respectively. We note that the confounder $\mathbf{Z}$ remains
unobserved. We concatenate the observational sample in a treatment matrix
$\mathbf{X}_{\textsc{o}}=(\mathbf{x}_{1},...,\mathbf{x}_{n})^{\top}\in\mathbb{R}^{n\times
p}$ and outcome vector
$\mathbf{y}_{\textsc{o}}=(y_{1},...,y_{n})^{\top}\in\mathbb{R}^{n}$, and
similarly with $\mathbf{X}_{\textsc{i}},\mathbf{y}_{\textsc{i}}$ for the
interventional sample. Finally, we denote the pooled data by
$\mathbf{X}_{\textsc{p}}=(\mathbf{X}_{\textsc{o}},\mathbf{X}_{\textsc{i}})\in\mathbb{R}^{(n+m)\times
p}$ and
$\mathbf{y}_{\textsc{p}}=(\mathbf{y}_{\textsc{o}},\mathbf{y}_{\textsc{i}})\in\mathbb{R}^{n+m}$.

#### Goal.

Our objective is to obtain an accurate estimate of the parameter vector
$\bm{\alpha}$, which characterizes the linear causal effect of $\mathbf{X}$ on
$Y$ in 3. Formally, it is given by

$\displaystyle\bm{\alpha}=\nabla_{\mathbf{x}}\mathbb{E}[Y|\mathrm{do}(\mathbf{X}\leftarrow\mathbf{x})],$

where the $\mathrm{do}(\cdot)$ operator denotes a manipulation of the
treatment assignment akin to 4, and the expectation is taken with respect to
the corresponding conditional distribution.

#### Confounding Issues.

In the general case with non-zero $\mathbf{B}$ and $\bm{\gamma}$, the
observational setting is confounded, meaning

$\mathbb{P}_{\text{obs}}(Y|\mathbf{X}=\mathbf{x})\neq\mathbb{P}(Y|\mathrm{do}(\mathbf{X}\leftarrow\mathbf{x}))=\mathbb{P}_{\text{int}}(Y|\mathbf{X}=\mathbf{x}),$

which complicates the use of observational data. Specifically, for our assumed
model 1–(3) the conditional expectation of $Y$ under $\mathbb{P}_{\text{obs}}$
is given by the following perturbed linear model (Ćevid et al., 2020):

$\mathbb{E}_{\text{obs}}[Y|\mathbf{X}=\mathbf{x}]=(\bm{\alpha}+\bm{\Delta})^{\top}\mathbf{x},$
(5)

where $\bm{\Delta}\in\mathbb{R}^{p}$ denotes the confounding bias, which is
given explicitly in terms of the model parameters as

$\bm{\Delta}=(\bm{\Sigma}_{\mathbf{N}_{\mathbf{X}}}+\mathbf{B}\bm{\Sigma}_{\mathbf{N}_{\mathbf{Z}}}\mathbf{B}^{\top})^{-1}\mathbf{B}\bm{\Sigma}_{\mathbf{N}_{\mathbf{Z}}}\bm{\gamma}.$
(6)

It can be seen from 6 that the confounding bias $\bm{\Delta}$ is zero if
$\mathbf{B}$ or $\bm{\gamma}$ are zero (i.e., $\mathbf{Z}$ only affects either
$\mathbf{X}$ or $Y$). Furthermore, we have that, in general,

$\text{Var}_{\text{obs}}(Y|\mathbf{X})=\sigma^{2}_{Y|\mathbf{X}}\neq\sigma^{2}_{Y|\text{do}(\mathbf{X})}=\text{Var}_{\text{int}}(Y|\mathbf{X})\,.$
(7)

#### Assessing Estimator Quality.

We rely on mean squared error with respect to the true parameter $\bm{\alpha}$
as a measure for comparing different estimators.

######  Definition 3.1 (MSE).

Let $\widehat{\bm{\alpha}}$ be any function of the pooled data
$(\mathbf{X}_{\textsc{p}},\mathbf{y}_{\textsc{p}})$ taking values in
$\mathbb{R}^{p}$. Then

$\mathrm{MSE}(\widehat{\bm{\alpha}})\coloneqq\mathbb{E}\left[\left\lVert\widehat{\bm{\alpha}}-\bm{\alpha}\right\rVert_{2}^{2}\right],$

where the expectation is taken over
$\mathbf{X}_{\textsc{p}},\mathbf{y}_{\textsc{p}}$.

We note that the mean squared error can also be written as follows:

$\mathrm{MSE}(\widehat{\bm{\alpha}})=\left\lVert\mathbf{Bias}(\widehat{\bm{\alpha}})\right\rVert_{2}^{2}+\mathrm{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}))\,,$
(8)

where

$\displaystyle\mathbf{Bias}(\widehat{\bm{\alpha}})$
$\displaystyle=\mathbb{E}[\widehat{\bm{\alpha}}]-\bm{\alpha}\,,$
$\displaystyle\mathbf{Cov}(\widehat{\bm{\alpha}})$
$\displaystyle=\mathbb{E}[(\widehat{\bm{\alpha}}-\mathbb{E}[\widehat{\bm{\alpha}}])(\widehat{\bm{\alpha}}-\mathbb{E}[\widehat{\bm{\alpha}}])^{\top}]\,.$

This decomposition highlights that biased estimators can dominate unbiased
ones through variance reduction.

#### Pure Estimators.

We study estimators for $\bm{\alpha}$ that are linear combinations of the
following ordinary least squares estimators obtained on the two data sets
individually.

######  Definition 3.2 (Pure Estimators).

For non-singular moment matrices
$\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}$ and
$\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}$, the _pure estimators_
based only on the observational/interventional sample are given by:

$\displaystyle\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$
$\displaystyle\coloneqq(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}})^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{y}_{\textsc{o}},$
$\displaystyle\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$
$\displaystyle\coloneqq(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{y}_{\textsc{i}}.$

Recall that $\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ is unbiased while
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ has bias $\bm{\Delta}$. Their
covariances conditionally on $\mathbf{X}_{\textsc{o}}$ and
$\mathbf{X}_{\textsc{i}}$ are given by

$\displaystyle\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$
$\displaystyle=(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}})^{-1}\sigma^{2}_{Y|\mathbf{X}},$
(9) $\displaystyle\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$
$\displaystyle=(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\sigma^{2}_{Y|\text{do}(\mathbf{X})}.$

Unlike previous work (see § 2), we do not make assumptions about the
covariance structure of either estimator.

#### Almost sure convergence.

To analyze the behavior of estimators in the infinite sample limit, we will
employ the following characterization known as almost sure convergence.

######  Definition 3.3 (Almost Sure Convergence).

Let $\mathbf{M}$ be a random matrix with realizations in $\mathbb{R}^{p\times
p}$. We say a sequence of random matrices $\widehat{\mathbf{M}}_{m}$ indexed
by $m\in\mathbb{N}$ _converges almost surely_ to $\mathbf{M}$, denoted
$\widehat{\mathbf{M}}_{m}\xrightarrow{\text{a.s.}}\mathbf{M}$, if and only if

$\lim_{m\to\infty}\text{P}\left(\widehat{\mathbf{M}}_{m}=\mathbf{M}\right)=1,$

where P denotes probability.

##  4 Matrix Weighted Linear Estimators

We now introduce our class of matrix weighted linear estimators, which combine
the two pure estimators from Def. 3.2 using a weight matrix $\mathbf{W}$ to
obtain a new (better) estimator.

######  Definition 4.1 ($\mathbf{W}$-weighted Linear Estimator).

Let $\mathbf{W}\in\mathbb{R}^{p\times p}$ (possibly random). The _$\mathbf{W}$
-weighted linear estimator_ for $\bm{\alpha}$ is given by

$\widehat{\bm{\alpha}}^{m}_{\mathbf{W}}\coloneqq\mathbf{W}\widehat{\bm{\alpha}}^{m}_{\textsc{i}}+(\mathbf{I}_{p}-\mathbf{W})\widehat{\bm{\alpha}}^{n}_{\textsc{o}}.$

We furthermore refer to $\mathbf{W}$ as a weight matrix.

We will generally think of $n$ as a function of $m$, where we sometimes even
explicitly write $n(m)$. However, to simplify notation we index estimators by
$m$ only, omitting the dependence $n(m)$.

Note that the purely interventional estimator is a special case of a
$\mathbf{W}$-weighted estimator with $\mathbf{W}=\mathbf{I}_{p}$. However,
while unbiased, it may be subject to high variance if $m$ is very
small.111E.g., consider a one-dimensional setting with $x_{i}=1$ if $i$ is
even and $-1$ otherwise. Then, for odd $m$,
$\text{Var}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}|\mathbf{X}_{\textsc{i}})\propto(\sum_{i}x_{i}^{2})^{-1}=\frac{1}{m}$.
Hence, we generally prefer to employ the observational data as well and choose
$\mathbf{W}\neq\mathbf{I}_{p}$.

###  4.1 Existing Methods as Special Cases

First, we show that several standard approaches can be viewed as special cases
of matrix-weighted estimators.

#### Data Pooling.

A straightforward approach for combining both data sets is to compute an
estimator on the pooled data. The resulting least-squares estimator
$\widehat{\bm{\alpha}}^{m}_{\textsc{p}}$ is:

$\displaystyle\widehat{\bm{\alpha}}^{m}_{\textsc{p}}:=$
$\displaystyle\,(\mathbf{X}_{\textsc{p}}^{\top}\mathbf{X}_{\textsc{p}})^{-1}\mathbf{X}_{\textsc{p}}^{\top}\mathbf{y}_{\textsc{p}}$
(10) $\displaystyle=$
$\displaystyle\,(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}+\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{y}_{\textsc{o}}+\mathbf{X}_{\textsc{i}}^{\top}\mathbf{y}_{\textsc{i}})$
$\displaystyle=$
$\displaystyle\,\mathbf{W}^{m}_{\textsc{p}}\widehat{\bm{\alpha}}^{m}_{\textsc{i}}+(\mathbf{I}-\mathbf{W}^{m}_{\textsc{p}})\widehat{\bm{\alpha}}^{n}_{\textsc{o}},$

where

$\mathbf{W}^{m}_{\textsc{p}}\coloneqq(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}+\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}.$
(11)

We see that $\widehat{\bm{\alpha}}^{m}_{\textsc{p}}$ indeed qualifies as a
valid matrix weighted estimator in the sense of Def. 4.1.

However, data pooling can lead to highly undesirable limiting behavior in
cases where the amount of observational data $n(m)$ does not vanish in the
limit of infinite interventional data $m\rightarrow\infty$. An example for
this is given in the following proposition.

######  Proposition 4.1.

Let $\lim_{m\to\infty}\frac{n(m)}{m}=c$ for some $c>0$ and
$\bm{\Delta}\neq\mathbf{0}$. Then, it holds that

$\lim_{m\to\infty}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{p}}\right)>0.$

The proof of Prop. 4.1 is provided in App. A.1. We note, however, that this
does not happen for a vanishing amount of observational data, that is
$\lim_{m\to\infty}\frac{n(m)}{m}=0$ (see Prop. 4.2 in App. A.2.

#### Ridge Regression.

The ridge regression estimator on the interventional data, which shrinks
$\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ towards the origin (see § 2 and Fig.
2), is given by

$\displaystyle\widehat{\bm{\alpha}}^{m}_{\text{ridge}}$
$\displaystyle=(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}+\lambda\mathbf{I}_{p})^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{y}_{\textsc{i}}$
$\displaystyle=\widehat{\mathbf{W}}^{m}_{\text{ridge}}\widehat{\bm{\alpha}}^{m}_{\textsc{i}}+(\mathbf{I}_{p}-\widehat{\mathbf{W}}^{m}_{\text{ridge}})\mathbf{0},$

where

$\widehat{\mathbf{W}}^{m}_{\text{ridge}}\coloneqq(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}+\lambda\mathbf{I}_{p})^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}.$
(12)

Hence, $\widehat{\bm{\alpha}}^{m}_{\text{ridge}}$ can also be seen as a
special case of a matrix weighted estimator with no observational data and
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}=\mathbf{0}$. Further, comparing 11 and
12 suggest an interpretation of ridge regression as a poor man’s data pooling
since access to observational data is replaced by a positive definite data
matrix $\lambda\mathbf{I}_{p}$. However, $\lambda$ is a constant, and
therefore
$\lim_{m\to\infty}\mathrm{MSE}(\widehat{\bm{\alpha}}^{m}_{\text{ridge}})=0$
even in the setting of Prop. 4.1, which contrasts data pooling.

###  4.2 Optimal Weighting Schemes

We now establish theoretically optimal weighting schemes that minimize the
mean squared error of $\mathbf{W}$-weighted linear estimators
$\widehat{\bm{\alpha}}^{m}_{\mathbf{W}}$ for different classes of weight
matrices $\mathbf{W}$ by exploiting the specific structure of our problem
setting (§ 3).

#### Optimal Scalar Weight.

First, we consider the special case of scalar estimator weighting by
considering weight matrices of the form $\mathbf{W}=w\mathbf{I}_{p}$ with
weight $w\in[0,1]$. The optimal scalar weight is then derived as follows:

$\displaystyle\frac{\partial}{\partial
w}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)\overset{!}{=}0$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:bias_variance}}}{{\iff}}$
$\displaystyle\frac{\partial}{\partial
w}\left(\left\lVert\mathbb{E}[\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}-\bm{\alpha}]\right\rVert_{2}^{2}+\text{Tr}\left(\mathbf{Cov}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)\right)\right)\overset{!}{=}0$
$\displaystyle\implies$ $\displaystyle
w^{m}_{*}=\frac{\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}))+\left\lVert\bm{\Delta}\right\rVert_{2}^{2}}{\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}))+\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}))+\left\lVert\bm{\Delta}\right\rVert_{2}^{2}}.$

#### Optimal Diagonal Weight Matrix.

A more general case is to weigh each dimension individually by different
scalars $w^{(k)}\in[0,1],\;k=1,...,p$, corresponding to a weight matrix of the
form $\mathbf{W}=\mathrm{diag}(\mathbf{w})$. The optimal diagonal weighting
$\text{diag}(\mathbf{w}^{m}_{*}$) is then given by

$w_{*}^{m(k)}=\frac{\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\Delta^{(k)\,2}}{\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\Delta^{(k)\,2}}\,,$

for $k=1,...,p$. The derivation is analogous to that for the optimal scalar
weight above, with the only difference being that we optimize over each
dimension separately.

#### Optimal Weight Matrix.

Finally, we can also determine the optimum weighting as follows:

$\displaystyle\begin{split}\mathbf{W}^{m}_{*}=&\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\bm{\Delta}\bm{\Delta}^{\top}\right)\\\
&\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\bm{\Delta}\bm{\Delta}^{\top}\right)^{-1}.\end{split}$
(13)

A thorough derivation of the proposed weighting schemes can be found in App.
C. In addition, we elaborate on how this weighting scheme handles sample
imbalance in App. E.

######  Remark 4.2.

If (i) $\bm{\Delta}=\mathbf{0}$ and (ii)
$\sigma^{2}_{Y|\mathbf{X}}=\sigma^{2}_{Y|\text{do}(\mathbf{X})}$, then
$\mathbf{W}^{m}_{*}=\mathbf{W}^{m}_{\textsc{p}}$, i.e., data pooling
corresponds to weighing with the optimal weight matrix under these two
assumptions.

Remark 4.2 can be verified by simplifying 13 with assumptions (i) and (ii) and
comparing to 11. It agrees with our intuition: Ordinary least squares relies
on the assumption that
$\mathbb{E}[Y|\mathbf{X}=\mathbf{x}_{i}]=\bm{\alpha}^{\top}\mathbf{x}_{i}$
with equal variance, for all $i$. Thus, data pooling recovers the optimal
estimator if these assumptions are true, i.e., the two conditional
distributions $\mathbb{P}_{\text{obs}}(Y|X)$ and
$\mathbb{P}_{\text{int}}(Y|\mathrm{do}(X))$ are identical. However, in
general, they will not be identical and data pooling then amounts to model
misspecification. This is likely to result in a non-vanishing mean squared
error for $m\rightarrow\infty$ as highlighted in Prop. 4.1.

###  4.3 Practical Estimators

Unfortunately, the optimal weighting derived in 13 cannot be implemented
directly, since the quantities $\bm{\Delta}$,
$\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$, and
$\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$ are unknown in
practice. To construct practical estimators informed by our theoretical
insights, one option is thus to rely on plug-in estimates of these unknown
quantities. For $\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$ and
$\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$, we use the standard
estimators

$\displaystyle\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$
$\displaystyle=(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\widehat{\sigma}^{2}_{Y|\mathrm{do}(\mathbf{X})},$
$\displaystyle\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$
$\displaystyle=(\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}})^{-1}\widehat{\sigma}^{2}_{Y|\mathbf{X}},$

which replace the conditional variances in 9 by

$\displaystyle\widehat{\sigma}^{2}_{Y|\mathrm{do}(\mathbf{X})}$
$\displaystyle=\frac{1}{m-1}\left\lVert\mathbf{y}_{\textsc{i}}-\mathbf{X}_{\textsc{i}}\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right\rVert_{2}^{2},$
$\displaystyle\widehat{\sigma}^{2}_{Y|\mathbf{X}}$
$\displaystyle=\frac{1}{n-1}\left\lVert\mathbf{y}_{\textsc{o}}-\mathbf{X}_{\textsc{o}}\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right\rVert_{2}^{2}.$

For $\bm{\Delta}$, one may consider using the unbiased estimator

$\widehat{\bm{\Delta}}_{m}=\widehat{\bm{\alpha}}^{n}_{\textsc{o}}-\widehat{\bm{\alpha}}^{m}_{\textsc{i}}.$
(14)

Substituting these into 13 then yields:

$\displaystyle\begin{split}\widehat{\mathbf{W}}^{m}_{*}=&\left(\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\widehat{\bm{\Delta}}_{m}\widehat{\bm{\Delta}}_{m}^{\top}+\epsilon\mathbf{I}_{p}\right)\\\
&\left(\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\widehat{\bm{\Delta}}_{m}\widehat{\bm{\Delta}}_{m}^{\top}+\epsilon\mathbf{I}_{p}\right)^{-1}.\end{split}$
(15)

The regularization with $\epsilon>0$ ensures that the inverse remains stable
even in the large sample limit where
$\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$ and
$\widehat{\mathbf{Cov}}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$ tend to zero.
The reason for instability without such regularization is that
$\mathbf{W}^{m}_{*}$ is not uniquely defined in the infinite sample limit.
With regularization, however, we can guarantee that
$\widehat{\mathbf{W}}^{m}_{*}$ converges to $\mathbf{I}_{p}$ almost surely.

######  Proposition 4.3 (Weight Matrix Convergence).

Let $\lim_{m\to\infty}\frac{n(m)}{m}=c$, for some constant $c>0$. Then,
$\widehat{\mathbf{W}}^{m}_{*}$ from 15 converges almost surely to
$\mathbf{I}_{p}$, i.e.,
$\widehat{\mathbf{W}}^{m}_{*}\xrightarrow{\text{a.s.}}\mathbf{I}_{p}$.

The proof for Prop. 4.3 is included in App. A.3. We can show that this
convergence implies that the mean squared error vanishes asymptotically.

######  Theorem 4.4 (Zero Mean Squared Error in the Sample Limit).

Let $\widehat{\mathbf{W}}^{m}$ be any sequence of random weight matrices such
that $\widehat{\mathbf{W}}^{m}\xrightarrow{\text{a.s.}}\mathbf{I}_{p}$ and
$\lim_{m\to\infty}\frac{n(m)}{m}=c$ for some constant $c>0$. Then,

$\underset{m\rightarrow\infty}{\text{lim}}\;\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\right)=0,$

where $\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}$ denotes the
matrix-weighted linear estimator with weight matrix
$\widehat{\mathbf{W}}^{m}$, as defined in Def. 4.1.

The proof of Thm. 4.4 is included in App. A.4.

Thm. 4.4 has the following relevant implication: we can incorporate an
arbitrarily large amount of biased observational data and are still guaranteed
that the bias (and also variance) of
$\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}_{*}}$ will vanish in the
infinite sample limit. Moreover, this guarantee is independent of
$\bm{\Delta}$ and
$|\sigma_{Y|\mathbf{X}}^{2}-\sigma_{Y|\mathrm{do}(\mathbf{X})}^{2}|$.

We also note that Thm. 4.4 does not imply unbiasedness of
$\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}_{*}}$ for any finite
sample size.

Further, we note that almost sure convergence of $\widehat{\mathbf{W}}^{m}$ to
$\mathbf{I}_{p}$ may generally not be the only option to achieve vanishing
mean squared error. For example, if $\bm{\Delta}=\mathbf{0}$ such that
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ is unbiased, we also obtain vanishing
mean squared error for almost sure convergence of $\widehat{\mathbf{W}}^{m}$
to $\mathbf{0}$.

###  4.4 Suitable Inductive Biases

Despite the desirable performance established in Thm. 4.4, the plug-in
estimates from § 4.3 will often not perform very well in finite sample
settings. The main issue is the estimation of $\bm{\Delta}$, which has a large
variance when done according to 14. To see this, we first note that

$\text{Tr}(\mathbf{Cov}(\widehat{\bm{\Delta}}_{m}))=\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}))+\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})),$
(16)

since the observational and interventional data are independent. Now, if we
only have a small interventional sample (as is typically the case),
$\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}))$ and hence
according to 16 also $\text{Tr}(\mathbf{Cov}(\widehat{\bm{\Delta}}_{m}))$ will
be large.

We therefore explore possible inductive biases in the form of additional
assumptions on the type of confounding that lead to reduced variance when
estimating $\widehat{\bm{\Delta}}_{m}$. These inductive biases can be
motivated from domain knowledge and validation techniques such as cross-
validation (Schaffer, 1993). Specifically, the application itself may provide
some prior knowledge about the nature of confounding, which can then be
confirmed by a better validation score compared to the other inductive
biases/methods proposed here.

To this end, we observe that 14 can be written as the solution of the
following two-step ordinary least squares procedure:

$\displaystyle\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$
$\displaystyle\leftarrow\;\mathrm{arg}\;\min_{\bm{\alpha}\in\mathbb{R}^{p}}\;\left\\{\left\lVert\mathbf{y}_{\textsc{o}}\;-\;\mathbf{X}_{\textsc{o}}\bm{\alpha}\right\rVert_{2}^{2}\right\\}$
$\displaystyle\mathbf{r}$
$\displaystyle\leftarrow\;\mathbf{y}_{\textsc{i}}-\mathbf{X}_{\textsc{i}}\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$
$\displaystyle\widehat{\bm{\Delta}}_{m}$
$\displaystyle\leftarrow\;\mathrm{arg}\;\min_{\bm{\Delta}\in\mathbb{R}^{p}}\;\left\\{\left\lVert\mathbf{r}\;+\;\mathbf{X}_{\textsc{i}}\bm{\Delta}\right\rVert_{2}^{2}\right\\}.$
(17)

#### Small $\left\lVert\bm{\Delta}\right\rVert_{2}$.

In some settings, we may be willing to assume that, despite the existence of
unobserved confounders, the resulting confounding bias is rather weak, i.e.,
that its Euclidean norm $\left\lVert\bm{\Delta}\right\rVert_{2}$ is small.
Since this is precisely the assumption underlying ridge regression, we
reformulate (17) using a regularizer $\lambda_{\ell^{2}}>0$ as

$\widehat{\bm{\Delta}}_{m}^{\ell^{2}}\;\leftarrow\;\mathrm{arg}\;\min_{\bm{\Delta}\in\mathbb{R}^{p}}\;\left\\{\left\lVert\mathbf{r}\;+\;\mathbf{X}_{\textsc{i}}\bm{\Delta}\right\rVert_{2}^{2}\;+\;\lambda_{\ell^{2}}\left\lVert\bm{\Delta}\right\rVert_{2}^{2}\right\\},$

for which a closed-form solution of the same computational complexity as least
squares exists. We refer to the weight matrix estimate obtained by using
$\widehat{\bm{\Delta}}_{m}^{\ell^{2}}$ in place of $\widehat{\bm{\Delta}}_{m}$
in 15 as $\widehat{\mathbf{W}}^{m}_{\ell^{2}}$. By Prop. 4.5, we still obtain
the same limiting guarantees of Thm. 4.4 for
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$, as long as $\lambda_{\ell^{2}}$ is
fixed ($\lambda_{\ell^{2}}$ is independent of $m$, $\mathbf{X}_{\textsc{p}}$,
$\mathbf{y}_{\textsc{p}}$).

######  Proposition 4.5.

Let $\lim_{m\to\infty}\frac{n(m)}{m}=c$ and $\lambda_{\ell^{2}}>0$ be fixed.
Then,

$\underset{m\rightarrow\infty}{\text{lim}}\;\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}_{\ell^{2}}}\right)=0.$

The proof for Proposition 4.5 is given in App. A.5.

#### Small $\left\lVert\bm{\Delta}\right\rVert_{0}$.

In other settings, we may have prior beliefs that only some treatment
variables $X_{i}$ are confounded, i.e., that the number of nonzero elements of
$\bm{\Delta}$, denoted by $\left\lVert\bm{\Delta}\right\rVert_{0}$, is small.
If we are unaware of which treatments are confounded, but $p$ is small, we can
simply fit all $2^{p}$ possible models or use best subset selection (James et
al., 2013, p. 205). For larger $p$, a more efficient technique known as the
LASSO employs $\ell^{1}$-regularization and has become a standard tool
(Tibshirani, 1996). For the LASSO, approximate optimization techniques exist
that have a computational complexity of $\mathcal{O}(p^{2}n)$ (Efron et al.,
2004), which is of the same order as ordinary least squares. In this case, we
reformulate (17) as

$\widehat{\bm{\Delta}}_{m}^{\ell^{1}}\;\leftarrow\;\mathrm{arg}\;\min_{\bm{\Delta}\in\mathbb{R}^{p}}\;\left\\{\left\lVert\mathbf{r}\;+\;\mathbf{X}_{\textsc{i}}\bm{\Delta}\right\rVert_{2}^{2}\;+\;\lambda_{\ell^{1}}\left\lVert\bm{\Delta}\right\rVert_{1}\right\\},$

for some $\lambda_{\ell^{1}}>0$, and where
$\left\lVert\;\cdot\;\right\rVert_{1}$ denotes the $\ell^{1}$-norm. We refer
to the weight matrix obtained by using $\widehat{\bm{\Delta}}_{m}^{\ell^{1}}$
in place of $\widehat{\bm{\Delta}}_{m}$ in 15 as
$\widehat{\mathbf{W}}^{m}_{\ell^{1}}$.

##  5 Experiments

Table 1: Mean squared error for the causal effect parameter $\bm{\alpha}$
using various weighting schemes for different types of confounding. The
standard plug-in optimal weight matrix estimator
$\widehat{\mathbf{W}}^{m}_{*}$ generally does not perform well, while
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ and
$\widehat{\mathbf{W}}^{m}_{\ell^{1}}$, which benefit from prior knowledge,
outperform prior work. Note that $\mathbf{W}^{m}_{*}$ is an oracle that is
generally not computable in practice. Numbers correspond to mean $\pm$ std.
dev. over 1000 runs; the best method is highlighted in bold.

|  | $\widehat{w}^{m}_{\text{rm}}$ | $\mathbf{W}^{m}_{\textsc{i}}$ | $\mathbf{W}^{m}_{\textsc{p}}$ | $\widehat{\mathbf{W}}^{m}_{*}$ | $\widehat{\mathbf{W}}^{m}_{\ell^{1}}$ | $\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ | ${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>
---|---|---|---|---|---|---|---|---
spread conf. | $\gamma=1$ | $\mathbf{0.07}\pm 0.02$ | $0.21\pm 0.06$ | $\mathbf{0.07}\pm 0.01$ | $0.21\pm 0.06$ | $0.10\pm 0.04$ | $0.08\pm 0.03$ | $0.04\pm 0.01$
$\gamma=5$ | $0.89\pm 0.20$ | $2.79\pm 0.78$ | $0.92\pm 0.14$ | $2.77\pm 0.77$ | $1.11\pm 0.42$ | $\mathbf{0.76}\pm 0.29$ | ${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.03}$
sparse conf. | $\gamma=1$ | $0.12\pm 0.02$ | $0.21\pm 0.06$ | $0.13\pm 0.02$ | $0.21\pm 0.06$ | $\mathbf{0.10}\pm 0.04$ | $0.16\pm 0.05$ | $0.05\pm 0.01$
$\gamma=5$ | $1.80\pm 0.37$ | $2.79\pm 0.78$ | $2.42\pm 0.24$ | $2.77\pm 0.77$ | $\mathbf{0.95}\pm 0.48$ | $2.28\pm 0.63$ | $0.30\pm 0.08$

0.31.03.010.0$\gamma=1$$\gamma=5$50010001500500100015000.030.100.301.00$m$MSE0.50.71.0$\gamma=1$$\gamma=5$0.31.03.00.31.03.00.030.050.10$\frac{n}{m}$$\mathbf{W}^{m}_{\textsc{i}}$$\mathbf{W}^{m}_{\textsc{p}}$$\widehat{\mathbf{W}}^{m}_{\ell^{1}}$$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$$\widehat{w}^{m}_{\mbox{rm}}$
Figure 3: Performance for varying dataset sizes and ratios. (Left) All methods
improve as the amount of data is increased. More sophisticated weighting
schemes outperform the purely interventional ($\mathbf{W}^{m}_{\textsc{i}}$)
and plug in estimators ($\widehat{\mathbf{W}}^{m}_{*}$) (not depicted due to
close performance overlap with pure interventional), whereas data pooling
($\mathbf{W}^{m}_{\textsc{p}}$) works well only for small $m$ and $\gamma$.
(Right) When keeping $m$ fixed and adding more observational data,
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ clearly works best in strongly
confounded ($\gamma=5$) settings. MSE and $\frac{n}{m}$ are plotted on a
$\text{log}_{10}$ scale. Shaded areas indicate $\pm 0.5$ standard deviations.

We investigate the empirical behavior of our proposed matrix weighted
estimators in a finite sample setting and compare them with baselines and
existing methods through simulations on synthetic data.222The source code for
all experiments is available at:
https://github.com/rudolfwilliam/matrix_weighted_linear_estimators To this
end, we consider different experimental settings in which we vary the strength
and sparsity of confounding, as well as the ratio and absolute quantity of
observational and interventional data.

#### Compared Methods.

We report the mean squared error attained by the theoretically optimal weight
matrix $\mathbf{W}^{m}_{*}$ from (13) as an oracle, as well as the plug-in
estimator $\widehat{\mathbf{W}}^{m}_{*}$ thereof from (15), and the
regularized regression-based $\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ and
$\widehat{\mathbf{W}}^{m}_{\ell^{1}}$ from § 4.4. For the latter two, we
choose the regularization hyperparameters $\lambda_{\ell^{2}}$ and
$\lambda_{\ell^{1}}$ by cross-validation on the interventional data. As
baselines, we consider only using interventional data
($\mathbf{W}^{m}_{\textsc{i}}=\mathbf{I}_{p}$) and data pooling according to
$\mathbf{W}^{m}_{\textsc{p}}$ from (11). We also compare to the Rosenman et
al. (2020) scalar weighting scheme which was proposed for vectors of binary
treatment effects and is given by
$\mathbf{W}=\widehat{w}^{m}_{\text{rm}}\mathbf{I}_{p}$ with

$\widehat{w}^{m}_{\text{rm}}\;:=\;\max\left\\{1\;-\;\frac{\text{Tr}\left(\widehat{\textbf{Cov}}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right)\right)}{\left\lVert\widehat{\bm{\alpha}}^{m}_{\textsc{i}}-\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right\rVert_{2}^{2}},\;0\right\\}.$

We emphasize that other commonly used methods for causal effect estimation
from observational data such as propensity score matching (Imai and Dyk, 2004)
are not applicable, because they require the relevant confounders to be
observed, which is not the case in our setting.

#### General Setup.

In all experiments, we use $p=30$ treatments, a one-dimensional ($d=1$)
confounder $Z$, and unit/isotropic (co)variances:
$\sigma_{N_{Y}}^{2}=\sigma_{N_{Z}}^{2}=1$,
$\bm{\Sigma}_{\mathbf{N}_{\mathbf{X}}}=\mathbf{I}_{p}$. We sample
$\tilde{\mathbf{N}}_{\mathbf{X}}\sim\mathcal{N}(\mathbf{0},\mathbf{Cov}(\mathbf{X}_{\textsc{o}}))$,
$\bm{\alpha}\sim\mathcal{N}(\mathbf{0},9\mathbf{I}_{p})$, and choose
$\mathbf{b}$ and $\gamma$ depending on the settings described below. Unless
otherwise specified, we then draw $m=300$ interventional and $n=600$
observational examples from $\mathbb{P}_{\text{int}}$ and
$\mathbb{P}_{\text{obs}}$, respectively, and compute estimates of
$\bm{\alpha}$ using the different weighting approaches. We repeat this
procedure $1000$ times and report the resulting mean and standard deviation of
the mean squared error.

#### Different Types of Confounding.

In our main experiment, we investigate how estimators perform under different
types of confounding encoded by (2) and (3), specifically by the parameters
$\mathbf{b}\in\mathbb{R}^{p}$ and $\gamma\in\mathbb{R}$ (for a scalar
confounder $Z$). For spread confounding, we sample
$\mathbf{b}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{p})$ such that the
confounder affects all treatment variables almost surely. For sparse
confounding, we sample $b^{(k)}\sim\mathcal{N}(0,1)$ for $k=1,...,5$, and
$b^{(k)}=0$ otherwise, such that only the first five treatments are
confounded. In both cases, we investigate $\gamma\in\\{1,5\\}$ which controls
the strength of $Z\to\leavevmode\nobreak\ Y$ and thus the extent to which
$\bm{\Delta}=0$ is violated.

#### Main Results.

The results are presented in Tab. 1. We find that our regularized estimators
generally perform well, particularly when the underlying assumptions are
satisfied: under sparse confounding $\widehat{\mathbf{W}}^{m}_{\ell^{1}}$
works best, and in the spread confounding case
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ is only narrowly outperformed by
$\widehat{w}^{m}_{\text{rm}}$ and $\mathbf{W}^{m}_{\textsc{p}}$ when
$\gamma=1$. Data pooling works relatively well when $\gamma=1$ (compared to
$\gamma=5$) where the violation of the identically distributed assumption is
weak and the variance from estimating unknown quantities is not compensated by
the bias reduction. In contrast, both the purely interventional approach
$\mathbf{W}^{m}_{\textsc{i}}$ and the plug-in estimator
$\widehat{\mathbf{W}}^{m}_{*}$ do not perform very well in this finite sample
setting due to high variance, as explained in § 4.4.

#### Varying Data Set Sizes and Ratios.

In Fig. 3, we investigate how the different estimators behave across different
data set sizes and ratios for the spread confounding setting. In the left two
plots, we vary the amount of interventional data $m$ while fixing the amount
of observational data to $n=3m$. The results confirm our theoretical results:
For small data set sizes, data pooling is a worthwhile alternative to more
sophisticated weights, in particular if the violation against the assumption
of identical distribution is minor ($\gamma=1$). However, for large enough
data set sizes, the approaches from both previous work and ours achieve a
better score. Particularly, we see that $\widehat{\mathbf{W}}^{m}_{\ell^{2}}$
outperforms all other weights in both scenarios for large enough data sets.

In the right two plots, we keep $m=500$ fixed and change $n$ and thus the
ratio of interventional to observational data. Unsurprisingly, we find that
the mean squared error of $\mathbf{W}^{m}_{\textsc{i}}$ remains constant. For
strong confounding ($\gamma=5$), we see that
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}$ adapts best with a considerable margin:
Unlike $\widehat{w}^{m}_{\text{rm}}$, it explicitly takes into account (an
estimate of) the covariance structure of
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ in constructing the weight matrix.

##  6 Discussion

#### Connection to Transfer Learning.

Our setting bears resemblance to transfer and multi-task learning (Thrun,
1995; Caruana, 1997), specifically to supervised domain adaptation, which aims
to leverage knowledge from a source domain to improve a model in a target
domain, for which typically much less data is available. In our case, we aim
to use the source model $\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$, learned by
estimating $\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]$ in the observational setting,
to improve our (high-variance) target model
$\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ of
$\mathbb{E}[Y|\text{do}(\mathbf{X}\leftarrow\mathbf{x})]$. Transfer learning
can only work if the domains are sufficiently similar, resulting in numerous
approaches leveraging different assumptions about shared components
(Quiñonero-Candela et al., 2008). These assumptions are often phrased in
causal terms (Schölkopf et al., 2012; Zhang et al., 2013; Gong et al., 2016;
Rojas-Carulla et al., 2018). Similarly, our observational (source) and
interventional (target) domains share the same causal model and only differ in
the treatment assignment mechanisms (2) and (4). Still, the bias in (6) can in
theory be arbitrary large, and our methods from § 4.4 implicitly rely on it
being small or sparse.

#### Beyond Linear Regression.

Some of our derivations and theoretical results rely on the fact that the
confounding bias in (5) is linear in $\mathbf{x}$. For the class of linear
SCMs (1)–(3), Gaussianity is necessary and sufficient333Note
$\mathbb{E}[Y|\mathbf{X}]=\bm{\gamma}^{\top}\mathbb{E}[\mathbf{Z}|\mathbf{X}]+\bm{\alpha}^{\top}\mathbf{X}$
and $\mathbb{E}[\mathbf{Z}|\mathbf{X}]$ is linear in $\mathbf{X}$ only in the
Gaussian case (Peters et al., 2017, Thm. 4.2). for this condition to hold, but
it may also hold for more general classes of SCMs. For binary treatments
$\mathbf{X}\in\\{0,1\\}^{p}$, in particular, it is always possible to write
the difference between the biased and unbiased average treatment effect
estimates using a constant offset $\bm{\Delta}$ akin to (14), irrespective of
the confounding relationship.444Specifically, we have
$\bm{\Delta}=\mathbb{E}[Y|\mathbf{X}=\mathbf{1}]-\mathbb{E}[Y|\mathbf{X}=\mathbf{0}]-\left(\mathbb{E}[Y|\text{do}(\mathbf{X}\leftarrow\mathbf{1})]-\mathbb{E}[Y|\text{do}(\mathbf{X}\leftarrow\mathbf{0})]\right)$.
Future work may thus investigate nonlinear extensions, e.g., by drawing
inspiration from semi-parametrics (Robins and Rotnitzky, 1995), doubly robust
estimation (Bang and Robins, 2005), and debiased machine learning
(Chernozhukov et al., 2018).

#### Incorporating Covariates.

Our current formulation does not explicitly account for observed confounders,
or pre-treatment covariates, which need to be adjusted for in the
observational setting to avoid introducing further bias. In principle, such
covariates can simply be included in $\mathbf{X}$, as different treatment
components $X_{i}$ are allowed to be dependent. However, this may result in
high-dimensional treatments and thus render full randomization in (4)
unrealistic. Other covariates, while unproblematic with regard to bias, may
help further reduce variance (Henckel et al., 2022). Extending our framework
to incorporate different types of covariates is thus a worthwhile future
direction.

##  7 Conclusion

In the present work, we have introduced a new class of matrix weighted linear
estimators for learning causal effects of continuous treatments from finite
observational and interventional data. Here, our focus has been on optimizing
statistical efficiency, which complements the vast causal inference literature
on identification from heterogeneous data. Our estimators are connected to
classical ideas from shrinkage estimation applied to causal learning and
provide a unifying account of data pooling and ridge regression, which emerge
as special cases. We show that our estimators are theoretically grounded and
compare favorably to baselines and prior work in simulations. While we
restricted our analysis to linear models for now, we hope that the insights
and methods developed here will also be useful for a broader class of causal
models and transfer learning problems.

###### Acknowledgements.

We thank the anonymous reviewers for useful comments and suggestions that
helped improve the manuscript. We thank the Branco Weiss Fellowship,
administered by ETH Zurich, for the support. This work was further supported
by the Tübingen AI Center and by the German Research Foundation (DFG) under
Germany’s excellence strategy – EXC number 2064/1 – project number 390727645.

## References

  * Angrist and Pischke (2009) J. D. Angrist and J.-S. Pischke.  _Mostly harmless econometrics: An empiricist’s companion_.  Princeton University Press, 2009.
  * Angrist et al. (1996) J. D. Angrist, G. W. Imbens, and D. B. Rubin.  Identification of causal effects using instrumental variables.  _Journal of the American Statistical Association_ , 91(434):444–455, 1996.
  * Bang and Robins (2005) H. Bang and J. M. Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61(4):962–973, 2005.
  * Bareinboim and Pearl (2012) E. Bareinboim and J. Pearl.  Causal Inference by Surrogate Experiments: z-Identifiability.  In _Proceedings of the 28th Conference on Uncertainty in Artificial Intelligence_ , pages 113–120, 2012.
  * Bareinboim and Pearl (2016) E. Bareinboim and J. Pearl.  Causal inference and the data-fusion problem.  _Proceedings of the National Academy of Sciences_ , 113(27):7345–7352, 2016.
  * Caruana (1997) R. Caruana.  Multitask Learning.  _Machine Learning_ , 28(1):41–75, 1997.
  * Ćevid et al. (2020) D. Ćevid, P. Bühlmann, and N. Meinshausen.  Spectral Deconfounding via Perturbed Sparse Linear Models.  _The Journal of Machine Learning Research_ , 21(1):9442–9482, 2020.
  * Cheng and Cai (2021) D. Cheng and T. Cai.  Adaptive Combination of Randomized and Observational Data.  _arXiv:2111.15012_ , 2021.
  * Chernozhukov et al. (2018) V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and structural parameters: Double/debiased machine learning.  _The Econometrics Journal_ , 21(1), 2018.
  * Colnet et al. (2020) B. Colnet, I. Mayer, G. Chen, A. Dieng, R. Li, G. Varoquaux, J.-P. Vert, J. Josse, and S. Yang.  Causal inference methods for combining randomized trials and observational studies: a review.  _arXiv:2011.08047_ , 2020.
  * Correa and Bareinboim (2020) J. Correa and E. Bareinboim.  General transportability of soft interventions: Completeness results.  _Advances in Neural Information Processing Systems_ , 33:10902–10912, 2020.
  * Eberhardt and Scheines (2007) F. Eberhardt and R. Scheines.  Interventions and causal inference.  _Philosophy of science_ , 74(5):981–995, 2007.
  * Efron (2012) B. Efron.  _Large-scale inference: empirical Bayes methods for estimation, testing, and prediction_ , volume 1.  Cambridge University Press, 2012.
  * Efron and Morris (1973) B. Efron and C. Morris.  Stein’s estimation rule and its competitors—an empirical Bayes approach.  _Journal of the American Statistical Association_ , 68(341):117–130, 1973.
  * Efron et al. (2004) B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani.  Least Angle Regression.  _The Annals of Statistics_ , 32(2), 2004.
  * Fisher (1936) R. A. Fisher.  Design of experiments.  _British Medical Journal_ , 1(3923):554, 1936.
  * Gong et al. (2016) M. Gong, K. Zhang, T. Liu, D. Tao, C. Glymour, and B. Schölkopf.  Domain adaptation with conditional transferable components.  In _International Conference on Machine Learning_ , pages 2839–2848, 2016.
  * Green and Strawderman (1991) E. J. Green and W. E. Strawderman.  A James-Stein Type Estimator for Combining Unbiased and Possibly Biased Estimators.  _Journal of the American Statistical Association_ , 86(416):1001–1006, 1991.
  * Green et al. (2005) E. J. Green, W. E. Strawderman, R. L. Amateis, and G. A. Reams.  Improved Estimation for Multiple Means with Heterogeneous Variances.  _Forest Science_ , 51(1):1–6, 2005.
  * Hastie et al. (2009) T. Hastie, R. Tibshirani, and J. Friedman.  _The Elements of Statistical Learning_.  Springer, 2009.
  * Hatt et al. (2022) T. Hatt, J. Berrevoets, A. Curth, S. Feuerriegel, and M. van der Schaar.  Combining observational and randomized data for estimating heterogeneous treatment effects.  _arXiv:2202.12891_ , 2022.
  * Henckel et al. (2022) L. Henckel, E. Perković, and M. H. Maathuis.  Graphical criteria for efficient total effect estimation via adjustment in causal linear models.  _Journal of the Royal Statistical Society Series B_ , 84(2):579–599, 2022.
  * Hernán and Robins (2020) M. A. Hernán and J. M. Robins.  _Causal inference: What if_.  Boca Raton: Chapman & Hall/CRC, 2020.
  * Hoerl (1970) A. E. Hoerl.  Ridge Regression: Biased Estimation for Nonorthogonal Problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Huang and Valtorta (2006) Y. Huang and M. Valtorta.  Identifiability in Causal Bayesian Networks: A Sound and Complete Algorithm.  In _Proceedings of the National Conference on Artificial Intelligence_ , volume 21, pages 1149–1154, 2006.
  * Ilse et al. (2021) M. Ilse, P. Forré, M. Welling, and J. M. Mooij.  Combining Interventional and Observational Data Using Causal Reductions.  _arXiv:2103.04786_ , pages 1–42, 2021.
  * Imai and Dyk (2004) K. Imai and D. A. V. Dyk.  Causal Inference with General Treatment Regimes: Generalizing the Propensity Score.  _Journal of the American Statistical Association_ , 99(467):854–866, 2004.
  * Imbens and Rubin (2015) G. W. Imbens and D. B. Rubin.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press, 2015.
  * James et al. (2013) G. James, D. Witten, T. Hastie, and R. Tibshirani.  _An Introduction to Statistical Learning_.  Springer, 2013.
  * James and Stein (1961) W. James and C. Stein.  Estimation with Quadratic Loss.  In _Proceedings of the 4th Berkeley Symposium on Probability and Statistics_. Berkeley, CA: University of California Press, 1961.
  * Kallus et al. (2018) N. Kallus, A. M. Puli, and U. Shalit.  Removing hidden confounding by experimental grounding.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Lee et al. (2020) S. Lee, J. D. Correa, and E. Bareinboim.  General Identifiability with Arbitrary Surrogate Experiments.  In _Proceedings of the 35th Uncertainty in Artificial Intelligence Conference_ , pages 389–398, 2020.
  * Morgan and Winship (2014) S. L. Morgan and C. Winship.  _Counterfactuals and Causal Inference: Methods and Principles for Social Research_.  Cambridge University Press, 2014.
  * Neyman (1923) J. Neyman.  On the application of probability theory to agricultural experiments: essay on principles.  _Statistical Science_ , 5:465–480, 1923.
  * Pearl (1995) J. Pearl.  Causal diagrams for empirical research.  _Biometrika_ , 82(4):669–688, 1995.
  * Pearl (2009) J. Pearl.  _Causality: models, reasoning, and inference_.  Cambridge University Press, 2nd edition, 2009.
  * Pearl and Bareinboim (2014) J. Pearl and E. Bareinboim.  External Validity: From Do-Calculus to Transportability Across Populations.  _Statistical Science_ , 29(4):579–595, 2014.
  * Peters et al. (2017) J. Peters, D. Janzing, and B. Schölkopf.  _Elements of Causal Inference: Foundations and Learning Algorithms_.  MIT Press, 2017.
  * Quiñonero-Candela et al. (2008) J. Quiñonero-Candela, M. Sugiyama, A. Schwaighofer, and N. D. Lawrence.  _Dataset Shift in Machine Learning_.  MIT Press, 2008.
  * Reichenbach (1956) H. Reichenbach.  _The Direction of Time_ , volume 65.  University of California Press, 1956.
  * Robbins (1964) H. Robbins.  The empirical Bayes approach to statistical decision problems.  _The Annals of Mathematical Statistics_ , 35(1):1–20, 1964.
  * Robins and Rotnitzky (1995) J. M. Robins and A. Rotnitzky.  Semiparametric efficiency in multivariate regression models with missing data.  _Journal of the American Statistical Association_ , 90(429):122–129, 1995.
  * Rojas-Carulla et al. (2018) M. Rojas-Carulla, B. Schölkopf, R. Turner, and J. Peters.  Invariant Models for Causal Transfer Learning.  _The Journal of Machine Learning Research_ , 19(1):1309–1342, 2018.
  * Rosenman et al. (2020) E. Rosenman, G. Basse, A. Owen, and M. Baiocchi.  Combining observational and experimental datasets using shrinkage estimators.  _Biometrics_ , 2020.
  * Rosenman et al. (2022) E. T. Rosenman, A. B. Owen, M. Baiocchi, and H. R. Banack.  Propensity score methods for merging observational and experimental datasets.  _Statistics in Medicine_ , 41(1):65–86, 2022.
  * Rubin (1974) D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688–701, 1974.
  * Schaffer (1993) C. Schaffer.  Selecting a Classification Method by Cross-Validation.  _Machine Learning_ , 13:135–143, 1993.
  * Schölkopf et al. (2012) B. Schölkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and J. M. Mooij.  On Causal and Anticausal Learning.  In _International Conference on Machine Learning_ , 2012.
  * Shpitser and Pearl (2006) I. Shpitser and J. Pearl.  Identification of Joint Interventional Distributions in Recursive Semi-Markovian Causal Models.  In _Proceedings of the National Conference on Artificial Intelligence_ , volume 21, pages 1219–1226, 2006.
  * Spirtes et al. (2000) P. Spirtes, C. Glymour, and R. Scheines.  _Causation, Prediction, and Search_.  MIT Press, 2000.
  * Stein (1956) C. Stein.  Inadmissibility of the Usual Estimator for the Mean of a Multivariate Normal Distribution.  In _Proceedings of the third Berkeley Symposium on Mathematical Statistics and Probability_ , volume 3, pages 197–207. University of California Press, 1956.
  * Thrun (1995) S. Thrun.  Is Learning The n-th Thing Any Easier Than Learning The First?  _Advances in Neural Information Processing Systems_ , 8, 1995.
  * Tian and Pearl (2002) J. Tian and J. Pearl.  A general identification condition for causal effects.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pages 567–573, 2002.
  * Tibshirani (1996) R. Tibshirani.  Regression Shrinkage and Selection via the Lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Wang et al. (2017) Y. Wang, L. Solus, K. Yang, and C. Uhler.  Permutation-based Causal Inference Algorithms with Interventions.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Wasserman (2006) L. Wasserman.  _All of Nonparametric Statistics_.  Springer, 2006.
  * Yang and Ding (2020) S. Yang and P. Ding.  Combining Multiple Observational Data Sources to Estimate Causal Effects.  _Journal of the American Statistical Association_ , 115(531):1540–1554, 2020.
  * Zhang et al. (2013) K. Zhang, B. Schölkopf, K. Muandet, and Z. Wang.  Domain Adaptation under Target and Conditional Shift.  In _International Conference on Machine Learning_ , pages 819–827, 2013.

Appendix

##  A Proofs

###  A.1 Proposition 4.1

###### Proof.

We begin by observing that we can write $\mathbf{W}^{m}_{\textsc{p}}$ as

$\mathbf{W}^{m}_{\textsc{p}}=\left(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\;+\;\frac{n}{m}n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\right)^{-1}\left(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\right).$
()

We apply the strong law of large numbers to obtain that

$m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\xrightarrow{\text{a.s.}}\mathbf{Cov}(\mathbf{X}_{\textsc{i}})\quad\text{and}\quad
n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\xrightarrow{\text{a.s.}}\mathbf{Cov}(\mathbf{X}_{\textsc{o}}).$

Due to the fact that
$\underset{m\rightarrow\infty}{\text{lim}}\frac{n(m)}{m}=c$ for some $c>0$, we
conclude

$\mathbf{W}^{m}_{\textsc{p}}\;\xrightarrow{\text{a.s.}}\;\mathbf{W}_{\infty}\;:=\;\left(\text{{Cov}}(\mathbf{X}_{\textsc{i}})\;+\;c\cdot\text{{Cov}}(\mathbf{X}_{\textsc{o}})\right)^{-1}\text{{Cov}}(\mathbf{X}_{\textsc{i}}).$

We observe that

$\left(\mathbf{I}-\mathbf{W}_{\infty}\right)=\left(\mathbf{Cov}(\mathbf{X}_{\textsc{i}})+c\cdot\mathbf{Cov}(\mathbf{X}_{\textsc{o}})\right)^{-1}c\cdot\mathbf{Cov}(\mathbf{X}_{\textsc{o}}).$

Since both covariance matrices are positive definite, so is
$\mathbf{Cov}(\mathbf{X}_{\textsc{i}})+c\cdot\mathbf{Cov}(\mathbf{X}_{\textsc{o}})$.
We conclude that the smallest singular value of
$\mathbf{I}-\mathbf{W}_{\infty}$ is strictly greater than 0. This means

$\big{|}\big{|}\mathbb{E}[\widehat{\bm{\alpha}}^{m}_{\mathbf{W}_{\infty}}]-\bm{\alpha}\big{|}\big{|}_{2}^{2}\;=\;||\left(\mathbf{I}_{p}-\mathbf{W}_{\infty}\right)\bm{\Delta}||_{2}^{2}\;\geq\;c^{\prime}||\bm{\Delta}||_{2}^{2},$

for some fixed constant $c^{\prime}>0$. We obtain therefore

$0<\underset{m\rightarrow\infty}{\text{lim}}\big{|}\big{|}\mathbb{E}[\widehat{\bm{\alpha}}^{m}_{\mathbf{W}_{\infty}}]-\bm{\alpha}\big{|}\big{|}_{2}^{2}\leq\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\,\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}_{\infty}}\right),$

where we invoked Jensen’s inequality. We see that $\mathbf{W}_{\infty}$ is
constant and bounded. We note that almost sure convergence implies convergence
in probability. We can thus apply Lemma B.1, which yields the desired result

$0<\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\,\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}_{\infty}}\right)\leq\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\,\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}_{\textsc{p}}}\right).$

∎

###  A.2 Proposition 4.2

######  Proposition 4.2.

Let $\lim_{m\to\infty}\frac{n(m)}{m}=0$. Then, it holds that

$\lim_{m\to\infty}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{p}}\right)=0.$

###### Proof.

Similar to the proof of Proposition 4.1, we employ the formulation of () and
consider the term

$\frac{n}{m}\cdot
n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}.$

We see that $\underset{m\rightarrow\infty}{\text{lim}}\frac{n(m)}{m}=0$ and by
the strong law of large numbers,
$n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\xrightarrow{\text{a.s.}}\mathbf{Cov}(\mathbf{X}_{\textsc{o}})$.
Hence, we obtain that

$\frac{n}{m}\cdot
n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\xrightarrow{\text{a.s.}}\mathbf{0}.$

By the continuous mapping theorem, we conclude that

$\mathbf{W}^{m}_{\textsc{p}}\xrightarrow{\text{a.s.}}\mathbf{I}_{p},$

and by Lemma B.2, this implies that

$\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}_{\textsc{p}}}\right)\;\leq\;\underset{m\rightarrow\infty}{\text{lim}}\text{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right)\;=\;0.$

∎

###  A.3 Proposition 4.3

###### Proof.

We rewrite $\widehat{\mathbf{W}}^{m}_{*}$ as follows:

$\displaystyle\widehat{\mathbf{W}}^{m}_{*}=$
$\displaystyle\left(n^{-1}\left(n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\right)^{-1}\hat{\sigma}^{2}_{Y|X}+\hat{\bm{\Delta}}\hat{\bm{\Delta}}^{\top}+\epsilon\mathbf{I}_{p}\right)$
$\displaystyle\left(n^{-1}\left(n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\right)^{-1}\hat{\sigma}^{2}_{Y|X}+m^{-1}\left(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\right)^{-1}\hat{\sigma}^{2}_{Y|\text{do}(X)}+\hat{\bm{\Delta}}\hat{\bm{\Delta}}^{\top}+\epsilon\mathbf{I}_{p}\right)^{-1},$

where we insert any almost surely converging estimators for $\bm{\Delta}$,
$\sigma_{Y|X}^{2}$ and $\sigma_{Y|\text{do}(X)}^{2}$ instead of their ground-
truth values. By almost sure convergence of linear estimators individually, we
see that this holds specifically for
$\hat{\bm{\Delta}}=\widehat{\bm{\alpha}}^{n}_{\textsc{o}}-\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$.
Also, we can use the strong law of large numbers to conclude almost sure
convergence of $\hat{\sigma}^{2}_{Y|X}$ and
$\hat{\sigma}^{2}_{Y|\text{do}(X)}$.

We now show
$\widehat{\mathbf{W}}^{m}_{*}\xrightarrow{\text{a.s.}}\mathbf{I}_{p}$: First,
we see that

$(cm)^{-1}\left(n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\right)^{-1}\hat{\sigma}^{2}_{Y|X}\xrightarrow{\text{a.s.}}\mathbf{0}\quad\text{and}\quad
m^{-1}\left(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\right)^{-1}\hat{\sigma}^{2}_{Y|\text{do}(X)}\xrightarrow{\text{a.s.}}\mathbf{0},$

since
$m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\;\hat{\sigma}^{2}_{Y|\text{do}(X)}$
and
$n^{-1}\mathbf{X}_{\textsc{o}}^{\top}\mathbf{X}_{\textsc{o}}\;\hat{\sigma}^{2}_{Y|X}$
converge almost surely to constants and $m^{-1}$ vanishes. Hence,

$\widehat{\mathbf{W}}^{m}_{*}\xrightarrow{\text{a.s.}}\left(\bm{\Delta}\bm{\Delta}^{\top}+\epsilon\mathbf{I}_{p}\right)\;\left(\bm{\Delta}\bm{\Delta}^{\top}+\epsilon\mathbf{I}_{p}\right)^{-1}=\mathbf{I}_{p}.$

∎

###  A.4 Theorem 4.4

###### Proof.

We have that $\mathbf{I}_{p}$ is bounded in norm, almost surely. So we can
apply Lemma B.2 to see that

$\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}_{*}}\big{)}\;\leq\;\underset{m\rightarrow\infty}{\text{lim}}\mathrm{MSE}\big{(}\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\big{)}\;=\;0.$

∎

###  A.5 Proposition 4.5

###### Proof.

By Theorem 4.4, it suffices to show that
$\widehat{\mathbf{W}}^{m}_{\ell^{2}}\xrightarrow{\text{a.s.}}\mathbf{I}_{p}$.
Since the other quantities
$\text{{Cov}}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})$,
$\text{{Cov}}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$ for estimating
$\mathbf{W}^{m}_{*}$ remain unchanged compared to
$\widehat{\mathbf{W}}^{m}_{*}$, it suffices to show that the modified
computation of $\widehat{\bm{\Delta}}_{m}$ we call
$\hat{\bm{\Delta}}_{m}^{\ell^{2}}$ converges almost surely to the true
$\bm{\Delta}=\bm{\alpha}_{\textsc{i}}-\bm{\alpha}_{\textsc{o}}$, where
$\bm{\alpha}_{\textsc{i}}$ and $\bm{\alpha}_{\textsc{o}}$ are short-hand for
$\mathbb{E}_{\text{int}}[Y|\mathbf{X}=\mathbf{x}]$ and
$\mathbb{E}_{\text{obs}}[Y|\mathbf{X}=\mathbf{x}]$, respectively. We observe
that $\hat{\bm{\Delta}}_{m}^{\ell^{2}}$ has a closed-form solution

$\displaystyle\hat{\bm{\Delta}}_{m}^{\ell^{2}}$
$\displaystyle=-(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}+\lambda_{\ell^{2}}\mathbf{I}_{p})^{-1}\mathbf{X}_{\textsc{i}}^{\top}(\mathbf{y}_{\textsc{i}}-\mathbf{X}_{\textsc{i}}\widehat{\bm{\alpha}}^{n}_{\textsc{o}})$
()

since $\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ is again a closed-form solution
to an ordinary least squares problem. Considering the first term in
(LABEL:equ:weak_bias_cons), we conclude almost sure convergence with respect
to $\bm{\alpha}_{\textsc{i}}$ (it is simply the ridge regression solution on
the interventional data, which is well-known to converge almost surely for
fixed $\lambda_{\ell^{2}}$). The second term satisfies

$(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}+\lambda_{\ell^{2}}\mathbf{I}_{p})^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}}\xrightarrow{\text{a.s.}}\mathbf{I}_{p}\quad\text{and}\quad\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\xrightarrow{\text{a.s.}}\bm{\alpha}_{\textsc{o}}.$

This leads to the desired conclusion. ∎

##  B Additional Lemmas

######  Lemma B.1.

Let $\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}\xrightarrow{\text{P}}\mathbf{0}$
555We note that $\mathbf{W}^{m}$ may be random. and let there exist $c>0$,
$m^{\prime}\in\mathbb{N}$, such that $||\mathbf{W}^{m}||_{2}\leq c,\;\text{for
all}\;m\geq m^{\prime}$, almost surely. Then, it holds that

$\underset{m\rightarrow\infty}{\text{lim}}\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\big{)}\;\leq\;\underset{m\rightarrow\infty}{\text{lim}}\;\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\big{)},$

where $\xrightarrow{\text{P}}$ denotes convergence in probability.

###### Proof.

We derive a lower bound on
$\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\big{)}$
by using the formulation

$\displaystyle\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\big{)}\;=$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}\;||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]\;+$
()
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right\\}\;||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right],\quad\forall\epsilon>0.$

We bound the second summand of () from below by zero. For the first summand,
we use reverse triangle inequality, which yields

$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]\;=\;\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-(\bm{\alpha}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}})||_{2}^{2}\right]$
() $\displaystyle\geq\;$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]-2\sqrt{\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}||_{2}^{2}\right]\;\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]}\;+$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}||_{2}^{2}\right]$
$\displaystyle\geq\;$
$\displaystyle\text{MSE}(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}})-\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]-$
$\displaystyle
2\sqrt{\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}||_{2}^{2}\right]\;\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]}.$

For any constant $\mathbf{W},\mathbf{W}^{\prime}\in\mathbb{R}^{p\times p}$, we
rewrite

$\displaystyle\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{\prime}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}}||_{2}^{2}\right]$
$\displaystyle=\;$
$\displaystyle\mathbb{E}\left[||(\mathbf{W}^{\prime}-\mathbf{W})\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\;+\;(\mathbf{W}-\mathbf{W}^{\prime})\widehat{\bm{\alpha}}^{n}_{\textsc{o}}||_{2}^{2}\right]$
$\displaystyle\leq$ $\displaystyle
2\left(||\mathbf{W}-\mathbf{W}^{\prime}||_{2}^{2}\text{Tr}\left(\mathbb{E}\left[\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\widehat{\bm{\alpha}}_{\textsc{i}}^{m\,\top}\right]\right)\;+\;||\mathbf{W}-\mathbf{W}^{\prime}||_{2}^{2}\text{Tr}\left(\mathbb{E}\left[\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\widehat{\bm{\alpha}}_{\textsc{o}}^{n\,\top}\right]\right)\right)$
$\displaystyle=$ $\displaystyle
2||\mathbf{W}-\mathbf{W}^{\prime}||_{2}^{2}\Bigg{[}\left(||\mathbb{E}\left[\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right]||_{2}^{2}\;+\;\text{Tr}\left(\text{{Cov}}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right)\right)\right)\;+\;\left(||\mathbb{E}\left[\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right]||_{2}^{2}\;+\;\text{Tr}\left(\text{{Cov}}\left(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right)\right)\right)\Bigg{]},$

where we have used Young’s inequality in the first step. We see that both
$||\mathbb{E}\left[\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right]||_{2}^{2}$
and
$||\mathbb{E}\left[\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right]||_{2}^{2}$
remain bounded $\forall m$, while
$\text{Tr}\left(\text{{Cov}}\left(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}\right)\right)$
and
$\text{Tr}\left(\text{{Cov}}\left(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}\right)\right)$
decrease monotonically in $m$. Hence, we conclude that for any
$\epsilon^{\prime}>0$, there exists an $\epsilon>0$ such that

$\displaystyle\mathbb{E}\left[\big{|}\big{|}\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{\prime}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}}\big{|}\big{|}_{2}^{2}\right]\leq\epsilon^{\prime},\;\forall
m\in\mathbb{N}\;\text{and}\;\forall\mathbf{W},\mathbf{W}^{\prime}\in\mathbb{R}^{p\times
p}\;\text{s.t.}\;||\mathbf{W}-\mathbf{W}^{\prime}||_{2}\leq\epsilon.$ ()

Since $||\mathbf{W}^{m}||_{2}\leq c$ for all $m\geq m^{\prime}$, we have that
$||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}$ is also
bounded by some constant $c^{\prime}>0$, for all $m\geq m^{\prime}$, almost
surely. We now fix an $\epsilon^{\prime}>0$ and choose a corresponding
$\epsilon$ such that () holds. We then conclude from () that

$\displaystyle\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\right)\geq$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]$
$\displaystyle\geq$
$\displaystyle\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)-2\sqrt{\epsilon^{\prime}\;\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]}-P\left(||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right)c^{\prime}$
$\displaystyle\geq$
$\displaystyle\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)-2\sqrt{\epsilon^{\prime}c^{\prime}}-P\left(||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right)c^{\prime},$

for all $m\geq m^{\prime}$. Thus, we conclude

$\underset{m\rightarrow\infty}{\text{lim}}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\right)\;\geq\;\underset{m\rightarrow\infty}{\text{lim}}\;\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)-2\sqrt{\epsilon^{\prime}c^{\prime}}.$

We can repeat this procedure for any $\epsilon^{\prime}>0$ and therefore
conclude

$\underset{m\rightarrow\infty}{\text{lim}}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\right)\;\geq\;\underset{m\rightarrow\infty}{\text{lim}}\;\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right),$

which is the desired result. ∎

######  Lemma B.2.

Let
$\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}\xrightarrow{\text{a.s.}}\mathbf{0}$
and let there exist some $c>0$, $m^{\prime}\in\mathbb{N}$, such that
$||\mathbf{W}^{m}||_{2}\leq c,\forall m\geq m^{\prime}$, almost surely. Then,
it holds that

$\underset{m\rightarrow\infty}{\text{lim}}\;\mathrm{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\big{)}\;\leq\;\underset{m\rightarrow\infty}{\text{lim}}\mathrm{MSE}\big{(}\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\big{)}.$

###### Proof.

We again employ the formulation from (), but this time to construct an upper
bound. For the first term of (), we see that

$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}\;||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]\;=\;\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}\;||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}+\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]$
() $\displaystyle\leq\;$
$\displaystyle\text{MSE}\,\big{(}\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\big{)}\;+\;2\sqrt{\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}||_{2}^{2}\right]\;\mathbb{E}[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}]}\;+$
$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}||_{2}^{2}\right],$

by triangle inequality and the Cauchy-Schwarz inequality. Since for $m\geq
m^{\prime}$ it holds that $||\mathbf{W}^{m}||_{2}\leq c$, almost surely, there
exists a constant $c^{\prime}>0$ such that
$\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]\;\leq\;c^{\prime}$,
for all $m\geq m^{\prime}$. This is true because the two estimators
$\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ and
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ have both bounded mean squared error
for any sample size $m$.

Analogously to the proof for Lemma B.1, we now fix an $\epsilon^{\prime}>0$
and choose a corresponding $\epsilon$ such that () holds. For $m\geq
m^{\prime}$, we then conclude from () that

$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}\leq\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]$
() $\displaystyle\leq$
$\displaystyle\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)+2\sqrt{\epsilon^{\prime}\;\mathbb{E}\left[||\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}-\bm{\alpha}||_{2}^{2}\right]}+\epsilon^{\prime}$
$\displaystyle\leq$
$\displaystyle\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)+2\sqrt{\epsilon^{\prime}c^{\prime}}+\epsilon^{\prime}.$

This bounds the first term of (). For the second term of (), we use almost
sure convergence of $\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}$. Since
$\mathbf{W}^{m}$ is bounded in the limit, almost surely, so is
$\widehat{\mathbf{W}}^{m}$. Formally, $||\widehat{\mathbf{W}}^{m}||_{2}\leq
c^{{}^{\prime\prime}},\forall m\geq m^{\prime}$ for some
$m^{\prime}\in\mathbb{N}$, almost surely.

We use this to bound
$||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}<c^{\prime\prime\prime}$
for all $m\geq m^{\prime}$, almost surely, for some
$c^{\prime\prime\prime}>0$. Now, we apply iterated expectations to the second
term of () to see that for all $m\geq m^{\prime}$

$\displaystyle\mathbb{E}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right\\}||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]\;=$
$\displaystyle\mathbb{E}_{\widehat{\mathbf{W}}^{m}}\left[\mathbbm{1}\left\\{||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right\\}\;\mathbb{E}_{\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}|\widehat{\mathbf{W}}^{m}}\left[||\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}-\bm{\alpha}||_{2}^{2}\right]\right]$
() $\displaystyle\leq$
$\displaystyle\text{P}\left(||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right)c^{\prime\prime\prime},$

almost surely. Now, we can combine the inequalities (LABEL:eq:first_term) and
() to obtain

$\text{MSE}\,\left(\widehat{\bm{\alpha}}^{m}_{\widehat{\mathbf{W}}^{m}}\right)\;\leq\;\text{MSE}\,\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}^{m}}\right)\;+\;2\sqrt{\epsilon^{\prime}c^{\prime}}\;+\;\epsilon^{\prime}\;+\;\text{P}\left(||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right)c^{\prime\prime\prime},$

for all $m\geq m^{\prime\prime}$. Almost sure convergence implies consistency
of $\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}$ with respect to $\mathbf{0}$, so
we see that
$\text{P}\left(||\widehat{\mathbf{W}}^{m}-\mathbf{W}^{m}||_{2}>\epsilon\right)$
vanishes in the limit $m\rightarrow\infty$, for all $\epsilon>0$. We can
repeat this procedure for any $\epsilon^{\prime}>0$. This implies the desired
result. ∎

##  C Detailed Derivation of Optimal Weighting Schemes

In general, we observe that

$\displaystyle{\mathbf{Bias}}(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}})$
$\displaystyle=\mathbf{W}\bm{\alpha}+(\mathbf{I}-\mathbf{W})(\bm{\alpha}+\bm{\Delta})-\bm{\alpha}=(\mathbf{I}-\mathbf{W})\bm{\Delta},$
$\displaystyle\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}})$
$\displaystyle=\mathbf{W}\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})\mathbf{W}^{\top}+(\mathbf{I}-\mathbf{W})\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})(\mathbf{I}-\mathbf{W})^{\top}.$

###  C.1 Optimal Scalar Weight

Here, we have

$\displaystyle\frac{\partial}{\partial
w}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)$
$\displaystyle=\quad$ $\displaystyle\frac{\partial}{\partial
w}\Big{|}\Big{|}{\mathbf{Bias}}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)\Big{|}\Big{|}_{2}^{2}\;+\;\frac{\partial}{\partial
w}\text{Tr}\left(\mathbf{Cov}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)\right)$
$\displaystyle=\quad$
$\displaystyle-2(1-w)||\bm{\Delta}||_{2}^{2}+2w\text{Tr}\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})\right)-2(1-w)\text{Tr}\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})\right)\overset{!}{=}0.$

By rearranging, we get

$w^{m}_{*}=\frac{\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}))+\left\lVert\bm{\Delta}\right\rVert_{2}^{2}}{\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}}))+\text{Tr}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}))+\left\lVert\bm{\Delta}\right\rVert_{2}^{2}}.$

###  C.2 Optimal Diagonal Weight Matrix

Here, we see that the objective decouples into a sum over the individual
dimensions

$\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{w\mathbf{I}_{p}}\right)=\sum_{k=1}^{p}\left(1-w^{(k)}\right)^{2}\bm{\Delta}^{(k)\,2}\;+\;w^{(k)\,2}\mathbf{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\left(1-w^{(k)}\right)^{2}\mathbf{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}}).$

Thus, we optimize for each dimension $k$ separately and obtain

$w_{*}^{m(k)}=\frac{\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\Delta^{(k)\,2}}{\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\mathrm{Cov}^{(k,k)}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\Delta^{(k)\,2}}.$

###  C.3 Optimal Weight Matrix

Using
$\frac{\partial}{\partial\mathbf{W}}\text{Tr}(\mathbf{W}\mathbf{A}\mathbf{W}^{\top})=2\mathbf{W}\mathbf{A}$,
since $\mathbf{A}$ is symmetric, we observe that

$\displaystyle\frac{\partial}{\partial\mathbf{W}}\mathrm{MSE}\left(\widehat{\bm{\alpha}}^{m}_{\mathbf{W}}\right)$
$\displaystyle=\quad$ $\displaystyle
2\mathbf{W}\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\bm{\Delta}\bm{\Delta}^{\top}\right)-2\left(\bm{\Delta}\bm{\Delta}^{\top}+\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})\right)$
$\displaystyle\overset{!}{=}\quad$ $\displaystyle\mathbf{0}.$

We see that this minimum is attained for

$\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\bm{\Delta}\bm{\Delta}^{\top}\right)\left(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\mathbf{Cov}(\widehat{\bm{\alpha}}^{n}_{\textsc{o}})+\bm{\Delta}\bm{\Delta}^{\top}\right)^{-1}.$

##  D Non Zero-Mean Exogenous Variables

All results established here can readily be extended to settings, where any of
the exogenous variables have non-zero mean, i.e.,
$\bm{\mu}_{\mathbf{N}_{\mathbf{X}}}$,
$\bm{\mu}_{\tilde{\mathbf{N}}_{\mathbf{X}}}\coloneqq\mathbb{E}[\tilde{\mathbf{N}}_{\mathbf{X}}]$,
$\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}}$, $\mu_{N_{Y}}$ (see (1)–(3)) may be non-
zero. In order to extend the practical estimators introduced here, one needs
to consider the following two pre-processing steps:

First, we center both treatment distributions separately, without scaling:

$\displaystyle\mathbf{x}^{\prime}_{i}\;\leftarrow\;\mathbf{x}_{i}-n^{-1}\sum_{j\in
1,...,n}\mathbf{x}_{j},\quad$ $\displaystyle\forall i\in 1,...,n,$ ()
$\displaystyle\mathbf{x}^{\prime}_{i}\;\leftarrow\;\mathbf{x}_{i}-m^{-1}\sum_{j\in
n+1,...,n+m}\mathbf{x}_{j},\quad$ $\displaystyle\forall i\in n+1,...,n+m.$ ()

In this manner, both treatment variables become zero-mean.

Furthermore, we add a dummy dimension with value one to all treatment vectors:

$\mathbf{x}^{\prime\prime}_{i}\;\leftarrow\;(\mathbf{x}^{\prime}_{i},\;1),\quad\forall
i\in 1,...,n+m.$

This naturally adds one more dimension also to $\bm{\alpha}$, which
corresponds to the intercept term. We then use the constructed
$\mathbf{x}^{\prime\prime}_{i}$ to compute the weight matrices proposed in
this work.

Finally, we see that the intercept term must be identical for both
distributions, interventional and observational:

$\mathbb{E}[Y\;|\;\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]\;=\;\bm{\gamma}^{\top}\mathbb{E}[\mathbf{Z}\;|\;\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]+\bm{\alpha}^{\top}\mathbf{x}^{\prime}+\mu_{N_{Y}}.$

We then have in the observational setting (data points $1,...,n$) that

$\displaystyle\bm{\gamma}^{\top}\mathbb{E}[\mathbf{Z}\;|\;\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]$
$\displaystyle=\bm{\gamma}^{\top}\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}}+\bm{\gamma}^{\top}\bm{\Sigma}_{\mathbf{N}_{\mathbf{Z}}}\mathbf{B}^{\top}(\bm{\Sigma}_{\mathbf{N}_{\mathbf{X}}}+\mathbf{B}\bm{\Sigma}_{\mathbf{N}_{\mathbf{Z}}}\mathbf{B}^{\top})^{-1}(\mathbf{x}^{\prime}-\mathbb{E}[\mathbf{X}^{\prime}])$
$\displaystyle=\bm{\gamma}^{\top}\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}}+\bm{\Delta}^{\top}\mathbf{x}^{\prime},$

where $\mathbb{E}[\mathbf{X}^{\prime}]=\mathbf{0}$ due to ().

For the interventional data, we have independence between
$\mathbf{X}^{\prime}$ and $\mathbf{Z}$ by definition and so we trivially get

$\bm{\gamma}^{\top}\mathbb{E}[\mathbf{Z}\;|\;\mathbf{X}^{\prime}=\mathbf{x}^{\prime}]=\bm{\gamma}^{\top}\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}}$

here. Thus, the intercept is
$\bm{\gamma}^{\top}\bm{\mu}_{\mathbf{N}_{\mathbf{Z}}}+\mu_{N_{Y}}$ for both
distributions and we fix $\hat{\Delta}^{(p+1)}=0$.

##  E Sample Imbalance

We see that the ground truth covariance matrices of
$\widehat{\bm{\alpha}}^{m}_{\textsc{i}}$ and
$\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ adapt to changes in the sample sizes,
keeping the distributions of all variables fixed. For instance, we see that

$\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})=(\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\sigma_{Y|\text{do}(X)}^{2}=m^{-1}(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\sigma_{Y|\text{do}(X)}^{2}.$

The term
$(m^{-1}\mathbf{X}_{\textsc{i}}^{\top}\mathbf{X}_{\textsc{i}})^{-1}\sigma_{Y|\text{do}(X)}^{2}$
is bounded in probability, for large enough $m$. Accordingly, this implies
that
$\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})\xrightarrow{\text{P}}\mathbf{0}$.
Thus, when keeping $n$ fixed, we obtain
$\mathbf{W}^{m}_{*}\xrightarrow{\text{P}}\mathbf{I}_{p}$, for
$m\rightarrow\infty$.

On the other hand, if we keep $m$ fixed and consider the limit
$n\rightarrow\infty$ instead, we observe that

$\mathbf{W}^{m}_{*}\xrightarrow{\text{P}}\bm{\Delta}\bm{\Delta}^{\top}(\mathbf{Cov}(\widehat{\bm{\alpha}}^{m}_{\textsc{i}})+\bm{\Delta}\bm{\Delta}^{\top})^{-1}.$

We note that we do not have
$\mathbf{W}^{m}_{*}\xrightarrow{\text{P}}\mathbf{0}$ here in general, because
the bias in $\widehat{\bm{\alpha}}^{n}_{\textsc{o}}$ remains, independent of
the sample size $n$.