Doubly Robust Proximal Causal Inference under Confounded Outcome-Dependent Sampling

Suppose the data contain a sample of $n$ identically and independently distributed observations drawn from the population of interest, referred to as the “target population”. For each observation in the target population, let $A$ and $Y$ be the treatment and outcome of primary interest, respectively, assumed to be binary for the time being. Let $S$ be a binary indicator for selection into the study sample, such that the available data only include observations with $S=1$. Let $X$ be a vector of measured pre-treatment covariates that may be associated with $A$, $Y$ and $S$. Suppose that in addition to $X$, there exist pre-treatment latent factors, denoted as $U$, that may cause $A$, $Y$ and $S$. Figure 1(a) shows the directed acyclic graph (DAG) for the causal relationships between the above variables. Similar to Li et al. (2022) Section 2.5, we assume the following model, where $\beta_{0}$, the conditional log odds ratio given $(U,X)$, encodes the treatment effect of primary interest: {model}[Homogeneous odds ratio model; Li et al. (2022) Assumption 3’]

where $\mathrm{logit}(x)=\log\left\{x/(1-x)\right\}$ denotes the logit function and $\eta$ is an unknown real-valued function. Model 2.1 describes a semiparametric logistic regression model that assumes a homogeneous association between $A$ and $Y$ on the odds ratio scale, across strata of all measured and latent factors $(U,X)$. Model 2.1 encodes the structural model

where $Y(a)$ ($a=0,1$) denotes the potential outcome were the individual given the treatment status $A=a$, under standard identifiability assumptions of (a) consistency, i.e. $Y(a)=Y$ if $A=a$, which requires there is only one version of treatment and an individual’s outcome is not affected by others’ treatment status (Cole & Frangakis, 2009); (b) latent ignorability, i.e. $Y(a)\perp\!\!\!\perp A\mid U,X$, which essentially requires $(U,X)$ to contain all confounders of A-Y association; and (c) positivity, i.e. $0<P(A=1\mid U,X)<1$ almost surely (Hernán MA, 2020). Under these conditions, $\beta_{0}$ encodes the log conditional causal odds ratio effect in every $(U,X)$ stratum.

In Li et al. (2022) Assumption 3, they also considered a homogeneous risk ratio model, by which the marginal causal risk ratio can be identified. We focus on Model 2.1 for two reasons: first, proximal identification of $\beta_{0}$ under Model 2.1 permits treatment-induced selection under Assumption 2.1 below, while identification under the homogeneous risk ratio model appears to require that $A$ has no causal effect on selection other than through $Y$; second, Li et al. (2022) invoked a rare outcome assumption in the target population for proximal identification under the homogeneous risk ratio model, which restricts the applicability of the method. Such a rare outcome assumption is not necessary for proximal identification of $\beta_{0}$ under Model 2.1. As discussed in Section S2 of the Supplementary Materials, however, under standard identifiability assumptions (a)-(c) above and a rare outcome assumption akin to that assumed by Li et al. (2022), $\beta_{0}$ in Model 2.1 still approximates the log marginal causal risk ratio $P(Y(1)=1)/P(Y(0)=1)$. Alternatively, as discussed in Section S3 of the Supplementary Materials, under settings where control subjects are selected to be positive for an outcome variable that is a priori known to be not affected by the treatment, $\beta_{0}$ may recover the marginal causal risk ratio without invoking the rare outcome condition. Such a setting may occur, for example, in test-negative design studies of vaccine effectiveness, where the control subjects are selected to have a pre-determined infection that is not affected by the vaccine (Jackson & Nelson, 2013; Sullivan et al., 2016; Chua et al., 2020; Schnitzer, 2022).

Although to simplify the presentation, results given in the main text focus primarily on Model 2.1, Section S11 of the Supplemental Materials extends these results to the more general model

which explicitly accounts for effect modification of odds ratio with respect to measured covariates $X$.

In general, the conditional log odds ratio $\beta_{0}$ is not identifiable: first, the latent factors $U$ defining the strata are not observed; second, the data only include selected subjects with $S=1$ whilst Model 2.1 is defined in the target population.

As in Li et al. (2022) Assumption 2’, we make an important assumption regarding the selection mechanism that the treatment $A$ does not modify the conditional risk ratio of selection against the outcome in every $(U,X)$ stratum. Formally, we assume: {assumption}[No effect modification by $A$ on the outcome-dependent selection]For $a=0,1$ and some positive-valued unknown function $\xi$,

