Introducing the b-value: combining unbiased and biased estimators from a sensitivity analysis perspective

Introducing the b-value: combining unbiased and biased estimators from a sensitivity analysis perspective

{5pt} {5pt} {5pt} {5pt} \hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}

{\bf Keywords}: data fusion, data integration, pretest, shrinkage estimator, soft-thresholding.

Introduction

In empirical research, it is common for researchers to employ different methods to estimate the same parameter of interest. These differences may arise from the use of distinct datasets or from imposing different model assumptions on the same dataset. We motivate our paper with the following two examples of combining estimators.

When the estimators are from different datasets, e.g., Example (ref), the estimators are independent as long as the datasets are independent. When the estimators are from the same dataset but with different model assumptions, e.g., Example (ref), the estimators are dependent in general. Given access to multiple potentially dependent estimators, some unbiased but less precise and others biased but more precise, a natural question is: How can we combine the unbiased and potentially biased estimators to improve the inference with unknown bias? From the point estimation perspective, this problem has been extensively studied bickel1984parametric,green1991james, giles1993pre,chen2015data,athey2020combining,de2020empirical, rosenman2023combining,gao2023pretest, yang2025cross. Many methods have been proposed for constructing combined estimators that perform well when the bias is small and have bounded risks when the bias is large. From the statistical inference perspective, this problem is less studied. In this paper, we answer the following question: how can we conduct valid statistical inference after combining the estimators?

This question has received considerably less attention. The primary difficulty lies in the impossibility of characterizing the distribution of the combined estimator with unknown bias armstrong2025adapting. Once the information about the bias is introduced, e.g., an upper bound on its magnitude, confidence intervals for the parameter of interest become possible. In the absence of such information, we focus on the following question:

The sequence of confidence intervals provides a way to quantify how the level of bias impacts the uncertainty in point estimation. One important application of this sequence is in hypothesis testing. Suppose the null hypothesis is not rejected based on the unbiased but less precise estimator. Now consider a scenario where the statistical test based on the biased but more precise estimator rejects the null hypothesis. If we have prior knowledge suggesting the bias is small, incorporating the biased estimator may yield a more precise estimator to reject the null hypothesis. In such cases, the sequence of confidence intervals enables us to address the following question:

The idea of constructing the sequence of confidence intervals and examining how conclusions change as the assumed level of bias varies is related to sensitivity analysis in causal inference with unmeasured confounding, e.g., cornfield1959smoking,rosenbaum1983assessing,vanderweele2017sensitivity. In observational studies, sensitivity analysis assesses how the causal conclusions change with respect to different degrees of unmeasured confounding by varying the sensitivity parameter rosenbaum2002observational,ding2016sensitivity. Our proposed framework has a similar flavor: by indexing inference results over a continuum of bias levels, we can assess the robustness of statistical inference.

We formalize two statistical inference questions, confidence interval and hypothesis testing, in the context of combining unbiased and biased estimators. Under regularity conditions, the estimators satisfy a joint central limit theorem. Consequently, we present our formulation in a finite-sample Gaussian setting, assuming exact normality for both the unbiased and biased estimators. This reduction to a Gaussian model is motivated by Le Cam’s classical asymptotic argument le1956asymptotic, and because of this, our Gaussian formulation should be viewed as an asymptotic idealization rather than a restrictive finite-sample assumption.

