A simple sensitivity analysis method for unmeasured confounders via linear programming with estimating equation constraints

In this paper, we consider to estimate the ATE for the overall mean over the population in an observational study with two treatment groups. Let $Z$ be the treatment assignment: $Z=1$ if the subject is in the treated (exposed) group and $Z=0$ if in the control group. Let $X$ be a vector of baseline covariates and $Y$ be the observed outcome. We follow Rubin’s causal model framework.²¹ Let $Y^{(1)}$ and $Y^{(0)}$ be the potential outcomes if the subjects were assigned to the treated group ($Z=1$) and the control group ($Z=0$), respectively. Suppose the observational study enrolls $n$ subjects, and the observed data $(Y_{i},Z_{i},X_{i})$ for subject $i$ ($i=1,2,\dots,n$) available, which are independent and identically distributed copies of $(Y,Z,X)$. Denote $\mu_{1}=E[Y^{(1)}]$ and $\mu_{0}=E[Y^{(0)}]$. The ATE, which is of our primary interest to estimate, is defined by

$$\psi=\mu_{1}-\mu_{0}=E[Y^{(1)}]-E[Y^{(0)}].$$

In observational studies, owing to the absence of randomization, the potential influence of confounders should be carefully handled in estimating the ATE. The propensity score is widely used to adjust the bias due to confounding. The propensity score is defined by $e(X_{i})=p(Z_{i}=1\mid X_{i})$. Various methods, such as stratification, matching, and IPW,^{3, 4, 22} can be employed to adjust for confounding with the propensity score. The standard propensity analysis is conducted under the following assumptions: {assumption} Consistency: $Y_{i}=Z_{i}Y_{i}^{(1)}+(1-Z_{i})Y_{i}^{(0)}$ . {assumption} Positivity: there exists a small positive parameter $\delta$ such that $0<\delta\leq e(X_{i})\leq 1-\delta$, for all subjects $i$. {assumption} SITA: $(Y^{(1)}_{i},Y^{(0)}_{i})\perp\!\!\!\perp Z_{i}\mid X_{i}$. Assumption 3 implies that the bias due to confounding can be adjusted by using $X$ in principle. The SITA assumption is corresponding to the Missing At Random (MAR) in the missing data analysis context. In this paper, we handle situations in which the SITA is violated, which is corresponding to the concept of the Missing Not At Random (MNAR) in the missing data problem. We use the terminologies SITA and MAR exchangeably. In practice, the propensity score is unknown, and then some parametric models such as the logistic regression model is usually assumed. Let $\mbox{logit}(e(X_{i};\theta,\alpha))=\theta+\alpha^{\top}X_{i}$. The unknown parameters are usually estimated by solving the following score equation:

$$\sum_{i=1}^{n}\left(\begin{array}[]{cc}1\\ X_{i}\end{array}\right)\left(Z_{i}-\frac{exp(\theta+\alpha^{\top}X_{i})}{1+exp(\theta+\alpha^{\top}X_{i})}\right)=0.$$

(1)

Let the solution to the score equation for $(\theta,\alpha)$ be denoted by $(\hat{\theta},\hat{\alpha})$, and $\hat{e}(X_{i})=e(X_{i};\hat{\theta},\hat{\alpha})$. Then, we can determine the unique set of propensity scores for all subjects. In this paper, the propensity score estimated under the SITA assumption is called the MAR-based propensity score to avoid confusion; another type of propensity score is introduced in a later section, which is called the outcome-dependent propensity score. The IPW estimator for $\mu_{1}$ is defined by

$$\hat{\mu}_{1}=\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{\hat{e}(X_{i})}.$$

(2)

Similarly, we can estimate $\mu_{0}$ with

$$\hat{\mu}_{0}=\frac{1}{n}\sum_{i=1}^{n}\frac{(1-Z_{i})Y_{i}}{1-\hat{e}(X_{i})},$$

(3)

and then the ATE $\psi$ is estimated with

$$\hat{\psi}=\hat{\mu}_{1}-\hat{\mu}_{0}.$$

The aforementioned IPW estimator has an unstabilized form, which may suffer from extremely large weights when some propensity scores are very close to one or zero, and then can cause instability in the estimation. The stabilized IPW (SIPW) estimator introduces a stabilization term to the weights, which helps mitigate the impact of extreme weights. Specifically, the SIPW estimator for $\mu_{1}$ is defined by

