Sensitivity analysis for inverse probability weighting estimators via the percentile bootstrap

Introduction

A common task in statistics is to estimate the average treatment effect in observational studies. In such problems, the estimand is not identifiable from the observed data without further assumptions. The most common identification assumption, namely the “no unmeasured confounder” (NUC) assumption, asserts that the confounding mechanism is completely determined by some observed covariates. Based on this assumption, many statistical procedures have been proposed and thoroughly studied in the past decades, including propensity score matching (Rosenbaum and Rubin 1983), inverse probability weighting (Horvitz and Thompson 1952), and doubly robust and machine learning estimators (Robins et al. 1994, Van Der Laan and Rubin 2006, Chernozhukov et al. 2017).

However, the underlying NUC assumption is not verifiable using empirical data, posing serious threats to the usefulness of the subsequent statistical inferences that crucially rely on this assumption. A prominent example is the antioxidant vitamin beta carotene. Willett 1990, after reviewing observational epidemiological data, concluded that “Available data thus strongly support the hypothesis that dietary carotenoids reduce the risk of lung cancer”. Quite unexpectedly, four years later a large-scale randomized controlled trial (The Alpha-Tocopherol Beta Carotene Cancer Prevention Study Group 1994) reached the opposite conclusion and found “a higher incidence of lung cancer among the men who received beta carotene than among those who did not (change in incidence, 18 percent; 95 percent confidence interval, 3 to 36 percent)”. The most probable reason for the disagreement between the observational studies and the randomized trial is insufficient control of confounding. People taking vitamin supplements tend to be healthier at baseline and have higher socioeconomic position, and the adjustment for social and environmental confounding may be insufficient in many observational studies (Lawlor et al. 2004). For the same reason, observational studies and randomized controlled trials often give different conclusions for other vitamin supplements (Lawlor et al. 2004, Rutter 2007). We refer the reader to Rutter 2007 for other more recent examples of observational studies with probably misleading causal claims due to unmeasured confounding.

Criticism of confounding bias in observational studies dates at least back to Fisher 1958 who suggested that the association between smoking and lung cancer may be due to genetic predisposition. In response to this, Cornfield et al. 1959 conducted the first formal sensitivity analysis in an observational study. They concluded that, in order to explain the apparent odds ratio of getting lung cancer when there is no real causal effect, those with the genetic predisposition (or any hypothetical unmeasured confounder) must be 9 times more prevalent in smokers than in non-smokers. Because a genetic predisposition was seen as unlikely to have such a strong effect, this strengthened the evidence that smoking had a causal harmful effect.

Cornfield et al. 1959’s initial sensitivity analysis only applies to binary outcomes and ignores sampling variability. These limitations were later removed in a series of pioneering work by Rosenbaum and his coauthors (Rosenbaum 1987, Gastwirth et al. 1998, Rosenbaum 2002b, Rosenbaum 2002c). Rosenbaum’s sensitivity model considers all possible violations of the NUC assumption as long as the violation is less than some degree. Using a matched cohort design, Rosenbaum attempts to quantify the strength of the unmeasured confounders needed to not reject the sharp null hypothesis of no treatment effect. However, Rosenbaum’s framework is only limited to matched observational studies and usually requires effect homogeneity to construct confidence intervals for the average treatment effect.

Besides matching, another widely used method in causal inference is inverse probability weighting (IPW), where the observations are weighted by the inverse of the probability of being treated (Horvitz and Thompson 1952). IPW estimators have good efficiency properties (Hirano et al. 2003) and can be augmented with outcome regression to become “doubly robust” (Robins et al. 1994). There is much recent work aiming to improve the efficiency of IPW-type estimators by using tools developed in machine learning and high dimensional statistics (Van der Laan and Rose 2011, Belloni et al. 2014, Athey et al. 2016, Chernozhukov et al. 2017). However, they all heavily rely on the NUC assumption. Robustness of the IPW-type estimators to unmeasured confounding bias are usually studied using pattern-mixture models (Birmingham et al. 2003) or selection models (Scharfstein et al. 1999), in which the unmeasured confounding is usually modeled parametrically. See Richardson et al. 2014 for a recent overview and Section 7.1 for more references and discussion.

In this paper we propose a new framework for sensitivity analysis of missing data and observational studies problems, which can be applied to “smooth” estimators such as the inverse probability weighting (IPW) estimator and the “doubly robust” augmented IPW estimators. We consider a marginal sensitivity model introduced in Tan 2006, which is a natural modification of Rosenbaum’s model. Compared to existing sensitivity analyses of IPW estimators, our sensitivity model is nonparametric in the sense that we do not require the unmeasured confounder to affect the treatment and the outcome in a parametric way (though the propensity score can be modeled in a parametric way). This is appealing because we can never observe the unmeasured confounder and thus cannot test any parametric assumption.

Our marginal sensitivity model measures the degree of violation of the NUC assumption by the odds ratio between the conditional probability of being treated given the measured confounders and conditional probability of being treated given the measured confounders and the outcome/potential outcome variable (see Section 3 for the precise definition). Given a user-specified magnitude for this odds-ratio, the goal is to obtain a confidence interval of the estimand (the mean response vector in missing data problems and the average treatment effect in observational studies), with asymptotically at least $(1-\alpha)$ coverage probability, for all data generating distributions which violate the MAR assumption within this threshold (see Definition 4.2). Additionally, we would like to also report an interval of point estimates which is the range of possible point estimates under a sensitivity model. Our main contribution in this paper is to show such an interval of point estimates can be computed very efficiently for IPW estimators, and a valid confidence interval for sensitivity analysis can then be obtained by bootstrapping the interval of point estimates. The relationship to Efron 1979’s Bootstrap is explained below:

Efron’s bootstrap for a point-identified parameter:

Proposed percentile bootstrap procedure for sensitivity analysis (partially identified parameter):

The obvious and direct approach to obtain confidence intervals in a sensitivity analysis involves using the asymptotic normal distribution of the IPW estimators. However, the asymptotic sandwich variance estimators of the IPW estimates are quite complicated. Finding extrema of the point estimates and variance estimates over a collection of sensitivity models is generally computationally intractable. This problem can be circumvented using the percentile bootstrap, reducing the problem to solving many linear fractional programs. The confidence interval our procedure constructs has the following desirable properties:

The interval covers the entire partially identified region of the estimand with asymptotic probability of the coverage at least $(1-\alpha)$, that is, it has asymptotically ‘strong nominal coverage’ (Vansteelandt et al. 2006, Definition 3). The proof involves establishing the limiting normal distribution of the bootstrapped IPW estimators and a generalized minimax/maximin inequality, which justifies the interchange of quantile and infimum/supremum (Theorem 4.4 and Corollary 5.1).