We develop a general framework that applies to any point estimator formed by combining such estimators. Within this framework, we construct a sequence of confidence intervals indexed by the bias level, and importantly, we introduce the notion of the b-value, a critical value of the unknown maximum relative bias at which combining estimators does not yield a significant result. We examine three canonical combined estimators: the precision-weighted estimator, the pretest estimator, and the soft-thresholding estimator. For each estimator, we derive either analytically or numerically the sequence of confidence intervals and the b-value. Among the three, we advocate for the soft-thresholding estimator, as it offers robustness to unknown bias compared with the precision-weighted estimator, and exhibits lower worst-case risk and more desirable properties for confidence interval construction than the pretest estimator. We provide a Python package for the proposed methods, available at \url{https://github.com/zhexiaolin/b-value}.

Notation. For a vector $a = (a_1,\ldots,a_d)^\top \in \mathbb{R}^d$, let $\lVert a \rVert_1 = \sum_{i=1}^d \lvert a_i \rvert$, $\lVert a \rVert_2 = ( \sum_{i=1}^d a_i^2 )^{1/2}$, and $\lvert a \rvert = (\lvert a_1 \rvert, \ldots, \lvert a_d \rvert)^\top$. For vectors $a = (a_1,\ldots,a_d)^\top$ and $b = (b_1,\ldots,b_d)^\top$, let $a \odot b = (a_1 b_1, \ldots, a_d b_d)^\top$ be the Hadamard (element-wise) product, and $a \le b$ denote $a_i \le b_i$ for all $i=1,2,\ldots,d$. For a scalar $b \in \mathbb{R}$ and a set $A \subset \mathbb{R}$, we write $b - A = \{b - a : a \in A\}$, which extends naturally to vectors and sets in $\mathbb{R}^d$. We use $\Phi(\cdot)$ and $\phi(\cdot)$ to denote the cumulative distribution function and density function of the standard normal distribution, respectively. We use $\phi_{\bm{\mu}, \bm{\Sigma}}(\cdot)$ to denote the density function of the multivariate normal distribution with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$. We use $\Psi_d(\cdot; \lambda)$ to denote the cumulative distribution function of the noncentral chi-squared distribution with noncentrality parameter $\lambda$ and degrees of freedom $d$. We use $c_{\alpha}$ to denote the $(1-\alpha)$ upper quantile of the standard normal distribution.

Problem setup and a review of point estimation

Problem Setup

We consider the following setting:

In practice, both the unbiased and biased estimators for the same parameter of interest are constructed from data. Under regularity conditions, these estimators are jointly asymptotically normal with an unknown covariance matrix. As long as this covariance matrix can be consistently estimated, the problem of combining estimators reduces to an exact normality framework under Le Cam’s asymptotic framework le1956asymptotic. le1956asymptotic showed that a wide class of estimators and test statistics can be approximated in large samples by Gaussian experiments with known covariances. Within this framework, we treat the estimators as arising from Gaussian experiments in which the variances are replaced by their consistent estimates, and the validity of the inference procedures is preserved asymptotically. For this reason, we present our analysis in the exact normality setting. Nevertheless, our results represent the asymptotic limits of a broad class of more general and practically relevant inference problems. The analysis can be generalized to dependent (Section (ref)), multivariate (Section (ref)) and multiple estimators (Section (ref)) cases.

Let $\gamma = \sigma_0^2/\sigma_1^2$ be the variance ratio between $\widehat\tau_0$ and $\widehat\tau_1$. In general, the unbiased estimator is less precise but the biased estimator is more precise. Therefore, we focus on the regime in which $\sigma_0^2$ is large and $\sigma_1^2$ is small, which indicates that $\gamma$ is large.

Under Assumption (ref), we address the following central question:

We focus on constructing two-sided symmetric confidence intervals centered at prespecified point estimators. In other words, we do not let the point estimator itself depend on the bias level $\Delta$. An alternative approach would be to use a bias-dependent point estimator, where the estimator incorporates $\Delta$. We will discuss the advantages of using prespecified point estimators later in Appendix (ref).

Since $\widehat\tau_0$ is an unbiased estimator of $\tau$, one natural approach is to construct a confidence interval based solely on $\widehat\tau_0$, which is invariant to the bias $\Delta$. Given a significance level $\zeta \in (0,1)$, we can construct a standard two-sided confidence interval \[ [\widehat\tau_0 - \sigma_0 c_{\zeta/2},\; \widehat\tau_0 + \sigma_0 c_{\zeta/2}], \] where $c_{\zeta/2}$ is the $(1-\zeta/2)$ upper quantile of the standard normal distribution. While being valid regardless of the bias $\Delta$, this confidence interval may be too wide when $\widehat\tau_0$ is not precise, whose length $2\sigma_0 c_{\zeta/2}$ scales proportionally with $\sigma_0$. As we have access to the additional biased but more precise estimator $\widehat\tau_1$, a natural question arises: can we shorten the confidence interval by incorporating information from $\widehat\tau_1$?

This motivates us to combine the two estimators $\widehat\tau_0$ and $\widehat\tau_1$ to construct a two-sided confidence interval with shorter length. Consider a generic combined estimator \[ \hat\tau = \hat\tau(\widehat\tau_0, \widehat\tau_1, \sigma_0^2, \sigma_1^2) \] of $\tau$, which depends only on observed data $(\widehat\tau_0, \widehat\tau_1)$ and known variances $(\sigma_0^2, \sigma_1^2)$, but not on the unknown bias $\Delta$. If $\Delta$ were known, constructing an exact two-sided confidence interval would be straightforward, as the distributions of both estimators are known, which lead to known distribution of $\hat\tau$. However, since $\Delta$ is unknown, the exact confidence interval depends on the magnitude of $\Delta$.

To analyze the problem, we impose bounds on the bias $\Delta$ and construct two-sided confidence intervals for different levels of bias. We assume that $\lvert \Delta/\sigma_0 \rvert \le b$ for some $b \ge 0$, and study how the confidence interval based on $\hat\tau$ changes as a function of the bias bound $b$. Here we focus on the relative bias $\Delta/\sigma_0$ instead of the absolute bias $\Delta$ to ensure that the parameterization is invariant to the scale of the estimators. For a given $b$, we aim to construct a two-sided confidence interval that achieves correct coverage uniformly over all $\Delta$ satisfying $\lvert \Delta/\sigma_0\rvert \le b$. We thus define the confidence interval below.

We next impose two natural monotonicity conditions on the interval ${\mathcal I}(b, \zeta)$ below.

The first condition in Assumption (ref) requires that as we allow greater bias in $\widehat\tau_1$, the confidence interval becomes wider and contains the previous intervals. The second condition in Assumption (ref) requires that to guarantee higher coverage rate, the confidence interval must widen. A common class of intervals satisfying both conditions is the class of symmetric fixed-length centered intervals, whose length does not depend on $(\widehat\tau_0, \widehat\tau_1)$, i.e., ${\mathcal I}(b,\zeta) = [-c(b,\zeta, \sigma_0^2, \sigma_1^2), c(b,\zeta, \sigma_0^2, \sigma_1^2)]$ for some $c(b,\zeta, \sigma_0^2, \sigma_1^2) \ge 0$ depending only on $b$, $\zeta$, $\sigma_0^2$, and $\sigma_1^2$. In this case, the confidence interval $\hat\tau - {\mathcal I}(b, \zeta)$ in Definition (ref) is given by $[\hat\tau - c(b,\zeta, \sigma_0^2, \sigma_1^2), \hat\tau + c(b,\zeta, \sigma_0^2, \sigma_1^2)]$.

Therefore, for given $(\widehat\tau_0, \widehat\tau_1, \sigma_0^2, \sigma_1^2)$ and $\zeta$, we can regard the confidence interval $\hat\tau - {\mathcal I}(b, \zeta)$ in Definition (ref) as a function of the bias bound $b$. Thus, to address the central question of constructing a sequence of two-sided confidence intervals for $\tau$ across different levels of bias, we compute this confidence interval for a range of values of $b$: $\{\hat\tau - {\mathcal I}(b, \zeta)\}_{b \ge 0}$.

Once we construct the sequence of confidence intervals, a natural application is hypothesis testing about the parameter $\tau$. We focus on the two-sided test:

We generalize the discussion to one-sided tests, such as testing $\tau = 0$ versus $\tau > 0$ or testing $\tau \le 0$ versus $\tau > 0$, in Section (ref). It is straightforward to extend the two-sided test to hypotheses of the form $\tau = \tau^*$ versus $\tau \neq \tau^*$ for any given $\tau^* \in \mathbb{R}$. A common practice for the two-sided test involves constructing confidence intervals for $\tau$ and checking whether the null value $0$ lies within these intervals.

Under Assumption (ref), the width of the confidence interval increases with the bias bound $b$. This leads to a central question when using the combined estimator $\hat\tau$ with the confidence interval $\hat\tau - {\mathcal I}(b, \zeta)$:

By the monotonicity of ${\mathcal I}(b, \zeta)$ in $b$ (holding $(\widehat\tau_0, \widehat\tau_1)$ and $\zeta$ fixed) from Assumption (ref), we define the limiting interval as ${\mathcal I}(\infty, \zeta) = \lim_{b \to \infty} {\mathcal I}(b, \zeta)$. Then there may exist a critical value of $b$ at which the confidence interval $\hat\tau - {\mathcal I}(b, \zeta)$ contains the null value tested in (ref). We focus on cases where the critical value exists. Specifically, we consider the case when $0 \notin \hat\tau - {\mathcal I}(0, \zeta)$ but $0 \in \hat\tau - {\mathcal I}(\infty, \zeta)$. In this scenario, we define the b-value below.

By the monotonicity conditions in Assumption (ref), we have $0 \notin \hat\tau - {\mathcal I}(b, \zeta)$ for $0 \le b < b^*$ and $0 \in \hat\tau - {\mathcal I}(b, \zeta)$ for $b > b^*$. Thus, we reject the null hypothesis $\tau = 0$ when $b < b^*$ and fail to reject it for $b > b^*$. The b-value $b^*$ thus represents the maximum relative bias beyond which the null hypothesis can no longer be rejected. We propose to report $b^*$ defined in (ref) for given estimators $(\widehat\tau_0,\widehat\tau_1)$ and significance level $\zeta$.

Compared with the sensitivity analysis literature, the bias bound $b$ in our framework plays a role analogous to the sensitivity parameter. The corresponding confidence interval $\hat\tau - {\mathcal I}(b, \zeta)$ serves as a sensitivity curve, illustrating how inference changes as the bias varies. In this context, the b-value $b^*$ serves a role similar to key robustness metrics in prior work: it parallels the design sensitivity by rosenbaum2004design, the E-value by vanderweele2017sensitivity, and the robustness value by cinelli2020making.

The b-value $b^*$ provides a criterion for comparing competing strategies of confidence interval construction. Different choices of the combined estimator $\hat\tau$ and different formulations of ${\mathcal I}(b, \zeta)$ can lead to different values of $b^*$. We prefer procedures that yield a larger b-value $b^*$, since this indicates that the resulting confidence interval is more robust to potential bias in the biased estimator. This role of $b^*$ is analogous to that use of the design sensitivity in rosenbaum2004design, where it serves as a criterion for comparing different test statistics and matched designs in observational studies.

Point estimation: a review

Before diving into the details of the method for constructing confidence intervals and obtaining the b-value, we first review the existing approaches to point estimation. Specifically, we review three point estimators: the precision-weighted estimator, the pretest estimator, and the soft-thresholding estimator.

First, we recall the precision-weighted estimator: \[ \widehat\tau_{\rm PW} := \frac{\sigma_1^2}{\sigma_0^2 + \sigma_1^2} \widehat\tau_0 + \frac{\sigma_0^2}{\sigma_0^2 + \sigma_1^2} \widehat\tau_1 = \widehat\tau_0 + \frac{\sigma_0^2}{\sigma_0^2 + \sigma_1^2} (\widehat\tau_1 - \widehat\tau_0) = \widehat\tau_0 + \frac{\gamma}{1+\gamma} (\widehat\tau_1 - \widehat\tau_0). \] When $\Delta$ is known to be 0, the precision-weighted estimator is the maximum likelihood estimator and best linear unbiased estimator of $\tau$. Moreover, by the classical Le Cam asymptotic decision theory le1956asymptotic and classical results for the Gaussian shift model, it is asymptotically admissible and minimax under $L_2$ risk in general, i.e., ${\mathrm E}[(\hat\tau - \tau)^2]$, among all regular estimators under standard regularity conditions.

However, $\widehat\tau_{\rm PW}$ is not robust to bias: its risk is large when $\lvert \Delta \rvert$ is large. This motivates us to consider a combined estimator that performs nearly as well as $\widehat\tau_{\rm PW}$ when the bias is small, and is robust to unknown bias $\Delta$, ensuring that the worst-case risk $\sup_{\Delta \in \bR} {\mathrm E}_\Delta[(\hat\tau - \tau)^2]$ remains bounded. This motivates the following two estimators.

Second, we recall the pretest estimator, which involves incorporating a pretest for $\Delta = 0$ versus $\Delta \neq 0$, a procedure that is commonly used bancroft1944biases,wallace1977pretest, bancroft1977inference,giles1993pre. Under $\Delta = 0$, given the independence between $\widehat\tau_1$ and $\widehat\tau_0$, their difference follows: $\widehat\tau_1 - \widehat\tau_0 \sim N(0, \sigma^2)$, where $\sigma^2 = \sigma_0^2 + \sigma_1^2$. For a fixed significance level $\alpha \in (0,1)$, we consider the test statistic $(\widehat\tau_1 - \widehat\tau_0)/\sigma$. Let $$\mathcal{A} = \{\lvert \widehat\tau_1 - \widehat\tau_0 \rvert \le \sigma c_{\alpha/2}\}$$ denote the event that the pretest fails to reject the null hypothesis ($\Delta = 0$). If the pretest fails to reject the null hypothesis, i.e., $\lvert \widehat\tau_1 - \widehat\tau_0 \rvert \le \sigma c_{\alpha/2}$, the pretest estimator uses the precision-weighted estimator: \[ \widehat\tau_0 + \frac{\gamma}{1+\gamma} (\widehat\tau_1 - \widehat\tau_0). \] If the pretest rejects the null hypothesis, i.e., $\lvert \widehat\tau_1 - \widehat\tau_0 \rvert > \sigma c_{\alpha/2}$, the pretest estimator employs hard-thresholding by reverting to the unbiased estimator $\widehat\tau_0$. Combining the two cases, the pretest estimator is: \[ \widehat\tau_{\rm PT} = \widehat\tau_0 + \frac{\gamma}{1+\gamma} (\widehat\tau_1 - \widehat\tau_0) \ind\left(\mathcal{A}\right). \]

Third, we recall the soft-thresholding estimator, which ensures continuity at the pretest boundary ($\lvert \widehat\tau_1 - \widehat\tau_0 \rvert = \sigma c_{\alpha/2}$). If the pretest fails to reject the null hypothesis, the soft-thresholding estimator also uses the precision-weighted estimator. If the pretest rejects the null hypothesis, the soft-thresholding estimator employs soft-thresholding by: \[ \widehat\tau_0 + \frac{\gamma}{1+\gamma} \sigma c_{\alpha/2}\,\sign(\widehat\tau_1 - \widehat\tau_0). \] Combining the two cases, the soft-thresholding estimator is: \[ \widehat\tau_{\rm ST}= \widehat\tau_0 + \frac{\gamma}{1+\gamma} (\widehat\tau_1 - \widehat\tau_0)\ind\left(\mathcal{A}\right) + \frac{\gamma}{1+\gamma} \sigma c_{\alpha/2}\,\sign(\widehat\tau_1 - \widehat\tau_0)\ind\left(\mathcal{A}^\textup{c}\right). \]

While $\widehat\tau_{\rm PT}$ reduces to $\widehat\tau_0$ with probability nearly one when the bias is large, the risk of $\widehat\tau_{\rm PT}$ is much higher than that of $\widehat\tau_{\rm ST}$ when the bias is moderate. See bickel1983minimax and armstrong2025adapting for a comparison of worst-case risk between $\widehat\tau_{\rm PT}$ and $\widehat\tau_{\rm ST}$ and the superior performance of $\widehat\tau_{\rm ST}$.

In both $\widehat\tau_{\rm PT}$ and $\widehat\tau_{\rm ST}$, the role of $\alpha$ is different from its usual interpretation in hypothesis testing. Rather than a Type I error rate, $\alpha$ here acts as a tuning parameter that balances the bias and variance of the point estimator. In this sense, $\alpha$ simply indexes a family of estimators, much like the penalty level in penalized regression methods (e.g., the Lasso). Choosing $\alpha$ optimally depends on $\Delta$, which is unknown in practice. In our empirical studies, we set $\alpha = 0.05$ as a default, following the standard convention, and examine how the estimator’s performance varies across different bias levels.

Confidence intervals, hypothesis testing, and the b-value

Given the canonical point estimators introduced in Section (ref), we now discuss the problem of constructing confidence intervals, hypothesis testing, and the b-value.

Confidence interval based on the precision-weighted estimator

As a warm-up, we use the precision-weighted estimator $\widehat\tau_{\rm PW}$ to construct the sequence of confidence intervals and to illustrate how the confidence interval changes with respect to $b$. Based on the theory below, we do not recommend using the precision-weighted estimator and its corresponding confidence intervals in practice. Under Assumption (ref), we have $\widehat\tau_{\rm PW} \sim N(\tau + (1+\gamma)^{-1} \gamma \Delta, (1+\gamma)^{-1} \sigma_0^2)$. The following theorem provides the sequence of confidence intervals based on $\widehat\tau_{\rm PW}$ depending on $b$.

The $\hat{L}_{\rm PW}$ in Theorem (ref) corresponds to the $(1-\zeta)$ quantile of the folded normal distribution $\lvert N(\frac{\gamma}{\sqrt{1+\gamma}} b, 1) \rvert$. This distribution also arises in the regression literature, where worst-case bias is incorporated into the construction of confidence intervals armstrong2020bias, armstrong2022robust.

In the special case when $b=0$, we have $\hat{L}_{\rm PW}(0, \zeta, \gamma) = c_{\zeta/2}$, which yields the length of the confidence interval based on $\widehat\tau_{\rm PW}$ to be $2 c_{\zeta/2} (1+\gamma)^{-1/2} \sigma_0$. For comparison, recall that the length of the confidence interval solely based on $\widehat\tau_0$ is $2c_{\zeta/2} \sigma_0$. Thus the length of the confidence interval reduces by a factor of $(1+\gamma)^{-1/2}$ when using $\widehat\tau_{\rm PW}$ instead of $\widehat\tau_0$, which is small when $\gamma$ is large, i.e., $\sigma_0^2$ is much larger than $\sigma_1^2$. This indicates that by using a more precise but potentially biased estimator, if we have prior knowledge that the bias is small, then we can achieve a much shorter confidence interval. Given $\zeta$ and $\gamma$, the $\hat{L}_{\rm PW}(b, \zeta, \gamma)$ in Theorem (ref) increases as $b$ increases, and goes to infinity as $b \to \infty$. Therefore, the precision-weighted estimator is not robust to unknown bias since the confidence interval is not bounded as the bias diverges.

One interesting observation is that $\hat{L}_{\rm PW}$ depends on $\sigma_0^2$ and $\sigma_1^2$ only through $\gamma$, the variance ratio. This explains why the bias $\Delta$ is scaled by $\sigma_0$ in the definition of $b$, and the confidence interval is represented as $[\widehat\tau_{\rm PW} - \hat{L}_{\rm PW} (1+\gamma)^{-1/2} \sigma_0, \widehat\tau_{\rm PW} + \hat{L}_{\rm PW} (1+\gamma)^{-1/2} \sigma_0]$.

Confidence interval based on the pretest estimator

We now construct the confidence interval for the pretest estimator $\widehat\tau_{\rm PT}$. In the following lemma, we present the distribution of $\hat\tau_{\rm PT} - \tau$.

By Lemma (ref), the distribution of $\hat\tau_{\rm PT}$ is a mixture of a normal distribution and a truncated normal distribution. We define $\hat{L}_{\rm PT} = \hat{L}_{\rm PT}(b, \zeta, \sigma_0^2, \sigma_1^2, \alpha)$ as the smallest length such that the confidence interval $[\widehat\tau_{\rm PT} - \hat{L}_{\rm PT} (1+\gamma)^{-1/2} \sigma_0, \widehat\tau_{\rm PT} + \hat{L}_{\rm PT} (1+\gamma)^{-1/2} \sigma_0]$ achieves correct coverage for all $\Delta$ satisfying $\lvert \Delta/\sigma_0 \rvert \le b$, where we choose the same scale as in Theorem (ref) for comparison. By Definition (ref), we can formulate $\hat{L}_{\rm PT}$ as the optimization problem

However, $\hat{L}_{\rm PT}$ generally does not admit a closed-form expression. Moreover, direct computation of $\hat{L}_{\rm PT}$ based on (ref) is computationally challenging since for each $L \ge 0$, the optimization problem (ref) involves finding the infimum of ${\mathrm P}_\Delta(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ over all $\Delta$ satisfying $\lvert \Delta/\sigma_0 \rvert \le b$. Nevertheless, we show that computing $\hat{L}_{\rm PT}$ is tractable due to the following theorem:

We have two main observations from Theorem (ref). First, for any $L>0$, the coverage probability ${\mathrm P}_\Delta(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ is symmetric about $\Delta=0$. As a result, instead of minimizing over all $\Delta$ satisfying $\lvert \Delta/\sigma_0 \rvert \le b$ as in (ref), we can restrict attention to the interval $0 \le \Delta/\sigma_0 \le b$. Second, Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ as a function of $t$ and $L$. The first term in (ref) corresponds to the event that the pretest fails to reject the null hypothesis ($\Delta = 0$), in which case $\hat\tau_{\rm PT}$ reduces to the precision-weighted estimator. The second and third terms in (ref) integrate over the regions where the pretest rejects the null hypothesis, in which case $\hat\tau_{\rm PT}$ reduces to the unbiased estimator.

Although the probability ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ admits a closed form for any $L > 0$, the minimum of ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ over $0 \le t \le b$ is not explicitly available in closed form, and thus must be computed numerically. Since the function ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ is monotonically increasing with respect to $L$ for $L \ge 0$, $L^*_{\rm PT}$ can be computed numerically using, for example, the bisection method, combined with numerical computation of the minimum of ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm PT} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ with respect to $t$ for each candidate $L$.

Confidence interval based on the soft-thresholding estimator

We now construct the confidence interval for the soft-thresholding estimator $\widehat\tau_{\rm ST}$. In the following lemma, we present the distribution of $\hat\tau_{\rm ST} - \tau$.

By Lemma (ref), the distribution of $\hat\tau_{\rm ST}$ is a mixture of a normal distribution and a truncated normal distribution. Suppose $\lvert \Delta/\sigma_0 \rvert \leq b$ for some $b > 0$. We seek the shortest length $\hat{L}_{\rm ST} = \hat{L}_{\rm ST}(b, \zeta, \sigma_0^2, \sigma_1^2, \alpha)$ such that the confidence interval $[\widehat\tau_{\rm ST}- \hat{L}_{\rm ST} (1+\gamma)^{-1/2} \sigma_0, \widehat\tau_{\rm ST}+ \hat{L}_{\rm ST} (1+\gamma)^{-1/2} \sigma_0]$ achieves correct coverage uniformly over all $\Delta$ with $\lvert \Delta/\sigma_0 \rvert \le b$. By Definition (ref), we can formulate $\hat{L}_{\rm ST}$ as the optimization problem:

However, $\hat{L}_{\rm ST}$ generally does not admit a closed-form expression. Moreover, direct computation of $\hat{L}_{\rm ST}$ based on (ref) is computationally challenging since for each $L \ge 0$, the optimization problem (ref) involves finding the infimum of ${\mathrm P}_\Delta(\lvert \widehat\tau_{\rm ST} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ over all $\Delta$ satisfying $\lvert \Delta/\sigma_0 \rvert \le b$. Nevertheless, we show that $\hat{L}_{\rm ST}$ can be computed efficiently due to the following theorem:

We have two main observations from Theorem (ref). First, for any $L>0$, the coverage probability ${\mathrm P}_\Delta(\lvert \widehat\tau_{\rm ST} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ is symmetric about $\Delta=0$, unlike the pretest estimator, monotonically decreasing in $\lvert \Delta \rvert$. This monotonicity implies that the worst-case coverage over the bias $\Delta$ satisfying $\lvert \Delta/\sigma_0 \rvert \le b$ as in (ref) is always attained at the boundary $\Delta/\sigma_0 = b$. The monotonicity of the coverage probability ${\mathrm P}_\Delta(\lvert \widehat\tau_{\rm ST} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ makes the computation of $\hat{L}_{\rm ST}$ more efficient than that of $\hat{L}_{\rm PT}$. Second, Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\Delta/\sigma_0 = t}(\lvert \widehat\tau_{\rm ST} - \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ as a function of $t$ and $L$. The first term in (ref) corresponds to the event that the pretest fails to reject the null hypothesis ($\Delta = 0$), in which case $\hat\tau_{\rm ST}$ reduces to the precision-weighted estimator. The second and third terms in (ref) integrate over the regions where the pretest rejects the null hypothesis, in which case $\hat\tau_{\rm ST}$ reduces to the unbiased estimator with a constant shift to ensure continuity at the pretest boundary.

Since ${\mathrm P}_{\Delta/\sigma_0 = b}(\lvert \widehat\tau_{\rm ST}- \tau \rvert \le L (1+\gamma)^{-1/2} \sigma_0)$ is monotonically increasing in $L$, $\hat{L}_{\rm ST}$ can be efficiently computed by, for example, the bisection method.

Comparison of the point estimators and confidence intervals

Compared with $\widehat\tau_{\rm PW}$, $\widehat\tau_{\rm ST}$ has bounded worst-case risk. The mean squared error of $\widehat\tau_{\rm ST}$ remains bounded regardless of the magnitude of the bias, whereas the mean squared error of $\widehat\tau_{\rm PW}$ grows without bound as the bias increases. Thus, $\widehat\tau_{\rm ST}$ offers a more balanced compromise between efficiency and robustness: it retains efficiency comparable to $\widehat\tau_{\rm PW}$ when the bias is small, yet its performance comparable to the unbiased estimator when the bias is large.

Compared with $\widehat\tau_{\rm PT}$, $\widehat\tau_{\rm ST}$ enjoys a desirable monotonicity property: for any $L > 0$, the probability ${\mathrm P}_{\Delta}(\lvert \widehat\tau_{\rm ST}- \tau \rvert \le L)$ decreases monotonically in $\lvert \Delta \rvert$ (Theorem (ref)). This behavior parallels the result in bickel1983minimax, which shows that the mean squared error of $\widehat\tau_{\rm ST}$ increases monotonically in $\lvert \Delta \rvert$. This monotonicity of $\widehat\tau_{\rm ST}$ allows us to compute the confidence interval based on $\widehat\tau_{\rm ST}$ more efficiently than $\widehat\tau_{\rm PT}$, as shown in Theorem (ref).

Figure (ref) compares the confidence intervals based on $\widehat\tau_{\rm ST}$, $\widehat\tau_{\rm PT}$, $\widehat\tau_{\rm PW}$, and the unbiased estimator $\widehat\tau_0$, as the bias bound $b$ varies. We fix $\sigma_0^2 = 1$ and $\widehat\tau_0 = 1$. Let $\zeta = 0.05$ and $\alpha = 0.05$. We set $\widehat\tau_1 = 2$ and examine two scenarios: $\gamma = 10$ and $\gamma = 100$. We observe that when the bias bound $b$ is small, the confidence intervals based on $\widehat\tau_{\rm ST}$, $\widehat\tau_{\rm PT}$, and $\widehat\tau_{\rm PW}$ are all shorter than that based on $\widehat\tau_0$. Importantly, when the bias bound $b$ is small, the confidence interval based on $\widehat\tau_{\rm ST}$ is shorter than $\widehat\tau_{\rm PT}$ and is comparable to $\widehat\tau_{\rm PW}$. Furthermore, the confidence interval length based on $\widehat\tau_{\rm ST}$ remains bounded, whereas that based on $\widehat\tau_{\rm PW}$ grows without bound as $b$ increases. Therefore, $\widehat\tau_{\rm ST}$ combines the efficiency gains of precision-weighted methods when bias is small with robustness to large biases, making it a superior choice for inference in practice.

Generalization to the dependent case

So far we assume independent unbiased and biased estimators. In this section, we generalize the above discussion to the case where $\widehat\tau_0$ and $\widehat\tau_1$ are jointly normal with known correlation $\rho$. We assume $\rho \sigma_1 \neq \sigma_0$, which trivially holds when $\sigma_0^2 > \sigma_1^2$, or equivalently $\gamma > 1$. The key is to construct a biased estimator which is independent of $\widehat\tau_0$ without losing information of $\widehat\tau_1$. To achieve that, we define the reparametrization

Then, $\widehat\tau_0$ and $\widehat\tau_1'$ are independent, with $\widehat\tau_1' \sim N(\tau + \Delta', \sigma_1'^2)$, where \[ \Delta' = \frac{\Delta}{1 - \rho \sigma_1/\sigma_0}, \quad \sigma_1'^2 = \frac{(1-\rho^2)\sigma_1^2}{(1-\rho \sigma_1/\sigma_0)^2}. \] Thus, the problem reduces to the independent case considered above.

In the reparametrization in (ref), we need to know how the transformation rescales the bias and how to interpret the relative bias. We can compute the confidence intervals and the b-value based on the reparametrization with independent unbiased estimator $\widehat\tau_0$ and transformed biased estimator $\widehat\tau_1'$. Let the confidence intervals and the b-value in the transformed problem be $\widehat\tau - {\mathcal I}'(b, \zeta)$ and $b^{*’}$, respectively. Then the confidence intervals for the original problem are $\widehat\tau - {\mathcal I}(b, \zeta) = \widehat\tau - {\mathcal I}'(b/\lvert 1 - \rho \sigma_1/\sigma_0 \rvert, \zeta)$ and the b-value is $b^* = \lvert 1 - \rho \sigma_1/\sigma_0 \rvert b^{*'}$, respectively. We relegate the technical details to Appendix (ref).

Generalization to the multivariate case

In this section, we extend our framework to the multivariate case. Many applications involve vector-valued parameters rather than scalars. Consider two leading examples in causal inference. First, when a treatment has multiple levels, the parameter of interest is a vector of treatment effects across those levels. This setting also includes factorial designs, where researchers aim to estimate factorial effects jointly. Second, when the population is partitioned into subgroups, the focus is often on subgroup treatment effects, yielding a vector of conditional average treatment effects indexed by the subgroup variable schwartz2025harmonized. Beyond causal inference, similar issues arise in regression settings. For instance, when combining OLS and IV estimators, the target parameter can be a vector of regression coefficients. Therefore, generalization to the multivariate case is essential for applications where researchers must make inferences about multiple parameters.

Setup

We consider the following multivariate setting:

As in the univariate case, we present our analysis in the exact normality setting which asymptotically captures the essential structure of inference as in wald1943tests and le1956asymptotic. Our goal is to construct confidence regions for $\bm{\tau}$ at different bias levels. Although we assume independent unbiased and biased estimators, the extension to the dependent case is straightforward; see Section (ref). Since $\bm{\widehat\tau}_0$ is unbiased for $\bm{\tau}$, given significance level $\zeta \in (0,1)$, a natural confidence region based solely on $\bm{\widehat\tau}_0$ is the ellipsoid:

where $\chi^2_{d,1-\zeta}$ is the $(1-\zeta)$ upper quantile of the chi-squared distribution with $d$ degrees of freedom.

Then we consider combining the two estimators $\bm{\widehat\tau}_0$ and $\bm{\widehat\tau}_1$, and we consider a generic combined estimator $\bm{\hat\tau} = \bm{\hat\tau}(\bm{\widehat\tau}_0, \bm{\widehat\tau}_1, \bm{\Sigma}_0, \bm{\Sigma}_1)$. We assume that $\lvert [\bm{\Sigma}^{-1/2} \bm{\Delta}]_j \rvert \le b_j$ for some $b_j \ge 0$ for all $j=1,2,\ldots,d$ and study how the confidence region changes with the maximum relative bias vector $\bm{b} = (b_1, b_2, \ldots, b_d)^\top$. Here $\bm{\Sigma} = \bm{\Sigma}(\bm{\Sigma}_0, \bm{\Sigma}_1) \in \bR^{d \times d}$ can be any fixed positive definite scaling matrix that depends on $\bm{\Sigma}_0$ and $\bm{\Sigma}_1$, which corresponds to $\sigma_0^2$ we used in the univariate case. For a given bias vector $\bm{b}$, different choices of $\bm{\Sigma}$ correspond to different regions for the bias $\Delta$. In this paper, we do not discuss how to choose the scaling matrix $\bm{\Sigma}$ optimally. For simplicity, one may take $\bm{\Sigma} = \bm{I}_d$, where $\bm{I}_d$ is the $d \times d$ identity matrix, or set $\bm{\Sigma} = \bm{\Sigma}_0$ as in the univariate case. Analogous to Definition (ref) for the univariate case, the confidence region is defined below.

Then we introduce the monotonicity conditions on the region $\bm{{\mathcal I}}(\bm{b},\zeta)$ below, which generalizes Assumption (ref) to the multivariate case.

Assume $\bm{{\mathcal I}}(\bm{b},\zeta)$ satisfies the monotonicity conditions of Assumption (ref). A common family of such regions is fixed-length centered ellipsoids: $\bm{{\mathcal I}}(\bm{b}, \zeta) = \{\bm{h} \in \bR^d: \bm{h}^\top \bm{A}^{-1} \bm{h} \le c(\bm{b}, \zeta, \bm{\Sigma}_0, \bm{\Sigma}_1)\}$ with some constant $c(\bm{b}, \zeta, \bm{\Sigma}_0, \bm{\Sigma}_1) \ge 0$ and a positive definite matrix $\bm{A}$. In this case, the confidence region $\bm{\hat\tau}-\bm{{\mathcal I}}(\bm{b},\zeta)$ in Definition (ref) is given by $\{\bm{\tau} \in \bR^d: (\bm{\widehat\tau} - \bm{\tau})^\top \bm{A}^{-1} (\bm{\widehat\tau} - \bm{\tau}) \le c(\bm{b}, \zeta, \bm{\Sigma}_0, \bm{\Sigma}_1)\}$.

Then we define the multivariate b-value below.

Definition (ref) extends the univariate notion of the b-value (Definition (ref)) to the multivariate setting. In higher dimensions $d > 1$, the set $\left\{\bm{b} \ge \bm{0}: \bm{0} \in \bm{\hat\tau} - \bm{{\mathcal I}}(\bm{b}, \zeta)\right\}$ is a convex region in the first quadrant. The b-value is defined as the boundary of this region, $\partial\{\bm{b} \ge \bm{0}: \bm{0} \in \bm{\hat\tau} - \bm{{\mathcal I}}(\bm{b}, \zeta)\}$, which is a $(d-1)$-dimensional surface, e.g., a curve when $d=2$, in the first quadrant. When $d=1$, this boundary reduces to the single left endpoint $\inf\left\{b \ge 0: 0 \in \hat\tau - {\mathcal I}(b, \zeta)\right\}$, which is exactly Definition (ref). By the monotonicity conditions in Assumption (ref), the geometry of this boundary yields a natural decision rule: we reject the null hypothesis $\bm{\tau} = \bm{0}$ if the bias bound vector $\bm{b}$ lies to the left of the b-value surface, and fail to reject it if $\bm{b}$ lies on or to the right of the b-value surface. To compute and visualize the multivariate b-value, we can plot the estimators and their confidence regions, and then we are ready to read the b-value for each estimator given the significance level $\zeta$.

Point estimation

First, we recall the precision-weighted estimator defined as \[ \bm{\widehat\tau}_{\rm PW} = (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1}(\bm{\Sigma}_0^{-1} \bm{\widehat\tau}_0 + \bm{\Sigma}_1^{-1} \bm{\widehat\tau}_1) = \bm{\widehat\tau}_0 + (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1} \bm{\Sigma}_1^{-1} (\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0). \] The precision-weighted estimator is the maximum likelihood estimator of $\bm{\tau}$ if $\bm{\Delta}$ is known to be all zeros. However, its risk is large when $\lVert \bm{\Delta} \rVert_2$ is large. Thus, we seek a combined estimator that performs nearly as well as $\bm{\widehat\tau}_{\rm PW}$ when the bias is small, and is robust to unknown bias $\bm{\Delta}$, ensuring that the maximum risk $\sup_{\bm{\Delta} \in \bR^d} {\mathrm E}_{\bm{\Delta}}[\lVert \bm{\widehat\tau}_{\rm PW} - \bm{\tau} \rVert_2^2]$ remains bounded. To develop estimators that perform well under zero bias yet remain robust to unknown bias, various estimators have been proposed in the empirical Bayes and shrinkage estimation literature j1981estimation,bickel1984parametric,green1991james,green2005improved,rosenman2023combining,rosenman2023empirical. In this paper, we focus on the generic pretest estimator and the generic soft-thresholding estimator, as discussed below.

Second, we recall the pretest estimator, which involves incorporating a pretest for $\bm{\Delta} = \bm{0}$ versus $\bm{\Delta} \neq \bm{0}$. Under $\bm{\Delta} = \bm{0}$, given the independence between $\bm{\widehat\tau}_1$ and $\bm{\widehat\tau}_0$, their difference follows: $\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0 \sim N(\bm{0}, \bm{\Sigma}_0 + \bm{\Sigma}_1)$. We consider the test statistic $\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2$ and the critical value $q \ge 0$. The critical value $q$ here plays a similar role as the significance level $\alpha$ in Section (ref). Let $$\mathcal{A} = \{\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2 \le q\}$$ denote the event that the pretest fails to reject the null hypothesis ($\bm{\Delta} = \bm{0}$). If the pretest fails to reject the null hypothesis, i.e., $\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2 \le q$, the pretest estimator uses the precision-weighted estimator: \[ \bm{\widehat\tau}_0 + (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1} \bm{\Sigma}_1^{-1} (\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0). \] If the pretest rejects the null hypothesis, i.e., $\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2 > q$, the pretest estimator employs hard-thresholding by reverting to the unbiased estimator $\bm{\widehat\tau}_0$. Combining the two cases, the pretest estimator is: \[ \bm{\widehat\tau}_{\rm PT} = \bm{\widehat\tau}_0 + (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1} \bm{\Sigma}_1^{-1} (\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \ind\left(\mathcal{A}\right). \]

Third, we recall the soft-thresholding estimator, which ensures continuity at the pretest boundary ($\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2 = q$). If the pretest fails to reject the null hypothesis, the soft-thresholding estimator also uses the precision-weighted estimator. If the pretest rejects the null hypothesis, the soft-thresholding estimator employs soft-thresholding by: \[ \bm{\widehat\tau}_0 + h_q^*(\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2)(\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1}\bm{\Sigma}_1^{-1}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0), \] where $h_q^*(\cdot): [q,\infty) \to [0,1]$ is a non-increasing function with $h_q^*(q)=1$. Combining the two cases, the soft-thresholding estimator is: \[ \bm{\hat\tau}_{\rm ST} = \bm{\widehat\tau}_0 + h_q(\lVert (\bm{\Sigma}_0+ \bm{\Sigma}_1)^{-1/2}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0) \rVert_2^2)(\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1}\bm{\Sigma}_1^{-1}(\bm{\widehat\tau}_1 - \bm{\widehat\tau}_0), \] where $h_q(\cdot): [0,\infty) \to [0,1]$ is defined as $h_q(r) = \ind(0 \le r \le q) + h_q^*(r) \ind(r > q)$. Here both $q$ and $h_q^*(\cdot)$ can depend on $\bm{\Sigma}_0$ and $\bm{\Sigma}_1$.

The generic soft-thresholding estimator $\bm{\hat\tau}_{\rm ST}$ contains many classical estimators as special cases, such as the estimator studied in j1981estimation and bickel1984parametric, which generalizes the univariate soft-thresholding estimator $\hat\tau_{\rm ST}$ to the multivariate setting. We provide a detailed discussion of the relationship between the generic soft-thresholding estimator and the estimator in j1981estimation, bickel1984parametric in Appendix (ref). In this paper, we focus on general $\bm{\Sigma}_0$, $\bm{\Sigma}_1$, and $h_q^*(\cdot)$.

Confidence intervals

First, we construct the confidence regions based on the precision-weighted estimator $\bm{\widehat\tau}_{\rm PW}$. Under Assumption (ref), we have $\bm{\widehat\tau}_{\rm PW} \sim N(\bm{\tau} + (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1} \bm{\Sigma}_1^{-1} \bm{\Delta}, (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})^{-1} )$. By using the ellipsoidal confidence region under the Mahalanobis distance and the covariance matrix of $\bm{\widehat\tau}_{\rm PW}$, we consider the confidence region based on $\bm{\widehat\tau}_{\rm PW}$ with the following form: \[ \{\bm{\tau}\in\mathbb{R}^d:(\bm{\widehat\tau}_{\rm PW}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm PW}-\bm{\tau})\le M\}, \] for some $M \ge 0$. The following theorem explicitly characterizes the confidence region based on $\bm{\widehat\tau}_{\rm PW}$ as a function of bias bound $\bm{b}$.

Theorem (ref) extends the univariate result in Theorem (ref) to the multivariate setting. In the univariate case, taking $\bm{\Sigma} = \sigma_0^2$ and $\bm{b} = b$ reduces the multivariate confidence region in Theorem (ref) exactly to the confidence interval in Theorem (ref). In the univariate case, the bias constraint $\{\Delta: \lvert \Delta/\sigma_0 \rvert \le b\}$ has only two boundary points, and the absolute bias of $\widehat\tau_{\rm PW}$ is the same at both endpoints. By contrast, the multivariate bias constraint $\{\bm{\Delta}: \lvert \bm{\Sigma}^{-1/2} \bm{\Delta}\rvert \le \bm{b}\}$ is a convex hyperrectangle with $2^d$ vertices, and the maximum $L_2$ norm of the bias of $\bm{\widehat\tau}_{\rm PW}$ may occur at any of these vertices. Theorem (ref) therefore characterizes the maximum possible bias of $\bm{\widehat\tau}_{\rm PW}$ over the entire bias constraint.

In the special case when $\bm{b}=0$, we have $\hat{M}_{\rm PW} = \chi^2_{d, 1-\zeta}$. Compared with the ellipsoid confidence region based on $\bm{\widehat\tau}_0$, the scaling matrix in the definition of ellipsoid confidence region reduces from $\bm{\Sigma}_0^{-1}$ to $\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1} = (\bm{\Sigma}_1^{-1} \bm{\Sigma}_0 + \bm{I}_d)\bm{\Sigma}_0^{-1}$, analogous to the efficiency gain in the univariate case, where $\bm{\Sigma}_1^{-1} \bm{\Sigma}_0$ can be regarded as the variance ratio $\gamma$ in the univariate case. As $\bm{b}$ increases, $\hat{M}_{\rm PW}$ increases, and goes to infinity as $\bm{b} \to \infty$.

