Efficient estimation of average treatment effects with unmeasured confounding and proxies

Efficient Estimation of Average Treatment Effects with Unmeasured Confounding and Proxies

\affil[1]{School of Management and Economics, The Chinese University of Hong Kong, Shenzhen} \affil[2]{Department of Biostatistics & Medical Informatics, University of Wisconsin-Madison}

\thispagestyle{empty}

{\bf Key Words:} Data-driven method; generalized method of moments; proximal causal inference; semiparametric efficiency.

Introduction

The unconfoundedness assumption is a key condition for identifying treatment effect parameters and establishing the consistency of many popular estimators. This condition may not hold if some confounders are unmeasured and omitted from the empirical analysis. One approach to address the omitted variables problem is proximal causal inference, which assumes the availability of outcome and treatment confounding proxies MiaoGengTchetgenTchetgen2018, TchetgenTchetgenYingCuiShiEtAl2024. Such an approach has seen a wide range of applications ShiMiaoNelsonTchetgenTchetgen2020, KallusMaoUehara2022,DukesShpitserTchetgenTchetgen2023, YingMiaoShiTchetgenTchetgen2023,EgamiTchetgenTchetgen2024, GhassamiYangShpitserTchetgenTchetgen2024,QiMiaoZhang2024,QiuShiMiaoDobribanEtAl2024,Ying2024.

With the outcome and treatment confounding proxies, MiaoGengTchetgenTchetgen2018 showed that an outcome bridge function, analogous to the outcome regression function one would use if all confounders were observed, identifies the average treatment effect. A common approach is to parameterize the bridge function, estimate it using standard methods, and then estimate the average treatment effect (ATE) via a plug-in. For instance, MiaoShiLiTchetgenTchetgen2024 suggested a recursive generalized method of moments (RGMM) by transforming the integral equation to fixed-dimensional, user-specified moment restrictions. TchetgenTchetgenYingCuiShiEtAl2024 introduced a proximal g-computation method, which results in a simple proximal two-stage least squares procedure in the special case of linear working models. The proximal g-computation method is typically more efficient than RGMM, benefiting from an additional parametric restriction on the joint distribution of covariates, but is also more prone to bias if the distribution model is misspecified. These methods suffer from efficiency losses because estimating the bridge function may be inefficient, and the sequential procedure may fail to utilize all information BrownNewey1998,AiChen2012. As a remedy, CuiPuShiMiaoEtAl2023 proposed a doubly robust (DR) locally efficient approach that incorporates an additional treatment bridge function. They proved that their proximal DR estimator achieves the semiparametric local efficiency bound if both bridge functions are correctly specified and consistently estimated, even if they are not efficiently estimated. If one of the bridge functions is misspecified, the DR estimator fails to achieve local efficiency, and it could benefit from more efficient estimation of the bridge functions. Moreover, the DR estimator may not be the most efficient, even if both bridge functions are consistently estimated, as illustrated in (ref).

This paper's primary contribution is to develop a simple, data-driven method for efficiently estimating the average treatment effect. The key step is to transform the conditional moment restriction that defines the bridge function into an expanding set of unconditional moment restrictions via a sieve basis, and to estimate the bridge function and the ATE jointly. As the number of moments increases, these unconditional moment restrictions provide a good approximation for the unknown conditional moment restriction. We show that: (i) the proposed estimator for the bridge function achieves the semiparametric efficiency bound established in CuiPuShiMiaoEtAl2023, without requiring any modeling of the data distribution beyond the bridge function; (ii) the proposed estimator for ATE has an asymptotic variance that is never larger than that of the DR locally efficient estimator; and (iii) the proposed estimator outperforms the theoretically optimal plug-in estimator. Moreover, the proposed method is readily implemented using widely available software for GMM estimation Hansen1982. We also propose a data-driven procedure for selecting the number of moments.