Maximizing/minimizing the IPW point estimate over the marginal sensitivity model is very efficient computationally. This is a linear fractional programming problem (ratio of two linear functions) which can be reduced to solving a linear program using the well-known Charnes-Cooper transformation (Section 4.4). In fact, by local perturbations, it can be shown that the solution of the resulting linear program has the same/opposite order as the outcome vectors, using which the confidence interval can be computed in time linear in sample size, for every bootstrap resample (Proposition 4.5).

Our method for constructing confidence intervals under the marginal sensitivity model, for the IPW estimator for the mean response with missing data, is described in Section 4. The extension to estimating the average treatment effect in observational studies is discussed in Section 5. The percentile bootstrap approach and the reduction to linear programming are very general and can be extended to sensitivity analysis of other smooth estimators, such as IPW estimates of the mean of the non-respondents and the average treatment effect on the treated (Section 6.1), the augmented inverse probability weighting estimator (Section 6.2), and inference for partially identified parameters (Section 7.3). We review some related sensitivity analysis methods and compare our framework with Rosenbaum’s sensitivity analysis in Section 7. In Section 8 we evaluate the performance of our method using some numerical examples, including a simulation study in Section 8.1 and a real data example in Section 8.2 where Rosenbaum’s sensitivity analysis is also applied. Theoretical proofs can be found in Appendix A and the supplementary file. The R code for the our method and the real data example in Section 8.2 can be found at https://github.com/qingyuanzhao/bootsens.

The Missing Data Problem: Background and Notation

Since some responses or potential outcomes are not observed, the estimand $\mu$ defined above is not identifiable without further assumptions. For the missing data problem, Rubin 1976 used the term “missing at random” (MAR) to describe data that are missing for reasons related to completely observed variables in the data set:

A\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y|\bm{X} under $F_{0}$.

With the additional assumption that no subject is missing or receives treatment/control with probability $1$, the parameter $\mu$ is identifiable from the data.

Since the seminal works of Rubin 1976 and Rosenbaum and Rubin 1983, many statistical methods have been developed for the missing data and observational studies problem based on first estimating the conditional probability $e_{0}(\bm{x})$ (often called the propensity score in observational studies). One distinguished example is the inverse probability weighting (IPW) estimator which dates back to Horvitz and Thompson 1952,

where $\hat{e}(\bm{X})$ is a sample-estimate of $e_{0}(\bm{X})$. Observe that under 2.1 and 2.2,

Sensitivity models

To rebut such criticism and make the statistical analysis more credible, a natural question is how sensitive the results are to the violation of the MAR assumption. This is clearly an important question, and, in fact, sensitivity analysis is often suggested or required for empirical publications. For example, sensitivity analysis is required in the Patient-Centered Outcome Research Institute (PCORI) methodology standards for handling missing data, see standard MD-4 in https://www.pcori.org/research-results/about-our-research/research-methodology/pcori-methodology-standards. See also Little et al. 2012. Many sensitivity analysis methods have been subsequently developed in the missing data and observational studies problems. In this paper we will consider a sensitivity model that is closely related to Rosenbaum’s sensitivity model (Rosenbaum 1987, Rosenbaum 2002c). We refer the reader to Robins 1999, Scharfstein et al. 1999, Imbens 2003, Altonji et al. 2005, Hudgens and Halloran 2006, Vansteelandt et al. 2006, McCandless et al. 2007, VanderWeele and Arah 2011, Richardson et al. 2014, Ding and VanderWeele 2016 for alternative sensitivity analysis methods and Section 7 for a more detailed discussion.

When the MAR assumption is violated, equation (3.1) is no longer valid. It is well well-known that $e_{0}(\bm{x},y)$ is generally not identifiable from the data without the MAR assumption (proof included in the supplement for completeness). For this reason we shall refer to a user-specified function $e_{0}(\bm{x},y)$ as a sensitivity model.

In this paper we consider the following collection of sensitivity models in which the degree of violation of the MAR assumption is quantified by the odds ratio of $e_{0}(\bm{x},y)$ and $e_{0}(\bm{x})$.

Fix a parameter $\Lambda\geq 1$. For the missing data problem, we assume $e(\bm{x},y)\in\mathcal{E}(\Lambda)$, where

The relationship between this sensitivity model and Rosenbaum’s sensitivity model is examined in Section 7.2.

The set $\mathcal{E}(\Lambda)$ becomes larger as $\Lambda$ increases. When $\Lambda=1$, $\mathcal{E}(1)$ contains the singleton $\{e_{0}(\bm{x})\}$ and corresponds to MAR. When $\Lambda=\infty$, $\mathcal{E}(\infty)$ contains all functions $e(\bm{x},y)$ that are bounded between $0$ and $1$. To understand 3.1, it can be conceptually easier to imagine that there is an unobserved variable $U$ that “summarizes” all unmeasured confounding, with the relation between $U$ and the potential outcomes being unconstrained in any way. This leads to an alternative “added variable” representation of sensitivity model (Rosenbaum 2002c). Robins 2002 pointed out that it suffices to consider the conditional probabilities when $U$ is either one of the potential outcomes. In the remainder of the paper we will follow Robins’ suggestion to simplify the notation.

It is often convenient to use the logistic representation of the marginal sensitivity model (3.2). Denote

and $h_{0}(\bm{x},y)=g_{0}(\bm{x})-g_{0}(\bm{x},y)$ be the logit-scale difference of the observed data selection probability and the complete data selection probability. If we further introduce the notation $e^{(h)}(\bm{x},y)=[1+\exp(h(\bm{x},y)-g_{0}(\bm{x}))]^{-1}$, then

and $\lambda=\log\Lambda$. This shows that the marginal sensitivity model puts bound on the $L_{\infty}$-norm of $h$. It is easy to verify that $e^{(h)}(\bm{x},y)\to 0$ or 1 if $h(\bm{x},y)\to\pm\infty$.

Confidence Interval for the Mean Response

In this section we construct confidence intervals for the mean response $\mu$ under the marginal sensitivity models.

which works seamlessly with our sensitivity model because the degree of violation of MAR is quantified by odds ratio. However our framework can be easily applied to other parametric models $e_{\bm{\beta}}(\bm{x})$ (e.g. different links).

As in (3.3), the constraint (4.3) is equivalent to $h_{\bm{\beta}_{0}}(\bm{x},y)\in\mathcal{H}(\lambda)$ for $h_{\bm{\beta}}(\bm{x},y)=g_{\bm{\beta}}(\bm{x})-g_{0}(\bm{x},y)$ and $\lambda=\log\Lambda$.

Next we define the shifted estimand under a specific sensitivity model $h$:

The definitions of $\mu^{(h)}$ and $e^{(h)}$ above can be easily extended to the non-parametric marginal sensitivity model by replacing $g_{\bm{\beta}_{0}}(\bm{x})$ by $g_{0}(\bm{x})$. Our statistical method can be applied regardless of the “baseline” choice of parametric or nonparametric model for $e_{0}(\bm{X})$, as long as the model is “smooth” enough so that the bootstrap is valid. For simplicity, hereafter we will often refer to $h(\bm{x},y)\in\mathcal{H}(\lambda)$ instead of $e^{(h)}(\bm{x},y)$ as the sensitivity model, since the former does not depends on the choice of the working model for $e_{0}(\bm{X})$. Consequently we will also call $\mathcal{H}(\lambda)$ a collection of sensitivity models without specifying which parametric/nonparametric model was used for $e_{0}(\bm{X})$.

Compared to 4.1, the only difference in (4.3) is that the postulated sensitivity model $e(\bm{x},y)$ is compared to the parametric model $e_{\bm{\beta}_{0}}(\bm{x})$ instead of the nonparametric probability $e_{0}(\bm{x})$. In other words, the parametric sensitivity model considers both

Model misspecification, that is, $e_{\bm{\beta}_{0}}(\bm{x})\neq e_{0}(\bm{x})$; and

Missing not at random, that is, $e_{0}(\bm{x})\neq e_{0}(\bm{x},y)$.

This is a desirable feature in the sense that the term “sensitivity analysis” is also widely used as the analysis of an empirical study’s robustness to parametric modeling assumptions. However, it might also make the choice of the sensitivity parameter $\lambda$ more difficult in practice. This issue is briefly explored in the simulation in 8.

With these notations, we can now define what is meant by a confidence interval under a collection of sensitivity models:

A data-dependent interval $[L,U]$ is called a confidence interval for the mean response $\mu$ with at least $(1-\alpha)$ coverage under the collection of sensitivity models $\mathcal{H}(\lambda)$ (may corresponds to $\mathcal{E}(\Lambda)$ or $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$, see Remark 4.1), if

2. The IPW Point Estimates

Intuitively, a confidence interval $[L,U]$ as in Definition 4.2 must at least include a point estimate of $\mu^{(h)}$ for every $h\in\mathcal{H}(\lambda)$. To this end, let $h$ be a postulated sensitivity model. The corresponding selection probability $e^{(h)}(\bm{x},y)$ can then be estimated by

where $\hat{e}^{(h)}$ is as in (4.5). However, it is well known that, even under the MAR assumption, the IPW estimator can be unstable when the selection probability $e_{0}(\bm{x})$ (or the parametric approximation $e_{\bm{\beta}_{0}}(\bm{x})$) is close to $0$ for some $\bm{x}\in\mathscr{X}$ (Kang and Schafer 2007). To alleviate this issue, the stabilized IPW (SIPW) estimator, obtained by normalizing the weights, is often used in practice:

It is easy to see that $\hat{\mu}^{(h)}$ estimates $\mu^{(h)}$ defined in (4.4).

Compared to IPW, the SIPW estimator is sample bounded (Robins et al. 2007), that is,

This property is even more desirable in sensitivity analysis because $e^{(h)}$ is almost always not the true selection probability, so the total unnormalized weights $\frac{1}{n}\sum_{i=1}^{n}A_{i}/\hat{e}^{(h)}(\bm{X}_{i},Y_{i})$ can be very different from $1$. For this reason, we will use the SIPW estimator for the remainder of this paper.

Heuristically, the confidence interval $[L,U]$ should at least contain the range of SIPW point estimates, $\big[\inf_{h\in\mathcal{H}(\lambda)}\hat{\mu}^{(h)},\sup_{h\in\mathcal{H}(\lambda)}\hat{\mu}^{(h)}\big]$. We defer the numerical computation of the extrema of SIPW point estimates till Section 4.4. For now we will focus on constructing the confidence interval $[L,U]$ assuming the range of point estimates can be efficiently computed.

3. Constructing the Confidence Interval

To construct a confidence interval for $\mu$, we need to consider the sampling variability of the SIPW estimator described above. In the MAR setting, the most common way to estimate the variance is the asymptotic sandwich formula or the bootstrap (Efron and Tibshirani 1994, Austin 2016). In sensitivity analysis, we also need to consider all possible violations of the MAR assumption in $\mathcal{H}(\lambda)$. In this case, optimizing the estimated asymptotic variance over $\mathcal{H}(\lambda)$ is generally computationally intractable, as explained below.

We begin by showing how individual confidence intervals of $\mu^{(h)}$ for $h\in H(\lambda)$ can be combined into a confidence interval in sensitivity analysis.

Suppose there exists data-dependent intervals $[L^{(h)},U^{(h)}]$ such that

holds for every $h\in\mathcal{H}(\lambda)$.

Let $L=\inf_{h\in\mathcal{H}(\lambda)}L^{(h)}$, $U=\sup_{h\in\mathcal{H}(\lambda)}U^{(h)}$. Then $[L,U]$ is an asymptotic confidence interval of $\mu$ with at least $(1-\alpha)$ coverage under the collection of sensitivity models $\mathcal{H}(\lambda)$.

Moreover, if there exists $\alpha^{\prime}\in[0,\alpha]$ (not depending on $h$) such that

for all $h\in\mathcal{H}(\lambda)$, then the union interval $[L,U]$ covers the partially identified region with probability at least $1-\alpha$,

4.1 suggests the following way to construct a confidence interval for $\mu$ in sensitivity analysis using the asymptotic distribution of $\hat{\mu}^{(h)}_{1}$. Using the general theory of $Z$-estimation, it is not difficult to establish that

See, for example, Lunceford and Davidian 2004 or C.2 in the supplement. Then using the sandwich variance estimator $(\hat{\sigma}^{(h)})^{2}$, an asymptotically confidence interval of $\mu^{(h)}$ at least $(1-\alpha)$ coverage is

However, the standard error $\hat{\sigma}^{(h)}$ is a very complicated function of $h$ (see C.2) and numerical optimization over $h\in\mathcal{H}(\lambda)$ is practically infeasible.

3.2. The Percentile Bootstrap

where $Q_{\alpha}(\hat{\hat{\mu}})$ is the $\alpha$-percentile of $\hat{\hat{\mu}}$ in the bootstrap distribution, that is,

We begin by showing that the percentile bootstrap interval $[L^{(h)},U^{(h)}]$ is an asymptotically valid confidence interval of $\mu^{(h)}$ for the parametric sensitivity model $e^{(h)}\in\mathcal{E}_{\beta_{0}}(\Lambda)$.

In the logistic model (4.2) and under regularity assumptions (see Assumption C.1 in the supplement), for every $e^{(h)}\in\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$ we have

The proof of this result, which invokes the general theory of bootstrap for $Z$-estimators (Wellner and Zhan 1996, Kosorok 2006, Chapter 10), is explained in Appendix C in the supplement.