Second, we construct the confidence regions based on the generic pretest estimator $\bm{\hat\tau}_{\rm PT}$. Similar to the confidence regions based on the precision-weighted estimator, we consider the confidence regions based on $\bm{\hat\tau}_{\rm PT}$ with the following form: \[ \{\bm{\tau}\in\mathbb{R}^d:(\bm{\widehat\tau}_{\rm PT}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm PT}-\bm{\tau})\le M\}, \] for some $M \ge 0$. Here, we focus on the ellipsoid confidence region defined using the same scaling matrix, $\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1}$, as that used in the precision-weighted estimator. This choice allows for a direct comparison with the confidence region based on the precision-weighted estimator through the upper bound $M$ in the confidence region. To ensure the confidence regions satisfy Definition (ref), we need to find the minimal value $\hat{M}_{\rm PT} = \hat{M}_{\rm PT}(\bm{b}, \zeta, \bm{\Sigma}_0, \bm{\Sigma}_1, \bm{\Sigma})$ such that the confidence region \[ \{\bm{\tau}\in\mathbb{R}^d:(\bm{\widehat\tau}_{\rm PT}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm PT}-\bm{\tau})\le \hat{M}_{\rm PT}\} \] achieves correct coverage for all $\bm{\Delta}$ satisfying $\lvert [\bm{\Sigma}^{-1/2} \bm{\Delta}]_j \rvert \le b_j$ for all $j=1,2,\ldots,d$. We can formulate $\hat{M}_{\rm PT}$ as the optimization problem