The theoretical results of this paper are closely related to established theories in causal inference under the assumption of no unmeasured confounding and in the missing data literature under the missing at random assumption, which demonstrate that regression-based estimators are more efficient than doubly robust estimators when the outcome model is correctly specified ScharfsteinRotnitzkyRobins1999,BangRobins2005. We extend this principle to the proximal causal inference. Our work is also related to several strands of the existing literature. First, it connects to research on instrumental variable estimation and partial identification strategies AiChen2003,NeweyPowell2003,Abadie2003,KlineTamer2023. Second, it is related to recent advances that incorporate high-dimensional and machine learning methods into proximal estimation frameworks MastouriZhuGultchinKorbaEtAl2021,KompaBellamyKolokotronesRobinsEtAl2022. Third, it is linked to developments in causal graphical models for identifying effects in the presence of unmeasured confounding, including work on the front-door criterion Pearl1995,RichardsonEvansRobinsShpitser2023,GuoBenkeserNabi2023,BhattacharyaNabiShpitser2022,GuoNabi2024, which also leverage proxy-like mediators to achieve nonparametric identification under structural assumptions.

The rest of the paper is organized as follows. (ref) reviews the identification results of proximal causal inference and challenges in constructing an efficient estimator. (ref) describes the proposed method for efficient estimation. (ref) establishes large-sample properties of the proposed estimators and compares them with the existing methods. (ref) discusses the choice of tuning parameters for practical implementation. (ref) provides a simulation study and applies our method to reanalyze the SUPPORT dataset. (ref) offers a brief discussion. Technical proofs and additional results are provided in the supplementary material.

Basic framework

Let $A\in\{0,1\}$ denote the treatment variable, and $Y$ denote the observed outcome. Let $Y(a)$ denote the potential outcome if treatment $A=a$ is assigned. The observed outcome is $Y=Y(A)$. Our objective is to estimate the average treatment effect, $\tau_0\triangleq \mathbb{E}\{Y(1)-Y(0)\}$, in settings with unmeasured confounders. Let $U$ denote the unmeasured confounders that may affect both the treatment $A$ and outcome $Y$. Suppose there are three types of observable covariates $(X,W,Z)$: $X$ affect both $A$ and $Y$; $W$ are outcome-inducing confounding proxies which are related to $A$ only through $(X, U)$; and $Z$ are treatment-inducing confounding proxies which are associated with $Y$ only through $(X, U)$. See (ref) for a causal directed acyclic graph (DAG). We formalize the relationships among these variables in the following assumption.

(ref)(i)-(ii) are standard in the literature of causal inference. (ref)(iii) describes the nature of the two types of proxies: treatment proxies $Z$ cannot affect the outcome $Y$, and the outcome proxies $W$ cannot be affected by either the treatment $A$ or the treatment proxies $Z$, upon conditioning on the confounders $X$ and $U$. These conditions are critical for identification but untestable as they involve conditional independence statements given the unmeasured variable $U$. Justifying these conditions and selecting proxy variables requires subject-matter knowledge; however, they can be satisfied in applications through careful study design. For instance, by incorporating post-outcome variables in $Z$ and pre-treatment measurements of the outcome in $W$, we can reasonably expect (ref)(iii) to hold as the future cannot affect the past. MiaoShiLiTchetgenTchetgen2024 provided an example of evaluating the short-term effect of air pollution on elderly hospitalization using time series data, in which air pollution measurements taken after hospitalization are included in $Z$, and hospitalization measurements taken before air pollution are included in $W$. (ref)(iv) can be intuitively interpreted as a requirement for $Z$ to exhibit sufficient variability relative to that of the unmeasured confounder $U$. In the case of categorical $Z$ and $U$, it requires $Z$ with at least as many categories as $U$. In the continuous case, ChenChernozhukovLeeNewey2014 and Andrews2017 showed that if the dimension of $Z$ exceeds that of $U$, then under mild conditions, completeness holds generically, except for distributions in negligible sets. TchetgenTchetgenYingCuiShiEtAl2024 also suggest measuring a rich set of baseline characteristics to make the completeness assumption more plausible; however, this benefit should be balanced against the efficiency loss incurred by including irrelevant proxies and the increased sample size required for estimating high-dimensional nuisance functions. Completeness is also shown to hold in many parametric and semiparametric models, such as exponential families NeweyPowell2003 and location-scale families HuShiu2018. For a detailed discussion of completeness, see DHaultfoeuille2011. We also refer readers to TchetgenTchetgenYingCuiShiEtAl2024 for a detailed discussion of these assumptions, additional examples of proxies, and approaches for sensitivity analysis in practice.

Under (ref), MiaoGengTchetgenTchetgen2018 established the identification of ATE through an outcome-confounding bridge function $h(w,a,x)$, which solves

It is worth noting that the bridge function solves a Fredholm integral equation of the first kind, which is known to be an ill-posed inverse problem Kress1989. Then, ATE is identified by