Our percentile bootstrap confidence interval under the collection of sensitivity models $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$ is given by $[L,U]$ where

The important thing to observe here is that the infimum/supremum is inside the quantile function in (4.10), which makes the computation especially efficient using linear programming (see Section 4.4). The interchange of quantile and infimum/supremum is justified in the following generalized (von Neumann’s) minimax/maximin inequalities (see Cohen 2013 for a similar result for finite sets).

Let $L,U$ be as defined in (4.10) and $L^{(h)},U^{(h)}$ be as defined in (4.9). Then

Asymptotic validity of the confidence interval $[L,U]$ in Equation 4.10 then immediately follows from the validity of the union method (4.1), the validity of the percentile bootstrap (4.2), and 4.3.

Under the same assumptions as in Theorem 4.2, $[L,U]$ is an asymptotic confidence interval of the mean response $\mu$ with at least $(1-\alpha)$ coverage, under the collection of sensitivity models $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$. Furthermore, $[L,U]$ covers the partially identified region $\{\mu^{(h)}:e^{(h)}\in\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)\}$ with probability at least $1-\alpha$.

4. Range of SIPW Point Estimates: Linear Fractional Programming

Theorem 4.4 transformed the sensitivity analysis problem to computing the extrema of the SIPW point estimates, $\inf_{h\in\mathcal{H}(\lambda)}\hat{\hat{\mu}}^{(h)}_{b}$ and $\sup_{h\in\mathcal{H}(\lambda)}\hat{\hat{\mu}}^{(h)}_{b}$. In practice, we only repeat this over $B~(\ll n)$ random resamples and compute the interval by

For notational simplicity, below we consider how to compute the extrema $\inf_{h\in\mathcal{H}(\lambda)}\hat{{\mu}}^{(h)}$ and $\sup_{h\in\mathcal{H}(\lambda)}\hat{{\mu}}^{(h)}$ using the full observed data instead of the resampled data. This also gives an interval of point estimates under the sensitivity model that is shorter than the confidence interval. We recommend reporting both intervals in any real data analysis, see Section 8 for some examples.

Recalling (4.6) and (4.7), computing the extrema is equivalent to solving

where the optimization variables are $z_{i}=e^{h(\bm{X}_{i},Y_{i})}$ for $i\in[n]$. All the other variables are observed or can be estimated from the data. Notice that $\hat{g}(\bm{x})$ needs to be re-estimated in every bootstrap resample. Without loss of generality, assume that the first $1\leq m<n$ responses are observed, that is $A_{1}=A_{2}=\cdots A_{m}=1$ and $A_{m+1}=\cdots=A_{n}=0$, and suppose that the observed responses are in decreasing order, $Y_{1}\geq Y_{2}\geq\cdots\geq Y_{m}$. Then (4.12) simplifies to

This optimization problem is the ratio of two linear functions of the decision variables $\bm{z}$, hence called a linear fractional programming. It can be transformed to linear programming by the Charnes-Cooper transformation (Charnes and Cooper 1962). Denote

This translates (4.13) into the following linear programming:

Therefore, the range of the SIPW point estimate can be computed efficiently by solving the above linear program. Furthermore, the following result shows that the solution of (4.4) must have the same or opposite order as the outcomes $\bm{Y}$, which enables even faster computation of (4.4).

Suppose $(z_{i})_{i=1}^{m}$ solves the maximization problem in (4.4). Then $(z_{i})_{i=1}^{m}$ has the same order as $(Y_{i})_{i=1}^{m}$, that is, if $Y_{s_{1}}>Y_{s_{2}}$ then $z_{s_{1}}>z_{s_{2}}$, for $1\leq s_{1}\neq s_{2}\leq m$. Furthermore, there exists a solution $(z_{i})_{i=1}^{m}$ and a threshold $M$ such that $z_{i}=\Lambda$, if $Y_{i}\geq M$, and $z_{i}=\frac{1}{\Lambda}$, if $Y_{i}<M$. The same conclusion holds for the minimizer in (4.4) with $(Y_{i})_{i=1}^{m}$ replaced by $(-Y_{i})_{s=1}^{m}$.

4.5 implies that we only need to compute the objective of (4.13) for at most $m$ choices of $(z_{i})_{i\in[m]}$ by enumerating the index where it changes from $\Lambda$ to $\Lambda^{-1}$:

It is easy to see that we only need $O(m)$ time to compute $\sum_{i=1}^{m}Y_{i}(1+z_{i}e^{-\hat{g}(\bm{X}_{i})})$ and $\sum_{i=1}^{m}(1+z_{i}e^{-\hat{g}(\bm{X}_{i})})$, for all the $m$ choices in (4.15). Hence, the computational complexity to solve the linear fractional programming (4.13) is $O(m)$. This is the best possible rate since it takes $O(m)$ time to just compute the objective once.

The above discussion suggests that the interval $[L_{B},U_{B}]$, which entails finding $\inf_{h\in\mathcal{H}(\lambda)}\hat{\hat{\mu}}^{(h)}_{b}$ and $\sup_{h\in\mathcal{H}(\lambda)}\hat{\hat{\mu}}^{(h)}_{b}$ for every $1\leq b\leq B$, can be computed extremely efficiently. The computational complexity is $O(nB+n\log n)$ ignoring the time spent to fit the logistic propensity score models, where the extra $n\log n$ is needed for sorting the data. Indeed, even under the MAR assumption, we need $O(nB)$ time to compute the Bootstrap confidence interval of $\mu$. In conclusion, our proposal requires almost no extra cost to conduct a sensitivity analysis for the IPW estimator than to obtain its bootstrap confidence interval under MAR.

Confidence Interval for the ATE in the Sensitivity Model

can consistently estimate the ATE $\Delta$, where $\hat{e}(\bm{X})$ is a sample-estimate of $e_{0}(\bm{X})$.

One can analogously define the difference between $e_{a}(\bm{x},y)$ and $e_{\bm{\beta}_{0}}(\bm{x})$ in the logit scale and rewrite (5.1) as an $L_{\infty}$-constraints on the difference. Similarly, one can define the shifted propensity score and the shifted ATE. We omit the details for brevity but would like to mention that the shifted estimand $\Delta^{(h_{0},h_{1})}$ now depends on two $h$ functions corresponding to the two potential outcomes. The SIPW estimator of $\Delta^{(h_{0},h_{1})}$ can be defined in the same way as (4.7).

Now, just as in Section 4.3.2, we can use the percentile bootstrap obtain a asymptotically valid interval for $\Delta^{(h_{0},h_{1})}$:

where $h_{0}$ and $h_{1}$ are held fixed and $Q_{\frac{\alpha}{2}}(\hat{\hat{\Delta}}^{(h_{0},h_{1})})$ is the $\alpha$-th bootstrap quantile of the SIPW estimates. Then, by interchanging the maximum/minimum and the quantile as in Theorem 5.1, we obtain a confidence interval for the ATE:

Under the same assumptions in Theorem 4.2, the confidence interval

covers the average treatment effect $\Delta$ is with probability at least $(1-\alpha)$ asymptotically under the collection of parametric sensitivity models (5.1).

The interval in (5.2) can be computed efficiently using linear fractional programming as in Section 4.4. To simplify notation, assume, without loss of generality, the first $m\leq n$ units are treated ($A=1$) and the rest are the control ($A=0$), and that the outcomes are ordered decreasingly among the first $m$ units and the other $n-m$ units. Then, as in (4.13), computing the interval (5.2) is equivalent to solving the following optimization problem:

where $z_{i}=e^{h_{1}(\bm{X}_{i},Y_{i})}$, for $1\leq i\leq m$ and $z_{i}=e^{-h_{0}(\bm{X}_{i},Y_{i})}$, for $m+1\leq i\leq n$. Note that the variables $(z_{i})_{i=1}^{m}$ and $(z_{i})_{i=m+1}^{n}$ are separable in (5.3), so we can solve the maximization/minimization problem in (5.3) by solving one maximization/minimization problem for $(z_{i})_{i=1}^{m}$ and one minimization/maximization problem for $(z_{i})_{i=m+1}^{n}$. Therefore, similar to the missing data problem, the time complexity to obtain the range of the SIPW estimates $\hat{\Delta}$, over a range of $B$ bootstrap resamples, is only $O(nB+n\log n)$.

Extensions

In this section we discuss three extensions of the general framework described in Section 4.

Then as in 4.4, under the collection of parametric sensitivity models $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$,

is an asymptotic confidence interval of the non-respondent mean $\mu_{0}$ with at least $(1-\alpha)$ coverage, where the SIPW estimate is

where $\hat{g}$ is as in (4.6) and $\hat{\hat{\mu}}^{(h)}_{01},\hat{\hat{\mu}}^{(h)}_{02},\ldots,\hat{\hat{\mu}}^{(h)}_{0B}$ are the $B$ bootstrap resamples of $\hat{\mu}_{0}^{(h)}$. As before, the interval (6.1) can be computed efficiently using linear fractional programming. In particular, it is easy to verify that 4.5 still holds and thus we only need to consider the $O(n)$ candidate solutions in Equation 4.15.

2. Augmented Inverse Probability Weighting

When MAR does not hold, by taking expectation conditioning on $\bm{X}$ and $Y$, (6.2) implies that

Now, as in Section 4, we consider the collection of parametric sensitivity models $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)$ so $\bar{e}(\bm{x})=e_{\bm{\beta}_{0}}(\bm{x})$. For $h\in\mathcal{H}(\lambda)$, define

where $\hat{e}^{(h)}(\bm{x},y)=\big[1+e^{h(\bm{x},y)-\hat{g}(\bm{x},y)}\big]^{-1}$ and where $\hat{g}$ is as in (4.6). The second term on the right hand side of (6.3) is not sample bounded, so it is often preferable to use the stabilized weights. This results in the following stabilized AIPW (SAIPW) estimator:

As before, this estimates $\mu$ when $h=h_{\bm{\beta}_{0}}$.

Compared to the SIPW estimator (4.7), in (6.4) we replace the response $Y_{i}$ by $Y_{i}-\hat{f}(\bm{X}_{i})$ and add an offset term $\frac{1}{n}\sum_{i=1}^{n}\hat{f}(\bm{X}_{i})$. Therefore, computing the extrema of (6.4) can still be formulated as linear fractional programming and the numerical computation is still efficient. To construct asymptotically valid confidence intervals, the outcome regression model $\hat{f}(\bm{X}_{i})$ must be parametric (for example, linear regression). The $Z$-estimation framework in Appendix C in the supplement can then be extended to show that 4.2 (validity of percentile bootstrap) still holds for SAIPW.

3. Lipschitz Constraints in the Sensitivity Model

So far we have focused on the marginal sensitivity models,

Although this model is very easy to interpret, some deviations $h\in\mathcal{H}(\lambda)$ may be deemed unlikely because the function $h$ is not smooth. Here, we consider an extension of our sensitivity model which assumes $h$ is also Lipschitz-continuous. Formally, define

As $\mathcal{H}_{\lambda,L}\subseteq\mathcal{H}(\lambda)$, the validity of the percentile bootstrap obviously holds for functions in $\mathcal{H}_{\lambda,L}$. To obtain a range of point estimates and confidence intervals under $\mathcal{E}_{\bm{\beta}_{0}}(\Lambda,L)$, we just need to add the following $n(n-1)$ constraints in the optimization problem (4.12):

These constraints are linear in $\{z_{i}\}_{i=1}^{m}$ and the resulting optimization problem is still a linear fractional program, which can be efficiently computed. However, 4.5 no longer holds, so we cannot use the algorithm described after 4.5 to solve the optimization problem corresponding to $\mathcal{H}_{\lambda,L}$.

Discussion

If there is no unmeasured confounder, the potential outcome $Y(a)$ is independent of the treatment $A$ given $\bm{X}$, Y(a)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}A|\bm{X}, for $a\in\{0,1\}$. Existing sensitivity analysis methods have considered at least three types of relaxations of this assumption:

Pattern-mixture models consider a specific difference between the conditional distribution $Y(a)|\bm{X},A$ and $Y(a)|\bm{X}$ (Robins 1999, Robins 2002, Birmingham et al. 2003, Vansteelandt et al. 2006, Daniels and Hogan 2008, e.g.).

Selection models consider a specific difference between the conditional distribution $A|\bm{X},Y(a)$ and $A|\bm{X}$ (Scharfstein et al. 1999, Gilbert et al. 2003, Gilbert et al. 2013, e.g.).

Rosenbaum’s sensitivity models consider a range of possible selection models so that within a matched set the probabilities of getting treated are no more different than a number $\Lambda$ in odds ratio. A worst case p-value of Fisher’s sharp null hypothesis is then reported (e.g. Rosenbaum 2002c, Chapter 4).

The first two approaches have the desirable property that, under a specified deviation, one can often use existing theory to derive an asymptotically normal (and sometimes efficient) estimator of the causal effect. However, they are arguably more difficult to interpret than Rosenbaum’s sensitivity model because it is impossible to exhaust all possible deviations in this way. Often, one considers just a few functional forms of the deviation and hopes the results of the sensitivity analysis can be extended to “similar” functional forms (Brumback et al. 2004, see e.g.).