We show how to compute $\hat{M}_{\rm PT}$ in the following theorem.

Theorem (ref) extends the univariate result in Theorem (ref) to the multivariate setting. Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\bm{\Sigma}^{-1/2}\bm{\Delta} = \bm{t}}((\bm{\widehat\tau}_{\rm PT}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm PT}-\bm{\tau})\le M)$ as a function of $M$ and $\bm{t}$. The first term in (ref) corresponds to the event that the pretest fails to reject the null hypothesis ($\bm{\Delta} = \bm{0}$), in which case $\bm{\widehat\tau}_{\rm PT}$ reduces to the precision-weighted estimator. The second term in (ref) integrates over the regions where the pretest rejects the null hypothesis, in which case $\bm{\widehat\tau}_{\rm PT}$ reduces to the unbiased estimator.

Third, we construct the confidence regions based on the generic soft-thresholding estimator $\bm{\hat\tau}_{\rm ST}$. Similar to the confidence regions based on the precision-weighted estimator, we consider the confidence regions based on $\bm{\hat\tau}_{\rm ST}$ with the following form: \[ \{\bm{\tau}\in\mathbb{R}^d:(\bm{\widehat\tau}_{\rm ST}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm ST}-\bm{\tau})\le M\}, \] for some $M \ge 0$. To ensure the confidence regions satisfy Definition (ref), we need to find the minimal value $\hat{M}_{\rm ST} = \hat{M}_{\rm ST}(\bm{b}, \zeta, \bm{\Sigma}_0, \bm{\Sigma}_1, \bm{\Sigma})$ such that the confidence region \[ \{\bm{\tau}\in\mathbb{R}^d:(\bm{\widehat\tau}_{\rm ST}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm ST}-\bm{\tau})\le \hat{M}_{\rm ST}\} \] achieves correct coverage for all $\bm{\Delta}$ satisfying $\lvert [\bm{\Sigma}^{-1/2} \bm{\Delta}]_j \rvert \le b_j$ for all $j=1,2,\ldots,d$. We can formulate $\hat{M}_{\rm ST}$ as the optimization problem