Equation (ref) is referred to as the proximal g-formula in TchetgenTchetgenYingCuiShiEtAl2024.

We consider the estimation of $\tau_0$ using $N$ independent and identically distributed observations drawn from the joint distribution of $\bm{{O}}=(A,Y,W,X,Z)$. We follow TchetgenTchetgenYingCuiShiEtAl2024 by parameterizing the bridge function $h(W,A,X;\bm\gamma)$ with the parameter value $\bm\gamma_0\in\Gamma\subset\mathbb{R}^p$. Then, one can estimate $\bm\gamma_0$ by solving

where ${\bm m}(\cdot)$ is a user-specified vector-valued function, with dimension equal to that of $\bm\gamma$, and satisfies that $\mathbb{E}\left\{\nabla_\gamma h(W,A,X;\bm\gamma_0){\bm m}(Z,A,X)^{\mathrm{\scriptscriptstyle T}}\right\}$ is nonsingular. For instance, if $h(W,A,X;\bm\gamma)=(1,W,A,X)\bm\gamma$, one can choose ${\bm m}(Z,A,X)=(1,Z,A,X)^{\mathrm{\scriptscriptstyle T}}$. Under mild regularity conditions White1982, $\widehat{\bm\gamma}({\bm m})$ that solves (ref) is consistent and asymptotically normally distributed for any ${\bm m}$. However, the choice of ${\bm m}$ may affect the efficiency, namely, the asymptotic variance of $\widehat{\bm\gamma}({\bm m})$. Let $\nabla_{\bm x}\bm{\ell}({\bm x})=\partial\bm{\ell}({\bm x})^{\mathrm{\scriptscriptstyle T}}/\partial{\bm x}$ denote the gradient of a (vector-valued) function $\bm{\ell}$. CuiPuShiMiaoEtAl2023 showed that among all score functions of $\bm\gamma$, the efficient score is the one with

that accounts for the potential heteroskedasticity. Unfortunately, such ${\bm m}_\text{\textit{eff}}$ is impractical because it requires modeling complex features of the observed data distribution, which are difficult to capture accurately. MiaoShiLiTchetgenTchetgen2024 suggested a recursive GMM approach, allowing the dimension of ${\bm m}$ to exceed that of $\bm\gamma$. While their approach may offer some efficiency gains for $\widehat{\bm\gamma}({\bm m})$, the specific low-dimensional ${\bm m}$ does not necessarily lead to an efficient estimator for $\bm\gamma$, in contrast to the increasing moment conditions proposed in this paper, as elaborated later. An inefficient $\widehat{\bm\gamma}({\bm m})$ may lead to inefficiency of the plug-in estimator:

Furthermore, the plug-in estimators, including $\widehat{\tau}^\textnormal{plug-in}({\bm m}_\text{\textit{eff}})$, fail to account for the correlation between the score functions of $\bm\gamma_0$ and $\tau_0$, thereby introducing an additional source of efficiency loss BrownNewey1998,AiChen2012.

Estimation

To improve existing approaches, we propose estimating $\bm\gamma_0$ and $\tau_0$ jointly and using an increasing number of moment restrictions. Specifically, let ${\bm u}_{K}(z,a,x) = \big(u_{K1}(z,a,x),\ldots, u_{KK}(z,a,x)\big)^{\mathrm{\scriptscriptstyle T}}$ denote a vector of known basis functions (such as power series, splines, Fourier series, etc.) with dimension $K\in\mathbb{N}$, which provides approximation sieves that can approximate a large class of smooth functions arbitrarily well as $K\rightarrow\infty$. The selection of ${\bm u}_{K}(z,a,x)$ is discussed in (ref). Model (ref) implies the following unconditional moment restrictions of $\bm\gamma_0$:

Denote the joint score function $${\bm g}_K(\bm{{O}};\bm\gamma,\tau)=\left(

\right),$$ and ${\bm G}_K(\bm\gamma,\tau)= N^{-1}\sum_{i=1}^N {\bm g}_K(\bm{{O}}_i;\bm\gamma,\tau)$. Since $K$ increases with the sample size, the number of moment restrictions typically exceeds that of unknown parameters. So, we apply the GMM method for estimation. For a user-specified $(K+1)\times(K+1)$ positive definite matrix $\bm\Omega$, the GMM estimator of $(\bm\gamma,\tau)$ is given by