Our marginal sensitivity model can be regarded as a hybrid of the selection model and Rosenbaum’s approach, in that we consider a range of possible differences between $A|\bm{X},Y(a)$ and $A|\bm{X}$. This model was considered first introduced by Tan 2006, who noticed that the range of IPW point estimates can be computed by linear programming. However, Tan 2006 did not consider sampling variation of the bounds, thus his method has limited applicability in practice.

2. Comparison with Rosenbaum’s sensitivity analysis

Rosenbaum 1987 proposed to quantify the degree of violation of the MAR/NUC assumption based on the largest odds ratio of $e_{0}(\bm{x},y_{1})$ and $e_{0}(\bm{x},y_{2})$:

Fix a parameter a $\Gamma\geq 1$ which will quantify the degree of violation from the MAR assumption.

For the missing data problem, assume $e(\cdot,\cdot)\in\mathcal{R}(\Gamma)$, where

For the observational studies problem, assume $e_{a}(\cdot,\cdot)\in\mathcal{R}(\Gamma)$, for $a=0,1$.

The proof of B.1 indicates that not all $e_{0}(\bm{X},Y)$ are compatible with the observed data. To compare the marginal sensitivity model $\mathcal{E}$ with Rosenbaum’s model $\mathcal{R}$, we introduce the following concept:

Then Rosenbaum’s sensitivity model and the marginal sensitivity model are related in the following way:

For any $\Lambda\geq 1$, $\mathcal{E}(\sqrt{\Lambda})\subseteq\mathcal{R}(\Lambda)$ and $\mathcal{R}(\Lambda)\cap\mathcal{C}\subseteq\mathcal{E}(\Lambda)\cap\mathcal{C}$.

As mentioned in the Introduction, Rosenbaum and his coauthors obtained point estimate and confidence interval of the causal effect under the collection of sensitivity models $\mathcal{R}(\Gamma)$. To this end, it is often assumed that the causal effect is additive and constant across the individuals, that is, $Y_{i}(1)-Y_{i}(0)\equiv\Delta$, for all $i\in[n]$. Then to determine if an effect $\Delta$ should be included in the $(1-\alpha)$-confidence interval, one just needs to test at level $\alpha$ the Fisher null $H_{0}:\,Y_{i}(0)=Y_{i}(1)$ for all $i\in[n]$, using $Y-\Delta A$ as the outcome (Hodges and Lehmann 1963) and under Rosenbaum’s sensitivity model. We refer the reader to Rosenbaum 2002c for an overview of this approach. Our approach is different from existing methods targeting Rosenbaum’s sensitivity model in many ways, sometimes markedly:

Population: Most if not all existing methods for Rosenbaum’s model treat the observed samples as the population, whereas we treat the observations as i.i.d. samples from a much larger super-population.

Design: Existing methods usually require the data are paired or grouped. Statistical theory assumes the matching is exact, which is usually not strictly enforced in practice. Our approach is based on the IPW estimator and does not require exact matching.

Sensitivity Model: We consider a different but closely related sensitivity model. Rosenbaum’s sensitivity model is most natural for matched designs, whereas the marginal model is most natural when using IPW estimators. We also consider a parametric extension of the marginal sensitivity model.

Statistical Inference: Most existing methods are based on randomization tests of Fisher’s sharp null hypothesis, utilizing the randomness in treatment assignment. Our approach takes a point estimation perspective by trying to estimate the average treatment effect directly. The distinction can be best understood by comparing to the distinction between hypothesis testing and point estimation, or in Ding 2017’s terminology, the subtle difference between Neyman’s null (the average causal effect is zero) and Fisher’s null (the individual causal effects are all zero).

Effect Heterogeneity: Constructing confidence intervals under Rosenbaum’s sensitivity model usually require the causal effect is homogeneous, apart from Rosenbaum 2002a who considered the “attributable effect” of a treatment. Some very recent advancements aim to remove this requirement in randomization inference. Fogarty et al. 2017 and Fogarty 2017 considered estimating the sample ATE in observational studies with a matched pairs design. Our approach inherently allows the causal effect to be heterogeneous.

Applicability to Missing Data Problems: Our approach can be easily applied to missing data problems.

3. Partially Identified Parameter

Our framework is also related to a literature in econometrics on partially identified parameters (Imbens and Manski 2004, Vansteelandt et al. 2006, Chernozhukov et al. 2007, Aronow and Lee 2012, Miratrix et al. 2017). See Richardson et al. 2014 for a recent review. The mean response $\mu$ or the ATE $\Delta$ can be regarded as partially identified under the marginal sensitivity model. In fact, we have adopted the terminology “partially identified region” for the set $\{\mu^{(h)}:e^{(h)}\in\mathcal{E}_{\bm{\beta}_{0}}(\Lambda)\}$.

The main distinction is that existing methods in this literature usually require estimates of the boundaries of the partially identified region with known asymptotic distributions. In Section 4.3.1 we have shown that this is inherently difficult for sensitivity analysis. Our work opens the door for inference of partially identified parameters when it is difficult to analyze the asymptotic behavior of the boundary estimates.

Numerical Examples

over $e^{(h)}(x,y)\in\mathcal{E}_{\beta_{0}}(\Lambda)$. Solving this population optimization problem is possible because $X$ and $Y$ are discrete. We use the percentile bootstrap procedure in Section 4 to construct the interval of point estimates and the confidence interval (nominal level is set to $90\%$). This is repeated for $1000$ times to obtain the non-coverage rate of the confidence intervals. We also report the median interval for the point estimates and the median confidence interval by taking the sample median of the corresponding extrema in the $1000$ repetitions.

The simulation results are reported in Table 1. The confidence intervals constructed by the percentile bootstrap have desired coverage (90%) when $\beta_{A}=0.5$ and do not appear to be conservative, even though we used the generalized minimax inequality in 4.3 to prove the main result. However, when $\beta_{A}=1.5$, the confidence intervals are anti-conservative and the non-coverage rate becomes larger as $\lambda$ increases. This is most likely due to the well known phenomenon that the IPW estimators are unstable when the selection probability is close to $0$ (Kang and Schafer 2007) and the relative small sample size being used ($n=200$). The instability problem is exacerbated as the sensitivity value increases.

Regardless of the value of $\beta_{A}$, the median interval of point estimates is always very close to the true partially identified interval. This shows that the range of SIPW estimates constructed by solving the optimization problem (4.4) is unbiased in the simulation.

2. Real data example