We show that $\hat{M}_{\rm ST}$ can be computed efficiently in the following theorem.

Theorem (ref) extends the univariate result in Theorem (ref) to the multivariate setting. Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\bm{\Sigma}^{-1/2}\bm{\Delta} = \bm{t}}((\bm{\widehat\tau}_{\rm ST}-\bm{\tau})^\top (\bm{\Sigma}_0^{-1} + \bm{\Sigma}_1^{-1})(\bm{\widehat\tau}_{\rm ST}-\bm{\tau})\le M)$ as a function of $M$ and $\bm{t}$. The first term in (ref) corresponds to the event that the pretest fails to reject the null hypothesis ($\bm{\Delta} = \bm{0}$), in which case $\bm{\widehat\tau}_{\rm ST}$ reduces to the precision-weighted estimator. The second term in (ref) integrates over the regions where the pretest rejects the null hypothesis, in which case $\bm{\widehat\tau}_{\rm ST}$ reduces to the unbiased estimator with a nonconstant shift to ensure continuity at the pretest boundary. As in the univariate case, the monotonicity of the coverage probability makes the computation of $\hat{M}_{\rm ST}$ more efficient than that of $\hat{M}_{\rm PT}$. For any $M>0$, the coverage probability is minimized at the boundary of the bias region, namely at the vertices $\bm{\Delta}$ satisfying $\bm{\Sigma}^{-1/2} \bm{\Delta} = \bm{b} \odot \bm{s}$ for some $\bm{s} \in \{-1,1\}^d$.