Hansen1982 showed that, with a fixed $K\ge p$, under some regularity conditions, $(\check{\bm\gamma},\check\tau)$ are consistent and asymptotically normally distributed but perform best only when $\bm\Omega$ is selected as the inverse of ${\bm\Upsilon}_{(K+1)\times(K+1)}= \mathbb{E}\{{\bm g}_K(\bm{{O}};\bm\gamma_0,\tau_0){\bm g}_K(\bm{{O}};\bm\gamma_0,\tau_0)^{\mathrm{\scriptscriptstyle T}}\}$. We use the initial estimator $(\check{\bm\gamma},\check\tau)$ to obtain an estimator of ${\bm\Upsilon}_{(K+1)\times(K+1)}$:

We then obtain the optimal GMM estimator

With a fixed $K$, Hansen1982 showed that under regularity conditions,

where ${\bm V}_K = \{{\bm B}_{(K+1)\times(p+1)}^{\mathrm{\scriptscriptstyle T}} {\bm\Upsilon}_{(K+1)\times(K+1)}^{-1} {\bm B}_{(K+1)\times(p+1)}\}^{-1}$ and

Large-sample properties

In this section, we establish the large-sample properties of the proposed estimator as the number of moment restrictions increases. We impose the following assumptions.

(ref)(i) rules out the degeneracy of moment restrictions. As indicated in (ref), we recommend orthonormalizing the basis functions such that the empirical second-moment matrix is the identity. (ref)(ii) requires sieve approximation error rates for the $p$-smooth function class, which have been well studied in the mathematical literature on approximation theory. For instance, suppose $\mathcal{X}$ and $\mathcal{Z}$ are compact subsets in $\mathbb{R}^{d_x}$ and $\mathbb{R}^{d_z}$, respectively, and ${\bm u}_K$ is a tensor product of polynomials or B-splines as introduced in (ref), Chen2007Handbook showed that (ref)(ii) holds with $\alpha=p/(d_x+d_z)$. (ref)(iii) restricts the number of moments to ensure the convergence and asymptotic normality of the proposed estimator. Newey1997 showed that if ${\bm u}_K$ is a power series, then $\zeta(K)=O(K)$, and it requires $K=o(N^{1/3})$. If ${\bm u}_K$ is a B-spline, then $\zeta(K)=O(\sqrt{K})$ and $K=o(\sqrt{N})$. These conditions guide the choice of ${\bm u}_K$, which is discussed in (ref). The asymptotic distribution of $\widehat{\bm\gamma}$ and $\widehat\tau$ is formally established in the following theorem.

(ref) demonstrates that $(\widehat{\bm\gamma},\widehat\tau)$ remains $\sqrt{N}$-consistent and asymptotically normally distributed when the number of moments increases slowly with the sample size. It follows intermediate lemmas in \Cref*{sec:app_intermediate_result} of the supplementary material. Specifically, we show that the initial estimators $(\check{\bm\gamma},\check\tau)$ obtained from (ref) are $\sqrt{N}$-consistent. We then demonstrate that (ref) still holds as $K\to\infty$ slowly. As a result, we obtain (ref) by calculating the limit of ${\bm V}_K$ in (ref). Notably, ${\bm V}_{\bm\gamma}$ is the semiparametric efficiency bound of $\bm\gamma_0$, derived in Theorem E.2 of CuiPuShiMiaoEtAl2023. Thus, the proposed estimator of $\bm\gamma_0$ is efficient. To see the efficiency of $\widehat\tau$, we consider the semiparametric local efficiency bound of $\tau_0$ derived in CuiPuShiMiaoEtAl2023, which requires the following assumption.

(ref)(i)-(ii) permits an alternative identification formula of ATE, given by $\tau=\mathbb{E}\{(-1)^{1-A}\allowbreak q(Z,A,X)Y\}$, using the treatment-confounding bridge function $q$ instead of $h$ in the proximal g-formula (ref) as established in CuiPuShiMiaoEtAl2023. (ref)(iii) ensures that $h$ and $q$ are uniquely identified by the integral equations (ref) and (ref), respectively. Combining these two identification results, CuiPuShiMiaoEtAl2023 showed that the efficient influence function of $\tau_0$, under the semiparametric model $\mathcal{M}_{sp}$, which does not restrict the observed data distribution other than the existence of a bridge function $h$ that solves (ref), evaluated at the submodel where (ref) holds, is

