SENSITIVITY ANALYSIS FOR INVERSE PROBABILITY WEIGHTING
ESTIMATORS VIA THE PERCENTILE BOOTSTRAP

arXiv:1711.11286v1 [stat.ME] 30 Nov 2017

QINGYUAN ZHAO, DYLAN S. SMALL AND BHASWAR B. BHATTACHARYA
Department of Statistics, The Wharton School, University of Pennsylvania

Abstract. To identify the estimand in missing data problems and observational studies, it
is common to base the statistical estimation on the “missing at random” and “no unmeasured confounder” assumptions. However, these assumptions are unverifiable using empirical
data and pose serious threats to the validity of the qualitative conclusions of the statistical
inference. A sensitivity analysis asks how the conclusions may change if the unverifiable
assumptions are violated to a certain degree. In this paper we consider a marginal sensitivity model which is a natural extension of Rosenbaum’s sensitivity model for matched
observational studies. We aim to construct confidence intervals based on inverse probability
weighting estimators, such that the intervals have asymptotically nominal coverage of the
estimand whenever the data generating distribution is in the collection of marginal sensitivity models. We use a percentile bootstrap and a generalized minimax/maximin inequality
to transform this intractable problem to a linear fractional programming problem, which
can be solved very efficiently. We illustrate our method using a real dataset to estimate the
causal effect of fish consumption on blood mercury level.

1. 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
E-mail address: {qyzhao,dsmall,bhaswar}@wharton.upenn.edu.
Date: December 1, 2017.
1

2

SENSITIVITY ANALYSIS FOR IPW

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.
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,c). 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 IPWtype 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, Section 4.2),
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. 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

SENSITIVITY ANALYSIS FOR IPW

3

the estimand (the mean response vector in missing data problems and the average treatment effect in observational studies), with asymptotically 1 − α coverage probability, for all
data generating distributions which violate the MAR assumption within this threshold (see
Definition 3).
The obvious approach to compute such a confidence interval 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. We circumvent this problem by using the percentile bootstrap, reducing the
problem to solving many linear fractional programs. The highlights of our strategy and the
summary of the results are given below:

(1) To obtain a confidence interval under the marginal sensitivity model, we simply maximize/minimize the IPW point estimates over the marginal sensitivity model over B
bootstrap resamples of the data. Then we report the α2 -sample quantile of the B minima and (1 − α2 )-sample quantile of the B maxima, as the lower and upper end-points
of the confidence interval, respectively (Section 4.3).
(2) The asymptotic (1 − α) coverage of this interval follows from the limiting normal
distribution of the IPW estimators and a generalized minimax/maximin inequality,
which justifies the interchange of quantile and infimum/supremum (Theorem 4 and
Corollary 6). In fact, our interval covers the entire partially identified region of the
estimand with probability tending to (1 − α).
(3) 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
5).
(4) 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).

The rest of this paper is organized as follows. We introduce the necessary notations in
Section 2 and the marginal sensitivity model in Section 3. Our method for constructing confidence 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. Other extensions, including the
use of augmented IPW estimators, are considered in Section 6. We review some related sensitivity analysis methods and compare our framework with Rosenbaum’s sensitivity analysis
in Section 7. In Section 8 we illustrate our method using a real data example and compare
it with Rosenbaum’s sensitivity analysis.

4

SENSITIVITY ANALYSIS FOR IPW

2. Background
In this paper, we consider sensitivity analysis in two closely related problems: estimation
of the mean response with missing data and estimation of the average treatment effect in
observational studies:
(1) In the missing data problem, we assume (A1 , X1 , Y1 ), (A2 , X2 , Y2 ), . . . (An , Xn , Yn )
are i.i.d. from a joint distribution F0 , where for each subject i ∈ [n] := {1, 2, . . . , n},
Ai is the indicator of non-missing response, Xi ∈ X ⊂ Rd is a vector of covariates,
and Yi ∈ R is the response (observed only if Ai = 1). In other words, we only observe
(Ai , Xi , Ai Yi ) for i ∈ [n]. Based on the observed data, our goal is to estimate the mean
response µ := E0 [Y ] in the complete data, where (A, X, Y ) ∼ F0 and E0 indicates
that the expectation is taken over the true data generating distribution F0 .
(2) In observational studies, we observe i.i.d. (A1 , X1 , Y1 ), (A2 , X2 , Y2 ), . . . , (An , Xn , Yn ),
where for each subject i ∈ [n], Ai is a binary treatment indicator (which is 1 if treated
and 0 in control), Xi ∈ X ⊂ Rd is a vector of measured confounders, and
Yi = Yi (Ai ) = Ai Yi (1) + (1 − Ai )Yi (0)
is the outcome. Furthermore we assume (Ai , Xi , Yi (0), Yi (1)) are i.i.d. from a joint
distribution F0 . Here, we are using Rubin (1974)’s potential outcome notation and
have assumed the stable unit treatment value assumption (Rubin, 1980). The goal
is to estimate the average treatment effect (ATE), ∆ := E0 [Y (1)] − E0 [Y (0)]. As
before, E0 means the expectation is taken over the data generating distribution F0 .
Notice that Yi (0) and Yi (1) are never observed together, which is often referred to as
the “fundamental problem of causal inference” (Holland, 1986). With the potential
outcome notation, the observational studies problem is essentially two missing data
problems: use (Ai , Xi , Ai Yi ) = (Ai , Xi , Ai Yi (1)) to estimate µ(1) := E0 [Y (1)], and
use (1 − Ai , Xi , (1 − Ai )Yi ) = (1 − Ai , Xi , (1 − Ai )Yi (0)) to estimate µ(0) := E0 [Y (0)].
Since some responses or potential outcomes are not observed, the estimands µ and ∆
defined above are 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:
Assumption 1 (Missing at random (MAR)). A ⊥
⊥ Y |X under F0 .
The corresponding assumption for observational studies is that the potential outcomes
(Y (0), Y (1)) are independent of the treatment A given the covariates X (Rosenbaum and
Rubin, 1983):
Assumption 2 (Strong ignorability or no unmeasured confounder (NUC)). A ⊥
⊥ (Y (0), Y (1))|X
under F0 .
With the additional assumption that no subject is missing or receives treatment/control
with probability 1, the parameters µ and ∆ are identifiable from the data.
Assumption 3 (Overlap). For x ∈ X ,
(1) In the missing data problem, e0 (x) := P0 (A = 1|X = x) ∈ (0, 1].
(2) In observational studies, e0 (x) ∈ (0, 1).
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 e0 (x) (often called the propensity score in

SENSITIVITY ANALYSIS FOR IPW

5

observational studies). One distinguished example is the inverse probability weighting (IPW)
estimator which dates back to Horvitz and Thompson (1952),
n

µ̂IPW

1 X Ai Yi
=
n
ê(Xi )
s=1

n

X Ai Yi
(1 − Ai )Yi
ˆ IPW = 1
and ∆
−
,
n
ê(Xi )
1 − ê(Xi )
s=1

where ê(X) is a sample-estimate of e0 (X). Observe that under Assumptions 1 and 3,






AY
AY
AY
µ = E0
= E0
= E0
,
P0 (A = 1|X, Y )
P0 (A = 1|X)
e0 (X)

(1)

which implies that µ̂IPW can consistently estimate the mean response µ if ê converges to e0 .
Notice that the first equality in (1) is due to the tower property of conditional expectation
and is always true if the denominator P0 (A = 1|X, Y ) is positive with probability 1. The
second equality in (1) uses P0 (A = 1|X, Y ) = P0 (A = 1|X) 6= 0 by Assumptions 1 and 3.
ˆ IPW can consistently estimate ∆ under Assumptions 2 and 3.
Similarly, ∆
3. Sensitivity models
A critical issue of the above approach is that the MAR and NUC assumptions are not verifiable using the data because some responses/potential outcomes are not observed. Therefore,
ˆ IPW is biased and
if the MAR or NUC assumption is violated, the IPW estimator µ̂IPW or ∆
the confidence interval based on it (and any other existing estimator that assumes MAR or
NUC) does not cover µ or ∆ at the nominal rate. In practice, statisticians often hope the data
analyst can use her subject knowledge to justify the substantive reasonableness of MAR/NUC
assumptions (Little and Rubin, 2014). However, speaking realistically, the MAR and NUC
assumptions are perhaps never strictly satisfied, so the statistical inference is always subject
to criticism of insufficient control of confounding.
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/NUC assumptions. This is clearly
an important question, and, in fact, sensitivity analysis is often suggested or required for
empirical publications.1 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,
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.
We begin with a description of the marginal sensitivity model used in this paper, which is
also considered by Tan (2006). To this end, with a slight abuse of notation, let e0 (x, y) =
P0 (A = 1|X = x, Y = y). Notice that, when showing the IPW estimator is consistent in (1),
the MAR assumption is useful because it implies that
e0 (x, y) = e0 (x).

(2)

Similarly, the NUC assumption in the observational studies problem implies that ea (x, y) =
ea (x), where ea (x, y) := P0 (A = 1|X = x, Y (a) = y), for a ∈ {0, 1}.
1For 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).

6

SENSITIVITY ANALYSIS FOR IPW

When the MAR assumption is violated, equation (2) is no longer valid. The following wellknown proposition (proof included in Appendix A.1 for completeness) shows that e0 (x, y) is
generally not identifiable from the data without the MAR assumption. For this reason we
shall refer to a user-specified function e0 (x, y) as a sensitivity model.
Proposition 1. In the missing data problem and assuming Assumption 3 holds, e0 (x, y) =
P0 (A = 1|X = x, Y = y) is not identifiable from the data.
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 odd ratio of the “complete data”
selection probability e0 (x, y) and the “observed data” selection probability e0 (x).
Definition 1 (Marginal Sensitivity Model). Fix a parameter Λ ≥ 1.
(1) For the missing data problem, we assume e(x, y) ∈ E(Λ), where
n
o
1
E(Λ) = e(x, y) :
≤ OR(e(x, y), e0 (x)) ≤ Λ, for all x ∈ X , y ∈ R .
Λ

(3)

(2) For the observational studies problem, we assume ea (x, y) ∈ E(Λ) for a = 0, 1.
The relationship between this sensitivity model and Rosenbaum’s sensitivity model is examined in Section 7.2.
Remark 1. To understand Definition 1, it can be conceptually easier to imagine that there
is an unobserved variable U that “summarizes” all unmeasured confounding. Definition 4
means that the odds ratio between P0 (A|X, U = u1 ) and P0 (A|X, U = u2 ) is always bounded
by Λ−1 and Λ, but the relation between U and the potential outcomes are not constrained
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 Robin’s suggestion to simplify the notation.
Remark 2. It is often convenient to use the logistic representation of the marginal sensitivity
model (3). For simplicity, let us consider the missing data problem. Denote
g0 (x) = logit(e0 (x)) = log