Using Theorems (ref)--(ref), we can construct confidence regions based on $\bm{\hat\tau}_{\rm PW}$, $\bm{\hat\tau}_{\rm PT}$, and $\bm{\hat\tau}_{\rm ST}$ for various bias bounds $\bm{b}$ and significance levels $\zeta$.

Generalization to multiple estimators

In this section, we extend our framework to the multiple estimators setting. In empirical research, analysts often face the challenge of integrating evidence from datasets of heterogeneous quality. A leading example in causal inference is to combine an RCT with multiple observational studies. While the RCT provides unbiased but often noisy estimates, observational studies can offer much larger sample sizes but are subject to hidden bias due to unmeasured confounding. Beyond causal inference, similar problems occur when synthesizing results from multiple registries, surveys, or administrative databases, where no single source is sufficient on its own. We present our theory below.

Setup

We consider the following multiple estimators setting with one unbiased estimator and $K$ potentially biased estimators. We can generalize to the case where there are multiple unbiased estimators estimating the same parameter. In that case, we can first combine all the unbiased estimators into a single unbiased estimator with smaller variance using the precision-weighted estimator. Then we can apply our current framework to combine this aggregated unbiased estimator with multiple biased estimators.

Let $\gamma_j = \sigma_0^2/\sigma_j^2$ be the variance ratio between the unbiased estimator and the $j$-th biased estimator. Let $\bm{\gamma} = (\gamma_1,\ldots,\gamma_K)^\top$ be the vector of variance ratios.