$$\hat{\mu}_{1,SIPW}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{\hat{e}(X_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{\hat{e}(X_{i})}.$$

(4)

Similarly, we can estimate $\mu_{0}$ with

$$\hat{\mu}_{0,SIPW}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{1-Z_{i}}{1-\hat{e}(X_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{(1-Z_{i})Y_{i}}{1-\hat{e}(X_{i})},$$

(5)

and then the ATE $\psi$ is estimated with

$$\hat{\psi}_{SIPW}=\hat{\mu}_{1,SIPW}-\hat{\mu}_{0,SIPW}.$$

In this paper, we focus on the SIPW estimator.

If the SITA assumption holds and the model of the propensity score is correctly specified, the ATE is consistently estimated. However, the SITA assumption does not hold in the presence of unmeasured confounders.

In this section, suppose that the SITA assumption does not necessarily hold in the presence of unmeasured confounder $U$. The estimation of the ATE using the method in Section 2.1 is no longer valid. To address the issue of unmeasured confounders, the outcome-dependent propensity score approach^{23, 24, 25} has been successfully introduced. We define the outcome-dependent propensity scores by $o^{1}(X_{i},Y_{i}^{(1)})=P(Z_{i}=1\mid X_{i},Y_{i}^{(1)})$ and $o^{0}(X_{i},Y_{i}^{(0)})=P(Z_{i}=1\mid X_{i},Y_{i}^{(0)})$ for subjects in the treated and control group, respectively.

One may consider the logistic regression models for $o^{1}(X_{i},Y_{i}^{(1)})$ and $o^{0}(X_{i},Y_{i}^{(0)})$. Let us consider the models $\mbox{logit}{(o^{1}(X_{i},Y_{i}^{(1)};\theta^{1},\alpha^{1},\beta^{1}))}=\theta^{1}+\alpha^{1\top}X_{i}+\beta^{1}Y_{i}^{(1)}$ and $\mbox{logit}{(o^{0}(X_{i},Y_{i}^{(0)};\theta^{0},\alpha^{0},\beta^{0}))}=\theta^{0}+\alpha^{0\top}X_{i}+\beta^{0}Y_{i}^{(0)}$. The score equation (1) does not work for estimation of the unknown parameters in these models, since $Y_{i}^{(z)}$ is observed only for subjects with $Z_{i}=z$. The unknown parameters in the model of $o^{1}(X_{i},Y_{i}^{(1)};\theta^{1},\alpha^{1},\beta^{1})$ can be estimated by solving the following estimating equation:

$$\sum_{i=1}^{n}g(X_{i})\left(1-\frac{Z_{i}}{o^{1}(X_{i},Y_{i};\theta^{1},\alpha^{1},\beta^{1})}\right)=0,$$

(6)

where $g(X)$ is a vector of the same dimensions as $(\theta^{1},\alpha^{1},\beta^{1})$ and the solution to the estimating equation (6) is denoted by $(\hat{\theta}^{1},\hat{\alpha}^{1},\hat{\beta}^{1})$. Similarly, the unknown parameters in the model of $o^{0}(X_{i},Y_{i}^{(0)};\theta^{0},\alpha^{0},\beta^{0})$ can be estimated by solving the following estimating equation:

$$\sum_{i=1}^{n}g(X_{i})\left(1-\frac{1-Z_{i}}{1-o^{0}(X_{i},Y_{i};\theta^{0},\alpha^{0},\beta^{0})}\right)=0.$$

(7)

The dimension of $g(X)$ should be equal to that of $(\theta^{0},\alpha^{0},\beta^{0})$ to obtain a solution. The solution to the estimating equation (7) is denoted by $(\hat{\theta}^{0},\hat{\alpha}^{0},\hat{\beta}^{0})$. Denote $\hat{o}^{1}(X_{i},Y^{(1)}_{i})=o^{1}(X_{i},Y_{i};\hat{\theta}^{1},\hat{\alpha}^{1},\hat{\beta}^{1})$ and $\hat{o}^{0}(X_{i},Y^{(0)}_{i})=o^{0}(X_{i},Y_{i};\hat{\theta}^{0},\hat{\alpha}^{0},\hat{\beta}^{0})$, respectively. We can then estimate $\mu_{1}$ under the MNAR with

$$\hat{\mu}^{MNAR}_{1,SIPW}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{\hat{o}^{1}(X_{i},Y^{(1)}_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{\hat{o}^{1}(X_{i},Y^{(1)}_{i})}.$$

Similarly, we can estimate $\mu_{0}$ under the MNAR with

$$\hat{\mu}^{MNAR}_{0,SIPW}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{1-Z_{i}}{1-\hat{o}^{0}(X_{i},Y^{(0)}_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{(1-Z_{i})Y_{i}}{1-\hat{o}(X_{i},Y^{(0)}_{i})},$$

and then the ATE is estimated with

$$\hat{\psi}^{MNAR}_{SIPW}=\hat{\mu}^{MNAR}_{1,SIPW}-\hat{\mu}^{MNAR}_{0,SIPW}.$$

The SIPW estimator with the outcome-dependent propensity score can consistently estimate the ATE without the SITA assumption as long as the parametric models for the outcome-dependent propensity score are correctly specified. However, estimations with (6) and (7) often encounter an unidentifiability issue, wherein the model coefficients obtained through solving the estimating equations may not be uniquely determined. Miao et al.⁹ pointed out that even if the model for the propensity score has a known parametric form, the model is not identifiable without specifying a parametric outcome distribution. A unique solution to the estimating equations is only achieved when both the outcome model and the propensity score model are appropriately specified. Specifically, without additional restrictions or assumptions, sorely solving the estimating equations (6) is not sufficient to determine the coefficients $(\hat{\theta}^{1},\hat{\alpha}^{1},\hat{\beta}^{1})$ uniquely. Therefore, the outcome-dependent propensity score cannot solve the issue of the unmeasured confounder completely.

Along with the lines of the work by Robins et al.,^{26, 27} Brumback et al. ²⁸ proposed to quantify the impact of the unmeasured confounders by modeling the mean between-group difference of the potential outcomes, conditional on all observed covariates. The sensitivity function is defined by $c(z,X)=E[Y^{(z)}\mid Z=1,X]-E[Y^{(z)}\mid Z=0,X]$. If the SITA assumption holds, $c(z,X)$ equals zero. Thus, the sensitivity function can describe the magnitude of the departure from SITA assumption or the impact of the unmeasured confounders. Once we specify the sensitivity function $c(z,X)$, one can predict the mean function of the counterfactual variables conditional on $X$ and then estimate the ATE without the SITA assumption. Li et al. ¹⁵ criticized a technical difficulty in defining the sensitivity function when covariates $X$ contain multiple dimensions. Of note, in practical sensitivity analysis, if $X$ is multi-dimensional, not only the functional form but also the specific coefficients for each covariate are required to be specified. Such specifications were criticized to be unlikely to accurately reflect the relationship between the departure from SITA assumption and the potential outcomes. Li et al. ¹⁵ proposed a refinement by defining the sensitivity function as a function of the MAR-based propensity score: $c(z,e(X))=E[Y^{(z)}\mid Z=1,e(X)]-E[Y^{(z)}\mid Z=0,e(X)]$. The MAR-based propensity score is a one-dimension summary of observed covariates, and this refinement made the specification of the sensitivity function much simpler.

However, in reality, even with the simplification by Li et al.¹⁵, it is not an easy task to define a plausible range of sensitivity functions. Furthermore, their method still relies on the estimation of the MAR-based propensity score. Misspecification of the parametric model for the MAR-based propensity score may result in difficulty in interpreting the results of the sensitivity analysis.

Tan ¹⁹ proposed the marginal sensitivity model, which describes a relaxation of the SITA assumption. The model assumes a single sensitivity parameter, which permits the presence of the unmeasured confounders $U$, but restricts the extent of selection bias that can be attributed to these confounders. One can specify a parameter $\lambda$, and then the following inequality is supposed to hold:

$$1/\lambda\leq\frac{e(X_{i},U_{i})/(1-e(X_{i},U_{i}))}{\hat{e}(X_{i})/(1-\hat{e}(X_{i}))}\leq\lambda,$$

where $e(X_{i},U_{i})$ refers to the true propensity score measuring all covariates, and $\hat{e}(X_{i})$ refers to the estimated MAR-based propensity score. The single parameter $\lambda$, that is, the odds ratio (OR) between true propensity score and estimated propensity score, can control degree of unconfoundedness. When $\lambda=1$, the inclusion of additional confounders has no effect on the treatment odds. This implies that the allocation of the treatment is not influenced by confounding factors. That is, the SITA assumption holds. Increasing $\lambda$ represents the allowance for stronger extent to which the SITA assumption is violated. Tan ¹⁹ proposed a sensitivity analysis method to assess how the estimates based on the nonparametric likelihood change under the violation of the SITA assumption.

By introducing the marginal sensitivity model, the sensitivity analysis for unmeasured confounders can be applied to the IPW estimator under the MNAR. If $U$ was observed, one can estimate $\mu_{1}$ with

$$\left(\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{e(X_{i},U_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{e(X_{i},U_{i})}.$$

(8)

In practice, ${U}$ is unobserved and $\hat{\mu}_{1}$ in equation (8) actually makes no sense. However, under the marginal sensitivity model, $\lambda$ can link the unobserved true propensity score and the estimated MAR-based propensity score, so that it is possible to evaluate bounds of (8) under some constraints. That is

$$\begin{split}\max\mbox{or}\min&\left(\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{e(X_{i},U_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{e(X_{i},U_{i})}\\ \mbox{subject to}\quad&1/\lambda^{1}\leq\frac{e(X_{i},U_{i})/(1-e(X_{i},U_{i}))}{\hat{e}(X_{i})/(1-\hat{e}(X_{i}))}\leq\lambda^{1},\end{split}$$

(9)

where $\lambda^{1}$ is the pre-specified constant, which describes the upper and lower bounds of the discrepancy of the true propensity score from the estimated MAR-based propensity score for the estimation of $\mu_{1}$. As long as the true propensity score for all the subjects satisfies the constraint, the true $\mu_{1}$ should be bounded by the minimum and maximum of (9) asymptotically. It is possible to have an interval for $\mu_{0}$ in a similar way. This method under the marginal sensitivity model was proposed firstly by Zhao.¹⁸ In this method, the sensitivity parameter $\lambda$ quantifies the extent to which the SITA assumption is violated. However, it still suffers from defining a plausible range for the sensitivity parameter and reliance on correct specification of the MAR-based propensity score model.

It was criticized that the interval obtained by (9) may not be tight and the interval was asymptotically conservative.²⁰ Dorn and Guo ²⁰ proposed the quantile balancing method, a refinement based on the marginal sensitivity model. Let $F(y\mid x,z)=P(Y\leq y\mid X=x,Z=z)$ and the quantile function is defined by $Q_{t}(x,z)=\mbox{inf}\{q:F(q\mid x,z)\geq t\}$. For bounding $\mu_{1}$, the quantile balancing method solves the following optimization problem:

	$\displaystyle\max\mbox{or}\min$	$\displaystyle\left(\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}}{e(X_{i},U_{i})}\right)^{-1}\frac{1}{n}\sum_{i=1}^{n}\frac{Z_{i}Y_{i}}{e(X_{i},U_{i})}$		(10)
	subject to	$\displaystyle\sum_{i=1}^{n}\left(\begin{array}[]{cc}1\\ \hat{Q}_{\tau}(X_{i},1)\end{array}\right)\left(\frac{Z_{i}}{e(X_{i},U_{i})}-\frac{Z_{i}}{\hat{e}(X_{i})}\right)=0$		(13)
		$\displaystyle 1/\lambda^{1}\leq\frac{e(X_{i},U_{i})/(1-e(X_{i},U_{i}))}{\hat{e}(X_{i})/(1-\hat{e}(X_{i}))}\leq\lambda^{1},$		(14)

where $\tau=\frac{\lambda^{1}}{1+\lambda^{1}}$ and $\hat{Q}_{\tau}(X_{i},1)$ is estimated with some quantile regression models.²⁰ Bounding $\mu_{0}$ and the ATE can be achieved in a similar way. The quantile balancing method refined Zhao’s sensitivity analysis method ¹⁸ by adding the quantile function to balance the treatment assignment $Z$ over the true propensity score at population level. This additional constraint based on the estimated quantile function ensured asymptotic optimality of the interval obtained by solving (10). Although it solves asymptotic conservativeness in Zhao’s method,¹⁸ it still suffers from the misspecification of the estimated MAR-based propensity score. Moreover, the quantile function also requires specifying some parametric models or machine learning-related methods.