e0 (x)
1 − e0 (x)

and g0 (x, y) = logit(e0 (x, y)),

(4)

and h0 (x, y) = g0 (x) − g0 (x, y) be the logit-scale difference of the observed data selection
probability and the complete data selection probability. If we write further introduce the
notation e(h) (x, y) = [1 + exp(h(x, y) − g0 (x))]−1 , then
E(Λ) = {e(h) (x, y) : h ∈ H(λ)}, where H(λ) = {h : X × R → R and ||h||∞ ≤ λ} ,

(5)

that is the marginal sensitivity model puts bound on the L∞ -norm of h.
4. Confidence Interval for the Mean Response
The primary goal of this paper is to construct confidence intervals for µ and ∆ under a
specified collection of sensitivity models. For simplicity, we focus on the missing data problem
in this section. The simple extension to observational studies is described in Section 5.

SENSITIVITY ANALYSIS FOR IPW

7

4.1. Confidence interval in sensitivity analysis. We begin by briefly extending the definition of sensitivity models. In Definition 1, a postulated sensitivity model e0 (x, y) is compared with the data-identifiable e0 (x) = P0 (A = 1|X = x). This model is most appropriate
when e0 (X) is estimated non-parametrically using the data, but this is only possible for
low-dimensional problems due to the curse of dimensionality. More often, for example, in
propensity score matching or IPW, e0 (x) is estimated by a parametric model. In this case
it is more appropriate to compare the postulated sensitivity model e0 (x, y) with the best
parametric approximation to e0 (x):

eβ0 (x) = arg minβ∈Θ KL P0 (A = 1|X = x) Pβ (A = 1|X = x)
= arg maxβ∈Θ E0 [A · log eβ (X) + (1 − A) · log(1 − eβ (X))|X = x] ,

(6)

where KL stands for the Kullback-Leibler divergence and {eβ (x) := Pβ (A = 1|X = x), β ∈
Θ ⊂ Rd } is a family of parametric models for the missingness probability. We recommend to
use the logistic regression model
0

eβ x
eβ (x) =
, or equivalently gβ (x) = logit(eβ (x)) = β 0 x,
(7)
1 + eβ0 x
which works seamlessly with our sensitivity model because the degree of violation of MAR is
quantified by odds ratios. However our framework can be easily applied to other parametric
models eβ (x) (e.g. different links).
Definition 2 (Parametric Marginal Sensitivity Model). Fix a parameter a Λ ≥ 1, this
collection of sensitivity models assumes the true missingness probability e0 (x, y) = P0 (A =
1|X = x, Y = y) satisfies
o
n
1
≤ OR(e(x, y), eβ0 (x)) ≤ Λ for all x ∈ X , y ∈ R . (8)
e0 (x, y) ∈ Eβ0 (Λ) := e(x, y) :
Λ
As in (5), the constraint (8) is equivalent to hβ0 (x, y) ∈ H(λ) for hβ (x, y) = gβ (x) − g0 (x, y)
and λ = log Λ.
Next we define the estimand (mean response) under a specific sensitity model h:
 
−1 h
i

A
(h)
E0 AY 1 + eh(X,Y )−gβ0 (X)
µ = E0 (h)
e (X, Y )
 
−1 

A
AY
= E0 (h)
E0 (h)
e (X, Y )
e (X, Y )

(9)

where e(h) (x, y) = (1 + eh(x,y)−gβ0 (x))−1 . By the definition of hβ and the first equality in (1),
it is easy to see that E0

A

e

(hβ )
0 (X,Y

)

= 1 and µ(hβ0 ) are equal to the true mean response µ.

The set {µ(h) : e(h) ∈ Eβ0 (λ)} will be referred to as the partially identifiable region under
Eβ0 (Λ); see Section 7.3 for more discussion.
Remark 3. In defining µ(h) and e(h) , we slightly abused the notation and did not indicate that
they also depend on the “baseline” choice of parametric or nonparametric model of e0 (X).
This is because the choice of the sensitivity model is a natural consequence of the choice of
the working model for e0 (X). If e0 (X) is modeled nonparametrically (or parametrically),
then we should use the collection of nonparametric (or parametric) sensitivity models E(Λ)
(or Eβ0 (Λ)). For notational simplicity, hereafter we will often refer to h(x, y) ∈ H(λ) instead
of e(h) (x, y) as the sensitivity model, since the former does not depends on the choice of

8

SENSITIVITY ANALYSIS FOR IPW

the working model for e0 (X). Consequently we will also call H(λ) a collection of sensitivity
models without specifying which parametric/nonparametric model was used for e0 (X).
Remark 4. Compared to Definition 4, the only difference in (8) is that the postulated
sensitivity model e0 (x, y) is compared to the parametric model eβ0 (x) instead of the nonparametric probability e0 (x). In other words, the parametric sensitivity model considers
both
(1) Model misspecification, that is, eβ0 (x) 6= e0 (x); and
(2) Missing not at random, that is, e0 (x) 6= e0 (x, y).
This is arguably a desirable feature, because the term “sensitivity analysis” is also widely
used as the analysis of an empirical study’s robustness to parametric modeling assumptions.
With these notations, we can now define what is meant by a confidence interval under a
collection of sensitivity models:
Definition 3. A data-dependent interval [L, U ] is called a (1 − α)-confidence interval for the
mean response µ under the collection of sensitivity models H(λ) (may corresponds to E(Λ)
or Eβ0 (Λ), see Remark 3), if
P0 (µ ∈ [L, U ]) ≥ 1 − α
(10)
is true for any data generating distribution F0 such that h0 ∈ H(λ) or hβ0 ∈ H(λ), depending
on whether E(Λ) or Eβ0 (Λ) is used. The interval is said to be an asymptotic (1−α)-confidence
interval, if lim inf n→∞ P0 (µ ∈ [L, U ]) ≥ 1 − α.
4.2. The IPW Point Estimates. Intuitively, a confidence interval [L, U ] as in Definition
3 must at least include a point estimate of µ(h) for every h ∈ H(λ). To this end, let h be a
postulated sensitivity model. The corresponding missingness probability e(h) (x, y) can then
be estimated by
1
,
(11)
ê(h) (x, y) =
h(x,y)−ĝ(x,y)
1+e
where
(
logit(P̂(A = 1|X = x))
if e(h) (x, y) ∈ E(Λ)
(12)
ĝ(x) =
logit(Pβ̂ (A = 1|X = x)) = β̂ 0 x if e(h) (x, y) ∈ Eβ0 (Λ).
We consider two point estimates of µ̃(h) or µ(h) , namely the IPW estimate and a stabilized
version of it:
(1) The IPW estimator of µ̃(h) is
n
1X
Ai Yi
(h)
µ̂IPW =
,
(h) (X , Y )
n
ê
i
i
i=1
where ê(h) is as in (11).
(2) It is well known that, even under the MAR assumption, the IPW estimator can be
unstable when the missingness probability e0 (x) (or the parametric approximation
eβ0 (x)) is close to 0 for some x ∈ X (Kang and Schafer, 2007). To alleviate this
issue, the stabilized IPW (SIPW) estimator, obtained by normalizing the weights, is
often used in practice:
 X
−1  X

n
n
1
Ai
1
Ai Yi
µ̂(h) =
.
(13)
(h) (X , Y )
(h) (X , Y )
n
n
ê
ê
i
i
i
i
i=1
i=1
It is easy to see that µ̂(h) estimates µ(h) .

SENSITIVITY ANALYSIS FOR IPW

9

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


(h)
µ̂ ∈ min Yi , max Yi .
i:Ai =1

i:Ai =1

This property is even more desirable in sensitivity analysis because e(h) P
is almost always not
1
the true missingness probability, so the total unnormalized weights n ni=1 Ai /ê(h) (Xi , Yi )
can be very different from 1. For this reason, we will use the SIPW estimator in the remainder
of this paper.
Heuristically, the confidence interval [L, U ] should at least contain the range of SIPW
point estimates, inf h∈H(λ) µ̂(h) , suph∈H(λ) µ̂(h) . 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 is given.
4.3. Constructing the Confidence Interval. To construct a confidence interval for µ, 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 H(λ). In this
case, optimizing the estimated asymptotic variance over H(λ) is generally computationally
intractable, as explained below.
4.3.1. The Union Method. We begin by showing how individual confidence intervals of µ(h)
for h ∈ H(λ) can be combined into a confidence interval in sensitivity analysis.
Proposition 2. Suppose there exists data-dependent intervals [L(h) , U (h) ] such that
lim inf P0 (µ(h) ∈ [L(h) , U (h) ]) ≥ 1 − α
n→∞

holds for every h ∈ H(λ).
(1) Let L = inf h∈H(λ) L(h) , U = suph∈H(λ) U (h) . Then [L, U ] is an asymptotic (1 − α)confidence interval of µ under the collection of sensitivity models H(λ).
(2) Moreover, if there exists α0 ∈ [0, α] (not depending on h) such that


lim sup P0 µ(h) < L(h) ≤ α0 and lim sup P0 µ(h) > U (h) ≤ α − α0 ,
(14)
n→∞

n→∞

for all h ∈ H(λ), then the union interval [L, U ] covers the partially identified region
with probability at least 1 − α,


lim inf P0 {µ(h) : h ∈ H(λ)} ⊆ [L, U ] ≥ 1 − α.
n→∞

Proposition 2 suggests the following way to construct a confidence interval for µ in sen(h)
sitivity analysis using the asymptotic distribution of µ̂1 . Using the general theory of Zestimation, it is not difficult to establish that

√  (h)
D
n µ̂ − µ(h) → N (0, (σ (h) )2 ),
(see, for example, Lunceford and Davidian (2004) or Corollary 9 in the Appendix). Then
using the sandwich variance estimator (σ̂ (h) )2 , an asymptotically (1 − α)-confidence interval
of µ(h) is
"
#
σ̂ (h) (h)
σ̂ (h)
(h)
(h)
(h)
,
[Lsand , Usand ] = µ̂ − z α2 · √ , µ̂ + z α2 · √
n
n

10

SENSITIVITY ANALYSIS FOR IPW

where z α2 is the upper α2 -quantile. Then, by Proposition 2, an asymptotically (1 − α)confidence interval under the collection of sensitivity models is [Lsand , Usand ],
Lsand =

(h)

inf Lsand ,

h∈H(λ)

(h)

Usand = sup Usand .
h∈H(λ)

σ̂ (h)

However, the standard error
is a very complicated function of h (see Corollary 9) and
numerical optimization over h ∈ H(λ) is practically infeasible.
4.3.2. The Percentile Bootstrap. The centerpiece of our proposal is to use the percentile
bootstrap to construct [L(h) , U (h) ]. Next we introduce the necessary notations to describe
the bootstrap procedure. Let Pn be the empirical measure of the sample T1 , T2 , . . . , Tn ,
where Ti = (Ai , Xi0 , Ai Yi ), and Tˆ1 , Tˆ2 , . . . , Tˆn be i.i.d. resamples from the empirical measure.
ˆ(h) the SIPW estimate (13) computing using the bootstrap resamples {T̂i }i∈[n] . For
Let µ̂
h ∈ H(λ), the percentile bootstrap confidence interval of µ(h) is given by
i
h
ˆ(h) ), Q1− α (µ̂
ˆ(h) ) ,
(15)
[L(h) , U (h) ] = Q α2 (µ̂
2
ˆ is the α-percentile of µ̂
ˆ in the bootstrap distribution, that is,
where Qα (µ̂)
ˆ := inf{t : P̂n (µ̂
ˆ ≤ t) ≥ α},
Qα (µ̂)
where P̂n is the bootstrap resampling distribution.
We begin by showing that the percentile bootstrap interval [L(h) , U (h) ] is an asymptotically
valid confidence interval of µ(h) for the parametric sensitivity model e(h) ∈ Eβ0 (Λ). Notice that
bootstrap is generally not valid if the missingness probability is modeled nonparametrically
(Abadie and Imbens, 2008).
Theorem 3 (Validity of the Percentile Bootstrap). In the logistic model (7) and under
Assumption 4 in Appendix B, for every e(h) ∈ Eβ0 (Λ) we have
 α
 α


and lim sup P0 µ(h) > U (h) ≤ .
lim sup P0 µ(h) < L(h) ≤
2
2
n→∞
n→∞
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 B.
Our percentile bootstrap confidence interval under the collection of sensitivity models
Eβ0 (Λ) is given by [L, U ] where




ˆ(h)
ˆ(h) .
and U = Q1− α2
L = Q α2
inf µ̂
sup µ̂
(16)
b
h∈H(λ)

h∈H(λ)

The important thing to observe here is that the infimum/supremum is inside the quantile
function in (16), 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).
Lemma 1. Let L, U be as defined in (17). Then
L≤