We consider a generic combined estimator $\hat\tau = \hat\tau(\widehat\tau_0, \widehat\tau_1, \ldots, \widehat\tau_K, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2)$. We assume that $\lvert \Delta_j/\sigma_0 \rvert \le b_j$ for some $b_j > 0$ for all $j = 1,\ldots,K$ and study how the confidence interval changes with the maximum relative bias vector $\bm{b} = (b_1,\ldots,b_K)$. Let $\bm{\Delta} = (\Delta_1,\ldots,\Delta_K)$ be the vector of unknown biases. Analogous to Definition (ref) for the univariate case, the confidence interval is defined below.

Then we introduce the monotonicity conditions on the interval ${\mathcal I}(\bm{b}, \zeta)$ below, which generalizes Assumption (ref) to the multiple estimators setting.

Assume ${\mathcal I}(\bm{b},\zeta)$ satisfies the monotonicity conditions of Assumption (ref). A common family of such intervals is fixed-length centered intervals: ${\mathcal I}(\bm{b}, \zeta) = [-c(\bm{b}, \zeta, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2), c(\bm{b}, \zeta, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2)]$ with some constant $c(\bm{b}, \zeta, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2) \ge 0$. In this case, the confidence interval $\hat\tau - {\mathcal I}(\bm{b}, \zeta)$ in Definition (ref) is given by $[\hat\tau - c(\bm{b}, \zeta, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2), \hat\tau + c(\bm{b}, \zeta, \sigma_0^2, \sigma_1^2, \ldots, \sigma_K^2)]$.

Then we define the b-value with multiple estimators below.

Definition (ref) extends the univariate notion of the b-value (Definition (ref)) to the multiple estimators setting. As in the multivariate b-value defined in Definition (ref), in higher dimensions $K > 1$, the set $\left\{\bm{b} \ge \bm{0}: \bm{0} \in \bm{\hat\tau} - \bm{{\mathcal I}}(\bm{b}, \zeta)\right\}$ is a convex region in the first quadrant, and the b-value is the boundary of this region, which is a $(K-1)$-dimensional surface. See Definition (ref) and the subsequent discussion for more details.

Point estimation

First, we recall the precision-weighted estimator defined as \[ \widehat\tau_{\rm PW} = \frac{\sigma_0^{-2}}{\sigma_0^{-2} + \sum_{\ell=1}^K \sigma_\ell^{-2}} \widehat\tau_0 + \sum_{j=1}^K \frac{\sigma_j^{-2}}{\sigma_0^{-2} + \sum_{\ell=1}^K \sigma_\ell^{-2}} \widehat\tau_j = \widehat\tau_0 + \sum_{j=1}^K \frac{\gamma_j}{1+\lVert \bm{\gamma} \rVert_1} (\widehat\tau_j - \widehat\tau_0). \] If $\Delta_j=0$ for all $j=1,\ldots,K$, the precision-weighted estimator is the maximum likelihood estimator and best linear unbiased estimator of $\tau$. Moreover, by the classical Wald–Le Cam asymptotic decision theory wald1949statistical, le1956asymptotic, it is asymptotically admissible and minimax under the $L_2$ risk, i.e., ${\mathrm E}[(\hat\tau - \tau)^2]$, among all regular estimators under standard regularity conditions.

Second, we introduce a pretest estimator. We incorporate the pretest for $\Delta_j = 0$ separately using $\widehat\tau_j$ and the unbiased estimator $\widehat\tau_0$, and combine those $\widehat\tau_j$'s that fail to reject the null. Here we assume the pretests share the same significance level $\alpha$ for simplicity. The generalization to different significance levels is straightforward. Since $\widehat\tau_j - \widehat\tau_0 \sim N(\Delta_j, \sigma_0^2 + \sigma_j^2)$ and $\sigma_0^2 + \sigma_j^2 = (1+\gamma_j^{-1}) \sigma_0^2$, let $$\mathcal{A}_j = \{\lvert \widehat\tau_j - \widehat\tau_0 \rvert \le (1+\gamma_j^{-1})^{1/2}\sigma_0 c_{\alpha/2}\}$$ denote the event that the pretest fails to reject the null hypothesis ($\Delta_j=0$). The pretest estimator is:

Third, we introduce a soft-thresholding estimator, which ensures continuity at the pretest boundary. Instead of dropping a biased estimator entirely when the pretest rejects, soft-thresholding shifts the unbiased estimator by making the estimator continuous at the pretest boundary. The soft-thresholding estimator is:

By the definitions of the pretest estimator and the soft-thresholding estimator, when the pretests fail to reject for all $j=1,\ldots,K$, both the pretest estimator and the soft-thresholding estimator reduce to the precision-weighted estimator. When $K=1$, i.e., there is only one biased estimator, both pretest estimator and soft-thresholding estimator reduce to the estimators introduced in Section (ref). Similar to the discussion in Section (ref), the choice of $\alpha$ is a tuning parameter.