Therefore, the semiparametric local efficiency bound of $\tau_0$ under $\mathcal{M}_{sp}$ is $V_{\tau,\text{\textit{eff}}}=\mathbb{E}\{\psi_{\text{\textit{eff}}}(\bm{{O}})^2\}$.

(ref) shows that the proposed estimator $\widehat{\tau}$ has an asymptotic variance that is always no larger than the semiparametric efficiency bound $V_{\tau,\text{\textit{eff}}}$. The efficiency gain arises from parameterizing the bridge function $h$. We note that the doubly robust estimator, constructed using the efficient score (ref), attains the semiparametric efficiency bound $V_{\tau,\text{\textit{eff}}}$ only if both bridge functions $h$ and $q$ are correctly specified and consistently estimated. We also note that the condition that guarantees $V_{\tau}= V_{\tau,\text{\textit{eff}}}$ is unlikely to hold, except under highly contrived data distributions. As demonstrated in the simulation studies in (ref), it does not hold even in the most straightforward linear data-generating processes. Therefore, the proposed estimator generally outperforms the doubly robust estimator equipped with consistently estimated $h$ and $q$.

To highlight the deficiency of the plug-in procedure, the following \namecref{coro:super_plug_in} compares the proposed estimator $\widehat\tau$ with the theoretically optimal plug-in estimator $\widehat{\tau}^\textnormal{plug-in}({\bm m}_\text{\textit{eff}})$.

The condition that guarantees $V_{\tau}=V^\textnormal{plug-in}_{\tau,\text{\textit{eff}}}$ is satisfied if the observed distribution concerning the correlation between the score functions of $\tau_0$ and $\bm\gamma_0$, characterized by $R(Z,A,X)$, meets a particular structure. A sufficient condition is that there is no additive interaction between $A$ and $W$ in the model of $h(W,A,X)$. In this case, equality holds for $\bm\alpha=\0$. In case of its violation, (ref) indicates that the proposed estimator of $\tau_0$ outperforms plug-in estimators, even though an efficient estimator of $\bm\gamma_0$ is employed.

Implementation details

The GMM algorithm can be easily implemented using routine software, such as $\mathsf{gmm}$ in R. Here, we briefly discuss the selection of tuning parameters. A large class of sieves ${\bm u}_K(z,a,x)$ is feasible for implementing the proposed approach. In this paper, we suggest a tensor-product linear sieve basis. To simplify the presentation, we suppose both $Z$ and $X$ are scalar variables. Let $\{\varphi_j(x) : j=1, \ldots, K_1\}$ denote a sieve basis for $\mathcal{L}_2(\mathcal{X})$, the space of square Lebesgue integrable functions on $\mathcal{X}$. For instance, if $\mathcal{X}=[0,1]$, one can choose the basis of polynomials $\varphi_j(x)=x^{j-1}$ or B-spline, and if $\mathcal{X}$ is unbounded, one can select the Hermite polynomial basis $\varphi_j(x)=\exp(-x^2)x^{j-1}$. Similarly, let $\{\phi_k(z):k=1,\ldots,K_2\}$ denote a basis for $\mathcal{L}_2(\mathcal{Z})$. Then, the tensor-product sieve basis for $\mathcal{Z}\times\{0,1\}\times\mathcal{X}$ is given by $\{I(a=l)\varphi_j(x)\phi_k(z):j=1,\ldots,K_1;k=1,\ldots,K_2;l=0,1\}$ with the number of terms $K=2K_1K_2$. However, the number of tensor-product basis terms grows exponentially with the dimension of the arguments. In practice, additive separable bases can be employed to mitigate the issue of a large $K$. We refer readers to Newey1997 and Chen2007Handbook for further details on sieve selection.

Another tuning parameter is the number of moments $K$. While the large-sample properties of the proposed estimator permit a wide range of values for $K$, practical guidance on selecting smoothing parameters is necessary for applied researchers who generally have only one finite sample at their disposal. Here, we propose an asymptotic mean-square-error-based criterion, as in DonaldImbensNewey2009, for the data-driven selection $K$. Specifically, we choose $K\in\{1,\ldots,\bar{K}\}$ to minimize $S_{\rm{GMM}}(K):=\sum_{j=1}^p\{\widehat\Pi(K;{\bm e}_j)^2/N+\widehat\Phi(K;{\bm e}_j)\}$ defined in (ref), where $\widehat\Pi(K;{\bm e}_j)^2/N$ is an estimate of the squared bias term of $\widehat\gamma_j$ derived in NeweySmith2004, and $\widehat\Phi(K;{\bm e}_j)$ is an estimate of the asymptotic variance term of $\widehat\gamma_j$. We refer interested readers to DonaldImbensNewey2009 for further insights.