Finally we illustrate the methods proposed in Sections 4 and 6 by an observational study, in which we are interested in estimating the causal effect of fish consumption on the blood mercury level. We obtained $2512$ survey responses from the National Health and Nutrition Examination Survey (NHANES) 2013-2014 who were at least $18$ years old, answered the questionnaire about seafood consumption, and had blood mercury measured. Among these individuals, $1$ has missing education, $7$ have missing smoking, and $175$ have missing income. We removed the individuals with missing education or smoking and imputed the missing income using the median income (we also added a binary indicator for missing income). Then we defined high fish consumption as more than 12 servings of fish or shellfish in the previous month, and low fish consumption as 0 or 1 servings of fish. In the end we were left with 234 treated individuals (high consumption), $873$ controls (low consumption), and $8$ covariates: gender, age, income, whether income is missing, race, education, ever smoked, and number of cigarettes smoked last month. The outcome variable is $\log_{2}$ of total blood mercury (in ug/L). This dataset was also analyzed by Zhao et al. 2017 and is publicly available in the R package CrossScreening on CRAN.

We compared the results with Rosenbaum’s sensitivity analysis as implemented in the senmwCI function (default options) in the R package sensitivitymw (Rosenbaum 2015). We used the $234$ matched pairs created by Zhao et al. 2017 as the basis of the sensitivity analysis. As mentioned previously, Rosenbaum’s sensitivity analysis assumes constant treatment effect (CTE) to construct confidence intervals.

The results of the five sensitivity analyses are reported in Table 2. Alternatively, one can report the results by plotting the confidence intervals against $\Lambda$ (Figure 1). Overall, the confidence intervals constructed by the percentile bootstrap were slightly wider than those constructed by Rosenbaum’s sensitivity analysis under the same $\Lambda$, while the percentile bootstrap intervals under $\sqrt{\Lambda}$ were shorter than Rosenbaum’s under $\Lambda$. This observation is not surprising given 7.1. Augmentation by outcome regressions (SAIPW estimators) helped to reduce the width of confidence intervals when the estimand is ATE, but did not reduce the width when the estimand is ATT. The IPW analyses suggested that the ATE/ATT is significantly positive for at least $\Lambda=2.72$, while the matching analysis found the effect is significantly positive for at least $\Gamma=7.39$. This means that, in order for the ATE/ATT to be non-significant, the estimated propensity score must be quite different from the actual propensity score given all confounders. We consider this as a fairly strong evidence that the qualitative conclusion “consuming fish increases blood mercury level” is somewhat insensitive to unmeasured confounding.

Lastly we want to comment on the computational costs reported in the last column of Table 2. The IPW analyses are slower than the matching analyses (assuming the matches are given), but the running time is quite acceptable given we have over a thousand observations. The reason for the apparent advantage of matching is that we do not include the time for generating the matches. The matching analysis uses analytic approximations to conduct sensitivity analysis, hence it is faster than the IPW analyses which use the bootstrap. Since the time complexity of the IPW analyses scales almost linearly with the sample size, we expect the running times for larger studies will still be acceptable.

Acknowledgement: The authors thank Colin Fogarty for pointing out the difference between Rosenbaum’s sensitivity model and the marginal sensitivity model.

References

Appendix A Proofs

Here we prove 4.1 and 4.3. Additional proofs can be found in the supplementary file.

By definition, under the sensitivity model $\mathcal{H}(\gamma)$, the true data generating distribution $F_{0}$ satisfies $h_{0}(\bm{x},y)\in\mathcal{H}_{\gamma}$. This implies that

The last inequality is true because $[L,U]\supseteq[L^{(h_{0})},U^{(h_{0})}]$. Now, taking limit on both sides gives

since $[L^{(h_{0})},U^{(h_{0})}]$ is an asymptotically $(1-\alpha)$-confidence interval for $\mu^{(h_{0})}$.

A.2. Proof of 4.3

For $1\leq b\leq N$, where $N=n^{n}$ is the total number of possible bootstrap resamples, denote by $\hat{\hat{\mu}}^{(h)}_{b}$ the estimate (4.7) in the $b$-th resample. Therefore, for every $h\in\mathcal{H}(\gamma)$,

Note that both sides of (A.1) above are sequences indexed by $b\in[N]$, and since the inequality is true entry by entry, it is also true for any order statistic of the sequences, namely for any $0<\alpha<1$,

by definition (4.11). Since (A.2) is true for any $h\in\mathcal{H}(\gamma)$, taking infimum on the LHS above gives

The lower bound on $U$ can be proved similarly.

A.3. Proof of 4.5

We prove the first claim by contradiction. Suppose there exists two indices $s_{1}<s_{2}$ such that $A_{s_{1}}=A_{s_{2}}=1$, $Y_{s_{1}}>Y_{s_{2}}$ but $z_{s_{1}}<z_{s_{2}}$. Then consider the following perturbation,

When $\varepsilon>0$ is sufficiently small, $(z_{i}^{\prime})_{i=1}^{m}$ is still feasible but the objective becomes larger, which contradicts the assumption that $(z_{i})_{i=1}^{m}$ is the maximizer.

Next, we prove the second claim. It is well known that if a linear programming is feasible and bounded, then there is at least one vertex (also called the basic feasible) solution (Dantzig 1951). Notice that in (4.4), there are $m+1$ optimization variables and $1$ equality constraint. It is also easy to verify that $t=0$ is not feasible. Therefore, among the $2m$ inequality constraints, $\frac{1}{\Gamma}t\leq\bar{z}_{i}\leq\Gamma t$, $1\leq s\leq m$, there exists a solution of (4.4) such that $m$ equalities hold. This implies that $z_{i}=\frac{1}{t}\bar{z}_{i}$ is either $\Gamma$ or $\frac{1}{\Gamma}$, for all $1\leq s\leq m$. Then, using the first part of the proposition, there exists a $M$ such that $z_{i}=\Gamma$, if $Y_{i}\geq M$, and $z_{i}=\frac{1}{\Gamma}$, if $Y_{i}<M$.

A.4. Proof of 7.1

Now, suppose the function $e\in\mathcal{R}(\Gamma)\cap\mathcal{C}$. By (7.1),

Appendix B Unidentifiability of e0​(x,y)e_{0}(\textit{{x}},y)

Proof: The complete data density can be factorized as

Appendix C Proof of Theorem 4.2

for $h\in\mathcal{H}(\gamma)$. We begin by showing that the estimates of the parameters $(\mu^{(h)},\kappa^{(h)},\bm{\beta}_{0}^{\prime})^{\prime}$ can be derived using the framework of $Z$-estimation.

It is easy to see that the $Z$ estimate $\hat{\mu}^{(h)}$ is exactly the SIPW estimate (4.7) for $\mu^{(h)}$ and $\hat{\bm{\beta}}$ is the MLE of $\bm{\beta}$ for the logistic regression model.