As mentioned above, a special case of Assumption 2.1 is if $A$ has no causal effect on $S$ other than through $Y$, a stronger assumption considered by Li et al. (2022) to identify the marginal causal risk ratio. As a result of Assumption 2.1, the conditional odds ratio given $(U,X,S=1)$ is identical to that given $(U,X)$, which leaves only the challenge of controlling for unmeasured confounding by $U$ in the selected sample.

To detect and correct for unmeasured confounding bias, the proximal inference framework proposes to leverage a pair of proxy measurements of unmeasured confounders (Miao et al., 2018a; Cui et al., 2020; Tchetgen Tchetgen et al., 2020). In this line of works, Li et al. (2022) developed proximal causal inference for $\beta_{0}$ under confounded outcome-dependent sampling as represented in Figure 1(a) (Jackson & Nelson, 2013; Chua et al., 2020). We similarly assume proxies of the latent factors are available. {assumption}[Proximal independence conditions] There exist a pair of proxies of the latent factor $U$: a treatment proxy (which Li et al. (2022) refers to as a “negative control exposure”), denoted as $Z$, and an outcome proxy (which Li et al. (2022) refers to as a “negative control outcome”), denoted as $W$, that satisfy

Assumption 2.2 requires the treatment proxy $Z$ and outcome proxy $W$ to satisfy certain conditional independences in the target population. Namely, the treatment proxy $Z$ has no direct effect on the primary outcome $Y$, selection indicator $S$, and the outcome proxy $W$; and the primary treatment $A$ has no direct effect on the outcome proxy $W$. It is crucial that both $Z$ and $W$ are both $U$-relevant, that is they carry information and therefore are associated with $U$ (Shi et al., 2020b; Tchetgen Tchetgen et al., 2020); this is formalized by Assumptions 2.2(a) and 2.3(a) below. As such, in the selected sample, an observed $A$-$W$ association after adjusting for $Y,X$, or an observed $Z$-$Y$ or $Z$-$W$ association after adjusting for $A,X$, may indicate the presence of bias due to the latent factor $U$. Figure 1 (b) and (c) show the DAGs of scenarios where Assumption 2.2 holds.

As in Li et al. (2022), we assumed the existence of a treatment confounding bridge function that connects the treatment proxy to the latent factor.

[Treatment confounding bridge function] There exists a treatment confounding bridge function $q(A,Z,X)$ that satisfies

almost surely for $a=0,1$.

Equation (1) formally defines a Fredholm integral equation of the first kind. To ensure the existence and identifiability of the solution to Equation (1), we make the following completeness assumptions, similar to those assumed by Li et al. (2022) and Cui et al. (2020):

[Completeness]    

1. (a)

   For any square-integrable function $g$, if $E[g(U)|A,Z,X,Y=0,S=1]=0$ almost surely, then $g(U)=0$ almost surely;

2. (b)

   For any square-integrable function $g$, if $E[g(Z)|A,W,X,Y=0,S=1]=0$ almost surely, then $g(Z)=0$ almost surely.

The completeness condition has been used to achieve identification for statistical functionals in econometrics and semiparametric statistics literature (Newey & Powell, 2003; D’Haultfoeuille, 2011). Essentially, Assumption 2.2(a) formalizes the $U$-relevance of $Z$ and requires that $Z$ has at least as much variability as $U$ in every $(A,X)$ stratum of the outcome-free subjects in the selected sample. When $U$ and $Z$ are both categorical variables, $Z$ necessarily has at least as many categories as $U$. Assumption 2.2(a) and the regularity conditions discussed in Section S4 of the Supplementary Materials, together constitute a sufficient condition for Assumption 2.2, i.e. the existence of $q$.

Similarly, Assumption 2.2(b) roughly requires that the outcome proxy $W$ is at least as variable as the treatment proxy $Z$, and is a sufficient condition under which $q(A,Z,X)$ can be uniquely identified from the selected sample, as stated in the theorem below:

Leveraging the proxies, Li et al. (2022) established the following identification result for $\beta_{0}$.

Equation (3) admits a closed-form solution for $\beta_{0}$:

Unique identification of $q$ by Assumption 2.2(b) is necessary. Although Equation (4) holds for any treatment confounding bridge function $q$ that satisfies Equation (1), $q$ can only be identified in the selected sample through solving Equation (2). If Equation (2) has multiple solutions, there is no guarantee that a solution to Equation (2) is also a solution to Equation (1).