inf L(h)
h∈H(λ)

and

U ≥ sup U (h) .
h∈H(λ)

Asymptotic validity of the confidence interval [L, U ] in Equation (16) then immediately
follows from the validity of the union method (Proposition 2), the validity of the percentile
bootstrap (Theorem 3), and Lemma 1.

SENSITIVITY ANALYSIS FOR IPW

11

Theorem 4. Under the same assumptions as in Theorem 3, [L, U ] is an asymptotic (1 − α)confidence interval of the mean response µ, under the collection of sensitivity models Eβ0 (Λ).
Furthermore, [L, U ] covers the partially identified region {µ(h) : e(h) ∈ Eβ0 (Λ)} with probability
at least 1 − α.
4.4. Range of SIPW Point Estimates: Linear Fractional Programming. Theorem
4 transformed the sensitivity analysis problem to computing the extrema of the SIPW point
ˆ(h) and sup
ˆ(h)
estimates, inf h∈H(λ) µ̂
h∈H(λ) µ̂b . In practice, we only repeat this over B ( N )
b
random resamples and compute the interval by






(h)
(h)
ˆ
ˆ
LB = Q α2
inf µ̂b
, UB = Q1− α2
sup µ̂b
.
(17)
h∈H(λ)

b∈[B]

h∈H(λ)

b∈[B]

For notational simplicity, below we consider how to compute the extrema inf h∈H(λ) µ̂(h) and
suph∈H(λ) µ̂(h) using the full observed data instead of the resampled data. Recalling (12) and
(13), this is equivalent to solving

Pn
−ĝ(Xi )


A
Y
1
+
z
e
i
i
i
 subject to zi ∈ Λ−1 , Λ , for all i ∈ [n],
max or min Pi=1
(18)
n
−ĝ(X
)
i
i=1 Ai 1 + zi e
where the optimization variables are zi = eh(Xi ,Yi ) for i ∈ [n]. All the other variables are
observed or can be estimated from the data. Notice that ĝ(x) needs to be re-estimated in
every bootstrap resample. Without loss of generality, assume that the first 1 ≤ m < n
responses are observed, that is A1 = A2 = · · · Am = 1 and Am+1 = · · · = An = 0, and
suppose that the observed responses are in decreasing order, Y1 ≥ Y2 ≥ · · · ≥ Ym . Then (18)
simplifies to

Pm
−ĝ(Xi )


i=1 Yi 1 + zi e
 , subject to zi ∈ Λ−1 , Λ , 1 ≤ i ≤ m.