Numerical studies

Simulation

We construct Monte Carlo simulations to examine the finite-sample performance of the proposed approach. We generate i.i.d samples from the following data-generating process:

with the parameters $\bm\beta_a=(-0.1,0.5,0.5)^{\mathrm{\scriptscriptstyle T}}$, $\bm\beta_z=(0.5,1,0.5,1)^{\mathrm{\scriptscriptstyle T}}$, $\bm\beta_w=(1,-1,1)^{\mathrm{\scriptscriptstyle T}}$, and $\bm\beta_y=(1,0.5,0.5,1,1)^{\mathrm{\scriptscriptstyle T}}$. We consider two scenarios for $\sigma_j(x)$. Scenario I is a simple setting with no heteroskedasticity, $\sigma_j(x)\equiv 1$ for $j=1,2,3$. In Scenario II, heteroskedasticity is present, with $\sigma_1(x)\equiv 1$, $\sigma_2(x)=(0.3+x^2)^{-1/2}$, $\sigma_3(x)=(0.5+0.8x^2)^{-1/2}$. The average treatment effect $\tau_0=0.5$. We show in \Cref*{sec:app_supp_simu} of the supplementary material that the above data-generating mechanism is compatible with the following models of $h$ and $q$:

Our proposed method, denoted GMM-div, is computed using a power series and a data-driven smoothing parameter $K$ with $\bar{K}=12$ in both scenarios. For comparison, we include five additional methods: (i) Naive estimation, which regards all three covariates $(W,Z,X)$ as confounders, and estimates ATE using classical g-formula GreenlandRobins1986. (ii) RGMM, recursive GMM estimation with a fixed number of moments introduced in MiaoShiLiTchetgenTchetgen2024. Following their suggestion, we choose the moment restrictions $\mathbb{E}[\{Y-h(W,A,X;\bm\gamma)\}(1,Z,A,X)^{\mathrm{\scriptscriptstyle T}}]=0$. In this case, it is equivalent to the proximal outcome regression estimator introduced in CuiPuShiMiaoEtAl2023. (iii) P2SLS, the proximal two-stage least squares introduced in TchetgenTchetgenYingCuiShiEtAl2024, which can be implemented using the command $\text{ivreg}(Y\!\sim\! A+X+W\mid A+X+Z)$ in R. (iv) PIPW and PDR, the proximal inverse probability weighting and the proximal doubly robust estimator introduced by CuiPuShiMiaoEtAl2023, which require modeling the treatment confounding bridge function $q(z,a,x;\bm\theta)$. Specifically, $\widehat\tau_{PIPW}=N^{-1}\sum_{i=1}^N (-1)^{1-A_i}q(Z_i,A_i,X_i;\widehat{\bm\theta})Y_i$, and $\widehat\tau_{DR}$ solves the empirical analogue of (ref) by replacing $h$ and $q$ with their estimates. Here, $\widehat{\bm\theta}$ is obtained by solving the estimating equation $\mathbb{E}\{(-1)^{1-A}q(Z,A,X;\bm\theta)(1,W,A,X)^{\mathrm{\scriptscriptstyle T}} - (0,0,1,0)^{\mathrm{\scriptscriptstyle T}}\}=0$. Across all methods, we replicate 500 simulations at sample sizes of 400 and 800.