The pretest estimator and the soft-thresholding estimator considered above are not the only way to combine $\hat\tau_0$ and $\hat\tau_1,\ldots,\hat\tau_K$. For example, one could apply shrinkage or soft-thresholding directly between the precision-weighted estimator $\hat\tau_{\rm PW}$ and the unbiased estimator $\hat\tau_0$, rather than at the level of individual biased components. We focus on thresholding at the component level because it admits a transparent interpretation in terms of testing and controlling each bias component $\Delta_j$ separately, which aligns naturally with our sensitivity analysis perspective. Our framework can accommodate other shrinkage schemes in principle, but we do not pursue them here.

Confidence intervals

First, we construct the confidence intervals based on the precision-weighted estimator $\widehat\tau_{\rm PW}$. Under Assumption (ref), we have $\widehat\tau_{\rm PW} \sim N(\tau + (1+\lVert \bm{\gamma} \rVert_1)^{-1} \langle \bm{\gamma}, \bm{\Delta} \rangle, (1+\lVert \bm{\gamma} \rVert_1)^{-1} \sigma_0^2)$. The following theorem provides the confidence intervals based on $\widehat\tau_{\rm PW}$ as a function of bias bound $\bm{b}$.

Theorem (ref) extends the single estimator result in Theorem (ref) to the multiple estimators setting. In the single estimator case, the bias constraint $\{\Delta: \lvert \Delta/\sigma_0 \rvert \le b\}$ has only two boundary points, and the absolute bias of $\widehat\tau_{\rm PW}$ is the same at both endpoints. By contrast, the multiple estimators bias constraint $\{\bm{\Delta}: \lvert \bm{\Delta}/\sigma_0 \rvert \le \bm{b}\}$ is a hyperrectangle with $2^K$ vertices, and the maximum absolute bias of $\widehat\tau_{\rm PW}$ occurs only at $\bm{\Delta}/\sigma_0 = \pm \bm{b}$. Theorem (ref) therefore characterizes the maximum possible bias of $\widehat\tau_{\rm PW}$ over the entire bias constraint.

Second, we construct the confidence intervals based on the pretest estimator $\widehat\tau_{\rm PT}$. We seek the shortest length $\hat{L}_{\rm PT} = \hat{L}_{\rm PT}(\bm{b}, \zeta, \bm{\gamma}, \alpha)$ such that the confidence interval $[\widehat\tau_{\rm PT}- \hat{L}_{\rm PT} (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0, \widehat\tau_{\rm PT}+ \hat{L}_{\rm PT} (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0]$ achieves correct coverage uniformly over all $\bm{\Delta}$ with $\lvert \Delta_j/\sigma_0 \rvert \le b_j$ for all $j=1,\ldots,K$. We can formulate $\hat{L}_{\rm PT}$ as the optimization problem:

We show how to compute $\hat{L}_{\rm PT}$ in the following theorem.

Theorem (ref) extends the single estimator result in Theorem (ref) to the multiple estimators setting. Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\bm{\Delta}/\sigma_0 = \bm{t}}(\lvert \widehat\tau_{\rm PT}- \tau \rvert \le L (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0)$ as a function of $\bm{t}$ and $L$.

Third, we construct the confidence intervals based on the soft-thresholding estimator $\widehat\tau_{\rm ST}$. We seek the shortest length $\hat{L}_{\rm ST} = \hat{L}_{\rm ST}(\bm{b}, \zeta, \bm{\gamma}, \alpha)$ such that the confidence interval $[\widehat\tau_{\rm ST}- \hat{L}_{\rm ST} (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0, \widehat\tau_{\rm ST}+ \hat{L}_{\rm ST} (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0]$ achieves correct coverage uniformly over all $\bm{\Delta}$ with $\lvert \Delta_j/\sigma_0 \rvert \le b_j$ for all $j=1,\ldots,K$. We can formulate $\hat{L}_{\rm ST}$ as the optimization problem:

We show that $\hat{L}_{\rm ST}$ can be computed efficiently in the following theorem.

Theorem (ref) extends the single estimator result in Theorem (ref) to the multiple estimators setting. Theorem (ref) provides an explicit expression for the coverage probability ${\mathrm P}_{\bm{\Delta}/\sigma_0 = \bm{t}}(\lvert \widehat\tau_{\rm ST}- \tau \rvert \le L (1+\lVert \bm{\gamma} \rVert_1)^{-1/2} \sigma_0)$ as a function of $\bm{t}$ and $L$. As in the single estimator case, the monotonicity makes the computation of $\hat{L}_{\rm ST}$ more efficient than that of $\hat{L}_{\rm PT}$. For any $L>0$, the coverage probability is minimized at the vertices $\bm{\Delta}$ satisfying $\bm{\Delta}/\sigma_0 = \pm \bm{b}$.

Empirical studies

In this section, we use the example from angrist1991does to demonstrate how our framework works in practice, where the authors studied the effect of years of schooling on earnings. angrist1991does used quarter of birth as an instrument to obtain the IV estimate. They also reported the OLS estimate, which would be biased in the presence of endogeneity.

We consider the OLS estimator and the IV estimator. The OLS estimator of the return to schooling is potentially biased because schooling decisions are endogenous: unobserved factors may influence both schooling and earnings. By contrast, the IV estimator uses quarter of birth as a source of exogenous variation, yielding a consistent estimate of the causal effect under standard IV assumptions. We take the IV estimator as the unbiased estimator $\widehat\tau_0$ and the OLS estimator as the biased estimator $\widehat\tau_1$. Using these two estimators, we construct the three combined estimators: the precision-weighted estimator $\widehat\tau_{\rm PW}$, the pretest estimator $\widehat\tau_{\rm PT}$, and the soft-thresholding estimator $\widehat\tau_{\rm ST}$.

Figure (ref) visualizes the confidence intervals of the IV, OLS, and combined estimators against the maximum relative bias $\lvert \Delta/\sigma_0 \rvert \le b$.

In both the full sample (black and white men) and the subsample (black men), we observe that the standard error of the OLS estimator is much smaller than that of the IV estimator. Moreover, although we compute the OLS estimator and the IV estimator using the same sample, their correlation is very low. Therefore, the new biased estimator constructed following Section (ref) is nearly identical to the original biased estimator (OLS). The combined estimators behave as expected. When the bias bound is small, $\widehat\tau_{\rm PW}$, $\widehat\tau_{\rm PT}$, and $\widehat\tau_{\rm ST}$ all yield confidence intervals substantially shorter than that of the unbiased estimator (IV), reflecting efficiency gains from incorporating the more precise OLS estimator. Notably, the confidence interval of $\widehat\tau_{\rm ST}$ is nearly identical to that of $\widehat\tau_{\rm PW}$ and much shorter than that of $\widehat\tau_{\rm PT}$ when the bias bound is small. Even when the bias bound becomes large, the confidence intervals of $\widehat\tau_{\rm ST}$ and $\widehat\tau_{\rm PT}$ remain nearly the same in length. These behaviors arise because the OLS estimator is far more precise than the IV estimator.

However, there are differences between the full sample (black men and white men) and the subsample (black men). In the full sample, the sample size is very large, and the IV estimator is sufficiently precise that the null hypothesis of zero returns to schooling can be rejected using IV alone. In contrast, for the subsample of black men, the sample size is much smaller, leading to a much larger standard error for the IV estimator. As a result, the IV-based confidence interval is wide enough that the null hypothesis cannot be rejected at the 5% significance level. This contrast illustrates a common empirical challenge: when the unbiased estimator is noisy, inference based solely on it can be inconclusive, even when the overall dataset is large. In such settings, the combined estimators, particularly the precision-weighted estimator $\hat\tau_{\rm PW}$ and the soft-thresholding estimator $\hat\tau_{\rm ST}$, yield much tighter confidence intervals under small biases. Under small biases, these estimators provide sufficient evidence to reject the null hypothesis of zero returns to schooling, while still allowing researchers to explicitly assess sensitivity to potential bias.

Discussion

This paper develops a strategy to combine unbiased and biased estimators from a sensitivity analysis perspective. In particular, we construct a sequence of confidence intervals indexed by the magnitude of bias and propose the notion of the b-value to quantify the maximum relative bias so that combining estimators yields an insignificant result.

Acknowledgment

We thank Avi Feller and Liyang Sun for helpful comments. Lin was partially supported by the Two Sigma PhD Fellowship. Ding was supported by the U.S. National Science Foundation (1945136, 2514234).