max or min Pm
(19)
−ĝ(X
)
i
i=1 1 + zi e
This optimization problem is the ratio of two linear functions of the decision variables z,
hence called a linear fractional programming. It can be transformed to linear programming
by the Charnes-Cooper transformation (Charnes and Cooper, 1962). Denote
zi
1
, for 1 ≤ s ≤ m, and t = Pm
.
−ĝ(X
)
−ĝ(Xi ) )
i )
i=1 (1 + zi e
i=1 (1 + zi e

z̄i = Pm

This translates (19) into the following linear programming:
X

m
m
X
−ĝ(Xi )
max or min
Yi e
z̄i + t
Yi
i=1

i=1

subject to t ≥ 0, Λ
m
X
i=1

−1

t ≤ z̄i ≤ Λt, for 1 ≤ i ≤ m,
X

m
e−ĝ(Xi ) z̄i + t
Ai = 1.

(20)

i=1

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 (20) must
have the same or opposite order as the outcomes Y , which enables even faster computation
of (20).

12

SENSITIVITY ANALYSIS FOR IPW

m
Proposition 5. Suppose (zi )m
i=1 solves the maximization problem in (20). Then (zi )i=1 has
m
the same order as (Yi )i=1 , that is, if Ys1 > Ys2 then zs1 > zs2 , for 1 ≤ s1 6= s2 ≤ m.
Furthermore, there exists a solution (zi )m
i=1 and a threshold M such that zi = Λ, if Yi ≥ M ,
and zi = Λ1 , if Yi < M . The same conclusion holds for the minimizer in (20) with (Yi )m
i=1
replaced by (−Yi )m
s=1 .

Proposition 5 implies that we only need to compute the objective of (19) for at most m
choices of (zi )i∈[m] by enumerating the index where it changes from Λ to Λ−1 :


1
(zi )i∈[m] : ∃ a ∈ [m] with zi = Λ, for 1 ≤ i ≤ a, and zi = , for a + 1 ≤ i ≤ m . (21)
Λ
P
Pm
It is easy to see that we only need O(m) time to compute i=1 Yi (1+zi e−ĝ(Xi ) ) and m
i=1 (1+
−ĝ(X
)
i
zi e
), for all the m choices in (21). Hence, the computational complexity to solve the
linear fractional programming (19) is O(m). This is the best possible rate since it takes O(m)
time to just compute the objective once.
ˆ(h)
The above discussion suggests that the interval [LB , UB ], which entails finding inf h∈H(λ) µ̂
b
ˆ(h) for every 1 ≤ b ≤ B, can be computed extremely efficiently. The compuand sup
µ̂
h∈H(λ)

b

tational 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 µ
, which is recommended in practice by Austin (2016). 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.
5. Confidence Interval for the ATE in the Sensitivity Model
As discussed in Section 2, an observational study is essentially two missing data problems,
so the framework developed above can be easily extended. Recall that ea (x, y) = P0 (A =
1|X = x, Y (a) = y). Suppose we use a parametric model eβ0 (x) to model the propensity
score P0 (A = 1|X = x) using the observed covariates. The parametric sensitivity model
assumes
1
≤ OR(ea (x, y), eβ0 (x)) ≤ Λ, for all x ∈ X , y ∈ R, a ∈ {0, 1}.
(22)
Λ
(β)
Let ga (x, y) = logit(ea (x, y)), gβ (x) = logit(eβ (x)) and ha (x, y) = gβ (x) − ga (x, y), for
a ∈ {0, 1}. Then (22) is equivalent to assuming ha,β0 ∈ H(λ) for a ∈ {0, 1}.
As in (9) define, for h ∈ H(λ),
"
#!−1 "
#
Aa (1 − A)1−a
Aa (1 − A)1−a Y
(h)
µ (a) = E
E0
,
(23)
(h)
(h)
ea (X, Y )
ea (X, Y )

−1
a+1
(h)
where ea (x, y) = 1 + e(−1) h(x,y)−logit(Pβ0 (A=a|X=x)) . Denote ∆(h0 ,h1 ) := µ(h1 ) (1) −
(β0 ) (β0 )
,h1 )

µ(h0 ) (0). It is straightforward to show that ∆(h0
(h1 )

Eβ0 (Λ) := {(e0

(h2 )

(x, y), e1

= ∆. For this reason, as before, let

(x, y)) : h1 , h2 ∈ H(λ)}

(24)

be the collection of parametric sensitivity models in the observational studies problem. These
can be estimated as in (11) by
1
(25)
ê(h)
a+1 [h(x,y)−ĝ(x)] ,
a (x, y) =
(−1)
1+e

SENSITIVITY ANALYSIS FOR IPW

13

where ĝ(x) = gβ̂ (x) is the estimated propensity score. Using this we can define the SIPW
estimate of the average treatment effect as follows:
ˆ (h0 ,h1 ) := µ̂(h1 ) (1) − µ̂(h0 ) (0),
∆

(26)

where
µ̂

(h1 )

−1  X

 X
n
n
1
1
Ai
Ai Yi
,
(1) =
(h )
(h )
n
n
ê 1 (Xi , Yi )
ê 1 (Xi , Yi )
i=1

i=1

1

1

µ̂(h1 ) (0)

and
defined analogously.
Now, as in Section 4.3.2, using the percentile bootstrap we can obtain a asymptotically
valid interval for ∆(h0 ,h1 ) :
h
i
ˆˆ (h0 ,h1 )
ˆˆ (h0 ,h1 )
Q α2 (∆
), Q1− α2 (∆
) ,
ˆ
ˆ (h0 ,h1 ) ) is the α-th bootstrap quantile of the SIPW estimates (26), as defined in
where Q α2 (∆
Section 4.3.2. Then, interchanging the maximum/minimum and the quantile as in Theorem
6, we get a confidence interval for the average-treatment effect for the collection of sensitivity
models (24).
Corollary 6. Under the same assumptions in Theorem 3,




ˆˆ (h0 ,h1 )
Q α2
inf
∆
, Q1− α2
inf
h0 ,h1 ∈H(λ)

h0 ,h1 ∈H(λ)

ˆˆ (h0 ,h1 )
∆


(27)

is an asymptotic (1 − α)-confidence interval of the average treatment effect ∆ under the
collection of parametric sensitivity models (24).
The interval in (27) can be computed efficiently using linear fractional programming as in
Section 4.4. To simplify notation, assume, without loss of generality, the first m ≤ 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 (19), computing
the interval (27) is equivalent to solving the following optimization problem:
Pn
Pm
ĝ(Xi ) )
−ĝ(Xi ) )
i=m+1 Yi (1 + zi e
i=1 Yi (1 + zi e
Pm
P
−
maximize or minimize
n
−ĝ(Xi ) )
ĝ(Xi ) )
i=1 (1 + zi e
i=m+1 (1 + zi e
(28)
1
subject to
≤ zi ≤ Λ, for 1 ≤ i ≤ n,
Λ
where zi = eh1 (Xi ,Yi ) , for 1 ≤ i ≤ m and zi = e−h0 (Xi ,Yi ) , for m + 1 ≤ i ≤ n. Note that
n
the variables (zi )m
i=1 and (zi )i=m+1 are separable in (28), so we can solve the maximization/minimization problem in (28) by solving one maximization/minimization problem for
n
(zi )m
i=1 and one minimization/maximization problem for (zi )i=m+1 . Therefore, similar to the
ˆ
missing data problem, the time complexity to obtain the range of the SIPW estimates ∆,
over a range of B bootstrap resamples, is only O(nB + n log n).
6. Extensions
In this section we discuss three extensions of the general framework described in Section 4.

14

SENSITIVITY ANALYSIS FOR IPW

6.1. Mean of the Non-Respondents and Average Treatment Effect on the Treated.
In many applications, it is also interesting to estimate the mean of the non-respondents µ0 =
E0 [Y |A = 0], and the average treatment effect on the treated (ATT) ∆1 = E0 [Y (1)−Y (0)|A =
1]. The method described in Section 4 can be easily applied to these estimands, as described
below.
To begin with note that
h
i
1−e0 (X,Y )
E
AY
0
e0 (X,Y )
E0 [Y · 1{A = 0}]
E0 [(1 − e0 (X, Y ))Y ]
µ0 = E0 [Y |A = 0] =
=
=
.
P0 (A = 0)
P0 (A = 0)
P0 (A = 0)
Then as in Theorem 4, under the collection of parametric sensitivity models Eβ0 (Λ),

 




(h)
ˆ(h)
ˆ
α
Q α2
inf µ̂
inf
,
Q
µ̂
,
1− 2
0b
0b
h∈H(λ)

b∈[B]

h∈H(λ)

b∈[B]

(29)

is an asymptotic (1 − α)-confidence interval of the non-respondent mean µ0 , where the SIPW
estimate is
Pn h(Xi ,Yi )−ĝ(Xi )
e
Ai Yi
(h)
,
µ̂0 = Pi=1
n
h(X
,Y
)−ĝ(X
)A
i
i
i
i
i=1 e
ˆ(h) , µ̂
ˆ(h) , . . . , µ̂
ˆ(h) are the B bootstrap resamples of µ̂(h) . As
where ĝ is as in (12) and µ̂
01
02
0
0B
before, the interval (29) can be computed efficiently using linear fractional programming. In
particular, it is easy to verify that Proposition 5 still holds and thus we only need to consider
the O(n) candidate solutions in Equation (21).
For the ATT, notice that ∆1 = E[Y (1)|A = 1]−E[Y (0)|A = 1]. The first term is identifiable
using the data and the second term is the non-respondent mean by treating Y (0) as the
response. We can use the same procedure in Section 5 to obtain confidence interval of ∆1
ˆ (h0 ) .
under Eβ0 (Λ), using the analogously defined SIPW estimates ∆
1
6.2. Augmented Inverse Probability Weighting. Besides the basic IPW and SIPW
estimators, another commonly used estimator in missing data and observational studies is
the augmented inverse probability weighting (AIPW) estimator which has a double robustness
property explained below (Robins et al., 1994). As usual, we consider the missing data
problem first for simplicity. Apart from the missingness probability, the AIPW estimator
also utilizes another nuisance parameter, f0 (x) = E0 [Y |A = 1, X = x]. Suppose this is
estimated by fˆ(x) from the sample, for example, by linear regression, and e0 (x) is estimated
by ê(x). Then the AIPW estimator (with weight stabilization) is given by
n

µ̂AIPW

Ai − ê(Xi ) ˆ
1 X Ai Yi
−
f (Xi ).
=
n
ê(Xi )
ê(Xi )
i=1

Let ē(x) and f¯(x) be the large-sample limits of ê(x) and fˆ(x). When e0 (x) and f0 (x)
are estimated non-parametrically, then ē(x) = e0 (x) and f¯(x) = f0 (x). When e0 (x) and
f0 (x) are estimated parametrically, ē(x) = eβ0 (x) and f¯(x) = fθ0 (X), the best parametric
approximations (see (6)). Then it is easy to show that µ̂AIPW always estimates (regardless
of the correctness of MAR)


A − ē(X)
P
µ̂AIPW → E0 [Y ] + E0
(Y − f¯(X)) .
(30)
ē(X)
Under MAR (Assumption 1), by first taking expectation conditioning on X it is straightforward to show that the second term is 0 if ē(x) = e0 (X) or f¯(x) = f0 (x). In other words,

SENSITIVITY ANALYSIS FOR IPW

15

µ̂AIPW is consistent for µ if MAR holds and at least one of ê(x) and fˆ(x) is consistent, a
property called double robustness.
When MAR does not hold, by taking expectation conditioning on X and Y , (30) implies
that


e0 (X, Y ) − ē(X)
P
¯
µ̂AIPW − µ → E0
(Y − f (X)) .
ē(X)
Now, as in Section 4, we consider the collection of parametric sensitivity models Eβ0 (Λ) so
ē(x) = eβ0 (x). For h ∈ H(λ), define
"
#
e0 (X, Y ) − e(h) (X, Y )
(h)
(Y − f¯(X)) ,
µ̃f¯ = µ + E0
e(h) (X, Y )

−1
(hβ )
where e(h) (x, y) = 1 + eh(x,y)−logit(P0 (A=1|X=x)) . As before, it is obvious that µ̃f¯ 0 = µ,
the true mean response, where hβ0 (x, y) = logit(eβ0 (x)) − logit(e0 (x, y)).
(h)
The AIPW estimator of µ̃f¯ is
n

(h)

µ̂AIPW =

Ai − ê(h) (Xi ) ˆ
1 X Ai Yi
−
f (Xi )
(h) (X )
(h) (X )
n
ê
ê
i
i
i=1

n
n
1Xˆ
1 X Ai (Yi − fˆ(Xi ))
=
,
f (Xi ) +
n
n
ê(h) (Xi )
i=1

(31)

i=1


−1
where ê(h) (x, y) = 1 + eh(x,y)−ĝ(x,y)
and where ĝ is as in (12). The second term on the
right hand side of (31) is not sample bounded, so it is often preferable to use the stabilized
weights. This results in the following stabilized AIPW (SAIPW) estimator:
" n
#
n
1 X Ai (Yi − fˆ(Xi ))
1
1Xˆ
(h)
f (Xi ) + 1 Pn
µ̂SAIPW =
(h) (X ) n
n
ê(h) (Xi )
i
i=1 Ai ê
n
i=1
i=1
(32)
Pn
n
ˆ(Xi ))(1 + eh(Xi ,Yi )−ĝ(Xi ) )
A
(Y
−
f
1Xˆ
i
i
.
f (Xi ) + i=1 Pn
=
h(Xi ,Yi )−ĝ(Xi ) )
n
i=1 Ai (1 + e
i=1
As before, this estimates µ when h = hβ0 .
Compared to the SIPW estimator (13), in (32) we replace the response Yi by Yi − fˆ(Xi )
P
and add an offset term n1 ni=1 fˆ(Xi ). Therefore, computing the extrema of (32) can still be
formulated as linear fractional programming and the numerical computation is efficient. To
construct asymptotically valid confidence intervals, the outcome regression model fˆ(Xi ) must
be parametric (for example, linear regression). The Z-estimation framework in Appendix B
can then be extended to show that Theorem 3 (validity of percentile bootstrap) still holds
for SAIPW.
Similar to Sections 5 and 6, the above procedures can be extended to observational studies
for estimating the ATE ∆, the non-respondent mean µ0 , or the ATT ∆1 , using analogously
ˆ SAIPW , µ̂0,SAIPW , and ∆
ˆ 1,SAIPW , respectively.
defined SAIPW estimates ∆
6.3. Lipschitz Constraints in the Sensitivity Model. So far we have focused on the
marginal sensitivity models,
E(Λ) or Eβ0 (Λ) = {e(h) (x, y) : h ∈ H(λ)}.

16

SENSITIVITY ANALYSIS FOR IPW

Although this model is very easy to interpret, some deviations h ∈ H(λ) 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
E(Λ) or Eβ0 (Λ) := {e(h) (x, y) : h ∈ Hλ,L },
where
Hλ,L


\
|h(x1 , y1 ) − h(x2 , y2 )|
h:
= H(λ)
≤ L, for all x1 , x2 ∈ X , and y1 , y2 ∈ R ,
d((x1 , y1 ), (x2 , y2 ))

and d is a distance metric defined on the space X × R and L > 0 is the Lipschitz constant.
As Hλ,L ⊆ H(λ), the validity of the percentile bootstrap obviously holds for functions in
Hλ,L . To obtain a range of point estimates and confidence intervals under Eβ0 (Λ, L), we just
need to add the following n(n − 1) constraints in the optimization problem (18):
zi ≤ eL·d((Xi ,Yi ),(Xj ,Yj )) zj , for all 1 ≤ i 6= j ≤ n,
These constraints are linear in {zi }m
i=1 and the resulting optimization problem is still a linear
fractional program, which can be efficiently computed. However, Proposition 5 no longer
holds, so we cannot use the algorithm described after Proposition 5 to solve the optimization
problem corresponding to Hλ,L .
7. Discussion
7.1. Related Frameworks of Sensitivity Analysis. If there is no unmeasured confounder,
the potential outcome Y (a) is independent of the treatment A given X, Y (a) ⊥
⊥ A|X, for
a ∈ {0, 1}. Existing sensitivity analysis methods have considered at least three types of
relaxations of this assumption:
(1) Pattern-mixture models consider a specific difference between the conditional distribution Y (a)|X, A and Y (a)|X (e.g. Robins, 1999, 2002, Birmingham et al., 2003,
Vansteelandt et al., 2006, Daniels and Hogan, 2008).
(2) Selection models consider a specific difference between the conditional distribution
A|X, Y (a) and A|X (e.g. Scharfstein et al., 1999, Gilbert et al., 2003, 2013).
(3) Rosenbaum’s sensitivity models consider a range of possible selection models so that
within a matched set the probabilies of getting treated are no more different than a
number Λ 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 (see e.g.
Brumback et al., 2004).
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|X, Y (a)
and A|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.

SENSITIVITY ANALYSIS FOR IPW

17

7.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 e0 (x, y1 ) and e0 (x, y2 ):
Definition 4 (Rosenbaum’s Sensitivity Models). Fix a parameter a Γ ≥ 1 which will quantify
the degree of violation from the MAR assumption.
(1) For the missing data problem, assume e(·, ·) ∈ R(Γ), where
o
n
1
≤ OR(e(x, y1 ), e(x, y2 )) ≤ Γ, for all x ∈ X , y1 , y2 ∈ R .
(33)
R(Γ) = e(·, ·) :
Γ
(2) For the observational studies problem, assume ea (·, ·) ∈ R(Γ), for a = 0, 1.
The proof of Proposition 1 indicates that not all e0 (X, Y ) are compatible with the observed
data. To compare the marginal sensitivity model E with Rosenbaum’s model R, we introduce
the following concept:
Definition 5 (Compatibility of Sensitivity Model). A sensitivity model e0 (x, y) = P0 (A =
1|X = x, Y = y) is called compatible if the RHS of (34) integrates to 1, for all x ∈ X , that
is, e0 (·, ·) ∈ C, where


ˆ
C = e(·, ·) : OR(e0 (x), e(x, y)) dP0 (y|A = 1, X = x) = 1, for all x ∈ X .
Then Rosenbaum’s sensitivity model and the marginal sensitivity model are related in the
following way:
√
Proposition 7. For any Λ ≥ 1, E( Λ) ⊆ R(Λ) and R(Λ) ∩ C ⊆ E(Λ) ∩ 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 R(Γ).
To this end, it is often assumed that the causal effect is additive and constant across the
individuals, that is, Yi (1) − Yi (0) ≡ ∆, for all i ∈ [n]. Then to determine if an effect ∆
should be included in the (1 − α)-confidence interval, one just needs to test the Fisher null
H0 : Yi (0) = Yi (1) for all i ∈ [n], using Y − ∆A as the outcome (Hodges and Lehmann,
1963) and under Rosenbaum’s sensitivity model. We refer the reader to Rosenbaum (2002c,
Chapter 4) 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

18

SENSITIVITY ANALYSIS FOR IPW

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.
7.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 µ or the ATE ∆ can be
regarded as partially identified under the marginal sensitivity model. In fact, we have adopted
the terminology “partially identified region” for the set {µ(h) : e(h) ∈ Eβ0 (Λ)}.
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.
8. Numerical Example
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 log2 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 used the percentile bootstrap to conduct sensitivity analyses for four estimators: ∆,
ˆ SAIPW , ∆
ˆ 1 , and ∆
ˆ 1,SAIPW . The propensity score P(A = 1|X) is estimated by a logistic
∆
regression and the outcome means E[Y (1)|X] and E[Y (0)|X] are estimated by linear regressions using all 8 covariates. We used B = 1000 bootstrap samples to obtain 90% confidence
intervals of ∆ and ∆1 under Eβ0 (Λ) for 5 values of λ = log Λ = 0, 0.5, 1, 2, 3.
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.

SENSITIVITY ANALYSIS FOR IPW

19

Table 1. Results of different sensitivity analyses under the collection of parametric
sensitivity models Eβ0 (Λ). Five methods were considered: the SIPW and SAIPW
estimators of ATE and ATT as described in Sections 4 and 6, and Rosenbaum’s
sensitivity analysis based on matched pairs. Eight sensitivity parameters were used:
λ = log Λ or γ = log Γ = 0, 0.5, 1, 2, 3. The running time of these five methods were
also reported. For the first four methods, we used B = 1000 bootstrap samples and
used 3 cores in parallel to obtain the confidence intervals. For matching, we do not
include the computational cost of generating the matches.
Estimand
ATE
ATT

Constant
Estimand
ATE
ATT

Constant

Method
SIPW
SAIPW
SIPW
SAIPW
Matching
Method
SIPW
SAIPW
SIPW
SAIPW
Matching

Point
90% CI
Λ = e0 = 1
(1.86,1.86)
(1.63,2.06)
(1.80,1.80)
(1.55,2.04)
(2.09,2.09)
(1.91,2.29)
(2.12,2.12)
(1.93,2.31)
Γ = e0 = 1
(2.08,2.08)
(1.90,2.25)

Point
Λ = e0.5
(1.33,2.37)
(1.34,2.28)
(1.59,2.55)
(1.64,2.56)
Γ = e0.5
(1.75,2.41)

90% CI
= 1.65
(1.11,2.55)
(1.16,2.53)
(1.38,2.72)
(1.45,2.74)
= 1.65
(1.57,2.59)

Point
90% CI
Λ = e2 = 7.39
(-0.10,3.78) (-0.30,4.01)
(0.12,3.61) (-0.02,3.83)
(-0.05,3.43) (-0.43,3.58)
(0.10,3.67) (-0.20,3.92)
Γ = e0 = 1
(0.87,3.36)
(0.58,3.65)

Point
90% CI
Λ = e3 = 20.09
(-0.91,4.78) (-1.15,4.99)
(-0.55,4.31) (-0.76,4.55)
(-1.07,3.53) (-1.36,3.68)
(-0.91,4.15) (-1.36,4.47)
Γ = e0.5 = 1.65
(0.28,3.97) (-0.23,4.48)

Point
Λ = e1
(0.83,2.84)
(0.89,2.76)
(1.04,2.95)
(1.15,2.95)
Γ = e1
(1.45,2.74)

90% CI
= 2.72
(0.61,2.99)
(0.73,2.99)
(0.80,3.12)
(0.95,3.16)
= 2.72
(1.25,2.94)

Running time (seconds)
32.8
47.7
18.0
29.9
1.8

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 1. Alternatively, one can
report the results by plotting the confidence intervals against Λ (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 Λ, while the percentile boot√
strap intervals under Λ were shorter than Rosenbaum’s under Λ. This observation is not
surprising given Proposition 7. 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 Λ = 2.72, while the matching analysis found the effect is
significantly positive for at least Γ = 7.39.
Lastly we want to comment on the computational costs reported in the last column of Table 1. 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.

20

SENSITIVITY ANALYSIS FOR IPW

Causal effect

4

2

●

●

●

●

●

●

Matching
SIPW (ATT)
SAIPW (ATT)

0

1.0

1.6

2.7

7.4

20.1

Λ or Γ

Figure 1. Graphical illustration of the sensitivity analysis results for three methods
(matching, SIPW for ATT, and SAIPW for ATT). The solid error bars are the range
of point estimates and the dashed error bars (together with the solid bars) are the
confidence intervals. The circles/triangles/squares are the mid-points of the solid
bars.

References
Alberto Abadie and Guido W Imbens. On the failure of the bootstrap for matching estimators.
Econometrica, 76(6):1537–1557, 2008.
Joseph G Altonji, Todd E Elder, and Christopher R Taber. Selection on observed and unobserved
variables: Assessing the effectiveness of catholic schools. Journal of Political Economy, 113(1):
151–184, 2005.
Peter M Aronow and Donald KK Lee. Interval estimation of population means under unknown but
bounded probabilities of sample selection. Biometrika, 100(1):235–240, 2012.
Susan Athey, Guido W Imbens, and Stefan Wager. Approximate residual balancing: De-biased inference of average treatment effects in high dimensions. arXiv preprint arXiv:1604.07125, 2016.
Peter C Austin. Variance estimation when using inverse probability of treatment weighting (IPTW)
with survival analysis. Statistics in Medicine, 35(30):5642–5655, 2016.
Alexandre Belloni, Victor Chernozhukov, and Christian Hansen. Inference on treatment effects after
selection among high-dimensional controls. The Review of Economic Studies, 81(2):608–650, 2014.
Jolene Birmingham, Andrea Rotnitzky, and Garrett M Fitzmaurice. Pattern-mixture and selection
models for analysing longitudinal data with monotone missing patterns. Journal of the Royal
Statistical Society: Series B (Statistical Methodology), 65(1):275–297, 2003.
Babette A Brumback, Miguel A Hernán, Sebastien JPA Haneuse, and James M Robins. Sensitivity
analyses for unmeasured confounding assuming a marginal structural model for repeated measures.
Statistics in medicine, 23(5):749–767, 2004.
Abraham Charnes and William W Cooper. Programming with linear fractional functionals. Naval
Research Logistics, 9(3-4):181–186, 1962.
Victor Chernozhukov, Han Hong, and Elie Tamer. Estimation and confidence regions for parameter
sets in econometric models. Econometrica, 75(5):1243–1284, 2007.
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney
Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, to appear, 2017.

SENSITIVITY ANALYSIS FOR IPW

21

Joel E Cohen. Generalized minimax and maximin inequalities for order statistics and quantile functions. Proceedings of the American Mathematical Society, 141(7), 2013.
Jerome Cornfield, William Haenszel, E Cuyler Hammond, Abraham M Lilienfeld, Michael B Shimkin,
and Ernst L Wynder. Smoking and lung cancer: recent evidence and a discussion of some questions.
Journal of the National Cancer Institute, 22(1):173–203, 1959.
Michael J Daniels and Joseph W Hogan. Missing data in longitudinal studies: Strategies for Bayesian
modeling and sensitivity analysis. CRC Press, 2008.
George B Dantzig. Maximization of a linear function of variables subject to linear inequalities. In
T Koopmans, editor, Activity Analysis of Production and Allocation, pages 339–347. Wiley, New
York, 1951.
Peng Ding. A paradox from randomization-based causal inference. Statistical Science, 32(3):331–345,
2017.
Peng Ding and Tyler J VanderWeele. Sensitivity analysis without assumptions. Epidemiology, 27(3):
368, 2016.
Bradley Efron and Robert J Tibshirani. An introduction to the bootstrap. CRC press, 1994.
Ronald A Fisher. Cigarettes, cancer and statistics. Centennial Review of Arts & Science, 2:151–166,
1958.
Colin B Fogarty. Studentized sensitivity analysis for the sample average treatment effect in paired
observational studies. arXiv preprint arXiv:1609.02112, 2017.
Colin B Fogarty, Pixu Shi, Mark E Mikkelsen, and Dylan S Small. Randomization inference and
sensitivity analysis for composite null hypotheses with binary outcomes in matched observational
studies. Journal of the American Statistical Association, 112(517):321–331, 2017.
Joseph L Gastwirth, Abba M Krieger, and Paul R Rosenbaum. Dual and simultaneous sensitivity
analysis for matched pairs. Biometrika, 85(4):907–920, 1998.
Peter B Gilbert, Ronald J Bosch, and Michael G Hudgens. Sensitivity analysis for the assessment of
causal vaccine effects on viral load in hiv vaccine trials. Biometrics, 59(3):531–541, 2003.
Peter B Gilbert, Bryan E Shepherd, and Michael G Hudgens. Sensitivity analysis of per-protocol
time-to-event treatment efficacy in randomized clinical trials. Journal of the American Statistical
Association, 108(503):789–800, 2013.
Keisuke Hirano, Guido W Imbens, and Geert Ridder. Efficient estimation of average treatment effects
using the estimated propensity score. Econometrica, 71(4):1161–1189, 2003.
Joseph L Jr Hodges and Erich L Lehmann. Estimates of location based on rank tests. The Annals of
Mathematical Statistics, pages 598–611, 1963.
Paul W Holland. Statistics and causal inference. Journal of American Statistical Association, 81:
945–960, 1986.
Daniel G Horvitz and Donovan J Thompson. A generalization of sampling without replacement from
a finite universe. Journal of the American statistical Association, 47(260):663–685, 1952.
Michael G Hudgens and M Elizabeth Halloran. Causal vaccine effects on binary postinfection outcomes. Journal of the American Statistical Association, 101(473):51–64, 2006.
Guido W Imbens. Sensitivity to exogeneity assumptions in program evaluation. American Economic
Review, 93(2):126–132, 2003.
Guido W Imbens and Charles F Manski. Confidence intervals for partially identified parameters.
Econometrica, 72(6):1845–1857, 2004.
Joseph DY Kang and Joseph L Schafer. Demystifying double robustness: A comparison of alternative
strategies for estimating a population mean from incomplete data. Statistical science, pages 523–
539, 2007.
Michael R. Kosorok. Introduction to Empirical Processes and Semiparametric Inference. Springer,
2006.
Roderick J Little, Ralph D’agostino, Michael L Cohen, Kay Dickersin, Scott S Emerson, John T
Farrar, Constantine Frangakis, Joseph W Hogan, Geert Molenberghs, Susan A Murphy, et al. The
prevention and treatment of missing data in clinical trials. New England Journal of Medicine, 367

22

SENSITIVITY ANALYSIS FOR IPW

(14):1355–1360, 2012.
Roderick JA Little and Donald B Rubin. Statistical analysis with missing data. John Wiley & Sons,
2014.
Jared K Lunceford and Marie Davidian. Stratification and weighting via the propensity score in
estimation of causal treatment effects: a comparative study. Statistics in medicine, 23(19):2937–
2960, 2004.
Lawrence C McCandless, Paul Gustafson, and Adrian Levy. Bayesian sensitivity analysis for unmeasured confounding in observational studies. Statistics in Medicine, 26(11):2331–2347, 2007.
Luke W Miratrix, Stefan Wager, and Jose R Zubizarreta. Shape-constrained partial identification of
a population mean under unknown probabilities of sample selection. Biometrika, to appear, 2017.
Amy Richardson, Michael G Hudgens, Peter B Gilbert, and Jason P Fine. Nonparametric bounds
and sensitivity analysis of treatment effects. Statistical Science, 29(4):596, 2014.
James Robins, Mariela Sued, Quanhong Lei-Gomez, and Andrea Rotnitzky. Comment: Performance of
double-robust estimators when” inverse probability” weights are highly variable. Statistical Science,
22(4):544–559, 2007.
James M Robins. Association, causation, and marginal structural models. Synthese, 121(1):151–179,
1999.
James M Robins. Comment on “covariance adjustment in randomized experiments and observational
studies”. Statistical Science, 17(3):309–321, 2002.
James M Robins, Andrea Rotnitzky, and Lue Ping Zhao. Estimation of regression coefficients when
some regressors are not always observed. Journal of the American Statistical Association, 89(427):
846–866, 1994.
Paul R Rosenbaum. Sensitivity analysis for certain permutation inferences in matched observational
studies. Biometrika, 74(1):13–26, 1987.
Paul R Rosenbaum. Attributing effects to treatment in matched observational studies. Journal of the
American statistical Association, 97(457):183–192, 2002a.
Paul R Rosenbaum. Covariance adjustment in randomized experiments and observational studies.
Statistical Science, 17(3):286–327, 2002b.
Paul R Rosenbaum. Observational Studies. Springer New York, 2002c.
Paul R Rosenbaum. Two R packages for sensitivity analysis in observational studies. Observational
Studies, 1(1):1–17, 2015.
Paul R Rosenbaum and Donald B Rubin. The central role of the propensity score in observational
studies for causal effects. Biometrika, 70(1):41–55, 1983.
Donald B Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies.
Journal of Educational Psychology, 66(5):688, 1974.
Donald B Rubin. Inference and missing data. Biometrika, 63(3):581–592, 1976.
Donald B Rubin. Comment on “Randomization analysis of experimental data: The fisher randomization test”. Journal of the American Statistical Association, 75(371):591–593, 1980.
Daniel O Scharfstein, Andrea Rotnitzky, and James M Robins. Adjusting for nonignorable dropout using semiparametric nonresponse models. Journal of the American Statistical Association, 94
(448):1096–1120, 1999.
Zhiqiang Tan. A distributional approach for causal inference using propensity scores. Journal of the
American Statistical Association, 101(476):1619–1637, 2006.
The Alpha-Tocopherol Beta Carotene Cancer Prevention Study Group. The effect of vitamin E and
beta carotene on the incidence of lung cancer and other cancers in male smokers. New England
Journal of Medicine, 330:1029–1035, 1994.
Mark J Van der Laan and Sherri Rose. Targeted learning: causal inference for observational and
experimental data. Springer, 2011.
Mark J Van Der Laan and Daniel Rubin. Targeted maximum likelihood learning. The International
Journal of Biostatistics, 2(1), 2006.

SENSITIVITY ANALYSIS FOR IPW

23

Tyler J VanderWeele and Onyebuchi A Arah. Bias formulas for sensitivity analysis of unmeasured
confounding for general outcomes, treatments, and confounders. Epidemiology, 22(1):42–52, 2011.
Stijn Vansteelandt, Els Goetghebeur, Michael G Kenward, and Geert Molenberghs. Ignorance and
uncertainty regions as inferential tools in a sensitivity analysis. Statistica Sinica, 16(3):953–979,
2006.
J. A. Wellner and Yihui Zhan. Bootstrapping z–estimators. manuscript, pages 1–37, 1996.
Jon A. Wellner. Empirical Processes: Theory and Applications. Delft Technical University Lecture
Notes, 2005.
Walter C Willett. Vitamin A and lung cancer. Nutrition Reviews, 48(5):201–211, 1990.
Qingyuan Zhao, Dylan S Small, and Paul R Rosenbaum. Cross-screening in observational studies that
test many hypotheses. Journal of American Statistical Association, 2017.

Appendix A. Proofs
A.1. Proof of Proposition 1. The complete data density can be factorized as
P0 (A, X, Y ) = P0 (A, X) · P0 (Y |A, X).
The first term in the RHS above is identifiable from the data, because A and X are always
observed. For the second term, note that P0 (Y |A = 1, X) is identifiable because Y is observed
if A = 1. For P0 (Y |A = 0, X), using Bayes rule, we get
P0 (Y |A = 0, X)
P0 (A = 0|Y, X) · P0 (Y |X)
P0 (A = 0|X)
P0 (A = 0|Y, X) · {P0 (Y |A = 0, X)P0 (A = 0|X) + P0 (Y |A = 1, X)P0 (A = 1|X)}
.
=
P0 (A = 0|X)
=

Notice that if the denominator P0 (A = 0|X = x) = 0, then we don’t need to consider the
conditional distribution of Y . By simple algebra, we have
P0 (Y |A = 0, X) = OR(e0 (X), e0 (X, Y )) · P0 (Y |A = 1, X),

(34)

where OR(p1 , p2 ) := [p1 /(1 − p1 )]/[p2 /(1 − p2 )] is the odds ratio of p1 , p2 ∈ (0, 1). The only
terms in the above equation that are not directly identifiable from the data are P0 (Y |A =
0, X) and e0 (X, Y ). We can arbitrarily specify e0 (x, y) as long as the RHS of (34) integrates
to 1. This gives us P0 (Y |A = 0, X) and hence the complete data distribution. It is easy
to see that different choices of e0 (x, y) can generate the same density for the observed data
(A, X, AY ) but different densities for the complete data (A, X, Y ).
A.2. Proof of Proposition 2. By definition, under the sensitivity model H(γ), the true
data generating distribution F0 satisfies h0 (x, y) ∈ Hγ . This implies that


P0 (µ ∈ [L, U ]) = P0 µ(h0 ) ∈ [L, U ] ≥ P0 µ(h0 ) ∈ [Lh0 , Uh0 ] .
The last inequality is true because [Lh0 , Uh0 ] is a (1 − α)-confidence interval for µ(h0 ) . Now,
taking limit on both sides gives

lim inf P0 (µ ∈ [L, U ]) ≥ lim inf P0 µ(h0 ) ∈ [Lh0 , Uh0 ] ≥ 1 − α.
n→∞

n→∞

For the second part of the proposition, let hmin = arg minh∈H(γ) µ(h) . By applying assump
tion (14) to the intervals [L(h) , ∞), for h ∈ H(γ), gives lim supn→∞ P0 µ(hmin ) < L ≤ α0 .

24

SENSITIVITY ANALYSIS FOR IPW


Similarly, defining hmax = arg maxh∈H(γ) µ(h) , we have lim supn→∞ P0 µ(hmax ) > U ≤ α − α0 .
Then using the union bound,


lim sup P0 {µ(h) : e(h) ∈ Eβ0 (Γ)} 6⊆ [L, U ]
n→∞
h
i

≤ lim sup P0 µ(hmin ) , µ(hmax ) 6⊆ [L, U ]
n→∞




≤ lim sup P0 µ(hmin ) < L + lim sup P0 µ(hmax ) > U
n→∞

n→∞

≤ α,
completing the proof.
A.3. Proof of Lemma 1. For 1 ≤ b ≤ N , where N = nn is the total number of possible
ˆ(h) the estimate (13) in the b-th resample. Therefore, for
bootstrap resamples, denote by µ̂
b
every h ∈ H(γ),
ˆ(h) ≥
µ̂
b

ˆ(h) , for all 1 ≤ b ≤ N.
inf µ̂
b

(35)

h∈H(γ)

Note that both sides of (35) above are sequences indexed by b ∈ [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 < α < 1,






(h)
(h)
ˆ
ˆ
α
α
Q2
inf µ̂b
= L,
(36)
µ̂b
≥ Q2
b∈[N ]

h∈H(γ)

b∈[N ]

by definition (17). Since (36) is true for any h ∈ H(γ), taking infimum on the LHS above
gives



ˆ(h)
µ̂
≥ L.
inf L(h) = inf Q α2
b
h∈H(γ)

h∈H(γ)

b∈[N ]

The lower bound on U can be proved similarly.
A.4. Proof of Proposition 5. We prove the first claim by contradiction. Suppose there
exists two indices s1 < s2 such that As1 = As2 = 1, Ys1 > Ys2 but zs1 < zs2 . Then consider
the following perturbation,


s 6= s1 , s2 ,
zi ,
0
ĝ(X
)
s
1 ,
zi = zs1 + εe
i = s1 ,


ĝ(X
)
s
2
zs2 − εe
, s = s2 .
When ε > 0 is sufficiently small, (zi0 )m
i=1 is still feasible but the objective becomes larger,
which contradicts the assumption that (zi )m
i=1 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 (20), 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, Γ1 t ≤ z̄i ≤ Γt, 1 ≤ s ≤ m, there exists a solution of (20) such that m equalities
hold. This implies that zi = 1t z̄i is either Γ or Γ1 , for all 1 ≤ s ≤ m. Then, using the first
part of the proposition, there exists a M such that zi = Γ, if Yi ≥ M , and zi = Γ1 , if Yi < M .

SENSITIVITY ANALYSIS FOR IPW

25

√
A.5. Proof of Proposition 7. To begin with suppose that e ∈ E( Γ). Then, for all x ∈ X ,
and y1 , y2 ∈ R,
| log OR(e(x, y1 ), e(x, y2 ))| ≤ | log OR(e(x, y1 ), e0 (x))| + | log OR(e0 (x), e(x, y2 ))|
≤ log Γ,
which implies e ∈ R(Γ).
Now, suppose the function e ∈ R(Γ) ∩ C. By (33),
!
1
≤ OR inf e(x, y), sup e(x, y) ≤ Γ, for all x ∈ X .
y∈R
Γ
y∈R
Notice that inf y e(x, y) ≤ e0 (x) ≤ supy e(x, y), because e(x, y) marginalizes to e0 (x) and
e(x, y) is compatible. Thus, OR(e(x, y), e0 (x)) must be between 1/Γ and Γ, which implies
e ∈ E(Γ) ∩ C.
Appendix B. Proof of Theorem 3
h 
i
To begin with, define κ(h) := E0 A 1 + eh(X,Y )−logit(Pβ0 (A=1|X)) . Recall from (9), that
i
h

1
µ(h) = (h) E0 AY 1 + eh(X,Y )−logit(Pβ0 (A=1|X)) ,
κ
for h ∈ H(γ). We begin by showing that the estimates of the parameters (µ(h) , κ(h) , β00 )0 can
be derived using the framework of Z-estimation.
B.1. The Z-Estimation Framework. Given a vector υ = (ν, κ, β 0 )0 ∈ Θ ⊂ R × R+ × Rd ,
where Θ is the compact parameter space, define the function Q : {0, 1} × R × Rd → Rd+2 as
follows: For t = (a, x0 , y)0 ∈ {0, 1} × R × Rd ,




0
eβ x


a − 1+e
x
β0 x
Q1 (t|υ)

 

0x


h(x,y)−β
Q(t|υ) =  Q2 (t|υ)  :=  κ − a 1 + e
(37)
.




Q3 (t|υ)
0
κν − ay 1 + eh(x,y)−β x
´
Next, define Φ(υ) = Q(t|υ)dP0 (t), where T = (A, X 0 , AY )0 ∼ P0 , the true distribution
generating the data. Note that Φ(υ0 ) = 0, where υ0 = (µ(h) , κ(h) , β00 ) is the true parameter
value. The Z-estimates υ̂ = (µ̂(h) , κ̂(h) , β̂ 0 )0 are obtained by solving the equations




0
1 Pn
eβ̂ Xi
X
A
−
i
i
0
i=1
n
n


 1+eβ̂ Xi

1X
0


1 Pn
(h)
h(X
,Y
)−
Q(Ti |υ̂) = 
(38)
Φn (υ̂) :=
κ̂ − n i=1 Ai 1 + e i i β̂ Xi
 = 0.
n




P
i=1
0
κ̂(h) µ̂(h) − n1 ni=1 Ai Yi 1 + eh(Xi ,Yi )−β̂ Xi
It is easy to see that the Z estimate µ̂(h) is exactly the SIPW estimate (13) for µ(h) and β̂ is
the MLE of β for the logistic regression model.
To formally define the bootstrap estimates, let Pn be the empirical measure of the sample
T1 , T2 , . . . , Tn , where Ti = (Ai , Xi0 , Ai Yi ), and Tˆ1 , Tˆ2 , . . . , Tˆn be i.i.d. samples from the
empirical measure. The bootstrap empirical distribution and the bootstrap empirical process
are
n
√
1X
P̂n =
δT̂i and Ĝn = n(P̂n − Pn ),
n
i=1

26

SENSITIVITY ANALYSIS FOR IPW

respectively. Noting that Φ̂n (υ) =
from the equations:




1
ˆ
ˆ
Φ̂n (υ̂) :=
Q(T̂i |υ̂) = 

n

i=1

n
X

´

ˆ are obtained
Q(t|υ)dP̂n (t), the bootstrap Z-estimates υ̂
1
n

Pn



ˆ0

eβ̂





X̂i

Âi −
X̂i
ˆ
β̂ 0 X̂

 1+e i
ˆ0
1 Pn
X̂
β̂
h(
X̂
,
Ŷ
)−
ˆ
i
i
i
κ̂ − n i=1 Âi 1 + e


P
ˆ0
n
1
(h)
h(
X̂
,
Ŷ
)−
β̂
X̂
ˆ µ̂
ˆ −
i i
i
κ̂
i=1 Âi Ŷi 1 + e
n
i=1




 = 0.




(39)

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. To this end, we need the following assumption.
Assumption 4. The parameter space Θ is compact and the true parameter υ0 is in the
interior of Θ. Moreover, the joint distribution of (Y, X) satisfies:
(1) E[Y 4 ] 
< ∞.

0

(2) det(E

eβ0 X
XX 0
0
(1+eβ0 X )2

) > 0.

h
i
0
(3) For every compact subset S ⊂ Rd , E supβ∈S eβ X < ∞.
Note that all the assumptions are trivially satisfied if we assume that the supports of
X and Y are bounded and E[XX 0 ] is positive definite. The assumptions allow for more
general distributions, for example, distributions with E[et||X|| ] < ∞, for all t ∈ R. Under
this assumption we can derive the asymptotic limiting distribution of the bootstrap estimates
(recall (38) and (39)).
Theorem 8. Suppose the joint distribution of (Y, X) satisfies Assumption 4. Then, for
h ∈ H(γ) be fixed, under P0 ,

√ ˆ
√
D
D
Σ
Φ̇
),
and
n (υ̂ − υ) → N (0, Φ̇−1
n
υ̂
−
υ
→ N (0, Φ̇−1
0
0
0 ΣΦ̇0 ),
where



0

Φ̇0 = E∇υ=υ0 Q(T |υ) = 


0
κ(h)

0
1
µ(h)

0



eβ0 X
−E
XX 0
0
(1+eβ0 X )2
0
EAX 0 eh(X,Y )−β0 X
0

EAY X 0 eh(X,Y )−β0 X



,


(40)

and Σ := E[Q(T , υ0 )Q(T , υ0 )0 ].
ˆ(h) (defined through (39)) is
The limiting distribution of µ̂(h) (recall (13) and (38)) and µ̂
an immediate consequence of the above theorem:
Corollary 9. Suppose the joint distribution of (Y, X) satisfies Assumption 4. Then, for
h ∈ H(γ) be fixed, under P0 ,


√  (h)
√  (h)
D
D
ˆ − µ̂(h) →
n µ̂ − µ(h) → N (0, (σ (h) )2 ), and
n µ̂
N (0, (σ (h) )2 ),
−1
where (σ (h) )2 = (Φ̇−1
0 ΣΦ̇0 )11 is the first diagonal element of Φ̇0 ΣΦ̇0 .

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

SENSITIVITY ANALYSIS FOR IPW

27

B.2. Proof of Theorem 8. We begin by verifying that the matrices Φ̇−1
0 and Σ in Theorem
8 are well-defined. To begin with, a direct computation gives




0
!
eβ0 X
0
XX
0
0
−E
0
β00 X


e
(1+eβ0 X )2
 | = κ(h) | det E
| det(Φ̇0 )| = | det 
XX 0 | > 0,
 EAX 0 eh(X,Y )−β00 X
β00 X )2
0
1 
(1
+
e
0
EAY X 0 eh(X,Y )−β0 X κ(h) µ(h)
by Assumption 4 (2). Therefore, Φ̇0 is invertible. Also, note that Σ < ∞, which follows by
direct multiplication and Assumption 4 (1).
We can now proceed to prove Theorem 8. 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:
(A) The class of functions {t → Q(t|υ) : υ ∈ Θ} is P0 -Glivenko-Cantelli (proved in
Section B.2.1).
(B) ||Φ(υ)||1 is strictly positive outside every open neighborhood of υ0 (proved in Section
B.2.2).
(C) The class of functions {t → Q(t|υ) : υ ∈ Θ} is P0 -Donsker and E[(Q(T |υn ) − Q(T |υ0 ))2 ] →
0, whenever ||υn − υ0 ||1 → 0 (proved in Section B.2.3).
B.2.1. Proof of (A). Define the envelope function B(t) := supυ∈Θ ||Q(t|υ)||1 . Then using
the compactness of Θ, ||h||∞ ≤ γ, and |a| ≤ 1,
||Q(t|υ)||1 ≤ ||Q1 (t|υ)||1 + |Q2 (t|υ)| + |Q3 (t|υ)|
0

≤ ||x||1 + |y| + eγ e−β x (1 + |y|) + M,

(41)

for some absolute constant M . Therefore, by Assumption 4, EB(T ) < ∞, and by (Wellner,
2005, Lemma 6.1), the class of functions {t → Q(t|υ) : υ ∈ Θ} is P0 -Glivenko-Cantelli.
B.2.2. Proof of (B). To begin with note that by Assumption 4 (2),
!
0
eβ X
X = 0,
E
(1 + eβ0 X )
has a unique root β = β0 , since its gradient is positive definite at β0 , and non-negative
definite everywhere. Then, fixing ε > 0, we have
"
#
0
0
eβ0 X
eβ X
X−
> 0,
(42)
||Φ(υ)||1 ≥ E
X
0
1 + eβ 0 X
1 + eβ0 X
1

ε
whenever ||β − β0 ||1 > M
, where M is a constant to be chosen later.
ε
ε
Next, assume that ||β − β0 ||1 ≤ M
. This implies ||β − β0 ||∞ ≤ M
and
h
i
h
i
0
0
0
0
E Aeh(X,Y )−β0 X − Aeh(X,Y )−β X ≤ E Aeh(X,Y ) · e−β0 X − e−β X
0

0

≤ eγ E e−β0 X − e−β X
≤ eγ E (β − β0 )0 Xe−t∗ (for some t∗ ∈ [β00 X, β 0 X])


≤ eγ ||β − β0 ||∞ E ||X||1 e−t∗
1
≤ eγ ||β − β0 ||∞ E[||X||21 E[e−2t∗ ] 2

28

SENSITIVITY ANALYSIS FOR IPW

#! 1

"
≤ ||β − β0 ||∞
≤ K1 (γ) ·

e2γ E[||X||21 ] · E

2

sup

e−2β∗ X

ε
β∗ :||β∗ −β0 ||1 ≤ M

ε
ε
≤
.
M
64K

(43)

by choosing M ≥ 64K · K1 (γ), where
#

"
2

2γ

K1 (γ) := e

E[||X||21 ]

·E

−2β∗ X

e

sup

< ∞,

ε
β∗ :||β∗ −β0 ||1 ≤ M

by Assumption 4, and K := supυ∈Θ |ν| ∈ (0, ∞) by the compactness of Θ. Therefore,
ε
ε
whenever ||β − β0 ||1 ≤ M
and |κ − κ(h) | > 4K
,
i
h
0
0
(44)
||Φ(υ)||1 ≥ κ − κ(h) + E Aeh(X,Y )−β0 X − Aeh(X,Y )−β X > 0.
ε
ε
Finally, assume that ||β − β0 ||1 ≤ M
and |κ − κ(h) | ≤ 4K
, but |ν − µ(h) | > 2κε(h) . Then, as
in (44),
i
h
ε
ε
0
0
E AY eh(X,Y )−β0 X − AY eh(X,Y )−β X ≤ K2 (γ) ·
≤
,
(45)
M
64K
by choosing M ≥ 64K · K2 (γ), where
#
"

K2 (γ)2 := e2γ E[||Y X||21 ] · E

e−2β∗ X < ∞,

sup
ε
β∗ :||β∗ −β0 ||1 ≤ M

by Assumption 4. Using (45) and |κ(h) ν − κν| ≤ 4ε now gives
i
h
0
0
||Φ(υ)||1 = κν − κ(h) µ(h) + E AY eh(X,Y )−β0 X − AY eh(X,Y )−β X > 0.

(46)

Combining (42), (44), and (46), we get, for all δ > 0, inf{||Φ(υ)||2 : ||υ − υ0 ||1 > δ} > 0,
as required.
B.2.3. Proof of (C). We begin by showing the class of functions {t → Q(t|υ) : υ ∈ Θ} is
P0 -Donsker. To this end, let υ1 = (ν1 , κ1 , β10 )0 and υ2 = (ν2 , κ2 , β20 )0 be two points in the
parameter space Θ, and recall the definition of Q(·|υ) from (37). Then, by the mean-value
theorem, there exists t∗ ∈ R such that t∗ ∈ [β10 x, β20 x] such that2
!


0
0
eβ2 x
eβ1 x
et∗
||Q1 (t|υ2 ) − Q1 (t|υ1 )||1 ≤
−
x
≤
(β2 − β1 )0 xx
0
1 + et∗
(1 + et∗ )2
1 + eβ1 x
1
1

≤ |(β2 − β1 )0 x| ||x||1
≤ ||β2 − β10 ||2 ||x||2 ||x||1
.d M1 (x)||β2 − β10 ||1 ,

(47)

where M1 (x) = ||x||21 . Next, observe that
0

0

|Q2 (t|υ2 ) − Q2 (t|υ1 )| ≤ |κ2 − κ1 | + eγ |e−β2 x − e−β1 x |
.γ |κ2 − κ1 | + (β2 − β1 )0 xe−t∗

(for some t∗ ∈ [β10 x, β20 x])

2For x, y ∈ R, x . y means x ≤ Cy, for some constant C > 0. We will use subscripts . , to denote that


the constant may depend on the subscripted parameters.

SENSITIVITY ANALYSIS FOR IPW

29
0

.γ |κ2 − κ1 | + ||β2 − β1 ||2 · ||x||2 sup e−β x
β∈Θ0

.γ,d M2 (x) (|κ2 − κ1 | + ||β2 − β1 ||1 ) ,

(48)

0

where M2 (x) = 1 + ||x||1 supβ∈Θ0 e−β x and Θ0 is the projection of the parameter space to
the last d coordinates. Finally,
0

0

|Q3 (t|υ2 ) − Q2 (t|υ1 )| ≤ |κ2 ν2 − κ1 ν1 | + eγ |ye−β2 x − ye−β1 x |
0

.γ |κ2 − κ1 | + |ν2 − ν1 | + ||β2 − β10 ||2 · ||yx||2 sup e−β x
β∈Θ0

.γ,d M3 (x, y) (|κ2 − κ1 | + |ν2 − ν1 | + ||β2 − β1 ||1 )
0

where M3 (x, y) = 1 + ||yx||1 supβ∈Θ0 e−β x .
Therefore, defining M (x, y) = M1 (x) + M2 (x) + M3 (x, y) and combining (47), (48), and
(49) gives
||Q(t|υ2 ) − Q(t|υ1 )||1 =

3
X

||Qb (t|υ2 ) − Qb (t|υ1 )||1 .γ,d M (x, y)||υ2 − υ1 ||1

(49)

b=1

Since E[M (X, Y )2 ] < ∞, by Assumption 4, this implies that the class {t → Q(t|υ) : υ ∈ Θ}
is P-Donsker.
The inequalities in (47), (48), and (49), combined with Assumption 4 also implies that
E[(Q(T |υn ) − Q(T |υ0 ))2 ] → 0, coordinate-wise, whenever ||υn − υ0 ||1 → 0.
B.3. Completing the Proof of Theorem 3. To begin with note that the α-th bootstrap
(h)
√ ˆ(h)
(h)
ˆ(h) ) = µ̂(h) + ẑ√α,n , where ẑα,n
quantile Qα (µ̂
− µ̂(h) ) ≤ t) ≥ α}. By
:= inf{t : P̂n ( n(µ̂
n
√
(h)
(h)
(h)
Corollary 9, as n → ∞, ẑα,n = zα + OP (1/ n), where zα is α-th quantile of N (0, (σ (h) )2 ).
Then, recalling (15), as n → ∞,



√ 

α
(h)
(h)
P0 µ(h) < LB = P0
n µ̂(h) − µ(h) ≤ z α ,n → ,
2
2


(h)
by Corollary 9. The limit of P0 µ(h) > UB
follows similarly.

Appendix C. Proof of Corollary 6
The proof of Corollary 6 is similar to the proof of Theorem 4. We begin by defining
"
#
"
#
1−A
A
(h)
(h)
κ0 := E (h)
, κ1 := E (h)
.
e0 (X, Y )
e1 (X, Y )
We now set up the SIPW estimate of the ATE (recall (26)) as a Z-estimation problem. To
(h) (h)
this end, given a vector υ = (ν0 , ν1 , κ0 , κ1 , β 0 )0 ∈ Θ ⊂ R × R × R+ × R+ × Rd , where
Θ is the parameter space, define the function Q : {0, 1} × R × Rd → Rd+4 as follows: For
t = (a, x0 , y)0 ∈ {0, 1} × R × Rd ,

30

SENSITIVITY ANALYSIS FOR IPW





Q(t|υ) = 



Q1 (t|υ)
Q2 (t|υ)
Q3 (t|υ)
Q4 (t|υ)
Q5 (t|υ)









 := 









0
eβ x
a − 1+e
x
β0 x


0
(h)
κ1 − a 1 + eh1 (x,y)−β x


0
(h)
κ1 ν1 − ay 1 + eh1 (x,y)−β x


0
(h)
κ0 − (1 − a) 1 + e−h0 (x,y)+β x


0
(h)
κ0 ν0 − (1 − a)y 1 + e−h0 (x,y)+β x








.





(50)

´
Next, define Φ(υ) = Q(t|υ)dP0 (t), where T = (A, X 0 , Y )0 ∼ P0 , the true distribution
(h) (h)
generating the data. Note that Φ(υ0 ) = 0, where υ0 = (µ(h0 ) (0), µ(h1 ) (1), κ0 , κ1 , β00 ) is the
(h ) (h )
true parameter value. Then, as in (38), Z-estimates υ̂ = (µ̂(h0 ) (0), µ̂(h1 ) (1), κ̂0 0 , κ̂1 1 , β̂ 0 )0
are obtained by solving the equations
n
1X
Q(Ti |υ̂) = 0.
Φn (υ̂) :=
n
i=1

µ̂(h1 ) (1) − µ̂(h0 ) (0)

It is easy to see that the Z estimate
is exactly the SIPW estimate (26) for
(h
,h
)
0
1
∆
and β̂ is the MLE of β for the logistic regression model. Moreover, as in (39), the
ˆ are obtained from the equations:
bootstrap Z-estimates υ̂
n
1X
ˆ = 0.
ˆ
Q(T̂i |υ̂)
Φ̂n (υ̂) :=
n
i=1

Then, as in Theorem 8, the limiting joint normality of (µ̂(h0 ) (0), µ̂(h1 ) (1)), and hence the
ˆ (h0 ,h1 ) = µ̂(h0 ) (0) − µ̂(h1 ) (1), can be derived. This would imply the
asymptotic normality of ∆
asymptotic validity of the confidence interval (27), as in Section B.3, completing the proof of
Corollary 6. Details are omitted.