(ref) reports the absolute bias, standard error, root mean squared error, coverage probability, average length of 95% confidence intervals, and the power of the hypothesis testing problem $H_0:\tau_0=0~~v.s.~~H_1:\tau_0\neq 0$ in both scenarios. The power is calculated as the proportion of rejected cases across 500 replications. It reveals that: (i) The naive estimator is severely biased due to unmeasured confounding, while the other five methods exhibit negligible bias as expected in both scenarios. (ii) In Scenario I without heteroskedasticity, the five methods perform similarly, and the proposed method mainly selects $K=4$ across replications as shown in (ref)(a). (iii) In Scenario II, the proposed method demonstrates significantly lower standard error, narrower confidence intervals, and, notably, higher power in detecting non-zero causal effects. It benefits from the increased efficiency of additional selected moments as shown in (ref)(b). (ref) compares the empirical distributions of different methods, illustrating that the proposed method performs more concentrated around the true value. (iv) Due to the tradeoff between type I and type II errors, the proposed estimator provides a relatively anti-conservative confidence interval. Given the small sample size in Scenario II, this yields a slightly smaller CP than other methods. Nevertheless, it approaches the nominal level as the sample size increases.

Real data application

We apply the proposed method to re-analyze the Study to Understand Prognoses and Preferences for Outcomes and Risks of Treatments (SUPPORT), as considered in CuiPuShiMiaoEtAl2023 and TchetgenTchetgenYingCuiShiEtAl2024. The study aims to evaluate the effectiveness of right heart catheterization (RHC) in the initial care of critically ill patients ConnorsSperoffDawsonThomasEtAl1996. This dataset has also been widely analyzed in the causal inference literature, assuming no unmeasured confounding HiranoImbens2001,Tan2006,VermeulenVansteelandt2015. The treatment $A$ indicates whether a patient received an RHC within 24 hours of admission. Among 5735 patients, 2184 received the treatment, while 3551 did not. The outcome $Y$ is the number of days between admission and death or censoring at 30 days. The data include 71 baseline covariates, comprising 21 continuous variables and 50 dummy variables derived from categorical variables. These covariates include demographics (age, sex, race, education, income, and insurance status), estimated probability of survival, comorbidity, vital signs, physiological status, and functional status. Following TchetgenTchetgenYingCuiShiEtAl2024, we select $Z = (\mathsf{pafi1}, \mathsf{paco21})$, $W = (\mathsf{ph1}, \mathsf{hema1})$, and use the remaining covariates as $X$.

To implement the proposed method, we select $(1,A,Z,X)$, along with quadratic and cubic polynomials of the continuous variables in $X$, as candidates for constructing moment equations. We then apply the proposed data-driven approach to select $K$. (ref) plots the loss curve against $K$, from which we choose $K=81$ for estimation. Our method yields a point estimate of $-1.610$ with a standard error of 0.272, and the corresponding 95% confidence interval is (-2.143,-1.077). As a comparison, OLS yields an estimate of -1.249 (SE = 0.275) with a 95% CI of (-1.789,-0.709), while the proximal 2SLS proposed by TchetgenTchetgenYingCuiShiEtAl2024 produces an estimate of -1.798 (SE = 0.431) with a 95% CI of (-2.643,-0.954). Our estimate is closely aligned with the proximal 2SLS estimate, suggesting that OLS, which assumes no unmeasured confounding, may underestimate the harmful effect of RHC on 30-day survival among critically ill patients. Moreover, our estimate has a smaller standard error than the proximal 2SLS estimate, resulting in a narrower confidence interval.

Discussion

This paper proposes a simple, data-driven method for obtaining fully efficient estimates of the bridge function and average treatment effect. A key feature of this method is that it requires model assumptions only for the bridge function without imposing any additional assumptions on the data distribution. The proposed estimator typically outperforms the locally efficient DR estimator, demonstrating substantial advantages in detecting significant causal effects. Its feasibility and ease of implementation make it highly useful in practical applications.

Although our method is sufficiently efficient, we acknowledge that it lacks the protective capacity against potential model misspecification that the DR estimator provides. Both reviewers suggested jointly estimating $h$ and $q$ efficiently, similar to our approach, using the DR estimating equation for the treatment effect. However, as noted by CuiPuShiMiaoEtAl2023, this would not improve the efficiency of the treatment effect estimator. The DR estimating equation is Neyman orthogonal and thus locally insensitive to perturbations in $h$ and $q$, so efficiency gains in estimating $h$ and $q$ do not translate into efficiency gains for the ATE estimator. This reflects a trade-off between efficiency and robustness. In practice, we therefore recommend using both methods and comparing their results. For example, when the proposed method and the DR method yield similar estimates, it is reasonable to assume the bridge function is correctly specified and to adopt the proposed method's narrower confidence intervals.

Supplementary Materials

The supplementary material includes intermediate lemmas, additional simulation studies, and all technical proofs.