Now, invoking the asymptotic theory of bootstrap for $Z$-estimators (Wellner and Zhan 1996, Kosorok 2006, Chapter 10), we can derive validity of the bootstrap confidence intervals discussed in Section 4.3.2. As a remark, the result of Wellner and Zhan (1996) extend to nonparametric bootstrap of infinite dimensional $Z$-estimators, under a collection of regularity conditions. Therefore, we can expect that the bootstrap procedure introduced above to also work if the missing probability is modeled non-parametrically, if the model satisfies these regularity conditions. However, it is important to note that bootstrap is generally not valid for general non-parametric models, as observed by Abadie and Imbens 2008.

To this end, we need the following assumption.

The parameter space $\Theta$ is compact and the true parameter $\bm{\upsilon}_{0}$ is in the interior of $\Theta$. Moreover, the joint distribution of $(Y,\bm{X})$ satisfies:

The limiting distribution of $\hat{\mu}^{(h)}$ (recall (4.7) and (C.1)) and $\hat{\hat{\mu}}^{(h)}$ (defined through (C.1)) is an immediate consequence of the above theorem:

where $(\sigma^{(h)})^{2}=(\dot{\Phi}_{0}^{-1}\Sigma\dot{\Phi}_{0})_{11}$ is the first diagonal element of $\dot{\Phi}_{0}^{-1}\Sigma\dot{\Phi}_{0}$.

The rest of the section is organized as follows: In Section C.2 we prove Theorem C.1, and in Section C.3 we complete the proof of Theorem 4.2 using Corollary C.2.

C.2. Proof of Theorem C.1

We begin by verifying that the matrices $\dot{\Phi}_{0}^{-1}$ and $\Sigma$ in Theorem C.1 are well-defined. To begin with, a direct computation gives

by Assumption C.1 (2). Therefore, $\dot{\Phi}_{0}$ is invertible. Also, note that $\Sigma<\infty$, which follows by direct multiplication and Assumption C.1 (1).

We can now proceed to prove Theorem C.1. This will be done by invoking (Kosorok 2006, Theorem 10.16), which gives conditions are asymptotic normality of bootstrapped $Z$-estimators. This entails verifying the following three conditions:

$||\Phi(\bm{\upsilon})||_{1}$ is strictly positive outside every open neighborhood of $\bm{\upsilon}_{0}$ (proved in Section C.2.2).

Define the envelope function $B(\bm{t}):=\sup_{\bm{\upsilon}\in\Theta}||Q(\bm{t}|\bm{\upsilon})||_{1}$. Then using the compactness of $\Theta$, $||h||_{\infty}\leq\gamma$, and $|a|\leq 1$,

C.2.2. Proof of (B)

To begin with note that by Assumption C.1 (2),

has a unique root $\bm{\beta}=\bm{\beta}_{0}$, since its gradient is positive definite at $\bm{\beta}_{0}$, and non-negative definite everywhere. Then, fixing $\varepsilon>0$, we have

whenever $||\bm{\beta}-\bm{\beta}_{0}||_{1}>\frac{\varepsilon}{M}$, where $M$ is a constant to be chosen later.

Next, assume that $||\bm{\beta}-\bm{\beta}_{0}||_{1}\leq\frac{\varepsilon}{M}$. This implies $||\bm{\beta}-\bm{\beta}_{0}||_{\infty}\leq\frac{\varepsilon}{M}$ and

by choosing $M\geq 64K\cdot K_{1}(\gamma)$, where

by Assumption C.1, and $K:=\sup_{\bm{\upsilon}\in\Theta}|\nu|\in(0,\infty)$ by the compactness of $\Theta$. Therefore, whenever $||\bm{\beta}-\bm{\beta}_{0}||_{1}\leq\frac{\varepsilon}{M}$ and $|\kappa-\kappa^{(h)}|>\frac{\varepsilon}{4K}$,

Finally, assume that $||\bm{\beta}-\bm{\beta}_{0}||_{1}\leq\frac{\varepsilon}{M}$ and $|\kappa-\kappa^{(h)}|\leq\frac{\varepsilon}{4K}$, but $|\nu-\mu^{(h)}|>\frac{\varepsilon}{2\kappa^{(h)}}$. Then, as in (C.19),

by choosing $M\geq 64K\cdot K_{2}(\gamma)$, where

by Assumption C.1. Using (C.20) and $|\kappa^{(h)}\nu-\kappa\nu|\leq\frac{\varepsilon}{4}$ now gives

Combining (C.17), (C.19), and (C.21), we get, for all $\delta>0$, $\inf\{||\Phi(\bm{\upsilon})||^{2}:||\bm{\upsilon}-\bm{\upsilon}_{0}||_{1}>\delta\}>0$, as required.

C.2.3. Proof of (C)

where $M_{1}(\bm{x})=||\bm{x}||_{1}^{2}$. Next, observe that

where $M_{2}(\bm{x})=1+||\bm{x}||_{1}\sup_{\bm{\beta}\in\Theta_{0}}e^{-\bm{\beta}^{\prime}\bm{x}}$ and $\Theta_{0}$ is the projection of the parameter space to the last $d$ coordinates. Finally,

where $M_{3}(\bm{x},y)=1+||y\bm{x}||_{1}\sup_{\bm{\beta}\in\Theta_{0}}e^{-\bm{\beta}^{\prime}\bm{x}}$.

Therefore, defining $M(\bm{x},y)=M_{1}(\bm{x})+M_{2}(\bm{x})+M_{3}(\bm{x},y)$ and combining (C.2.3), (C.23), and (C.2.3) gives

The inequalities in (C.2.3), (C.23), and (C.2.3), combined with Assumption C.1 also implies that

holds coordinate-wise, whenever $||\bm{\upsilon}_{n}-\bm{\upsilon}_{0}||_{1}\rightarrow 0$.

C.3. Completing the Proof of Theorem 4.2

Appendix D Proof of Corollary 5.1

The proof of Corollary 5.1 is similar to the proof of Theorem 4.4. We begin by defining

It is easy to see that the $Z$ estimate $\hat{\mu}^{(h_{1})}(1)-\hat{\mu}^{(h_{0})}(0)$ is exactly the SIPW estimate for $\Delta^{(h_{0},h_{1})}$ and $\hat{\bm{\beta}}$ is the MLE of $\bm{\beta}$ for the logistic regression model. Moreover, as in (C.1), the bootstrap $Z$-estimates $\hat{\hat{\bm{\upsilon}}}$ are obtained from the equations:

Then, as in Theorem C.1, the limiting joint normality of $(\hat{\mu}^{(h_{0})}(0),\hat{\mu}^{(h_{1})}(1))$, and hence the asymptotic normality of $\hat{\Delta}^{(h_{0},h_{1})}=\hat{\mu}^{(h_{0})}(0)-\hat{\mu}^{(h_{1})}(1)$, can be derived. This would imply the asymptotic validity of the confidence interval (5.2), as in Section C.3, completing the proof of Corollary 5.1. Details are omitted.