From Theorem 2.2(b), it is clear that the treatment confounding bridge function $q(A,Z,X)$ plays a crucial role for identifying $\beta_{0}$. However, in case Assumption 2.2 does not hold, or Assumption 2.2 holds but Assumption 2.2(b) does not hold, then the solution to Equation (2) in the selected sample may not be unique and satisfy Equation (1) that defines a treatment confounding bridge function. As a result, the parameter $\beta_{0}$ identified by Equation (4), where $q$ is a solution to Equation (2), may not have a causal interpretation. In this section, we establish a new identification result for $\beta_{0}$ that does not rely on $q(A,Z,X)$. We instead define an outcome confounding bridge function as below: {assumption}[Outcome confounding bridge function]There exists an outcome confounding bridge function $h(A,W,X)$ that satisfies

almost surely. Miao et al. (2018b) and Cui et al. (2020) previously proposed proximal identification of average treatment effect using an outcome confounding bridge function. Our definition of $h(A,W,X)$ by Equation (5) is different from those in previous works of proximal causal inference, in that the conditional expectation of our $h(A,W,X)$ is defined to equal the conditional odds of the outcome, whilst the conditional expectation of the outcome confounding bridge function in previous works was defined to equal the conditional mean outcome. Theorem 2.4 below suggests that the standard outcome bridge function of prior literature would indeed fail to identify $\beta_{0}$ while our alternative definition yields the desired identification result.

In contrast with Assumption 2.2, we make the following completeness assumption for the existence and identifiability of the outcome confounding bridge function $h$.

[Completeness]    

1. (a)

   For any square-integrable function $g$, if $E[g(U)|A,W,X,Y=0,S=1]=0$ almost surely, then $g(U)=0$ almost surely;

2. (b)

   For any square-integrable function $g$, if $E[g(W)|A,Z,X,Y=0,S=1]=0$ almost surely, then $g(W)=0$ almost surely.

Similar to Assumption 2.2(a), Assumption 2.3(a) essentially requires that $W$ has sufficient variability relative to $U$ in every $(A,X)$ stratum of the outcome-free subjects in the selected sample, and constitutes a sufficient condition for Assumption 2.3, i.e. the existence of $h$, together with other regularity conditions, as discussed in Section S4 of the Supplementary Materials. Assumption 2.3(b) requires the opposite of Assumption 2.2(b), that the treatment proxy $Z$ is sufficiently variable relative to the outcome proxy $W$, and is a sufficient condition under which $h(A,W,X)$ can be uniquely identified from the selected sample, as stated in Theorem 2.3 below.

We introduce the new identification result in Theorem 2.4 below:

The integral equation (7) also admits a closed-form solution

Until now, identification of $\beta_{0}$ has required the existence and identification of either a treatment confounding bridge function or an outcome confounding bridge function. In this section, we present a closed-form expression for $\beta_{0}$, which has a desirable doubly-robustness property in the sense that it only requires one of the two confounding bridge functions to exist and be identified.

In Equation (9), setting $h(A,Z,X)=0$ recovers the closed-form solution (4), while setting $q(A,Z,X)=0$ gives the closed-form solution (8), illustrating the double robustness property of the closed-form expression (9).

In Section S5 of the Supplementary Materials, we compare our assumptions and results with those in Li et al. (2022), and summarise the additional conditions under which $\beta_{0}$ identifies the marginal causal risk ratio.

Denote data for the $i$th subject in the selected sample as $O_{i}=(A_{i},Y_{i},Z_{i},W_{i},X_{i})$, $i=1,\dots,n$. As suggested by Equations (4), (8) or (9), consistent estimators for $q(A,Z,X)$ and $h(A,W,X)$ can be used to construct estimators for $\beta_{0}$ following the standard plug-in principle (Bickel & Ritov, 2003). It remains to estimate the two confounding bridge functions. Below, we present crucial moment conditions for the identification of the confounding bridge functions.

Equations (4), (8) and (9) and Theorem 2.6 together suggest a parametric approach to estimate $\beta_{0}$: one may postulate suitable parametric working models $q(A,Z,X;\tau)$ and $h(A,W,X;\psi)$ where $\tau$ and $\psi$ are unknown parameters of finite dimensions. One can then estimate the nuisance parameters $\tau$ and $\psi$ by solving the estimating equations

