Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
196,710 characters · 57 sections · 55 citation commands
Cross-Validated Causal Inference: a Modern Method to Combine Experimental and Observational Data
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 {
} \fi \if01 {
} \fi
{\it Keywords:} Causal inference; treatment effect; data integration; unobserved confounders, randomized experiments, observational data.
\spacingset{1.9}
We focus on the problem of estimating the average treatment effect (ATE) in causal inference. Over the past decades, a wide range of statistical methods have been developed to draw causal conclusions from either experimental or observational studies. Experimental data, collected from randomized controlled trials (RCTs), offer high internal validity. However, such data can be costly to obtain. In contrast, observational data are often cheaper, but their internal validity is suspect. Specifically, ATE estimates based on observational data, assuming unconfoundedness, may suffer from biases due to unobserved confounders.
In this paper we consider the combination of experimental and observational data, with the goal of producing robust (to the presence of unobserved confounders) and precise (by including observational data) causal conclusions. We propose a framework that minimizes a weighted combination of losses: the experimental loss, which assesses the causal parameter’s validity; the observational loss, which measures the full model’s fit; and their relative weighting, chosen adaptively via cross-validating the causal parameter.
To illustrate the basic ideas, consider a setting with no covariates. We have an experimental sample where we observe both treated and control units, and an observational sample where we observe only control units, based on the widely used LaLonde data lalonde1986evaluating, dehejia1999causal. Because in this setting there is no question about estimating the average outcome for the treated, for which we only have the experimental data, the question is how to estimate the average control outcome for the experimental population, $\mathbb{E}[Y_i^{\mathrm{exp}}(C)]$. The average control outcome in the experimental sample, $\overline{Y}_C^{\mathrm{exp}}$, is unbiased for this expectation (but possibly imprecise due to limited data size). The average of the control outcome in the observational sample, $\overline{Y}_C^{\mathrm{obs}}$, may be biased for the experimental population’s average control outcome. We consider a weighted average of the average control outcome in the observational sample and the average of the control outcome in the experimental sample, with weights $\lambda\in[0,1]$ and $1-\lambda$ respectively:
What properties would we like $\lambda$ to have? If the experimental sample is large, then even if the bias in the observational sample is very small, as long as there is some bias we would like $\lambda$ to be close to zero. If on the other hand the bias in the observational sample is negligible, then we would like to choose $\lambda$ close to one. In other words, we would like to shrink our experimental estimate towards the observational data, but do so in a data-adaptive fashion, that is, with a data-driven $\lambda$. In this simple no-covariate case where the focus is on the expected control outcome in the experimental population, we implement this objective by selecting $\lambda$ through cross-validating on the experimental data: \[ \widehat\lambda= \arg\min_{\lambda\in[0,1]} \underbrace{ \frac{1}{K} \sum_{k = 1}^{K} \Big( \overline{Y}_{C,B_k}^{\mathrm{exp}} - \Big( (1 - \lambda) \overline{Y}^{\mathrm{exp}}_{C,-B_k} + \lambda \overline{Y}_C^{\mathrm{obs}}\Big) \Big)^2}_{\mathsf{CV}(\lambda), \text{ the cross-validation objective}}, \] where the subscripts $\{B_k, -B_{k}\}_{ k\in[K]}$ denote the complementary subsets in $K$-fold cross-validation. In the paper we extend this to the case with more general models for the observational data involving covariates.
In Figure (ref), we present some results for this example based on the LaLonde data lalonde1986evaluating, dehejia1999causal. In the bottom two panels we present two sets of three estimates of the ATE. First, in both panels, results based on the experimental data alone (corresponding to $\lambda=0$). Second, again in both panels, results based on the observational data alone (corresponding to $\lambda=1$). Both are intended to set the stage for our preferred results based on the cross-validated $\hat\lambda$. The cross-validation is based on five fold splits, leading to a unique $\hat\lambda$. We repeat this many times to get a distribution of selected $\hat\lambda.$ In the case without covariates, we find that the selected $\hat\lambda$ is always close to or exactly equal to 0, corresponding to the experimental estimates. The cross-validation makes clear that the the data can tell us that the observational data are of little value in this case. For a covariate-adjusted version of the observational data estimator, the cross-validated $\hat\lambda$ is much closer to 1, with the average value for $\hat\lambda$ over many choices of five folds equal to 0.77. Here the data imply that the observational data are valuable. The combination of the two sets of results shows that in this case our proposed method can detect when the observational sample is valuable, and when it is not, in a fully data-driven way.
Our contributions are three-fold: First, we introduce a novel method to systematically combine experimental and observational data. The methodological advantages include: $(i)$ we do not require additional model specifications or identification assumptions; $(ii)$ the method allows for the setting where observational and experimental samples have different covariates; $(iii)$ the method allows the treatment to have a different effect in two populations. Second, we conduct experiments on both synthetic and real data to show the effectiveness of our method. For synthetic data, we address the most common cases (no-covariate and covariate-adjusted linear setting). For real data, we use the LaLonde-Dehejia-Wahba dataset lalonde1986evaluating, dehejia1999causal. Third, we develop supporting non-asymptotic theories for the robustness of our method. Under regularity conditions, we show that our method achieves an $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$ error rate regardless of the level of bias in the observational data, where ${N^{{\mathrm{exp}}}}$ is the experimental sample size (Corollary (ref)). This is known to be optimal for estimators that are based solely on experimental data. Moreover, in the no-covariate setting, we show that the $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$-rate is minimax optimal (when the observational data are unbiased) over a class of robust estimators that combine experimental and observational data (Theorem (ref)).
We choose the weight for the experimental and observational estimates, denoted by $\lambda$, through cross-validation of the causal parameter. This is both inspired by a broader cross-validation-based statistical learning family that includes stacking wolpert1992stacked, breiman1996stacked, aggregation tsybakov2003optimal, tsybakov2004optimal, and super learner van2007super. We adapt these tools for causal inference by addressing issues such as identification, confounding, and distributional shifts. We design our cross-validation criterion to be explicitly tailored to causal estimands, rather than predictive objectives. Specifically, in each split, we fit on $K{-}1$ experimental folds and all observational data, and evaluate the causal parameter on the held-out experimental fold. The $\lambda$ that optimizes for the average experimental loss is then used to refit on all data. Intuitively, when the observational sample exhibits low bias, our method assigns more weight to the observational loss, exploiting the additional sample size.
A systematic and unified framework to combine experimental and observational data remains largely absent---existing literature is often {\it ad hoc} in nature and hinges on auxiliary assumptions, such as extrapolatable bias kallus2018removing, additional model specifications yang2020combining, prespecified study structures rosenman2023combining, or covariate similarity gui2024combining.
The state of the existing literature is summarized at a high level in Table (ref). As the table indicates, the three major lines of work—pooling, shrinkage, and error-prone estimators—each have their limitations. Pooling methods treat all data as coming from a single source, breaking experimental randomization and requiring unconfoundedness assumptions to incorporate observational data ross2009pooled, gao2023pretest, yang2023elastic, xiong2023federated. Our method could be seen as a “soft” version of pooling, dynamically adjusting the weighting of each data source rather than making an all-or-nothing decision. Shrinkage methods are most similar in spirit to our proposed method. They tolerate bias from observational data but depend on predefined strata, often assuming that the average effects are equal across data sources within each stratum---a condition that may not hold in practice stein1956inadmissibility, green1991james, green2005improved, rosenman2023combining. Although both these methods and ours involve weighting, they adjust stratum-level estimators, while we bypass the need for discrete stratification and instead optimize weights at the loss level. Error-prone estimators carefully balance biased components for each source to cancel out confounding effects yang2020combining. This is an appealing idea, though in practice it can be challenging to construct such estimators. Both our method and error-prone approaches exploit the consistency of experimental estimates, but the mechanisms differ fundamentally. Instead of relying on delicate bias-cancellation conditions, we directly cross-validate on experimental data to prevent incorporating observational bias. In Section (ref), we provide an extended discussion of related methods, with a focus on unmeasured confounding in observational data and a broader discussion on cross-validation in machine learning.
Suppose we have access to two datasets $X^{\mathrm{exp}}$ and $X^{\mathrm{obs}}$. The former is comprised of ${N^{{\mathrm{exp}}}}$ experimental samples, $X^{{\mathrm{exp}}}_i=(Y^{{\mathrm{exp}}}_i,W^{{\mathrm{exp}}}_i,Z^{{\mathrm{exp}}}_i)\in\mathcal{X}^{{\mathrm{exp}}}$, where $Y^{{\mathrm{exp}}}_i\in{\mathbb R}$, $W^{{\mathrm{exp}}}_i\in\{0,1\}$, $Z^{{\mathrm{exp}}}_i\in{\mathbb R}^{d_{{\mathrm{exp}}}}$ are the observed outcome, binary treatment ($0$ for control, $1$ for treated), and covariate/pre-treatment vector, respectively. The latter consists of ${N^{{\mathrm{obs}}}}$ observational samples, $X^{{\mathrm{obs}}}_i=(Y^{{\mathrm{obs}}}_i,W^{{\mathrm{obs}}}_i,Z^{{\mathrm{obs}}}_i)\in\mathcal{X}^{{\mathrm{obs}}}$, defined analogously.
We adopt the classical potential outcome framework rubin1974estimating, rubin1977assignment, rubin1978bayesian, where the potential outcomes are denoted by $(Y^{\sf s}_i(1),Y^{\sf s}_i(0))$, and $Y^{\sf s}_i = Y^{\sf s}_i(W^{\sf s}_i)$ for $\sf s\in\{{\mathrm{exp}},{\mathrm{obs}}\}$. For the experimental data, we make standard assumptions: (1) $(Y^{{\mathrm{exp}}}_i(1),Y^{{\mathrm{exp}}}_i(0),W^{{\mathrm{exp}}}_i,Z^{{\mathrm{exp}}}_i)\simiidP^{{\mathrm{exp}}}$ for some distribution $P^{{\mathrm{exp}}}$; (2) there is no unobserved confounder, {\it i.e.}, $(Y^{{\mathrm{exp}}}_i(1),Y^{{\mathrm{exp}}}_i(0))\mathbin{\perp\!\!\!\perp} W^{{\mathrm{exp}}}_i|Z^{{\mathrm{exp}}}_i$; (3) the overlap condition is satisfied, {\it i.e.}, the propensity score $\ensuremath{{\mathbb{P}}}(W^{{\mathrm{exp}}}_i=1|Z^{{\mathrm{exp}}}_i)$ lies in the open interval $(0,1)$.
For the observational data, we impose no distributional assumptions. In particular, we do not assume that the two data sources share the same covariate distributions, thus allowing for covariate shift; we also do not require their outcome models to be the same, permitting label shift and differing response mechanisms. Additionally, we allow the observational data to be non–independent and non-identically distributed (non-i.i.d.), and we allow both unmeasured confounders and outcome model misspecification—conditions under which standard doubly robust estimators will fail to provide valid inference.
We want to estimate the ATE on the population associated with the experimental data:
where the expectation is over the distribution in the experimental population $P^{{\mathrm{exp}}}$. This estimand can be easily extended to targeting other populations ({\it e.g.,} observational or mixed) by modifying the cross-validation objective in Section (ref).
Let $\theta$ denote the parameter of the full model, with $\beta \coloneqq \beta(\theta)$ being its causal estimand, which can be characterized in terms of this full parameter. For example, when there are no covariates, $\beta=\theta$; in a linear model, $\beta$ is the coefficient for the treatment. More generally, $\beta$ indexes the counterfactuals implied by $\theta$ we estimate from the data.
We start with the standard no-covariate setting where only response and treatment are observed in both sources. For a random experimental sample $X=(Y, W) \sim P^{\mathrm{exp}}$, we are interested in the ATE \[\tau^{\star}= {\mathbb E} (Y \mid W = 1) - {\mathbb E} (Y \mid W = 0 ), \] which is estimated by the difference in means: \[ \widehat\tau^{\mathrm{exp}} = \frac{1}{\sum_i \mathbbm{1}\{W^{{\mathrm{exp}}}_i=1\}} \sum_{W^{{\mathrm{exp}}}_i=1} Y^{{\mathrm{exp}}}_i - \frac{1}{\sum_i \mathbbm{1}\{W^{{\mathrm{exp}}}_i=0\}} \sum_{W^{{\mathrm{exp}}}_i=0} Y^{{\mathrm{exp}}}_i . \]
Consider, for example, the LaLonde dataset. In the experimental data, the treatment group has an average outcome of $\$ 6.3$k, and the control group has an average outcome of $\$4.6$k, yielding an ATE estimate of $\widehat\tau^{\mathrm{exp}}=\$ 1.8$k. The observational data share the same treatment group, but the control group's average outcome is $\$ 21.6$k, yielding an estimate of $\widehat\tau^{\mathrm{obs}}=\$-15.2$k. Notably, the observational control group is much larger (2,490 vs.\ 260 samples), offering potential efficiency gains despite its bias. How can we systematically combine them to improve estimation? When there are no covariates, our method utilizes a weighted average of the means of the two control group, where the relative weighting $\lambda$ is selected through cross-validation.
In the LaLonde example, the treatment mean is the same across data sources and we focus on estimating the control mean. In the general case where we need to estimate the treatment mean (or the control mean, analogously), our estimate is $\beta(\widehat\theta(\widehat\lambda)) = \widehat\theta(\widehat\lambda) \in {\mathbb R}$, where its closed-form expression given $\lambda \in [0, 1]$ is as follows:
where the overline denotes sample mean. Intuitively, it shrinks the experimental estimate towards the observational one. See Section (ref) for its derivation and additional discussions.
We select $\lambda$ by cross-validating on the experimental data. The overall mean estimator is \[ \widehat\theta(\widehat\lambda), \quad \widehat\lambda= \arg\min_{\lambda\in[0,1]} \underbrace{ \frac{1}{K} \sum_{k = 1}^{K} \Big( \overline{Y}_{B_k}^{\mathrm{exp}} - \Big( (1 - \lambda) \overline{Y}^{\mathrm{exp}}_{-B_k} + \lambda \overline{Y}^{\mathrm{obs}}\Big) \Big)^2}_{\mathsf{CV}(\lambda), \text{ the cross-validation objective}}, \] where the subscripts $\{B_k, -B_{k}\}_{ k\in[K]}$ denote the complementary subsets in $K$-fold cross-validation.
We now consider the case where the full model is linear, ${\mathbb E}(Y_i^{\mathrm{obs}}|W^{\mathrm{obs}}_i,Z^{\mathrm{obs}}_i)=\theta_W W_i+\theta_Z ^\top Z_i$ with $\theta^\top=(\theta_W\ \theta_Z^\top)$. This setting does not require the experimental data source to include covariates, and if covariates are present in the experimental data, they may differ entirely from those present in the observational data. We define each component for the overall estimator $\widehat\theta(\widehat\lambda)$ as follows: first, $\theta$ represents the parameter vector of a linear outcome model fit on observational data. The first entry of $\theta$ corresponds to the treatment effect $\beta$. The observational loss is \[ L^{\mathrm{obs}}(\theta; X^{\mathrm{obs}}) \coloneqq \frac{1}{N^{\mathrm{obs}}}\sum_{i=1}^{N^{\mathrm{obs}}} \Big(Y_i^{\mathrm{obs}} - \Big( W_i^{\mathrm{obs}}\quad {Z_i^{\mathrm{obs}}}^\top \Big) \theta \Big)^2 . \] Second, for the causal parameter $\beta$, we define the experimental loss $L^{\mathrm{exp}}$: \[ L^{\mathrm{exp}}(\beta; X^{\mathrm{exp}}_\mathcal{J}) \coloneqq \Big( \beta - \widehat\tau^{\mathrm{exp}} \Big)^2, \] where $\widehat\tau^{\mathrm{exp}}$ is obtained from a subset of experimental data $X^{{\mathrm{exp}}}_{\mathcal{J}}$ indexed by $\mathcal{J}$. This could be the simple difference in means based on the experimental data, $\overline{Y}_T^{\mathrm{exp}}-\overline{Y}^{\mathrm{exp}}_C$, or a more complex estimator that involves some covariate adjustment. Here, we use the standard $\ell_2$ loss as it is strongly convex in $\beta$ (which facilitates the theoretical analysis) and admits a desirable additive structure (formalized as Lemma (ref)) when the experimental estimate $\widehat\tau^{\mathrm{exp}}$ can be expressed as an average over individual units. This structure applies to common estimators including the difference-in-means, plug-in, and the AIPW estimators, {\it i.e.,}
where, for example,
with an outcome model $\widehat\mu:\{0, 1\} \times {\mathbb R}^{d_{{\mathrm{exp}}}} \rightarrow {\mathbb R} $ and propensity score $\widehat\pi: {\mathbb R}^{d_{{\mathrm{exp}}}} \rightarrow (0, 1) $.
For example, if we use the plug-in estimator for $\widehat\tau^{\mathrm{exp}}$ with a linear experimental outcome model, then for a vector $\theta^{\mathrm{exp}}$ with its first entry as the treatment coefficient,
Knowing how to evaluate the causal parameter, we now provide a closed-form solution for the full model $\widehat\theta(\lambda)$. We denote $W^{\mathrm{obs}}, Z^{\mathrm{obs}}, Y^{\mathrm{obs}}$ as the respective matrices containing all observational samples, where each column corresponds to one sample. We append a 1 to each $Z_i^{\mathrm{obs}}$ to include an intercept term in the linear model. For $\lambda\in[0, 1]$, the full model \[ \widehat\theta(\lambda) = \arg\min_{\theta} (1-\lambda) \underbrace{(\theta^\top e_1- \overbrace{\widehat\tau^{\mathrm{exp}}}^{\text{from Eq.}\eqref{eq:m-estimate_all}})^2}_{\text{experimental loss}} + \lambda \underbrace{\frac{1}{N^{\mathrm{obs}}}\sum_{i=1}^{N^{\mathrm{obs}}} \Big(Y_i^{\mathrm{obs}} - \Big( W_i^{\mathrm{obs}} \quad {Z_i^{\mathrm{obs}}}^\top \Big) \theta \Big)^2}_{\text{observational loss}} \] is given by the solution to
where $e_1 = (1 \quad 0 \cdots 0 )^\top . $ Intuitively, the first term on both sides regularizes the treatment coefficient toward the experimental estimate, while the second term on both sides fits the full model to the observational data. The derivation is provided in Section (ref). Similar to the no-covariate setting, we select a $\widehat\lambda$ by cross-validation to provide the final estimate $\beta(\widehat\theta(\widehat \lambda))$.
To estimate the values of $\theta$ in the general parametric setting, we formulate the problem as an empirical risk minimization (ERM) problem. Suppose the loss on the experimental and observational data is denoted by $L^{\mathrm{exp}} (\beta; X^{\mathrm{exp}})$ and $L^{\mathrm{obs}}(\theta; X^{\mathrm{obs}})$, respectively. The experimental loss quantifies validity of the causal parameter on experimental data. Since the experimental data are assumed to be unconfounded, this loss serves as a benchmark for consistent causal estimation. The observational loss evaluates how well the full model (including its causal parameter) explains the observational dataset under its data-generating process. Intuitively, when the experimental sample size goes to infinity, we would expect to converge to the true ATE $\tau^\star$ by minimizing the loss: \[ \tau^\star = \lim_{|X^{\mathrm{exp}}| \to \infty} \arg\min_{\beta} L^{\mathrm{exp}}(\beta; X^{\mathrm{exp}}). \] Meanwhile, the observational data could give a biased estimate even in the limit: \[ \tau^\star + \varepsilon = \lim_{|X^{\mathrm{obs}}| \to \infty} \beta( \arg\min_{\theta} L^{\mathrm{obs}}(\theta; X^{\mathrm{obs}})), \] where $\varepsilon$ is unobserved and unestimable. We do not impose structural or source-specific assumptions on $\varepsilon$, allowing it to capture diverse real-world scenarios. For instance, $\varepsilon$ can be interpreted as the effect of an unobserved binary confounder that aligns with the treatment assignment, or more generally, as the combined effect of multiple unobserved confounders. It could also arise from both unmeasured confounders and model misspecification (in the case of AIPW), or treatments having a different effect on the observational population.
We now present our method for the general case. Our overall estimate is \[ \beta(\widehat{\theta} (\widehat{\lambda}; X^{\mathrm{exp}}, X^{\mathrm{obs}})), \quad \widehat{\lambda} = {\arg\min}_{\lambda\in[0,1]} \mathsf{CV}(\lambda; X^{\mathrm{exp}}, X^{\mathrm{obs}}), \] where each component is defined as follows:
\paragraph{Learning $\widehat{\theta}(\lambda)$.} Given $\lambda$, the full model fitted on $X^{\mathrm{exp}}, X^{\mathrm{obs}}$ is obtained by \[ \widehat{\theta} (\lambda) \coloneqq \widehat{\theta} (\lambda; X^{\mathrm{exp}}, X^{\mathrm{obs}}) = \arg \min_{\theta} \Big\{ (1 - \lambda) \underbrace{L^{\mathrm{exp}} (\beta (\theta); X^{\mathrm{exp}})}_{\text{causal parameter}} + \lambda \underbrace{L^{\mathrm{obs}} (\theta; X^{\mathrm{obs}})}_{\text{full model}} \Big\}. \] We have provided closed-form solutions for the most common cases (no-covariate and linear setting). For other cases, we can employ gradient-based, (quasi-)Newton, or other optimization techniques suited to the structure of the objective function. \paragraph{Selecting $\widehat{\lambda}$ by cross-validating the causal parameter.} We use $\{X^{\mathrm{exp}}_{B_k}, X^{\mathrm{exp}}_{-B_k}\}_{ k\in[K]}$ to denote complementary subsets in the $K$-fold splitting in cross-validation. Denote $D\coloneqq(X^{\mathrm{exp}}, X^{\mathrm{obs}})$, $D_{k}\coloneqq(X^{\mathrm{exp}}_{B_k} , X^{\mathrm{obs}})$, and $D_{-k}\coloneqq (X^{\mathrm{exp}}_{-{B_k}} , X^{\mathrm{obs}})$, as we only split experimental data and always reuse observational data. For each fold $k$, fit a model on $D_{-k}$: \[ \widehat{\theta} (\lambda; D_{-k}) = \arg \min_{\theta} \Big\{ (1 - \lambda) L^{\mathrm{exp}} (\beta (\theta); X_{-{B_k}}^{\mathrm{exp}}) + \lambda L^{\mathrm{obs}} (\theta; X^{\mathrm{obs}} ) \Big\}. \] Then evaluate the causal parameter on $D_{-k}$ for the cross-validation objective $\mathsf{CV}$:
$\mathsf{CV}$ quantifies how well the estimated treatment effect aligns with experimental evidence. See Section (ref) for pseudo-code and analysis of the computational complexity of our procedure.
To summarize the motivation, the loss-based objective explicitly encodes the trade-off between bias and variance in a unified optimization framework. Specifically, the observational data are leveraged as a source of potential efficiency gains to aid fitting the full model that contains the causal parameter. We employ cross-validation to safeguard for causal validity. When the causal estimate from the combined data are well aligned the experimental evidence, cross-validation favors models that leverage this alignment. Otherwise, it reverts toward the experimental data to control for potential bias.
We now discuss the choices of $L^{\mathrm{exp}}$ and $L^{\mathrm{obs}}$. For simplicity of presentation, we demonstrate using the squared error loss for both $L^{\mathrm{exp}}$ and $L^{\mathrm{obs}}$. Since $L^{\mathrm{exp}}$ evaluates a scalar $\beta$, the squared error is a natural choice. For $L^{\mathrm{obs}}$, we assume strong convexity and smoothness conditions ({\it i.e.}, three times differentiable with bounded second and third derivatives), as formalized later in Assumption (ref). This class includes squared loss, L2 regularization (Ridge loss), and L$_p$ loss ({\it i.e.}, $|y-y’|^p/p, p\ge3$). On the other hand, this class excludes L1 regularization (LASSO, due to the non-differentiability at zero), elastic net (a combination of L1 and L2 regularization), and Huber loss (because it is not twice differentiable at the threshold). These requirements are imposed to facilitate the theoretical analysis in Section (ref). Violations in practice are unlikely to result in catastrophic failure.
In this section, we present empirical evidence on the following questions: How does the bias $\varepsilon$ affect the performance of our method? How does $N^{\mathrm{obs}}$ affect the estimation error? Can our cross-validation procedure reliably select a “good” value of $\widehat\lambda$?
Without loss of generality, we estimate the treatment mean and take our samples to be $Y^{\mathrm{exp}}_1, \ldots, Y^{\mathrm{exp}}_{N^{\mathrm{exp}}} \overset{iid}{\sim} {\mathcal N}(\tau^\star, {\sigma}^2)$ and $Y^{\mathrm{obs}}_1, \ldots, Y^{\mathrm{obs}}_{N^{\mathrm{obs}}} \overset{iid}{\sim} {\mathcal N}(\tau^\star + \varepsilon, {\sigma}^2)$. We compare the proposed method with the empirical risk minimizer ({\it i.e.}, sample means) on either data source, and an additional baseline to determine the value of $\lambda$ via a t-test. We use (empirical) Mean Squared Error (MSE) for assessment. For implementation details see Section (ref).
\thispagestyle{empty}
The results demonstrates a clear advantage for our method. As shown in Figures (ref) and (ref), it reliably adapts to varying, unknown values of $\varepsilon$ and outperforms at least one of the single-source methods. When $\varepsilon$ is small, it improves upon the $X^{\mathrm{exp}}$-only approach; for intermediate $\varepsilon$, it yields the lowest error among all baselines; for large $\varepsilon$, it outperforms the $X^{\mathrm{obs}}$-only and t-test approaches while remaining comparable to using $X^{\mathrm{exp}}$ alone. As shown in Figures (ref) and (ref), our estimator increasingly resembles the experimental estimate as $\varepsilon$ grows, with only minor fluctuations observed before $\widehat{\lambda}$ approaches zero. This adaptivity underscores a key strength of our method: its ability to dynamically adjust the reliance on two data sources via cross-validation, which is implicitly governed by the finite-sample error and observational bias. When observational data are scarce or less reliable, cross-validation leans more heavily on experimental data. This flexibility enables robust performance across diverse data regimes without requiring prior knowledge of $\varepsilon$, making the method suitable for practical applications where the experimental-observational trade-off is unknown or context-dependent.
We note that our method’s performance improves as the number of observational samples increases (Figure (ref)). For a fixed $N^{\mathrm{obs}}$, it consistently outperforms the $X^{\mathrm{exp}}$-only baseline (Figure (ref)), demonstrating the benefit brought by incorporating observational data.
The above observations hold in both low ($\sigma^2 = 1$) and high ($\sigma^2 = 100$, in Figures (ref)) noise settings. Additional results on the impact of noise level are provided in Section (ref).
Assume each data sample consists of a tuple of response $Y$, covariates $Z$, and binary treatment $W$. For experimental data, generate the response as a linear combination of the covariates plus an exogenous noise: $Y = Z^\top \theta^{\mathrm{exp}} + W \tau^\star + \xi$, where $\xi \sim {\mathcal N} (0, \sigma^2) \perp Z, W . $ For observational data, we incorporate a bias $\varepsilon$ to capture unmeasured confounders associated with the treatment: $Y = Z^\top \theta^{\mathrm{obs}} + W (\tau^\star + {\varepsilon}) + \xi$, where $\xi \sim {\mathcal N} (0, \sigma^2) \perp Z, W.$ Here, $\theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}}$ denote the parameters of the respective linear outcome models. The two parameter vectors can differ entirely in both values and dimensions. We consider two scenarios, $\theta^{\mathrm{obs}} = \theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}} \neq \theta^{\mathrm{exp}}$. For implementation details see Section (ref).
We observe trends similar to those in the no-covariate case: Figure (ref) shows that our method consistently outperforms at least one of the baselines relying on one data source alone. This advantage holds regardless of whether the two data sources share the same covariates (Figure (ref)) or not (Figure (ref)), and whether the experimental dataset is small ($N^{\mathrm{exp}} = 50$ in Figures (ref), (ref), (ref), and (ref)) or large ($N^{\mathrm{exp}} = 1000$ in Figures (ref) and (ref)). When the bias is moderately low, our method achieves the most accurate causal estimates. Such low-bias regime corresponds roughly to $\varepsilon \leq 0.5$ when $N^{\mathrm{exp}}=50$ (Figures (ref) and (ref)) and narrows to $\varepsilon \leq 0.1$ when $N^{\mathrm{exp}}=1000$ (Figure (ref)). Incorporating more observational samples, even when they contain minor bias, can enhance estimation accuracy (Figures (ref) and (ref)).
\thispagestyle{empty}
\thispagestyle{empty}
We now turn to Figure (ref) to examine the behavior when $\varepsilon$ is large. When the experimental data are abundant ($N^{\mathrm{exp}}=1000$ in Figures (ref) and (ref)), our method reliably converges to the experimental estimates, as cross-validation consistently selects $\widehat\lambda = 0$. This outcome aligns with our expectations: our method appropriately downweights the observational component when the experimental estimates are sufficiently more reliable. In contrast, when the experimental data are limited ($N^{\mathrm{exp}} = 50$ in Figures (ref) and (ref)), the selected $\widehat\lambda$ generally remains close to zero for large $\varepsilon$, but occasionally small non-zero values are selected (Figure (ref)). When this happens, the resulting squared error can be large because a large $\varepsilon$ would amplify the error for very small $\widehat\lambda$ (Figure (ref)). Consequently, the overall MSE suffers from these rare but high-error instances. To address these results, we make the following comments: first, a bias beyond the order of several thousand percent is highly unlikely in real-world settings when data are collected by trained professionals. Nonetheless, under extreme bias, for example, $\varepsilon = 200$, Figure (ref) shows that among 5000 simulation runs, the selection of $\widehat\lambda > 0.1$ occurred only a handful of times. We argue that this rarity can be interpreted as a form of high-probability safeguard: while the method is not immune to error under severe confounding for small $N^{\mathrm{exp}}$, it exhibits robust behavior in the majority of cases. We supplement Figures (ref) and (ref) for additional evidence.
Interestingly, the results remain unaffected by whether $\theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}}$ are equal or different (comparing Figures (ref) to (ref), and (ref)–(ref) to (ref)–(ref)). This reflects our framework’s design: only estimated treatment effects, not raw covariates or outcome models, are shared across data sources. As a result, it naturally accommodates entirely distinct outcome models, including differences in functional forms, learned weights, sets of covariates, and their underlying distributions. Such flexibility is typically not supported by existing methods, which often require stronger assumptions about model alignment or covariate overlap across data sources.
In 1986, Robert LaLonde published a seminal paper that compared the results of a field experiment with the range of estimates an econometrician might have produced using only nonexperimental data, concluding that the nonexperimental methods at that time failed to systematically replicate the trial results lalonde1986evaluating. The original study examined the effect on trainee earnings of an employment program implemented through a field experiment, wherein participants were randomly allocated to either treatment or control groups. Discussion and analysis of the LaLonde dataset has led to significant methodological advances in causal inference heckman1987we, smith2005does, athey2017state, imbens2024lalonde. In our paper, we compare the ATE estimates on the LaLonde dataset from our method with various baselines, using the widely adopted data-selection process outlined by Dehejia and Wahba in dehejia1999causal.
The National Supported Work Demonstration (NSW) was the randomized trial where the treatment is to receive a job training between 1975 and 1977. LaLonde and later Dehejia and Wahba analyzed its impact on real earnings (RE) in 1978, with the latter restricts on a smaller subgroup. The resulting NSW dataset contains 185 treated and 260 control samples. We detail the data selection process in Section (ref)
The observational control data comes from the Panel Study of Income Dynamics (PSID) and Westat's Matched Current Population Survey (CPS). They are control-only datasets. We term them observational control group. They are partitioned by pre-intervention variables into subgroups PSID-2, PSID-3, CPS-2, and CPS-3, with the full datasets denoted PSID-1 and CPS-1. We detail the partition procedure in Section (ref).
To estimate ATE, we apply various linear models on full NSW data to produce the “gold standard" experimental estimates. Same linear models are applied on the NSW treatment group, but with different observational control groups instead.
We will show the following four sets of methods in each panel: First, the approach using experimental data alone (corresponds to $\lambda=0$). Second, our proposed method, which selects $\lambda$ via cross-validation. Third, the approach using observational data alone (corresponds to $\lambda=1$), which uses NSW treatment group and observational controls dehejia1999causal. Lastly, pooling all data together ross2009pooled. This can be interpreted as treating the NSW treatment group, NSW control group, and observational control groups collectively as observational data, and setting $\lambda=1$ .
Under these settings, we produce Tables (ref), (ref), and (ref): the first table highlights selected configurations in the main text, while the latter two with full configurations are deferred to Section (ref). Specifically, Table (ref) focuses on two major observational control group (PSID-1 and CPS-1) and three covariate settings (matching columns 1, 3, and 8 of Tables (ref) and (ref)). It integrates the point estimates from Table (ref) and bootstrapped standard deviations from Table (ref). We detail their setup in Section (ref).
As a starting point, it is encouraging to see that results from nearly 26 years ago can still be largely replicated precisely today. Our reproduced point estimates using single data source (first and third panels in Tables (ref) and (ref)) match exactly with those of Dehejia and Wahba's (columns 1-4 of panels B and C in Table 2, which originally correspond to LaLonde's Table 5 without data selection). We note, however, that column 5 of their panels B and C were described as controlling for all pre-intervention variables, but simply including all such variables did not allow us to exactly replicate their results. We provide additional discussions of these reproduced results in Section (ref).
Importantly, our method gives more accurate and reliable ATE estimates compared to other methods. By “more accurate” we mean estimates closer to those obtained from experimental data only (the first panel). We note that this may not be ground-truth effect because it still contains finite-sample error. As shown in Table (ref), our method (the second panel) consistently outperforms the approach relies solely on observational data (the third panel). In fact, this holds in all 48 configurations in Table (ref). In such a case, we confirm what LaLonde found: there is inherent difficulty in accurate nonexperimental modeling due to extreme inter-model variability, even after choosing more suitable subsets. Comparing to pooling (the fourth panel), our method is more accurate in majority of cases (30 out of 48 configurations). While pooling occasionally performs better (CPS-1, column 8, last row in Table (ref)), such gains are usually marginal. In contrast, when pooling fails, it could produce drastically biased estimates ({\it e.g.}, $-13,598$ with a p-value $< 0.0001$).
Moreover, we identify two trends that align with our intuition: First, inclusion of additional informative covariates leads to a greater weight on the observational component. In Table (ref), $\widehat\lambda$ values are generally small in column 1, where only the treatment is used, and increase in columns 3 and 8, where additional covariates are included. Such trend is also presented in Table (ref)--columns 6 and 7, which corresponds to columns 2 and 4 with the addition of RE74, exhibit noticeably larger $\widehat\lambda$ values. Second, for observational control subgroups that are more similar to the NSW control group, the selected $\widehat\lambda$ are generally larger. This agrees with LaLonde’s assertion that subgroups such as PSID-2, PSID-3, CPS-2, and CPS-3 are more comparable to the NSW control group in distributions of pre-intervention variables.
Finally, the bootstrapped standard deviations in Table (ref) show that our estimates have variability comparable to other methods on the LaLonde dataset. Such uncertainty reflects both data re-sampling and cross-validation splitting, though the latter contributes little (Table (ref)). The same pattern holds in full configurations (Table (ref)), indicating that our method matches the stability of existing approaches while offering greater flexibility.
When using the LaLonde dataset, experimental estimates are treated as the ground-truth effect. How to determine whether our method offers gains compared to using experimental data alone? We conduct experiments on synthetic data derived from the LaLonde dataset. To generate synthetic $X^{\mathrm{exp}}$ and $X^{\mathrm{obs}}$, we fit linear models on respective real data sets and re-sample the residuals from Gaussian distributions under sample mean and variance. This ensures the experimental estimate to be unbiased for the ground-truth effect in expectation.
The presented results are generally a callback to the analysis in Section (ref). Our method achieves the lowest RMSE on CPS, column 8 (Table (ref)), corresponding to the regions where it has the lowest error in Figure (ref). In other cases, we refer back to discussions of Figure (ref) where our method underperforms the experimental estimate, given the substantial bias present in the LaLonde observational data and the small experimental sample size.
Recall the setup in Section (ref), where we are given ${N^{{\mathrm{exp}}}}$ i.i.d. experimental samples, $X^{{\mathrm{exp}}}_i=(Y^{{\mathrm{exp}}}_i,W^{{\mathrm{exp}}}_i,Z^{{\mathrm{exp}}}_i)\in\mathcal{X}^{{\mathrm{exp}}}, i\in[{N^{{\mathrm{exp}}}}]$, and ${N^{{\mathrm{obs}}}}$ observational samples, $X^{{\mathrm{obs}}}_i=(Y^{{\mathrm{obs}}}_i,W^{{\mathrm{obs}}}_i,Z^{{\mathrm{obs}}}_i)\\\in\mathcal{X}^{{\mathrm{obs}}},i\in[{N^{{\mathrm{obs}}}}]$. We consider the quadratic experimental loss, {\it i.e.}, ${L^{\mathrm{exp}}}(\beta;X^{{\mathrm{exp}}}_{\mathcal{J}}) =(\beta-\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}}))^2$ for any set of experimental samples $X^{{\mathrm{exp}}}_{\mathcal{J}}$ indexed by $\mathcal{J}\subseteq[{N^{{\mathrm{exp}}}}]$. In addition, we make the following two assumptions for our analysis:
Assumption (ref) assumes that the ATE estimator $\widehat\tau^{\mathrm{exp}}$ based on experimental samples is $\sqrt{{N^{{\mathrm{exp}}}}}$-consistent and admits a linear approximation. For example, this is satisfied in our linear setting. A sufficient condition for Assumption (ref) is that the $\widehat\tau^{\mathrm{exp}}$ is derived from some $Z$-estimation problem (Assumption (ref)). We refer to Section (ref) for more details. In Assumption (ref), we require the observational loss ${L^{\mathrm{obs}}}$ to be strongly convex and have smooth higher-order derivatives. These are standard regularity conditions for analyzing empirical risk minimization. Moreover, the assumption $\beta(\theta)= \theta_1$ can be generalized to $\beta$ being a linear function of $\theta$, as they are equivalent up to a linear transformation on $\theta$. We choose the former for simplicity of presentation.
Throughout our presentation, we use ${C},{C'}>0$ to denote constants that depend polynomially on the parameters in the assumptions (See Section (ref) for details). Our main result holds for any experimental sample size exceeding a threshold determined by a user-specified parameter $\delta \in (0,1/2)$, which controls the probability of failure. Namely, we assume
for some constant ${C} = C(B)>0$.
See the proof in Section (ref). A direct consequence of Theorem (ref) is
The proof is presented in Section (ref). Theorem (ref) and Corollary (ref) indicate that our estimator $\beta(\widehat\theta(\widehat{\lambda}))$ is robust to the choice of observational samples---it achieves an $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$ error rate regardless of the level of bias in observational data. Notably, this $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$ rate is known to be optimal and can be attained, for instance, by the AIPW estimator hahn1998role using ${N^{{\mathrm{exp}}}}$ experimental samples and no observational data. Moreover, even with a sufficiently large number of observational samples, one cannot achieve a rate faster than $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$ without imposing additional assumptions on the observational data. We demonstrate the following matching minimax lower bound on the estimation error over a class of robust estimators $\widehat{\mu}$ in the no-covariate setting:
The proof can be found Section (ref). Theorem (ref) shows that when taking both experimental and observational data as input, no robust estimator ({\it i.e.}, one with an error rate of order $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$ uniformly over $\varepsilon\in[-1,1]$) can achieve an error rate better than $\mathcal{O}(1/{N^{{\mathrm{exp}}}})$, even when $\varepsilon$ is zero.
We have proposed a simple, general method for integrating experimental and observational data, leveraging cross-validation to adaptively tune their relative contribution. Our approach requires no additional specification or identification assumptions and accommodates a broad range of scenarios beyond the scope of existing methods. We demonstrated its efficacy, adaptivity, and robustness through experiments on both real-world and synthetic datasets. Furthermore, we provided theoretical analysis showing that it is robust to the bias in observational data and achieves the minimax optimal rate over a class of robust estimators.
We focus on ATE in this paper, which allows broad applicability for transformed outcomes such as logarithms of the original outcome. Future work could extend our framework to other causal estimands, such as the conditional average treatment effect. Another direction is to explore extensions involving instrumental variables. On one hand, these tools may help reduce bias in observational components to improve upon experimental estimates, as demonstrated in our experiments. On the other hand, the generality of our framework opens opportunities to exploit problem-specific structure, such as the relationship between experimental and observational models, for tailored adaptations in case-by-case applications.
We thank Fan Chen, Zeyu Jia, Ian Waudby-Smith, and Shu Yang for helpful discussions.
The supplementary material is organized as follows: Section (ref) details the setup for Figure (ref). Section (ref) extends the discussion of prior work. Section (ref) presents the pseudocode of our proposed method along with its analysis. For experiments, implementation details and additional results are provided for the no-covariate setting (Section (ref)), the linear setting (Section (ref)), and the LaLonde dataset (Section (ref)), along with their reproducibility (Section (ref)). Finally, Sections (ref) and (ref) contain proofs organized by section.
We supplement the Figure (ref) for the CPS control group (where Figure (ref) uses the PSID group) and then provide a detailed explanation of both. \paragraph{Left versus right} In both Figures (ref) and (ref), the left subfigures estimate the control mean in the no-covariate case. The right subfigures adjust for all available covariates (that correspond to column 8 in Table (ref)). We employ five-fold cross-validation. Error bars reflect $\pm$1 bootstrapped standard deviations.
\paragraph{Top panel} The top panels illustrate the process of cross-validating $\lambda$. In each run, experimental data are split into $K$ folds to perform cross-validation, and the $\lambda$ that minimizes the cross-validation objective $\mathsf{CV}(\lambda)$ is selected. The curves in the top panels show $\mathsf{CV}(\lambda)$ averaged over 5000 runs, and the average selected $\widehat\lambda$ is marked by the blue dashed vertical line. \paragraph{(Top panel) Why does the average $\widehat\lambda$ not minimize the average $\mathsf{CV}(\lambda)$ in Figures (ref) and (ref)?} As shown by the density plots ({\it i.e.}, how often we select a particular $\lambda$ over 5000 runs), there is a small tail around 0 that skews the average of $\widehat \lambda$ leftwards.
\paragraph{Bottom panel} The bottom panels illustrates ATE estimates across different $\lambda$: $\lambda=0$ corresponds to using experimental data alone; $\lambda=1$ to using observational data alone; and the light blue curve in between show the estimates for $\lambda\in(0,1)$.
\paragraph{(Bottom panel) Why does our method's estimate (the red square) not exactly align with the light blue curve in the back?} We note that our method’s ATE point estimate (the red square) does not necessarily coincide with this curve at the average $\widehat\lambda$. This discrepancy arises because we may select different $\lambda$ in different runs (due to varying $K$-fold splits), whereas the light blue curve represents the average estimate at a fixed $\lambda$. Similarly, the bootstrapped standard deviation for our method incorporates uncertainty from both data resampling and cross-validation splitting, while the light blue shaded area reflects data resampling alone at a fixed $\lambda$.
Unconfoundedness in observational data is inherently untestable, but there have been efforts to assess it indirectly. With just the observational data, sensitivity analysis was proposed to measure the impact of potential unmeasured confounders on estimated causal effects by sensitivity parameters rosenbaum1983assessing, imbens2003sensitivity. To mitigate bias from unmeasured confounding, a major advance is the development of doubly robust estimators robins1994estimation, rotnitzky1998semiparametric, bang2005doubly, van2006targeted. These estimators remain consistent provided that either the treatment assignment mechanism (propensity score) or the outcome model is correctly specified. Among them, the AIPW estimator combines regression-based outcome modeling with inverse probability weighting (IPW) to achieve double robustness robins1994estimation. We use it as our experimental component in Section (ref).
We review the following three lines of work that are popular and most relevant in the space:
First, a widely adopted and methodologically straightforward line of work is pooling, which aggregates all the samples together and treat the pooled data as if it comes from a single study ross2009pooled. One drawback is that it breaks the randomization in experimental data, possibly resulting in a biased overall estimate. Follow-up work introduces a test-then-pool strategy to mitigate this aspect by conducting hypothesis testing to decide whether to include observational data gao2023pretest, yang2023elastic. Specifically, yang2023elastic performs hypothesis tests on transportability (whether the observational estimate aligns with the experimental data) and internal validity (to check for unmeasured confounding in the observational data). If the test passes, the method derives the efficient estimate on both data sources; otherwise, it relies solely on the experimental data. However, this approach generally requires common support between the datasets—when covariate overlap is insufficient, the transportability assumption likely fails, and the test excludes observational data. Nonetheless, when there is no common support, the test automatically fails. In contrast, our method can adapt to the scenario where covariates of both data sources are completely different. Moreover, rather than making an all-or-nothing decision, our method offers flexibility by adjusting the weight assigned to each data source, allowing it to adapt to a wider range of scenarios. In fact, this line of work can be viewed as a specific case in our framework (assigning weight $1$ to the pooled source). We compare our proposed method with the pooling approach in Section (ref).
Second, from a statistical perspective, combining biased and unbiased estimators has been studied through techniques in Stein Shrinkage and Empirical Bayes stein1956inadmissibility, green1991james, green2005improved. For causal setting, rosenman2023combining uses James-Stein type shrinkage estimator on the strata of samples based on (stabilized) IPW estimators that do not require an outcome model. This approach operates on fixed, predefined strata and allows residual bias of unknown magnitude to remain. A key limitation is its reliance on stratification and the strong assumption that ATEs are equal—or differ by at most $\mathcal{O}(1/n)$—across data sources within each stratum, which may not hold in practice. Furthermore, like classical Stein shrinkage, it requires at least four strata to ensure risk reduction. While their method and ours share the high-level idea of weighting, we are conceptually different: rosenman2023combining uses stratum-level weighting to directly combine estimates, whereas our approach performs loss-level weighting within a model-agnostic empirical risk minimization framework. This allows us to avoid assumptions about per-stratum ATE equality and to flexibly incorporate different types of models beyond (stabilized) IPW.
Third, a relatively assumption-light approach is error-prone estimators—estimators derived from two data sources that are individually biased for the ATE but share the same expected bias yang2020combining. They first construct an asymptotically normal estimator from experimental data, and then adds and subtracts two such error-prone estimators—one from each source—to cancel out the bias in expectation. It is assumption-light in a way that it permits different outcome models and scenarios without covariate overlap—provided that the error-prone estimators can be constructed using only treatment and outcome. While their approach and ours both leverage the consistency of the estimator derived from experimental data, the way we incorporate this consistency differs. While their method perform on bias-cancellation through algebraic manipulation of two error-prone estimators assumed to share the same expected bias, our method is a joint optimization over experimental and observational loss functions with a tunable trade-off to anchor for consistency. One limitation of their approach lies in the additional specification of the multi-dimensional error-prone estimators, which, as noted in their Remark 3, can significantly affect the efficiency of the overall estimator. Finally, their theoretical guarantees are primarily asymptotic, while ours are non-asymptotic, providing bounds that hold in finite-sample regimes, which is especially desirable when experimental data are limited.
Techniques for combining multiple statistical or machine learning estimators via data-driven weighting have a rich history, offering improvements over single-model selection. Well-established methods including stacking wolpert1992stacked, breiman1996stacked, aggregation tsybakov2003optimal, tsybakov2004optimal, and super learner van2007super leverage cross-validation to determine weights that optimally blend different estimators. The goal is typically to enhance predictive performance and robustness by integrating the strengths of diverse models, avoiding the brittleness of relying on a single “best” model. While developed primarily for general-purposed prediction tasks, the underlying principle of using cross-validation to build robust, data-driven combinations of estimators holds significant potential for causal inference. However, adapting these powerful tools for causal inference requires careful methodological design due to non-trivial challenges, such as: ensuring adherence to identification assumptions, appropriately incorporating the information on treatment assignment mechanism, and selecting cross-validation criteria specifically targeted at the causal objective rather than just predictive accuracy. Our work presents a principled way to conduct cross-validated causal inference to combine experimental and observational data.
The first panel represents whether the method can give a consistent estimate in the presence of outcome model misspecification (outcome model misspecification) for experimental data. The second panel represents whether the overall method can give a consistent estimate when observational data has unmeasured confounders (unmeasured confounders), outcome model misspecification (outcome model misspecification), or both (both). \textbf{The third panel} represents whether each model allows an inconsistent observational estimate to be included in the final result (\texttt{inconsistent observational estimate}), common covariates having different distributions (\texttt{shift in common covariates}), completely non-overlapping covariates (\texttt{no covariate overlap}), different experimental and observational outcome models (\texttt{different outcome models}), no additional model specifications (\texttt{allow no extra model specification}), and treatment acting differently on either data sources after eliminating observational biases (\texttt{allow different ATE across sources}).
Algorithm (ref) proceeds as follows: Line 1-4 define a subroutine that fits a model by minimizing a combination of the experimental and observational losses, where the weight is given by $\lambda$. Line 5-14 evaluate the performance of the models fit using each candidate $\lambda$ via $K$-fold cross-validation. Importantly, only the experimental dataset is partitioned for training and evaluation during cross-validation. The value $\widehat\lambda$ that yields the lowest average cross-validation loss is then selected. A final model $\widehat\theta(\widehat\lambda)$ is trained using the full dataset.
Our method involves training models $\mathcal{O}(K|\Lambda|)$ times, with the overall complexity depending on the cost of each individual training. For example, in the no-covariate case, each training reduces to computing sample means, which takes $\mathcal{O}(N^{\mathrm{exp}}+N^{\mathrm{obs}})$ time. Under the linear setting, each training requires solving linear systems. For an observational linear model with $d^{\mathrm{obs}}$ covariates, the closed-form solution can be computed in up to $\mathcal{O}((d^{\mathrm{obs}})^2 N^{\mathrm{obs}} + (d^{\mathrm{obs}})^3)$, depending on the solver. To compute the experimental estimate, using a linear outcome model with $d^{\mathrm{exp}}$ covariates for the plug-in estimator or AIPW estimator requires up to $\mathcal{O}((d^{\mathrm{exp}})^2 N^{\mathrm{exp}} + (d^{\mathrm{exp}})^3)$ time. In practice, the cross-validation step (Lines 7–11 in Algorithm (ref)) could be implemented efficiently by batching computations for multiple $\lambda$ values in parallel.
We utilize the closed-form solution in Eq. (ref). We set $\tau^\star=0.5$, where the specific value chosen would not affect the qualitative results. For cross-validation, we set $K=N^{\mathrm{exp}}$ and conduct a grid search over for candidate values of $\lambda \in [0, 1]$ in 50 linearly spaced bins. The t-test baseline is as follows: the null hypothesis is that the two populations have the same mean, while the alternative is that their means differ. If it fails to reject the null hypothesis, we set $\lambda = 0$ to rely solely on experimental samples. Otherwise, we set $\lambda = N^{\mathrm{obs}} / (N^{\mathrm{exp}}+N^{\mathrm{obs}})$ to incorporate both sources. For experiments varying $N^{\mathrm{obs}}$ (or $N^{\mathrm{exp}}$), we generate a large observational (or experimental) dataset and draw random subsets of the desired size for each run. We repeat 5000 runs for each experiment.
For figure production, the insets in Figure (ref), (ref), (ref), and (ref) display zoomed-in views of the plots over $\varepsilon \in [0, 2]$ to highlight the performance gains in that region, and over $\varepsilon \in [0.53, 1.47]$ to provide a closer examination of the model's behavior. In Figures (ref), (ref), (ref), (ref), we apply a continuous piecewise transformation to the vertical axis to improve visual clarity. Specifically, values below a threshold $b$ are scaled linearly, while values above $b$ are log-transformed relative to the threshold. This transformation takes the form \[ stretch(y) =
\] where $a$ controls the intensity of the stretch and ensures continuity at the transition point $y = b$. This approach preserves detail for small values while compressing the dynamic range of larger values, making trends and comparisons more visually accessible. The transformation is invertible, allowing us to recover the original values on the vertical axis. We set $b$ to be the maximum of our method's empirical MSE, and $a$ to be $5$.
Raising the noise level $\sigma^2$ from 1 to 100, we observe that each sub-figure in Figure (ref) mirrors its counterpart in Figure (ref). While the overall behaviors remain qualitatively unchanged, the MSEs scale up by a factor of roughly 100. This is due to the bias-variance decomposition of MSE, where the variance component dominates as the noise level increases. The scaling also shifts the threshold of $\varepsilon$ beyond which biased observational data lose its utility: from $\varepsilon \approx 0.125$ in Figures (ref) and (ref) to $\varepsilon \approx 1.25$ in Figures (ref) and (ref). \thispagestyle{empty}
We utilize the closed-form solution in Eq. (ref). For each experimental and observational sample, we independently generate the covariates $ Z \sim {\mathcal N} (0, \sigma^2 {\mathbf I})$, the response $W \sim \mathrm{Bern}(0.5)$ for experimental data and $\mathrm{Bern}(0.2)$ for observational data, and an exogenous noise $\xi \sim {\mathcal N} (0, \sigma^2) \perp Z, W$. We set $\sigma^2=1$. For experimental samples, the response is generated as $Y = Z^\top \theta^{\mathrm{exp}} + W \tau^\star + \xi$. We set $\tau^\star=0.5$. For observational samples, we introduce the bias via $Y = Z^\top \theta^{\mathrm{exp}} + W (\tau^\star + \varepsilon) + \xi$. The weights of $\theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}}$ are sampled from a multivariate normal distribution ${\mathcal N}(0, {\mathbf I})$. We then append a 1 to each $Z$ and $0$ to $\theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}}$ to account for the intercept term. The dimensions of $\theta^{\mathrm{exp}}$ and $\theta^{\mathrm{obs}}$ are set to 6 (including the intercept). For the experiments varying $\varepsilon$, weights are sampled independently for each simulation. For experiments varying $N^{\mathrm{obs}}$, weights are sampled once to generate a large observational dataset, from which random subsets of the desired size are drawn in each run. For cross-validation, we set $K=5$ and use 50 linearly spaced bins for candidate values of $\lambda$. To calculate the experimental estimate $\widehat\tau^{\mathrm{exp}}$, we employ the average of AIPW estimates with a known propensity score (0.5). A linear outcome model is fit on half of the experimental data, and the AIPW estimates are computed using the remaining half. When splitting the data either for computing the AIPW estimate or for cross-validation, we stratify by treatment assignment, resulting in each fold containing approximately 50% treated and 50% control samples. We repeat 5000 runs for each experiment.
For figure production, the insets in Figures (ref), (ref), (ref), (ref), (ref), and (ref) provide zoomed-in views: over small values of $\varepsilon$ to highlight performance gains, and over the range $\varepsilon \in [24.49, 118.37]$ to enable a closer examination of the model's behavior. We apply the same linear-log transformation described in Section (ref) to figures involving varying $N^{\mathrm{obs}}$. The threshold $b$ is set to the maximum MSE of our method. The transformation intensity parameter $a$ is set to 3 in Figures (ref) and (ref), and to 5 in Figures (ref) and (ref).
Dehejia and Wahba's paper narrowed the focus to male participants who were under 55 years of age at the time of the program's initiation. This specific subgroup was chosen because it allows for a consistent examination of labor outcomes and eliminates the potential impact of retirement. The selection criteria were further refined to those individuals who had earnings data available for both 1975 (pre-intervention) and 1978 (post-intervention). This subset comprised 297 treated and 425 control participants. However, to enhance the analysis's robustness and to focus on those with more complete data profiles, the dataset was further narrowed to those who also had earnings data available for 1974, reducing the sample to 185 treated and 260 control participants. As illustrated in Figure (ref), the subset selection is based on only pre-intervention variables. Table (ref) shows the proportion of each pre-intervention variables in treatment and control group, for the sake of assessing internal validity.
We now discuss our reproduced results from Dehejia and Wahba in the first and third panels in Table (ref). The correspondence is as follows: point estimates in our column 1 align exactly with their panel B(1) and C(1), 2 with B(2), 3 with B(3) and C(3), 4 with B(4), 6 with C(2), and 7 with C(4). Interestingly, our column 5 (which regresses on all covariates excluding RE74 and 1974 employment status) yields the most accurate statistically significant ATE estimate (\$1167). This result was not highlighted by either Dehejia and Wahba or LaLonde, suggesting a potentially overlooked finding. Notably, this improvement is achieved using CPS-1, which is not considered a specially selected subgroup. This outcome contradicts LaLonde’s assertion that subgroups such as PSID-2, PSID-3, CPS-2, and CPS-3 are more comparable to the NSW control group.
We then critically assess Dehejia and Wahba’s claim that 1974 earnings (RE74) are a valuable predictor in estimating treatment effects. In our reproduced results, the impact of RE74 could be mixed. For example, we compare results in columns 2 and 6, where they share the same set of covariates but column 6 additionally includes RE74. Column 2 has four statistically significant estimates of treatment effect, none of which are statistically significant in column 6. The column 2 estimates deviate a large amount from those of the NSW data, implying that the training had a negative impact on future earnings. In column 6, incorporating RE74 gives less negative results, but at the cost of losing statistical significance.
Interestingly, RE74 appears to add value in a specific setting: the linear model in column 6, applied to the CPS-3 subgroup, produces the best near-significant estimate ($1,326$) across all configurations with a p-value of 0.09. Although this does not meet our predetermined 0.05 significance threshold, it still indicates marginal significance. This suggests that CPS-3 may serve as a promising comparison group. Additionally, we note that caution is necessary when interpreting the CPS-3 results. The selection criteria for CPS-3 excluded individuals with 1975 incomes below the poverty line, whereas NSW participants were not restricted in this way. Specifically, the ATE may differ depending on whether individuals' incomes in 1975 were above or below the poverty threshold. The result may be attributable to good luck rather than a meaningful underlying effect.
We now present the full configurations. The observational data are partitioned into six subgroups:
Each column represents the estimated effect of treatment according to a specific linear setting as follows:
In the following Table (ref), each panel contains result for different methods detailed in Section (ref). Each row represents the data configuration with (T) for treatment group and (C) for control group. For the second panel, we report the estimated treatment effect with $\pm 1$ standard deviation over 5000 runs, followed by $\widehat\lambda$ in parentheses selected by five-fold cross-validation. For the other panels, the p-values (in parentheses) comes from testing the null hypothesis that the treatment coefficient is zero. Statistically significant results (under $0.05$) are in bold.
We note that in Table (ref), the reported standard deviations of our method come from random $K$-fold splits in each run. In contrast, Table (ref) presents bootstrap standard deviations: for our method, this captures uncertainty from both data resampling and cross-validation splitting, while for the other methods, it reflects uncertainty from data resampling alone.
\thispagestyle{empty}
\thispagestyle{empty}
Each applicable setting in this work are repeated 5000 times. Given the large number of replications, we expect our results to be robust to the choice of random seed, as random fluctuations introduced by any particular seed are likely to average out. Codes are available in \url{https://github.com/xyang23/cross_validated_causal}.
We now derive Eq. (ref), which is \[ \widehat\theta(\lambda) = \arg \min_\theta \sum_{i = 1}^{N^{\mathrm{exp}}} (1 - \lambda) (Y_i^{\mathrm{exp}} - \theta)^2 + \lambda (\overline{Y}^{\mathrm{obs}} - \theta)^2 = (1 - \lambda) \overline{Y}^{\mathrm{exp}} + \lambda \overline{Y}^{\mathrm{obs}}. \] The result follows from the following calculation:
Moreover, we note that the following four minimizers are equivalent:
The equivalence of these formulations follows directly from Lemma (ref). Specifically, the first and second terms in each formulation resemble $ (1 - \lambda) \theta^2 - 2 (1 - \lambda) \overline{Y}^{\mathrm{exp}} \theta$ and $\lambda \theta^2 - 2 \lambda \overline{Y}^{\mathrm{obs}} \theta$, respectively, up to additive constants that do not affect the minimizer. This equivalence implies that aggregate- and unit-level losses yield the same minimizer, reflecting their alignment in the underlying principle across granularity.
Informally, treating the experimental data as fixed, the squared error between a given scalar (e.g., the causal parameter) and the average experimental estimate is proportional to the average squared error between that scalar and each individual estimate, up to constants.
We prove the following additive property for squared loss:
We now derive the closed-form solution for the linear setting. We are interested in
We can write $\beta(\theta) = e_1^\top \theta$, where $e_1^\top = ( 1 \quad 0 \cdots 0 ) $. Then we have
We take the gradient with respect to $\theta$ and set it to 0:
Solving this linear system gives the desired minimizer. When $\lambda = 0$, the minimizer may not be unique, but every solution must satisfy $\beta(\theta)=\widehat\tau^{\mathrm{exp}}$, thereby matching the experimental estimate. When $\lambda=1$, the objective reduces to ordinary least squares on observational data, yielding the observational estimate.
By Lemma (ref), a sufficient condition for Assumption (ref) is the following condition assuming that $\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}})$ is derived from some $Z$-estimation problem. In this case, $(B_{\widetilde{h}},B_{\tau^{\star}},B_{\tau^{\star},0},B_{\tau^{\star},1},B_{\tau^{\star},\mathrm{num}})$ can be chosen as constants that depend polynomially on the parameters $(d_{\eta},1/\gamma,B_{{\sf H}},B_{h,0},B_{h,1},B_{h,2})$ in Assumption (ref).
In Assumption (ref), we posit that the ATE estimator $\widehat\tau^{\mathrm{exp}}$ is given by the first coordinate of some $Z$-estimator. Specifically, Assumption (ref) assumes that the true ATE $\tau^{\star}$ equals the first coordinate of the true parameter $\eta^{\star}$ of the $Z$-estimation problem. Note that this can be generalized to any linear function of $\eta^{\star}$ by a simple change of variables. Assumption (ref) imposes standard smoothness conditions on the $Z$-function and its derivatives. Assumption (ref) assumes $\sqrt{{N^{{\mathrm{exp}}}}}$-convergence of the $Z$-estimator. This is satisfied when e.g., the $Z$-function is the gradient of some convex loss. In fact, a sufficient condition for Assumption (ref) is the following convexity condition (ref). We refer to Lemma (ref) for more details.
It is readily verified that the ordinary least squares (OLS) estimator satisfies Assumption (ref) when the observed outcome $Y^{{\mathrm{exp}}}_i$ is linear in the covariates $Z^{{\mathrm{exp}}}_i$ and the treatment assignment $W^{{\mathrm{exp}}}_i$. Additionally, under proper conditions, the inverse propensity weighted (IPW) estimator horvitz1952generalization satisfies Assumption (ref) when the true propensity score $p(Z^{{\mathrm{exp}}}_i)\coloneqq \ensuremath{{\mathbb{P}}}(W^{{\mathrm{exp}}}_i=1|Z^{{\mathrm{exp}}}_i)$ follows a logistic model, {\it i.e.}, $p(Z^{{\mathrm{exp}}}_i) = \exp(Z_i^{{\mathrm{exp}}\top }\omega^\star)/(1+\exp(Z_i^{{\mathrm{exp}}\top } \omega^\star))$ for some $\omega^\star\in{\mathbb R}^{d_{{\mathrm{exp}}}}$, and is estimated via logistic regression (see Example 3 in lin2024worthwhile).
We now restate and clarify the notation. For any set $\mathcal{J} \subseteq [{N^{{\mathrm{exp}}}}]$, we define $X^{{\mathrm{exp}}}_\mathcal{J} \coloneqq (X^{{\mathrm{exp}}}_i)_{i \in \mathcal{J}}$ as the subset of experimental samples indexed by $\mathcal{J}$. In particular, recall that $X^{{\mathrm{exp}}}_{B_i}$ denote the set of experimental samples in the $i$-th fold, for $i \in [K]$. We write $X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}]} = X^{{\mathrm{exp}}}$ and $X^{{\mathrm{obs}}}_{[{N^{{\mathrm{obs}}}}]} = X^{{\mathrm{obs}}}$ to denote the full set of experimental and observational samples, respectively. With this notation, the full dataset is $D = (X^{{\mathrm{exp}}}, X^{{\mathrm{obs}}}) = (X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}]}, X^{{\mathrm{obs}}}_{[{N^{{\mathrm{obs}}}}]})$, and the dataset excluding the $i$-th experimental fold is $D_{-i} = (X^{{\mathrm{exp}}}_{-B_i}, X^{{\mathrm{obs}}}) = (X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}] \setminus B_i} , X^{{\mathrm{obs}}}_{[{N^{{\mathrm{obs}}}}]})$, for $i \in [K]$. We write $\widehat\theta(\lambda)=\widehat\theta(\lambda;D)$ to specify the dependence of $\widehat\theta(\lambda)$ on $D.$ We also define $D^{{\mathrm{obs}}} \coloneqq X^{{\mathrm{obs}}} = X^{{\mathrm{obs}}}_{[{N^{{\mathrm{obs}}}}]}$.
For each subset of experimental samples $X^{{\mathrm{exp}}}_{\mathcal{J}}=(X^{{\mathrm{exp}}}_j)_{j\in[\mathcal{J}]}$, we write the experimental loss ${L^{\mathrm{exp}}}(\beta(\theta); X^{{\mathrm{exp}}}_{\mathcal{J}}) = (\beta(\theta) - \widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}}))^2$, where $\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}})$ denotes an estimate of the average treatment effect (ATE) based on the samples indexed by $\mathcal{J}$. We also write ${L^{\mathrm{exp}}}(\beta(\theta); P^{{\mathrm{exp}}}) = (\beta(\theta) - \tau^{\star})^2$ for the population loss. In addition, for any function $f$, with slight abuse of notation, we let $\widehat{{\mathbb E}}_{\mathcal{J}}[f(X^{{\mathrm{exp}}})] \coloneqq \frac{1}{|\mathcal{J}|} \sum_{j \in \mathcal{J}} f(X^{{\mathrm{exp}}}_j)$ denote the empirical average over a subset $\mathcal{J}$ of the experimental samples.
We use $\| \cdot\|_{2}$ to denote the Euclidean norm for vectors and $\ensuremath{|\!|\!| \cdot | \! | \!|_{{\tiny{\mbox{op}}}}}$ to denote the spectral norm (or operator norm) for matrices and third-order tensors. Concretely, for a third-order tensor $\mathcal{T} \in \mathbb{R}^{d_1 \times d_2 \times d_3}$, its spectral norm (or operator norm) is defined as \[ \ensuremath{|\!|\!| \mathcal{T} | \! | \!|_{{\tiny{\mbox{op}}}}} \coloneqq \sup_{\| x\|_{2} = \| y\|_{2} = \| z\|_{2} = 1} \sum_{i=1}^{d_1} \sum_{j=1}^{d_2} \sum_{k=1}^{d_3} \mathcal{T}_{ijk} x_i y_j z_k. \]
Throughout the proofs, we use ${C},{C'}>0$ to denote constants that depend polynomially on the parameters in the assumptions. We allow their values to change from place to place. More specifically, when Assumption (ref) and (ref) hold, the constants ${C}=C(B)$ (or ${C'}=C'(B)$) depends polynomially on the parameters $(B_{\Theta},1/b_{\mathrm{obs},2},B_{\mathrm{obs},2},B_{\mathrm{obs},3}; B_{\widetilde{h}},B_{\tau^{\star}},\\B_{\tau^{\star},0},B_{\tau^{\star},1},B_{\tau^{\star},\mathrm{num}})$. Alternatively, when Assumption (ref) and (ref) hold, the constants ${C}={C(d,\gamma,B)}$ (or ${C'}={C'(d,\gamma,B)}$) depends polynomially on the parameters $(B_{\Theta},1/b_{\mathrm{obs},2},B_{\mathrm{obs},2},B_{\mathrm{obs},3};d_{\eta},1/\gamma,B_{{\sf H}},B_{h,0},B_{h,1},B_{h,2})$. The set of parameters the constants ${C}, {C'}$ depend on should be clear from context, as it only depends on what assumptions are made. We therefore omit the explicit dependence in the notation.
Under Assumption (ref) and (ref) and the sample size condition $\sqrt{{N^{{\mathrm{exp}}}}}\geq {C'} {K}(\log^{1.5}K\\+\log^{0.5}(1/\delta))$ in Eq. (ref), we will show that
See the proof in Section (ref).
See the proof in Section (ref).
With the two lemmas at hand, we are ready to prove Theorem (ref). Let
be the optimal regularization parameter that minimizes the estimation error given the dataset $D$. Since $\widehat\theta(0;D)=\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}]})$ satisfies $|\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}]})-\tau^{\star}|\leq B_{\tau^{\star},0}\sqrt{\log(1/\delta)}/\sqrt{{N^{{\mathrm{exp}}}}}$ with probability at least $1-\delta$ by Assumption (ref), we have
Consequently, on the event where Eq. (ref) and (ref) hold, we have
Adopt the shorthands $\widehat\tau^{\mathrm{exp}}_i=\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{B_i}),\widehat\tau^{\mathrm{exp}}_{-i}=\widehat\tau^{\mathrm{exp}}(\{X^{{\mathrm{exp}}}_{B_j},j\neq i \}), \widehat\tau^{\mathrm{exp}}=\widehat\tau^{\mathrm{exp}}(\{X^{{\mathrm{exp}}}_{B_j},j\in[K] \})$. Also define
By some basic algebra, we have
We make the following claims which will be shown at the end of the proof:
By claim (ref), when $\lambda>0$, we have
for some $R_1$ such that $|R_1|\leq {C}\log^{1.5}(K/\delta)/{N^{{\mathrm{exp}}}}/\sqrt{{N^{{\mathrm{exp}}}}/K}\leq{C}\log(1/\delta)/{N^{{\mathrm{exp}}}}$ with probability at least $1-\delta$. Moreover, we have by Eq. (ref) in Lemma (ref) that
with probability at least $1-\delta$. Combining the last two bounds and using claim (ref) yields
for all $\lambda\in(0,1]$ for some ${C}=C(B)>0$ with probability at least $1-\delta$. The bounds on $R_1$ for the case $\lambda=0$ is similar and we thus omit the details.
Moreover, for ${\mathcal R}_2$, we have with probability at least $1-\delta$ that, for all $\lambda\in[0,1]$,
where the second inequality follows from Eq. (ref) in Lemma (ref). Finally, note that
where the first inequality uses $(a+b)^2\leq 2(a^2+b^2)$, and the second inequality follows from Lemma (ref). Combining the bounds on ${\mathcal R}_1$, ${\mathcal R}_2$ and ${L^{\mathrm{exp}}}(\beta(\widehat\theta(\lambda;D^{{\mathrm{obs}}})); P^{{\mathrm{exp}}})$ yields the desired result.
\paragraph{Proof of claim (ref).} Note that $\widehat\theta(\lambda;D^{{\mathrm{obs}}}),\widehat\theta(\lambda;D_{-i})$ are empirical risk minimizers. Taking the derivatives with respect to $\theta$, we have
Introduce the shorthand $\widetilde\Delta_i\coloneqq \widehat\theta(\lambda;D_{-i})-\widehat\theta(\lambda;D^{{\mathrm{obs}}})$. Taking the difference and performing a Taylor expansion yields
where
with $e_1\coloneqq (1,0,\cdots,0)^\top\in{\mathbb R}^{d_{\theta}}$ and $E_{11}\in{\mathbb R}^{d_{\theta}\timesd_{\theta}}$ being the matrix where the $(1,1)$-th entry is one and all other entries are zero. Recall that $T(\lambda)= \lambda \nabla_\theta^2 {L^{\mathrm{obs}}}(\widehat\theta(\lambda;D^{{\mathrm{obs}}});D^{{\mathrm{obs}}}) + (1-\lambda)\nabla_\theta^2 {L^{\mathrm{exp}}}(\beta(\widehat\theta(\lambda;D^{{\mathrm{obs}}}));P^{{\mathrm{exp}}})$. By Lemma (ref), we have $\| \widetilde\Delta_i\|_{2}\leq {C} |\widehat\tau^{\mathrm{exp}}_{-i}-\tau^{\star}|$. Therefore,
for some ${C'}=C'(B)>0$ with probability at least $1-\delta$, where the first inequality uses Eq. (ref) and the last inequality follows from Assumption (ref). Finally, applying an union bound over all $i\in[K]$ yields the desired result in Eq. (ref). The case $\lambda=0$ follows immediately from Lemma (ref).
\paragraph{Proof of claim (ref).} The proof of claim (ref) follows from the same arguments as in the proof of claim (ref) since $\widehat\theta(\lambda;D_{-i})$ in the proof of claim (ref) can be replaced by $\widehat\theta(\lambda;D)$ without loss of generality.
\paragraph{Proof of claim (ref).} By the expression of Schur's complement, we have
Thus, we have $2+\frac{\lambda}{1-\lambda}\Big(b_{\mathrm{obs},2} -\frac{B_{\mathrm{obs},2}^2}{b_{\mathrm{obs},2}} \Big)\geq 1$ (and therefore $|(1-\lambda) e_1^\topT(\lambda)^{-1}e_1|\leq 1$) when $\lambda\leq 1/{C}_1$ for some ${C}_1=C(B)>0$ sufficiently large. On the other hand, when $\lambda\geq {C}_1$, we have
Combining the two cases completes the proof.
By defintion of ${L^{\mathrm{exp}}}$, we have
Similarly to the proof of Lemma (ref), we claim that
The proof of this claim will be given momentarily.
With these two claims at hand and using Eq. (ref) and Assumption (ref), we have
for some $R_3$ such that $|R_3|\leq \log^2(K/\delta)/({N^{{\mathrm{exp}}}})^2$ for all $\lambda\in(0,1]$ with probability at least $1-\delta$. The bound on ${\mathcal R}_3$ when $\lambda=0$ follows similarly. Thus we have
for all $\lambda\in[0,1]$ with probability at least $1-\delta$ since $\sqrt{{N^{{\mathrm{exp}}}}}\geq \polyshortK(\log^{1.5}(K)+\log^{0.5}(1/\delta))$.
Moreover, for ${\mathcal R}_4$, we have by the Cauchy-Schwarz inequality that
for all $\lambda\in(0,1]$ with probability at least $1-\delta$. Combining the bounds on ${\mathcal R}_3$ and ${\mathcal R}_4$ completes the proof.
\paragraph{Proof of claim (ref).} Since $\widehat\theta(\lambda;D),\widehat\theta(\lambda;D)$ are both assumed to be the empirical risk minimizer on the respective datasets, we have
Let $\widebar{\Delta}_i\coloneqq \widehat\theta(\lambda;D_{-i})-\widehat\theta(\lambda;D)$. Taking the difference and performing a Taylor expansion yields
where
with $e_1\coloneqq (1,0,\cdots,0)^\top\in{\mathbb R}^{d_{\theta}}$ and $E_{11}\in{\mathbb R}^{d_{\theta}\timesd_{\theta}}$ being the matrix where the $(1,1)$-th entry is one and all other entries are zero. Recall that $\widetilde{T}(\lambda)\coloneqq \lambda \nabla_\theta^2 {L^{\mathrm{obs}}}(\widehat\theta(\lambda;D);D^{{\mathrm{obs}}}) + (1-\lambda)\nabla_\theta^2 {L^{\mathrm{exp}}}(\beta(\widehat\theta(\lambda;D));X^{{\mathrm{exp}}}_{[{N^{{\mathrm{exp}}}}]})$. By Lemma (ref), we have $\| \widetilde\Delta_i\|_{2}\leq {C} |\widehat\tau^{\mathrm{exp}}_{-i}-\widehat\tau^{\mathrm{exp}}|$. Therefore, similar to the proof of Lemma (ref),
for some ${C'}=C'(B)>0$ with probability at least $1-\delta$, where the first inequality uses Eq. (ref) and the last inequality follows from Assumption (ref) and a triangle inequality. Finally, applying an union bound over all $i\in[K]$ gives Eq. (ref). The case $\lambda=0$ follows immediately from Lemma (ref).
\paragraph{Proof of claim (ref).} The proof follows from the same argument as the proof of claim (ref) in the proof of Lemma (ref). We thus omit the details here.
Note that the experimental sample size condition (ref) is satisfied with $\delta =1/{N^{{\mathrm{exp}}}}$ when $K\leq {C}\sqrt{{N^{{\mathrm{exp}}}}}/\log^{1.5}{N^{{\mathrm{exp}}}}$. Therefore, we have by Theorem (ref) that
with probability at least $1-\delta$ for any $\delta\geq 1/{N^{{\mathrm{exp}}}}$. Let $\mathcal{E}\coloneqq\{(\beta(\widehat\theta(\widehat{\lambda};D))-\tau^{\star})^2\geq {C'}\log({N^{{\mathrm{exp}}}})/{N^{{\mathrm{exp}}}}\}$. Then we have $\ensuremath{{\mathbb{P}}}(\mathcal{E})\leq 1/{N^{{\mathrm{exp}}}}$. Thus,
This completes the proof.
We prove the theorem by contradiction. Let $\Delta\in[0,1/2]$ be some value which will be specified later. It there exists some $\widehat{\mu}\in{\mathcal{M}}_{c_1}$ such that
for some value $L>0$, then we have by Chebyshev's inequality that, for $\tau^{\star} = 0$,
Suppose we have chosen $\Delta$ such that
Then it follows from the triangle inequality that
and therefore
We will show at the end of the proof that there exist some absolute constants $c_3,c_4>0$ such that, when $L\leq c_3$, one can choose $\Delta = \min\{\sqrt{c_4\log(1/L)}/\sqrt{{N^{{\mathrm{exp}}}}},{1}/{2}\}$ such that Eq. (ref) (and therefore Eq. (ref)) holds.
As a consequence of Eq. (ref), $\widehat{\mu}$ does not belong to the class ${\mathcal{M}}_{c_1}$ for $c_1\leq \min\{c_4\log(1/L)/2,\\{{N^{{\mathrm{exp}}}}}/8\}$ when Eq. (ref) holds. Therefore, conversely, for the absolute constant $\widetilde{c}_1\coloneqq c_4\log(1/c_3)/2$ and any $c_1\in[\widetilde{c}_1,{N^{{\mathrm{exp}}}}/8]$, Eq. (ref) is not satisfied for any $\widehat{\mu}\in{\mathcal{M}}_{c_1} $ with any $L< \exp(-2c_1/c_4)\eqqcolon c_2$. This completes the proof.
\paragraph{Verification of Eq. (ref).} Let $\tau^{\star}=0.$ Denote the event $\{|\widehat{\mu}((Y^{{\mathrm{exp}}}_i)_{i=1}^{{N^{{\mathrm{exp}}}}};(Y^{{\mathrm{obs}}}_i)_{i=1}^{{N^{{\mathrm{obs}}}}})-\tau^{\star}|\geq \Delta\}$ by $\mathcal{E}$. Introduce the shorthand notations $\ensuremath{{\mathbb{P}}}_{\tau^{\star}}$ and $\ensuremath{{\mathbb{P}}}_{\tau^{\star}+2\Delta}$ to denote the joint distribution $(Y^{{\mathrm{exp}}}_i)_{i=1}^{{N^{{\mathrm{exp}}}}}\overset{iid}{\sim}{\mathcal N}(\tau^{\star},1), (Y^{{\mathrm{obs}}}_i)_{i=1}^{{N^{{\mathrm{obs}}}}}\overset{iid}{\sim} {\mathcal N}(\tau^{\star},1)$ and $(Y^{{\mathrm{exp}}}_i)_{i=1}^{{N^{{\mathrm{exp}}}}}\overset{iid}{\sim}{\mathcal N}(\tau^{\star}+2\Delta,1), (Y^{{\mathrm{obs}}}_i)_{i=1}^{{N^{{\mathrm{obs}}}}}\overset{iid}{\sim} {\mathcal N}(\tau^{\star},1)$, respectively. When $\ensuremath{{\mathbb{P}}}_{\tau^{\star}}(\mathcal{E}) \leq \frac{1}{8}$, we have
where step (i) follows from the formula of KL divergence between two Gaussian distributions, and step (ii) follows from data-processing inequality. Therefore, to ensure that $\ensuremath{{\mathbb{P}}}_{\tau^{\star}+2\Delta}(\mathcal{E}) \leq \frac{1}{2}$, it suffices to choose $\Delta$ such that
It can be verified that there exist absolute constants $c_3,c_4>0$ sufficiently small such that when $L\leq c_3$, $\Delta = \sqrt{c_4\log(1/L)}/\sqrt{{N^{{\mathrm{exp}}}}}$ satisfies the conditions in Eq. (ref).
See the proof in Section (ref).
See the proof in Section (ref).
Consequently, Assumption (ref) is satisfies with $\widetilde{h}(X^{{\mathrm{exp}}})=h(X^{{\mathrm{exp}}};\eta^{\star})$, $\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}})=\widehat\eta_1$ and some $(B_{\widetilde{h}},B_{\tau^{\star}},B_{\tau^{\star},0},B_{\tau^{\star},1},B_{\tau^{\star},\mathrm{num}})$ depending polynomially on the parameters $(d_{\eta},1/\gamma,\\B_{{\sf H}},B_{h,0},B_{h,1},B_{h,2})$ in Assumption (ref).
See the proof in Section (ref).
See the proof in Section (ref).
We only prove the bounds on $ \widehat\theta(\lambda;D_{-i})-\widehat\theta(\lambda;D^{{\mathrm{obs}}})$ (and $\beta(\widehat\theta(\lambda;D_{-i}))-\beta(\widehat\theta(\lambda;D^{{\mathrm{obs}}}))$ when $\lambda=0$). The bounds on $\widehat\theta(\lambda;D_{-i})-\widehat\theta(\lambda;D)$ and $\widehat\theta(\lambda;D)-\widehat\theta(\lambda;D^{{\mathrm{obs}}})$ follow from similar arguments. \\
{Case 1: $\lambda>0$:} from Eq. (ref), we have
It suffices to show $\| (1-\lambda) T_i(\lambda)^{-1} e_1\|_{2}\leq (1+B_{\mathrm{obs},2}/b_{\mathrm{obs},2})/2$. Adopt the shorthands $\mathbf{T}^1_i(\lambda)\coloneqq \int_0^1 \nabla_\theta^2 {L^{\mathrm{obs}}}(\widehat\theta(\lambda;D^{{\mathrm{obs}}})+t\widetilde\Delta_i;D^{{\mathrm{obs}}})dt$ and $\mathbf{T}^2_{i,{\varepsilon}}(\lambda)\coloneqq {\varepsilon}{\mathbf I}_{d_{\theta}}+(2-{\varepsilon})E_{11}$ for ${\varepsilon}\geq 0$. Then we have
where the first equality follows from Woodbury's matrix identity. Since $ \| \mathbf{T}^2_{i,{\varepsilon}}(\lambda)^{-1}e_1\|_{2}=1/2$, $\ensuremath{|\!|\!| \mathbf{T}^1(\lambda) | \! | \!|_{{\tiny{\mbox{op}}}}}\leq B_{\mathrm{obs},2}$ and $\ensuremath{|\!|\!| (\mathbf{T}^1_i(\lambda)+(1-\lambda)\mathbf{T}^2_{i,{\varepsilon}}(\lambda)/\lambda)^{-1} | \! | \!|_{{\tiny{\mbox{op}}}}}\leq \ensuremath{|\!|\!| \mathbf{T}^1_i(\lambda)^{-1} | \! | \!|_{{\tiny{\mbox{op}}}}}\leq b_{\mathrm{obs},2}^{-1}$, it follows that
for any ${\varepsilon}\geq 0$. When $\lambda>0$, since $T_i(\lambda)=\lambda\mathbf{T}^1_i(\lambda)+(1-\lambda)\mathbf{T}^2_{i,0}(\lambda)$ is non-singular, taking ${\varepsilon}\to 0$ in the bound above yields the desired result.\\
{Case 2: $\lambda=0$:} we have $\beta(\widehat\theta(\lambda;D_{-i}))=\widehat\tau^{\mathrm{exp}}_{-i}$ and $\beta(\widehat\theta(\lambda;D^{{\mathrm{obs}}}))=\tau^{\star}$. The result follows immediately.
\paragraph{Proof of Eq. (ref).} Eq. (ref) follows by noting that $\widehat\tau^{\mathrm{exp}}_i-\tau^{\star},i\in[K]$ are i.i.d. ${C}/\sqrt{{N^{{\mathrm{exp}}}}/K}$ sub-Gaussian random variables by Assumption (ref), and applying Hoeffding's inequality.
\paragraph{Proof of Eq. (ref).} By Assumption (ref), it can be verified that
for some $R_a$ such that $|R_a|\leq {C}\log^{1.5}(K/\delta)/({N^{{\mathrm{exp}}}}^{1.5}/K)$ with probability at least $1-\delta$. Moreover,
with probability at least $1-\delta$ for some constant ${C}=C(B)>0$, where the second line follows from Hoeffding's inequality and Bernstein's inequality (noting that $(\widehat{{\mathbb E}}_{B_i}[\widetilde{h}(X^{{\mathrm{exp}}})])^2$ is ${C}/({N^{{\mathrm{exp}}}}/K)$ sub-Exponential). Putting the pieces together and using the fact that $\sqrt{{N^{{\mathrm{exp}}}}}\geq \polyshortK(\log^{1.5}(K)+\log^{0.5}(1/\delta))$ yields Eq. (ref).
\paragraph{Proof of Eq. (ref).} Similarly, we have by Lemma (ref) that
for some $R_b$ such that $|R_b|\leq {C}\log^{2}(K/\delta)/({N^{{\mathrm{exp}}}})^2$ with probability at least $1-\delta$. Moreover, basic algebra gives
where the second inequality follows from the same argument as in the proof of Eq. (ref). Putting the pieces together and using the fact that $\sqrt{{N^{{\mathrm{exp}}}}}\geq \polyshortK(\log^{1.5}(K)+\log^{0.5}(1/\delta))\geq{C} (\log K+\log^{0.5}(1/\delta))$ yields Eq. (ref).
Adopt the shorthand notations $\widehat{\Delta}\coloneqq \widehat\eta-\eta^{\star}$. By a Taylor expansion on $\sum_{j\in\mathcal{J}} h(X^{{\mathrm{exp}}}_j;\widehat\eta)- \sum_{j\in\mathcal{J}} h(X^{{\mathrm{exp}}}_j;\eta^{\star})$, we have
Thus,
Recall that $\| \widehat{\Delta}\|_{2}\leq \frac{{C}\sqrt{\log(1/\delta)}}{\sqrt{|\mathcal{J}|}}$ with probability at least $1-\delta$ by Assumption (ref). We claim that there exist some constants ${C},{C'}>0$ such that when $|\mathcal{J}|\geq {C'} \log(1/\delta)$, with probability at least $1-\delta$,
Putting the claims together, noting that $\sigma_{\min}(\nablaH(\eta^{\star}))\geq\gamma$ and applying the triangle inequality and a union bound, we have
with probability at least $1-\delta$ when $|\mathcal{J}|\geq {C'}\log(1/\delta)$ for some constant ${C'} = {C'(d,\gamma,B)}>0$ sufficiently large. The bound on $\widehat\tau^{\mathrm{exp}}(X^{{\mathrm{exp}}}_{\mathcal{J}})-\tau^{\star}$ follows immediately from taking the first coordinate of $\widehat\eta-\eta^{\star}$.
\paragraph{Proof of the claims.} Claim (ref) follows from applying Hoeffding's inequality to each element of the matrix and a union bound; claim (ref) again follows from Hoeffding's inequality for each element of the vector and a union bound; claim (ref) uses the assumption that $\ensuremath{|\!|\!| \nabla^2h(X^{{\mathrm{exp}}};\eta) | \! | \!|_{{\tiny{\mbox{op}}}}}\leq B_{h,2}$.
First, $\sigma_{\min}(\nablaH(\eta^{\star}))\geq \gamma$ since condition (ref) assumes $\nablaH(\eta^{\star})\succeq \gamma {\mathbf I}$. The proof of the second part of this lemma follows from standard nonasymptotic analysis of the maximum likelihood estimator (MLE) (see e.g., Lemma 9 in lin2024worthwhile). Namely, we will show the following claims:
From claim (ref), we have
Combining this with claim (ref) and noting $|\mathcal{J}|\geq {C}\log(1/\delta)$ for some ${C}>0$ sufficiently large yields Lemma (ref).
\paragraph{Proof of claim (ref).} Let $H_i(\eta)$ denote the $i$-th element of the vector $H(\eta)$. For each $i\in[d_{\eta}]$, let $\widetilde{H}^i_{\eta}\coloneqq H_i(\eta)-\widehat{{\mathbb E}}_{\mathcal{J}}[ h_i(X^{{\mathrm{exp}}};\eta)],~\eta\in{\sf H}$ is a sub-Gaussian process with respect to the metric $\rho(\eta_a,\eta_b)\coloneqq \frac{B_{h,0}\cdot\| \eta_a-\eta_b\|_{2}}{\sqrt{|\mathcal{J}|}}$. Thus, by Dudley's entropy integral bound (see e.g., Theorem 5.22 in wainwright2019high), we have
where step (i) follows from the fact that ${\mathcal N}( t;\| \cdot\|_{2},{\sf H})\leq (1+2B_{{\sf H}}/t)^{d_{\eta}}$ (see e.g., example 5.8 in wainwright2019high). Combining this with a concentration inequality for functions with bounded differences (see e.g., Corollary 2.21 in wainwright2019high), we arrive at
with probability at least $1-\delta$ for some constant ${C}={C(d,\gamma,B)}>0$.
\paragraph{Proof of claim (ref).} When $\| \eta-\eta^{\star}\|_{2}\leq \frac{\gamma}{2B_{h,2}}$, we have
This proves the first case. Introduce the unit-norm vector $\widebar{\Delta}\coloneqq \frac{\eta-\eta^{\star}}{\| \eta-\eta^{\star}\|_{2}}$. Similarly, when $\| \eta-\eta^{\star}\|_{2}> \frac{\gamma}{2B_{h,2}}$, we have
where step (i) uses $\nablaH(\eta)\succeq {\boldsymbol 0}$ for all $\eta\in{\sf H}$ and step (ii) follows from
when $t\leq \gamma/(2B_{h,2})$.