with user-specified functions $\kappa_{1}$ and $\kappa_{2}$ of dimensions equal to that of $\tau$ and $\psi$ respectively. Closed-form parametric models for the confounding bridge functions can be derived in some settings as described below.

After obtaining the estimators $\hat{\tau}$ and $\hat{\psi}$ by solving Equations (12) and (13), $\beta_{0}$ can then be estimated with the following plug-in estimators:

Alternatively, one can jointly estimate $\beta_{0}$ and nuisance parameters $\tau$ and $\psi$ by solving estimating equations with the following estimating functions:

The resulting estimators are consistent and asymptotically normal, following standard estimating equation theory (Van der Vaart, 2000). Similar to Cui et al. (2020), we refer to $\hat{\beta}_{c,\text{PIPW}}$, $\hat{\beta}_{c,\text{POR}}$ and $\hat{\beta}_{c,\text{PDR}}$ as the proximal inverse probability weighted (PIPW) estimator, the proximal outcome regression (POR) estimator, and the proximal doubly robust (PDR) estimator, respectively. As discussed in Sections S2 and S3 in the Supplementary Materials, the PIPW estimator can estimate the marginal causal risk ratio assuming either (a) a homogeneous risk model, a treatment-independent selection mechanism, and a rare outcome condition in the target population, or (b) a test-negative design where the controls are selected to have an irrelevant outcome to the treatment. Because they estimate the same functional, the new POR and PDR estimators can be viewed as estimators of the marginal causal risk ratio under the same conditions.

As implied by Theorem 2.5, the estimator $\hat{\beta}_{c,\text{PDR}}$ also enjoys the desirable double robustness property– it is consistent for $\beta_{0}$ if either $q(A,Z,X)$ or $h(A,W,X)$ can be consistently estimated. We stated this result formally in the theorem below.

Theorem 2.9 indicates that the proposed PDR estimator has improved robustness over the PIPW and POR estimator against model misspecification of the two confounding bridge functions. If candidate proxies $W$ and $Z$ have different cardinalities, for example, when they are categorical variables of unequal dimensions, then only one of the completeness Assumptions 2.2 and 2.3 may hold. Say $W$ is of higher dimension, then Assumption 2.3(b) cannot hold, the outcome bridge function $h(A,W,X)$ may not be not uniquely identified, and $\hat{\beta}_{c,\text{POR}}$ may be biased for $\beta_{0}$. However, it may be possible to coarsen levels of $W$ to match those of $Z$ (Shi et al., 2020a). We highlight that the consistency of $\hat{\beta}_{c,\text{PDR}}$ only requires one of the two completeness assumptions to hold.

Our identification and estimation strategies involve user-specified functions $c(\cdot)$, $\kappa_{1}(\cdot)$ and $\kappa_{2}(\cdot)$. In principle, these nuisance functions can be chosen to construct a potentially more efficient estimator for $\beta_{0}$, $q$, and $h$ respectively, but the optimal choice of these functions may require modeling complex features of the observed data distribution. As discussed in Cui et al. (2020) and Stephens et al. (2014), such features can be considerably difficult to model correctly, and unsuccessful attempts to do so may lead to loss of efficiency, thus are not always worthwhile. In our simulation and real data analysis, we found that our generic choices for $c(\cdot)$, $\kappa_{1}(\cdot)$ and $\kappa_{2}(\cdot)$ deliver reasonable efficiency.

Until now, our estimation has required specification of parametric working models for $q(A,Z,X;\tau)$ and/or $h(A,W,X;\psi)$. When $X$ is continuous or high-dimensional, however, neither of these working models may contain the true confounding bridge functions. Ghassami et al. (2021) and Kallus et al. (2021) concurrently proposed a cross-fitting minimax kernel learning method for proximal learning of average treatment effect. Their method allows nonparametric estimation of the nuisance functions and $\sqrt{n}$-consistent estimation for the average treatment effect, even when both estimated confounding bridge functions are consistent at slower than $\sqrt{n}$ rates. In Section S7 of the Supplementary Materials, we demonstrate that their method is also applicable to our setting, and describe a semiparametric kernel estimator for the odds ratio parameter with flexibly estimated confounding bridge functions, $\sqrt{n}$-consistency, and asymptotic normality.

To prove Equation (2), we have