# Multistage Estimators for Missing Covariates and Incomplete Outcomes

Daniel Suen
Department of Statistics, University of Washington
and
Yen-Chi Chen
Department of Statistics, University of Washington

###### Abstract

We study problems with multiple missing covariates and partially observed
responses. We develop a new framework to handle complex missing covariate
scenarios via inverse probability weighting, regression adjustment, and a
multiply-robust procedure. We apply our framework to three classical problems:
the Cox model from survival analysis, missing response, and binary treatment
from causal inference. We also discuss how to handle missing covariates in
these scenarios, and develop associated identifying theories and asymptotic
theories. We apply our procedure to simulations and an Alzheimer’s disease
dataset and obtain meaningful results.

_Keywords:_ Missing covariates, missing not at random, multiple imputation,
inverse probability weighting, multiply-robustness

##  1 Introduction

Missing data occurs across many different fields, and as the size of datasets
becomes ever larger, it becomes increasingly unrealistic to expect fully
observed datasets. In many statistical problems, we have outcome variables of
interest and a set of associated covariates. However, these outcome variables
may not be completely observed. For example, in a survival problem, either the
censoring time or the time-to-event are observed. In a missing response
problem, a response may be subject to missingness and only be observed for a
subset of the sample. In the binary treatment problem, only one of the two
potential outcomes will be observed. In a traditional setting, all the
covariates are assumed to be completely measured and there has been a vast
amount of literature on how to estimate a parameter of interest. However,
little is known when the covariates could be missing.

Some early attempts to handle missing covariates rely on ignorability
assumptions, utilize estimation procedures in which the assumptions are not
well-specified, or assume that there is only one missing pattern (one
covariate subject to missing). For example, many standard assumptions to
identify an effect of interest take the form of a conditional independence
statement, where the conditioning is done on the covariates $X$. In the case
of missing covariates, strong ignorability assumptions have been proposed by
conditioning on the observed covariates $X_{R}$ and the missing covariate
pattern $R$ and using generalized propensity scores with propensity score
matching (Rosenbaum and Rubin, 1984, 1985). Other approaches include
augmenting the support of the covariates with a new category for missing data
$\\{\texttt{NA}\\}$ (Mayer et al., 2020, 2021). One downside to both of these
methods is that in some situations, the interpretation of these missing data
assumptions may not be so natural as there are effectively different sets of
covariates being conditioned on for different individuals. Various imputation
methods have also been proposed in the Cox model such as those in Paik (1997);
White and Royston (2009); Beesley et al. (2016), but those proposed methods
rely on specific assumptions that may not extend to more general settings. The
multiple imputation approach discussed in Qu and Lipkovich (2009); Crowe et
al. (2010) rely on missing data assumptions that are not clearly stated, which
can lead to difficulty in interpretation. Berg and Yu (2021) also recently
examined estimation under missing-not-at-random (MNAR) but they only consider
a single covariate that are subject to missingness.

In the case of multiple missing variables, the missing-at-random (MAR)
assumption becomes hard to interpret (Robins and Gill, 1997), so a
nonignorable/MNAR missing data assumption may be preferred. There have been
previous developments using nonignorable missing data assumptions in the
literature. For example, Robins and Gill (1997) proposed the permutation
missingness process as a generalization of MAR processes. Also, Tchetgen
Tchetgen et al. (2018) used discrete-choice models as a way to generate a
class of MNAR assumptions. Shpitser (2016); Sadinle and Reiter (2017);
Malinsky et al. (2020) considered a “no self-censoring” or “itemwise
conditionally independent nonresponse” assumption. Chen (2020) introduced the
idea of a pattern graph to generate further MNAR assumptions. Other approaches
include various graphical ones: Mohan et al. (2013); Mohan and Pearl (2014);
Nabi et al. (2020); Mohan and Pearl (2021).

In this paper, we build on this previous work by introducing a general
framework for constructing a set of compatible assumptions on problems with
missing covariates and incomplete outcomes. This framework has several
advantages. First, one can construct identifying restrictions on the missing
covariates that are compatible with common identifying restrictions for the
outcome variables. This ensures that the final outcome of interest is
nonparametrically identifiable (Robins et al., 2000). Also, because we start
with a clear set of assumptions and then develop an estimation method based on
those assumptions, this framework provides a principled approach for
estimating the final outcome of interest. Previous work on MNAR mechanisms can
be utilized in our estimation procedure. Lastly, our framework can be applied
to various settings including but not limited to: survival analysis with
censoring, binary treatment from causal inference, and missing response.

_Outline._ In Section 2, we introduce our framework and several problems of
interest. In Section 3, we outline how to construct a multistage estimator
with the Cox model as a motivating example. In Sections 4 and 5, we explore
simulation results and include a specific data analysis. As a matter of
theoretical interest, we include details on a multistage estimator for the
missing response problem in Appendix A. We discuss the asymptotic theory of
our estimators and demonstrate the $\sqrt{n}$-consistency of our estimators in
Appendix B. Further simulations, data analysis results, and theoretical
details such as sensitivity analyses are provided in Appendices C, D, E, F,
and G.

##  2 Background and framework

In this section, we discuss the problem setup and introduce notation that will
be helpful in describing our framework. We also discuss the setting of three
classical problems from the statistics literature: the Cox model from survival
analysis, missing response, and binary treatment from causal inference.

###  2.1 Problem Setup

In our problems, there are both covariate and outcome-related variables. Let
$Z=(X,R)$ represent the covariate-related variables, where
$X\in\mathbb{R}^{d}$ is a set of covariate variables. Since $X$ is
multivariate, it can have an arbitrary missing pattern. We use the binary
vector $R\in\\{0,1\\}^{d}$ to represent the missing response for $X$. If
$R_{j}=1$, then the variable $X_{j}$ is observed; otherwise, it is missing.
Let $1_{d}=(1,1,\ldots,1)^{\top}$ denote the case where $X_{j}$ is observed
for all $j$. We denote $\bar{r}=1_{d}-r$ as the flipped version of pattern
$r$. Thus, the notation $X_{r}=(X_{j}:r_{j}=1)$ and
$X_{\bar{r}}=(X_{j}:r_{j}=0)$ refer to observed and missing covariate
variables, respectively.

Next, we define $S=(W,A)$ as the outcome-related variables. Here
$W\in\mathbb{R}^{k}$ is a set of outcome variables, and $A\in\mathcal{A}$ is a
random variable over a discrete support $\mathcal{A}$ that represents the
observed pattern of the outcome variables. Note that $A$ has a more general
interpretation than the missing pattern $R$ for $X$ because $A$ does not
necessarily correspond to a missing pattern for $W$ in the traditional sense.
This idea will be made more precise in the following subsections. Similar to
the previous case with the covariates, we will use the notation $W_{A}$ to
denote the set of outcome variables that are observed under the observed
pattern $A$.

Therefore, with a slight abuse of notation, the observed dataset consists of
$n$ IID observations of the form
$\\{(X_{i,R_{i}},R_{i},W_{i,A_{i}},A_{i})\\}_{i=1}^{n}$. Here $R_{i}$ is the
covariate missingness pattern for the $i$th individual. As $X_{i}$ consists of
the covariates for the $i$th individual, some of its entries may be subject to
missingness; thus, we use $X_{i,R_{i}}$ to denote the observed covariates for
the $i$th individual. Similarly, the observed outcome variables may differ
between individuals, so we use $W_{i,A_{i}}$ to denote the observed outcome
variables for the $i$th individual, where $A_{i}$ is the observed pattern of
the outcome variables for the $i$th individual.

To further describe our problem formulation, we provide a brief description of
each of the three problems here (Cox model, missing response, and binary
treatment).

###  2.2 Cox model

In survival analysis, we typically wish to determine the effect of covariates
on a time-to-event variable. However, the event time is often not known for
some individuals due to censoring. For example, people participating in a
medical trial may dropout prematurely, which may mean that the event of
interest is not observed for those individuals. Thus, the observed data is of
the form

$(X_{1,R_{1}},R_{1},Y_{1},\Delta_{1}),\ldots,(X_{n,R_{n}},R_{n},Y_{n},\Delta_{n}),$

where $Y_{i}=\min\\{T_{i},C_{i}\\}\in\mathbb{R}$ is the observed time, and
$\Delta_{i}=1(Y_{i}=T_{i})=1(C_{i}\geq T_{i})$ is the status indicator for a
time-to-event $T_{i}$ and a censoring time $C_{i}$. Thus, the outcome
variables are $W=(C,T)$ and the observed outcome pattern corresponds to
$A=\Delta$. We see that $W_{0}=C$ and $W_{1}=T$, as the observed outcome
pattern denotes whether or not the time-to-event variable was censored. One
can see that the observed outcome pattern $A$ is not associated with a
missingness pattern. Importantly, the outcome variables $W=(C,T)$ are never
completely observed for any individuals.

The Cox proportional hazards model assumes that the effect of the covariates
is multiplicative on $\lambda(t|x)$ through a baseline hazard function
$\lambda(t|x)=\lambda_{0}(t)\exp(\beta^{\top}x)$. Under the Cox model and an
independent censoring assumption, the parameter of interest $\beta^{*}$
satisfies the population estimating equation:

$\displaystyle
U(\beta)=\mathbb{E}\left[\Delta\left(X-\frac{s^{(1)}(Y;\beta)}{s^{(0)}(Y;\beta)}\right)\right]=0,$
(1) $\displaystyle s^{(1)}(t;\beta)=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X],$ $\displaystyle s^{(0)}(t;\beta)=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)].$

Our goal is to estimate $\beta^{*}$ in this semiparametric setting. When the
covariates $X$ are fully observed, under suitable regularity conditions, it is
known that finding a near-solution $\widehat{\beta}$ to the estimating
equation

$\displaystyle
U_{n}(\beta)=\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\left(X_{i}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right)=o_{p}(1/\sqrt{n}),$
$\displaystyle\widehat{s}^{(1)}(t;\beta)=\frac{1}{n}\sum_{i=1}^{n}1(Y_{i}\geq
t)\exp(\beta^{\top}X_{i})X_{i},$
$\displaystyle\widehat{s}^{(0)}(t;\beta)=\frac{1}{n}\sum_{i=1}^{n}1(Y_{i}\geq
t)\exp(\beta^{\top}X_{i})$

provides a consistent and asymptotically normal estimator of $\beta^{*}$.
Since $X$ can now be subject to missingness, we need a way to estimate these
nuisance functions.

###  2.3 Missing response

In the missing response problem, we assume that we have a single outcome
variable of interest $Y\in\mathbb{R}$ that is possibly subject to missingness.
We use the binary random variable $A\in\\{0,1\\}$, which is $1$ if and only if
$Y$ is observed. The observed data is either in the form of
$(X_{i,R_{i}},R_{i},Y_{i},A_{i}=1)$ or $(X_{i,R_{i}},R_{i},A_{i}=0)$. In the
second case ($Y_{i}$ is missing), we can think of the data as
$(X_{R_{i}},R_{i},\texttt{NA},A_{i}=0)$. Using our aforementioned notations,
$W_{0}=\texttt{NA}$ and $W_{1}=Y$. Since $Y$ can be subject to missingness,
$A$ takes a traditional interpretation as a missing response pattern. Then,
$Y_{A_{i}}$ is the observed outcome variable for the $i$th individual. In a
typical case, the parameter of interest is the marginal mean $\mathbb{E}[Y]$
and we often assume $Y\perp A|X$ to identify the marginal mean. Under this
assumption, we have

$\displaystyle\mathbb{E}[Y]$
$\displaystyle=\mathbb{E}\left[\frac{YA}{\pi_{A}(X)}\right]$ (2)
$\displaystyle=\mathbb{E}[m_{1}(X)],$ (3)

where $\pi_{A}(X)=P(A=1|X)$ is the propensity score and
$m_{1}(X)=\mathbb{E}[Y|A=1,X]$ is the outcome regression. The expressions (2)
and (3) are associated with the inverse probability weighting and regression
adjustment estimators, respectively. When the covariates $X$ are always
observed, these quantities can be easily estimated from the observed data, but
the situation becomes more complicated when the covariates are multivariate
and subject to nonmonotone missingness.

###  2.4 Binary treatment

Under the Neyman-Rubin causal model (Rubin, 1978; Splawa-Neyman et al., 1990),
we have a binary treatment $A\in\\{0,1\\}$ that corresponds to potential
outcomes $Y(0)$ and $Y(1)$, respectively. Here the outcome variables are
$W=(Y(0),Y(1))$. Similar to the survival analysis problem with censoring, only
one of the outcome variables is observed for a given individual. The observed
outcome variables are $W_{0}=Y(0)$ and $W_{1}=Y(1)$.

A common parameter of interest is the average treatment effect (ATE)
$\mathbb{E}[Y(1)]-\mathbb{E}[Y(0)]$ and under the classical assumptions
(consistency, SUTVA, weak ignorability; (Rubin, 1974, 1978)) the ATE admits
the following expression

$\displaystyle\mathbb{E}[Y(1)]-\mathbb{E}[Y(0)]$
$\displaystyle=\mathbb{E}\left[\frac{YA}{\pi_{A}(X)}\right]-\mathbb{E}\left[\frac{Y(1-A)}{1-\pi_{A}(X)}\right]$
$\displaystyle=\mathbb{E}[m_{1}(X)]-\mathbb{E}[m_{0}(X)],$

where $\pi_{A}(X)=P(A=1|X)$ is the propensity score and
$m_{a}(X)=\mathbb{E}[Y|A=a,X]$ is the outcome regression. Thus, this leads to
the inverse probability weighting or regression adjustment estimators for
inferring the ATE.

##  3 A multistage estimator

The intermediate quantities that we described in the previous section such as
$s^{(1)}(t;\beta)$ and $\pi_{A}(X)$ are calculated using the covariates $X$
and the outcome-related variables $(W,A)$. Thus, due to this shared property,
before we construct our multistage estimator, we first consider a general
real-valued function $h$ of the covariate-related variables $X$ and the
outcome-related variables $(W,A)$, as defined below.

######  Definition 1 (Observed-decomposable functions).

Let $h:\mathbb{R}^{d}\times\mathbb{R}^{j}\times\mathcal{A}\mapsto\mathbb{R}$
be a function of interest. Call $h$ an observed-decomposable function (ODF) if

$h(X,W,A)=\sum_{a\in\mathcal{A}}f_{a}(X,W_{a})1(A=a),$

where
$f_{a}:\mathbb{R}^{d}\times\mathbb{R}^{\text{dim}(W_{a})}\mapsto\mathbb{R}$ is
a fixed (non-random) real-valued function for all $a\in\mathcal{A}$.

The interpretation of an ODF is that for each observed outcome pattern $a$, we
have a separate real-valued function $f_{a}$ that is dependent on all of the
covariates $X$ and the observed outcome variables $W_{a}$ under pattern $a$.
Since $h$ has this decomposition, it is computable in the traditional setting
with incomplete outcomes and completely observed covariates. As we will see in
the construction of a multistage estimator, many quantities of interest can be
written in this way, so this is a very natural decomposition. We will either
typically be interested in estimating $\theta=\mathbb{E}[h(X,W,A)]$ or
$\mathbb{E}[h(X,W,A)]=0$ will be an estimating equation. We provide an
illustrative example of $h$ using the missing response problem.

######  Example 1 (Missing response).

In the traditional missing response problem without missing covariates where
we commonly wish to estimate $\mathbb{E}[Y]$, there are 3 traditional
estimators: regression adjustment on $Y$, inverse probability weighting on
$Y$, and doubly robust on $Y$ estimators. Because our goal is to estimate
$\mathbb{E}[Y]$, the corresponding $h$, $f_{0}$, and $f_{1}$ functions take
the following forms in these three situations. Note that in the case $A=0$
($Y$ is missing), we simplify notation such that
$f_{0}(X,W_{0})=f_{0}(X,\texttt{NA})=f_{0}(X)$ since $W_{0}$ can only take one
possible outcome NA.

  1. 1.

Regression adjustment on $Y$. We have $h(X,W,A)=\mathbb{E}[Y|A=1,X],$ which
implies that $f_{0}(X)=f_{1}(X,Y)=\mathbb{E}[Y|A=1,X]$.

  2. 2.

Inverse probability weighting on $Y$. We have
$h(X,W,A)=\frac{1(A=1)Y}{P(A=1|X)},$ which implies that $f_{0}(X)=0$ and
$f_{1}(X,Y)=\dfrac{Y}{P(A=1|X)}$.

  3. 3.

Doubly-robust on $Y$. We have
$h(X,W,A)=\frac{1(A=1)Y}{P(A=1|X)}+\mathbb{E}[Y|A=1,X]\left(1-\frac{1(A=1)}{P(A=1|X)}\right),$
which implies that $f_{0}(X)=\mathbb{E}[Y|A=1,X]$ and
$f_{1}(X,Y)=\dfrac{Y-\mathbb{E}[Y|A=1,X]}{P(A=1|X)}+\mathbb{E}[Y|A=1,X]$.

######  Example 2 (Cox model).

In the Cox model, $\Delta=A$ and $W_{0}=W_{1}=Y$ and there are three key
equations as described in equation (1):
$U(\beta),s^{(1)}(t;\beta),s^{(0)}(t;\beta)$.

  1. 1.

$U(\beta)$. For this quantity, $U(\beta)=\mathbb{E}[h(X,W,A)]$, so
$f_{0}(X,W_{0})=0$ and
$f_{1}(X,W_{1})=\left(X-\frac{s^{(1)}(Y;\beta)}{s^{(0)}(Y;\beta)}\right)$.

  2. 2.

$s^{(1)}(t;\beta)$. In this case, by setting
$s^{(1)}(t;\beta)=\mathbb{E}[h(X,W,A)]$, we obtain $f_{0}(X,W_{0})=1(Y\geq
t)\exp(\beta^{T}X)X$ and $f_{1}(X,W_{1})=1(Y\geq t)\exp(\beta^{T}X)X$.

  3. 3.

$s^{(0)}(t;\beta)$. This case is similar to $s^{(1)}(t;\beta)$, so that
$f_{0}(X,W_{0})=1(Y\geq t)\exp(\beta^{T}X)$ and $f_{1}(X,W_{1})=1(Y\geq
t)\exp(\beta^{T}X)$.

###  3.1 The IPW approach for missing covariates

To handle missing covariates, we start with an IPW procedure based on using
the complete odds as defined below.

######  Definition 2 (Complete odds).

For all missingness patterns $r$ and observed outcome patterns $a$, the
complete odds is

$O_{r,a}(X,W_{a})=\frac{P(R=r|X,W_{a},A=a)}{P(R=1_{d}|X,W_{a},A=a)}.$

Thus, these complete odds $O_{r,a}(X,W_{a})$ are a function of all of the
outcome variables completely observed under outcome pattern $a$ and the
covariates. With the complete odds, $\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$ can be
decomposed via the following lemma.

######  Lemma 1 (Complete odds decomposition).

The complete odds can be used to decompose the following two expectations as
follows

$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$
$\displaystyle=\mathbb{E}\left[f_{a}(X,W_{a})1(A=a)O_{r,a}(X,W_{a})1(R=1_{d})\right],$
$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$
$\displaystyle=\sum_{r}\mathbb{E}\left[f_{a}(X,W_{a})1(A=a)O_{r,a}(X,W_{a})1(R=1_{d})\right].$
(4)

This lemma shows that the expectation can be estimated using the complete
cases since we have the indicator $1(R=1_{d})$ in the expression on the RHS.
Thus, we see that $\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$ is identifiable, provided
that the complete odds $O_{r,a}(X,W_{a})$ are identifiable for all $r$. So, it
is clear what we need to show before estimation. Additionally, note that

$\sum_{r}O_{r,a}(X,W_{a})=\sum_{r}\frac{P(R=r|X,W_{a},A=a)}{P(R=1_{d}|X,W_{a},A=a)}=\frac{1}{P(R=1_{d}|X,W_{a},A=a)}$
(5)

is the inverse propensity score, so we see that (4) recovers the form of the
IPW estimator. In light of this, one can view the complete odds as an
alternate version of a propensity score. In this situation, it is not required
to use the complete odds approach since one may use the propensity score
directly. However, modeling via complete odds has several nice properties; we
discuss this in Remark 1.

In the following proposition, we show that a set of identifiable complete odds
is equivalent to a set of identifiable selection probabilities. This
connection offers one key insight: modeling the covariate missingness
mechanism through the complete odds is equivalent to modeling the covariate
missingness mechanism using the selection model (Little and Rubin, 2002).
Specific missing data assumptions may make the former approach easier than the
latter, so this can be a useful result.

######  Proposition 2 (Equivalence of complete odds and selection
probabilities).

Fix an observed outcome pattern $a$. The complete odds $O_{r,a}(X,W_{a})$ are
identifiable for all $r$ if and only if the selection probabilities
$P(R=r|X,W_{a},a)$ are identifiable for all $r$.

Furthermore, another consequence of this proposition is that a set of
identifiable complete odds implies that the distribution $p(x,r,w_{a},a)$ is
identifiable for all $a$. Marginalizing out $r$ implies that the distribution
$p(x,w_{a},a)$ is nonparametrically identified (Robins et al., 2000), which is
sufficient for showing the effect of interest is identifiable under the
traditional assumptions that assume fully observed covariates. This implies
that with a set of identifiable complete odds, we can use existing problem-
specific assumptions to identify the parameter of interest.

######  Remark 1 (Advantages to modeling the complete odds).

_A modeling approach using the complete odds can be preferable to one using
the selection probabilities directly. Suppose $r$ and $a$ are fixed. Then,
$O_{r,a}(X,W_{a})$ is variation independent of any $O_{r^{\prime},a}(X,W_{a})$
for $r^{\prime}\neq r$. Indeed, $O_{r,a}(X,W_{a})$ can take any value in
$[0,\infty)$ and is not constrained by any other odds function
$O_{r^{\prime},a}(X,W_{a})$. On the other hand, the selection probabilities
$P(R=r|x,w_{a},a)$ have an explicit sum to 1 constraint, so the selection
probabilities for all missing patterns $r$ need to be modeled at the same
time. In principle, one can use different methods to model each odds
$O_{r,a}(X,W_{a})$, provided it is identifiable under a missing data
assumption for $X$. These allows for some flexibility because one can choose
to have different binary classifiers for each pattern $r$. Due to the
variation independence, they can all be compatible with each other._

Because of Proposition 2, we only need to identify the complete odds. Here we
consider a (conditional) complete case missing value (CCMV) assumption
(Little, 1993; Tchetgen Tchetgen et al., 2018).

######  Definition 3 (CCMV missing data assumption).

Using a pattern mixture model factorization, the complete case missing value
(CCMV) assumption equates the following two distributions

$p(x_{\bar{r}}|x_{r},r,w_{a},a)\stackrel{{\scriptstyle\text{CCMV}}}{{=}}p(x_{\bar{r}}|x_{r},R=1_{d},w_{a},a)$
(6)

for all covariate missingness patterns $r$ and observed outcome patterns $a$.

Note that the above condition is a conditional CCMV in the sense that we are
using the CCMV for the missing $X$ conditioned on the observed outcomes. We
call it CCMV for abbreviation. Under Proposition 3, this particular missing
data assumption implies a simple closed form expression $Q_{r,a}(X_{r},W_{a})$
for the odds.

######  Proposition 3 (Odds under CCMV).

Under the CCMV assumption for the covariates $X$, the odds $O_{r,a}(X,W_{a})$
has the following simplified form in terms of observable quantities $X_{r}$
and $W_{a}$, namely

$O_{r,a}(X,W_{a}):=\frac{P(R=r|X,W_{a},a)}{P(R=1_{d}|X,W_{a},a)}\stackrel{{\scriptstyle\text{CCMV}}}{{=}}\frac{P(R=r|X_{r},W_{a},a)}{P(R=1_{d}|X_{r},W_{a},a)}=:Q_{r,a}(X_{r},W_{a}).$
(7)

Because each $Q_{r,a}(X_{r},W_{a})$ is identifiable from the data, Proposition
3 shows that the CCMV assumption allows us to estimate the parameter of
interest under missing covariates. The modeling $Q_{r,a}(X_{r},W_{a})$ is easy
because it is comparing the proportions between two missing patterns.

We provide some comments on generalizing the CCMV assumption to other
assumptions in Appendix C. Other missing data assumptions can yield other
identifiable and estimable expressions for the odds, but one would simply
obtain a different equality than the one in (7). For example, Chen (2020)
exploits a pattern graph identifying restriction to present a recursive form
for computing the odds.

######  Remark 2 (Relation to density ratio estimation).

_Wang and Kim (2021) have examined the problem of estimating an inverse
propensity score using a density ratio representation under a missing at
random assumption. Interestingly, because they assume the variables can be
partitioned into one set of always missing variables and one set of always
observed variables, the missing at random assumption is equivalent to CCMV.
Additionally, our complete odds can be expressed as a density ratio as well,
so these problems are related. We elect to model our complete odds directly
rather than reexpress it as a ratio of densities, but note that one can borrow
methods from the density ratio literature to also tackle this problem. See
Sugiyama et al. (2012) for various approaches._

######  Theorem 4 (Nonparametric identification of CCMV).

The CCMV assumption in Definition 6 nonparametrically identifies the full-
covariate distribution $p(x,r,w_{a},a)$.

With Theorem 4, we are able to identify the full-covariate distribution
$p(x,r,w_{a},a)$. This implies that the problem reduces to scenarios where the
covariates are observed. Thus, we can estimate the parameter of interest as if
all covariates are observed. In what follows, we give an example of Cox model,
and in Appendix A, we provide examples of missing response problems.

####  3.1.1 Example: Cox model and survival problems

Recall that in the survival analysis problem, $\Delta=A$ and $W_{0}=W_{1}=Y$.
In the case of $s^{(1)}(t;\beta)$, we have

$f_{0}(X,W_{0})=f_{1}(X,W_{1})=1(Y\geq t)\exp(\beta^{\top}X)X.$ (8)

For any fixed time point $t$, a direct application of Lemma 1 yields

$\displaystyle s^{(1)}(t;\beta)$ $\displaystyle=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X]$ $\displaystyle=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X1(\Delta=1)]+\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X1(\Delta=0)]$ $\displaystyle=\sum_{r}\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X1(\Delta=1)O_{r,1}(X,Y)1(R=1_{d})]$
$\displaystyle\quad+\sum_{r}\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X1(\Delta=0)O_{r,0}(X,Y)1(R=1_{d})]$
$\displaystyle=\sum_{r,\delta}\mathbb{E}[1(Y\geq
t)Xe^{\beta^{\top}X}O_{r,\delta}(X,Y)1(\Delta=\delta)1(R=1_{d})],$

where $O_{r,\delta}(X,Y)=P(R=r|X,Y,\delta)/P(R=1_{d}|X,Y,\delta)$ are the
complete odds. A similar expression can be obtained for $s^{(0)}(t;\beta)$.
Next, we consider the estimating equation, which rewrites as

$\displaystyle U(\beta)$
$\displaystyle=\mathbb{E}\left[\Delta\left(X-\frac{s^{(1)}(Y;\beta)}{s^{(0)}(Y;\beta)}\right)\right]$
$\displaystyle=\sum_{r}\mathbb{E}\left[\Delta\left(X-\frac{s^{(1)}(Y;\beta)}{s^{(0)}(Y;\beta)}\right)O_{r,1}(X,Y)1(R=1_{d})\right]$
(9)

using Lemma 1.

With Proposition 3, $O_{r,\delta}(X,Y)=Q_{r,\delta}(X_{r},Y)$, which provides
an elegant way of computing the odds because the quantity
$Q_{r,\delta}(X_{r},Y)$ is a ratio of two conditional probabilities that no
longer depend on all of the covariates but rather observed covariates $X_{r}$.
We now discuss the construction of a multistage estimator under CCMV using the
complete odds. Henceforth, we will use $Q_{r,\delta}(X_{r},Y)$ rather than
$O_{r,\delta}(X,Y)$ because they are equal by Proposition 3. This choice of an
identifying restriction allows for us to compute the inverse propensity score
in a fairly straightforward manner.

There are three steps to our multistage estimator, which we outline below.

  1. 1.

Stage 1 (Modeling the missingness of $X$ using the complete odds). We start by
modeling $Q_{r,\delta}(X_{r},Y)$ for every $r$ and $\delta$. For example, we
can impose a parametric model $Q_{r,\delta}(X_{r},Y;\alpha_{r,\delta})$ such
as a logistic regression. Because this essentially becomes a binary
classification problem (comparing the classes $R=r$ and $R=1_{d}$ based on the
observed data $X_{r}$ and $Y$), one may also appeal to various machine
learning methods such as decision trees and support vector machines. Other
parametric binary classification tools are also available.

  2. 2.

Stage 2 (Estimating intermediate quantities). Next, we estimate the
intermediate quantities of interest that are dependent on the potentially
missing covariates $X$. Because of the model described in Stage 1, this
estimation becomes possible. For example, assuming that the model for
$Q_{r,\delta}(X_{r},Y;\alpha_{r,\delta})$ is correctly specified and
$\widehat{\alpha}_{r,\delta}$ is a consistent estimator of
$\alpha_{r,\delta}^{*}$, under the suitable regularity conditions, the
following is a consistent estimator for $s^{(1)}(t;\beta)$

$\displaystyle\widehat{s}^{(1)}(t;\beta)$
$\displaystyle=\sum_{r,\delta}\left(\frac{1}{n}\sum_{i=1}^{n}1(Y_{i}\geq
t)X_{i}e^{\beta^{\top}X_{i}}Q_{r,\delta}(X_{i,r},Y_{i};\widehat{\alpha}_{r,\delta})1(\Delta_{i}=\delta)1(R_{i}=1_{d})\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\biggr{(}1(Y_{i}\geq
t)X_{i}e^{\beta^{\top}X_{i}}\cdot\underbrace{\biggr{[}\sum_{r,\delta}Q_{r,\delta}(X_{i,r},Y_{i};\widehat{\alpha}_{r,\delta})1(\Delta_{i}=\delta)1(R_{i}=1_{d})\biggr{]}}_{\lambda_{i}}\biggr{)}$
(10) $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot 1(Y_{i}\geq
t)X_{i}e^{\beta^{\top}X_{i}},$

where $\lambda_{i}$ is a weight for the $i$th individual, as defined in (10).
Similarly, we can mimic the same construction to estimate $s^{(0)}(t;\beta)$,
yielding

$\widehat{s}^{(0)}(t;\beta)=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot
1(Y_{i}\geq t)e^{\beta^{\top}X_{i}}.$

  3. 3.

Stage 3 (Constructing the final estimator). Lastly, we assemble the
intermediate quantities of interest together to create the final estimator.
Following the decomposition in (9), we have

$\displaystyle U_{n}(\beta)$
$\displaystyle=\sum_{r,\delta}\left(\frac{1}{n}\sum_{i=1}^{n}\Delta_{i}\left(X_{i}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right)Q_{r,\delta}(X_{i,r},Y_{i};\widehat{\alpha}_{r,1})1(R_{i}=1_{d})\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\biggr{(}\Delta_{i}\left(X_{i}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right)\cdot\underbrace{\sum_{r}Q_{r,1}(X_{i,r},Y_{i};\widehat{\alpha}_{r,1})1(\Delta_{i}=1)1(R_{i}=1_{d})}_{\lambda_{i,1}}\biggr{)}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot\Delta_{i}\left(X_{i}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right),$

where the second to last equality follows because $\Delta_{i}\in\\{0,1\\}$.
The last equality follows because
$\lambda_{i}\Delta_{i}=\lambda_{i,1}\Delta_{i}$. Thus, $U_{n}(\beta)=0$ is a
reweighted estimating equation and can be solved to provide a consistent
estimator $\widehat{\beta}$ for $\beta^{*}$.

A nice property is that the weights $\lambda_{i}$ and $\lambda_{i,1}$ need
only to be computed from the data once and can be reused thereafter. We
summarize the implementation of this reweighted Cox model in Algorithm 1.

Algorithm 1 Fitting a weighted Cox model under CCMV for $X$

  1. 1.

Input: $\\{(X_{i,R_{i}},R_{i},Y_{i},\Delta_{i})\\}_{i=1}^{n}$

  2. 2.

For each pair $(r,\delta)$, fit the parametric model
$Q_{r,\delta}(X_{r},Y;\alpha_{r,\delta}).$

  3. 3.

For $i=1,\ldots,n$, do

    1. 3-1.

Form the weights

$\lambda_{i}=\sum_{\delta}\sum_{r}Q_{r,\delta}(X_{i,r},Y_{i};\widehat{\alpha}_{r,\delta})1(\Delta_{i}=\delta)1(R_{i}=1_{d}).$

  4. 4.

Estimate the intermediate quantities

$\widehat{s}^{(0)}(t;\beta)=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot
1(Y_{i}\geq t)e^{\beta^{\top}X_{i}},$
$\widehat{s}^{(1)}(t;\beta)=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot
1(Y_{i}\geq t)X_{i}e^{\beta^{\top}X_{i}}.$

  5. 5.

Find $\widehat{\beta}$ by finding a near-solution to the finite sample
estimating equation

$U_{n}(\beta)=\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\cdot\Delta_{i}\left(X_{i}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right)=0.$

  6. 6.

Return: $\widehat{\beta}$.

###  3.2 Regression adjustment and imputation for missing covariates

In addition to the IPW approach for handling missing covariates, we can use a
regression adjustment approach as well. Moreover, this regression adjustment
naturally leads to an imputation procedure. The following lemma is a key step
to the regression adjustment (RA).

######  Lemma 5 (RA decomposition).

Regression adjustment can be used to decompose the following two expectations

$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$
$\displaystyle=\mathbb{E}[m_{r,a}(X_{r},W_{a})1(R=r)1(A=a)],$
$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$
$\displaystyle=\sum_{r}\mathbb{E}[m_{r,a}(X_{r},W_{a})1(R=r)1(A=a)],$

where $m_{r,a}(X_{r},W_{a})=\mathbb{E}[f_{a}(X,W_{a})|X_{r},R=r,W_{a},A=a]$ is
a regression function.

Thus, in order to estimate the expectation, it suffices to model the
regression function $m_{r,a}(X_{r},W_{a})$, which is a conditional
expectation.

In practice, placing a parametric model for each regression function via
$m_{r,a}(x_{r},w_{a};\beta_{r,a})$ may be problematic because naively fitting
models can lead to incompatible regression functions and uncongenial sources
of input (Meng et al., 1994). Instead, we recommend placing models on the
extrapolation densities $p(x_{\bar{r}}|x_{r},r,a,w_{a})$ for all $r$. This
ensures that the modeled regression functions will be compatible with each
other.

Under the CCMV assumption (equation 6), we have

$p(x_{\bar{r}}|x_{r},r,w_{a},a)=p(x_{\bar{r}}|x_{r},R=1_{d},w_{a},a),$

which further implies

$m_{r,a}(X_{r},W_{a})=\mathbb{E}[f_{a}(X,W_{a})|X_{r},R=r,W_{a},A=a]\stackrel{{\scriptstyle\text{CCMV}}}{{=}}\mathbb{E}[f_{a}(X,W_{a})|X_{r},R=1_{d},W_{a},A=a].$

The RHS is identifiable using observations with complete covariates. In this
case, there are two ways of estimating the regression function:

  1. 1.

Method 1: Closed-form computation. In some good scenarios, we may be able to
derive the the closed-form of the expectation
$m_{r,a}(X_{r},W_{a})=\mathbb{E}[f_{a}(X,W_{a})|X_{r},R=r,W_{a},A=a]$ from a
parametric model of the extrapolation densities
$p(x_{\bar{r}}|x_{r},r,a,w_{a})$. For instance, a multivariate Gaussian model
implies that the conditional distributions will be conditional Gaussians.

  2. 2.

Method 2: Imputation (Monte Carlo approximation). Alternatively, we can
approximate the conditional expectation by a Monte Carlo method. This is
equivalent to the frequentist multiple imputation. We generate from the
extrapolation density $p(x_{\bar{r}}|x_{r},r,a,w_{a})$ and then evaluate the
conditional expectation $\mathbb{E}[f_{a}(X,W_{a})|X_{r},R=r,W_{a},A=a]$.
Under CCMV (6), the extrapolation density $p(x_{\bar{r}}|x_{r},r,y,\delta)$ is
equated to the identifiable distribution
$p(x_{\bar{r}}|x_{r},R=1_{d},y,\delta)$. The latter distribution can be
estimated using the complete cases.

####  3.2.1 Example: Cox model and survival problems

Using Lemma 5 with the functions defined in (8), we have

$\displaystyle s^{(1)}(t;\beta)$ $\displaystyle=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X]$
$\displaystyle=\sum_{r,\delta}\mathbb{E}[m^{(1)}_{r,\delta}(X_{r},Y;t,\beta)1(R=r)1(\Delta=\delta)].$

Thus, to estimate $s^{(1)}(t;\beta)$, it suffices to model
$m^{(1)}_{r,\delta}(X_{r},Y;t,\beta):=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)X|X_{r},R=r,Y,\Delta=\delta]$ for all $(r,\delta)$ pairs.
We obtain a similar expression for $s^{(0)}(t;\beta)$ via

$s^{(0)}(t;\beta)=\sum_{r,\delta}\mathbb{E}[m^{(0)}_{r,\delta}(X_{r},Y;t,\beta)1(R=r)1(\Delta=\delta)]$

for $m^{(0)}_{r,\delta}(X_{r},Y;t,\beta):=\mathbb{E}[1(Y\geq
t)\exp(\beta^{\top}X)|X_{r},R=r,Y,\Delta=\delta]$.

Applying Lemma 5 on the overall estimating equation, we obtain

$\displaystyle
U(\beta)=\sum_{r}\mathbb{E}\left[\mathbb{E}\left[X-\frac{s^{(1)}(Y;\beta)}{s^{(0)}(Y;\beta)}\biggr{|}X_{r},R,Y,\Delta\right]\cdot
1(R=r)\cdot 1(\Delta=1)\right].$

Therefore, the regression adjustment method is equivalent to finding a near-
solution to the following finite sample estimating equation

$\displaystyle
U_{n}(\beta)=\sum_{r}\left[\frac{1}{n}\sum_{i=1}^{n}\left(\widehat{\mathbb{E}}\left[X-\frac{\widehat{s}^{(1)}(Y;\beta)}{\widehat{s}^{(0)}(Y;\beta)}\biggr{|}X_{i,r},R_{i},Y_{i},\Delta_{i}\right]\right)\cdot
1(R_{i}=r)\cdot 1(\Delta_{i}=1)\right]=o_{p}(1/\sqrt{n}),$
$\displaystyle\widehat{s}^{(1)}(t;\beta)=\sum_{r,\delta}\left[\frac{1}{n}\sum_{i=1}^{n}\widehat{m}^{(1)}_{r,\delta}(X_{i,r},Y_{i};t,\beta)1(R_{i}=r)1(\Delta_{i}=\delta)\right],$
$\displaystyle\widehat{s}^{(0)}(t;\beta)=\sum_{r,\delta}\left[\frac{1}{n}\sum_{i=1}^{n}\widehat{m}^{(0)}_{r,\delta}(X_{i,r},Y_{i};t,\beta)1(R_{i}=r)1(\Delta_{i}=\delta)\right].$
(11)

This can be achieved using a Monte Carlo procedure. We describe the multistage
process below, and we summarize our sampling method in Algorithm 2. This
algorithm is the same as performing a frequentist (type B) multiple imputation
approach described in Wang and Robins (1998) and Little and Rubin (2002).

  1. 1.

Stage 1 (Modeling the missingness of $X$ with extrapolation densities). As
previously discussed, instead of modeling the regression functions
$m^{(1)}_{r,\delta}(X_{r},Y;t,\beta)$ and
$m^{(0)}_{r,\delta}(X_{r},Y;t,\beta)$ for all $r$ and $\delta$ directly, we
will instead focus on the extrapolation densities. Under CCMV, we simply need
to model the joint distribution $p(x,R=1_{d},y,\delta)$. This will induce a
model $p(x_{\bar{r}}|x_{r},R=1_{d},y,\delta)$ for any covariate missing
pattern $r$. By CCMV, that conditional distribution equates exactly to the
extrapolation distribution for the missing data.

  2. 2.

Stage 2 (Estimating intermediate quantities). Next, we again outline how to
estimate the nuisance functions $s^{(1)}(t;\beta)$ and $s^{(0)}(t;\beta)$ that
are dependent on the potentially missing covariates $X$. We estimate these
using a Monte Carlo procedure through iterative imputation of the missing
covariates.

By estimating the regression functions using $M$ draws from the extrapolation
distribution, the following is a consistent estimator for $s^{(1)}(t;\beta)$

$\displaystyle\widehat{s}^{(1)}(t;\beta)$
$\displaystyle=\sum_{r,\delta}\left[\frac{1}{n}\sum_{i=1}^{n}\widehat{m}^{(1)}_{r,\delta}(X_{i,r},Y_{i};t,\beta)1(R_{i}=r)1(\Delta_{i}=\delta)\right]$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{r,\delta}\widehat{m}^{(1)}_{r,\delta}(X_{i,r},Y_{i};t,\beta)1(R_{i}=r)1(\Delta_{i}=\delta)\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{1}{M}\sum_{j=1}^{M}1(Y_{i}\geq
t)\widetilde{X}_{i}^{(j)}e^{\beta^{\top}\widetilde{X}_{i}^{(j)}}\right),$
$\displaystyle=\frac{1}{Mn}\sum_{j=1}^{M}\sum_{i=1}^{n}1(Y_{i}\geq
t)\widetilde{X}_{i}^{(j)}e^{\beta^{\top}\widetilde{X}_{i}^{(j)}},$

where
$\widetilde{X}_{i}^{(j)}=\left(X_{i,R_{i}},\widetilde{X}_{i,\bar{R}_{i}}^{(j)}\right)$.
Here $\widetilde{X}_{i,\bar{R}_{i}}^{(j)}$ denotes the $j$th imputed tuple of
missing covariates for the $i$th observation, and $\widetilde{X}_{i}^{(j)}$ is
the result of attaching the imputed variables to the observed covariates
$X_{i,R_{i}}$. A similar construction applied to $s^{(0)}(t;\beta)$ yields

$\widehat{s}^{(0)}(t;\beta)=\frac{1}{Mn}\sum_{j=1}^{M}\sum_{i=1}^{n}1(Y_{i}\geq
t)e^{\beta^{\top}\widetilde{X}_{i}^{(j)}}.$

  3. 3.

Stage 3 (Constructing the final estimator). Lastly, we again merge everything
together in the final stage. Following the decomposition in (11), we have

$\displaystyle U_{n}(\beta)$
$\displaystyle=\frac{1}{Mn}\sum_{j=1}^{M}\sum_{i=1}^{n}\Delta_{i}\left(\widetilde{X}_{i}^{(j)}-\frac{\widehat{s}^{(1)}(Y_{i};\beta)}{\widehat{s}^{(0)}(Y_{i};\beta)}\right).$

We solve the above estimating equation to obtain the final estimate
$\widehat{\beta}$.

The terms we have derived suggest a simple procedure for the estimation. These
estimators $\widehat{s}^{(1)}(t;\beta),\widehat{s}^{(0)}(t;\beta),$ and
$U_{n}(\beta)$ are all formed by merging the $M$ imputed datasets together to
form one large completed dataset with $Mn$ entries. Thus, this regression
adjustment method with a Monte Carlo approximation is equivalent to a stacked
multiple imputation approach utilized by Robins and Wang (2000); Beesley and
Taylor (2020). Note that this approach differs from a traditional imputation
approach where $M$ imputed datasets are generated separately, analysis is
performed on each one, and each analysis is pooled to form a final estimator.
In contrast, in the stacked imputation approach, the pooling operation is done
before the analysis such that $M$ imputed datasets are pooled first and
analysis is performed on all $Mn$ observations to form a final point
estimator.

An advantage of the traditional imputation approach is that the separate $M$
analyses allows for parallelization. On the other hand, the stacked imputation
does not exhibit this nice property, but a single analysis yields a more
straightforward implementation, which can be helpful in situations where we
perform the bootstrap.

Algorithm 2 Cox model using multiple imputation under CCMV for $X$

  1. 1.

Input: $\\{(X_{i,R_{i}},R_{i},Y_{i},\Delta_{i})\\}_{i=1}^{n}$

  2. 2.

For $\delta=0,1$, fit the parametric model $p(x|R=1_{d},y,\delta).$

  3. 3.

For $j=1,\ldots,M$, do

    1. 3-1.

For $i=1,\ldots,n$, do

      1. i.

Impute the missing covariates $\widetilde{X}_{i,\bar{R}_{i}}^{(j)}$ for
observation $(X_{i,R_{i}},R_{i},Y_{i},\Delta_{i})$ using the estimated
extrapolation distribution $p(x_{\bar{r}}|x_{r},R=1_{d},y,a)$.

      2. ii.

Set the $j$th imputation for the $i$th observation to be
$(\widetilde{X}_{i}^{(j)},R_{i},Y_{i},\Delta_{i})$, where
$\widetilde{X}_{i}^{(j)}=\left(X_{i,R_{i}},\widetilde{X}_{i,\bar{R}_{i}}^{(j)}\right)$.

    2. 3-2.

Let
$\widetilde{D}^{(j)}=\\{(\widetilde{X}_{i}^{(j)},R_{i},Y_{i},\Delta_{i})\\}_{i=1}^{n}$
denote the $j$th imputed data set.

  4. 4.

Pool all $\widetilde{D}^{(j)},j=1,\ldots,M$ data sets together into a single
large one $\check{D}$.

  5. 5.

Fit a Cox model using the completed data set $\check{D}$ to get an estimator
$\widehat{\beta}$.

  6. 6.

Return: $\widehat{\beta}$.

In a frequentist multiple imputation approach, we first use the observed data
to estimate the extrapolation densities and thereby, perform the imputation.
Therefore, all further analysis is based on the estimated extrapolation
densities, and it is important to incorporate this uncertainty in the final
variance calculation. Rubin’s estimator for the variance is commonly used in
practice, but it will underestimate the variance in this situation.
Specifically, there will be an additional term in the asymptotic variance
related to the estimator for the extrapolation distribution that will not be
included. See Sections 14.4 and 14.5 of Tsiatis (2007) for more details.

To avoid this issue in practice, we recommend using the nonparametric
bootstrap to estimate the variance. We bootstrap the data (including the
missing values) many times and repeat the same analysis on each bootstrapped
sample to obtain many point estimates. Each analysis on a bootstrapped sample
will utilize a different estimated extrapolation distribution, so this
uncertainty will be taken into account. Chen and Sadinle (2019) provide
additional related details in their work.

###  3.3 A multiply-robust estimator

In the previous two sections, we have described two approaches based on an
inverse probability weighting and a regression adjustment. In the seminal
paper by Robins et al. (1994), the authors introduced a general framework of
an augmented inverse probability weighted (AIPW) estimating function under a
MAR assumption for the covariates, and thus, the selection probability did not
depend on any missing covariates. Wang and Chen (2001) discussed further
details on constructing an AIPW estimator with augmented functions for the Cox
model under MAR. Additionally, in their paper, the covariates took the form of
$(X_{\text{miss}},X_{\text{obs}})$, where $X_{\text{miss}}$ was a covariate
vector possibly missing for some individuals and $X_{\text{obs}}$ was a
covariate vector that was always observed, thereby inducing only two possible
missingness patterns for the covariates. These AIPW methods rely on two
parametric models for the selection probabilities and the missing covariates.
To relax this parametric assumption, Qi et al. (2005) outlined a fully
augmented weighted estimator using a kernel regression approach for a one-
dimensional covariate. However, it is unclear how to handle the case with
multiple covariates that may all be missing.

We extend the previous work by describing a multiply-robust estimator that
allows missing not at random (MNAR) for a multidimensional covariate and
arbitrary missing patterns. The following lemma shows that there exists an
estimator that is robust to the misspecification of at most one of
$Q_{r,a}(X_{r},W_{a})$ or $m_{r,a}(X_{r},W_{a})$ for each $r$.

######  Lemma 6 (Doubly robust linear form under CCMV).

Fix a missing pattern $r$ and an observed outcome pattern $a$. Define

$\displaystyle g_{r}(X,W_{a})$
$\displaystyle=f_{a}(X,W_{a})1(A=a)1(R=1_{d})Q_{r,a}(X_{r},W_{a})$ (12)
$\displaystyle\quad+m_{r,a}(X_{r},W_{a})1(R=r)1(A=a)$ (13)
$\displaystyle\quad-m_{r,a}(X_{r},W_{a})1(A=a)1(R=1_{d})Q_{r,a}(X_{r},W_{a}).$
(14)

If $Q_{r,a}(X_{r},W_{a})$ or $m_{r,a}(X_{r},W_{a})$ is correctly specified,
then under a CCMV assumption of equation (6), we have

$\mathbb{E}[f_{a}(X,W_{a})1(A=a)]=\sum_{r}\mathbb{E}[g_{r}(X,W_{a})].$

The expression of Lemma 6 shows an interesting connection to equation (10.95)
of Tsiatis (2007). In the case where the there is only two missing patterns of
the covariates, i.e., the covariates can be expressed as
$(X_{\text{miss}},X_{\text{obs}})$, where $X_{\text{miss}}$ are possibly
missing covariates and $X_{\text{obs}}$ are completely observed covariates,
the MAR and CCMV assumptions on the covariates can be shown to be equivalent.
Our expression and equation (10.95) in Tsiatis (2007) will be the same.

######  Theorem 7 (Multiply-robust Cox model).

Under CCMV, the estimator $\widehat{\beta}_{\text{mr}}$ that solves the
augmented estimating function is consistent provided either
$Q_{r,\delta}(X_{r},Y)$ or $m_{r,\delta}(X_{r},Y)$ is correctly specified for
each $r$.

Thus, this implies that we can construct a multistage estimator based on Lemma
6 that is multiply-robust since we have doubly robustness for each covariate
pattern $r\neq 1_{d}$. However, this comes at a cost. We see that the doubly
robust expression above contains the terms from Lemmas 1 and 5 along with an
additional term. Term (14) is of interest because it shows that in order to
compute it, we will have to perform imputation even for cases with completely
observed covariates. Therefore, we can see explicitly that this doubly robust
estimator is even more computationally intensive than the regression
adjustment estimator in Lemma 5.

##  4 Simulation

We examine two scenarios of the Cox model and the binary treatment problem to
observe the consistency of our proposed estimators as well as behavior of the
estimators under perturbations of the assumed identifying restriction.

###  4.1 Cox model simulation

We consider a Cox model with two binary covariates
$X=(X_{1},X_{2})\in\\{0,1\\}^{2}$ from the following generative model:

  1. 1.

Sample $X=(X_{1},X_{2})$ such that $P(X_{1}=1)=1/2$, $P(X_{2}=1)=3/10$, and
$\text{Corr}(X_{1},X_{2})=0.3$.

  2. 2.

For $T$ and $C$ (and $Y$ and $\Delta$), we do the following:

     * •

Generate $T|X$ from an Exponential distribution with rate
$\exp(x^{\top}\beta)$ and $\beta=\begin{bmatrix}-0.5\\\ 2\end{bmatrix}$. This
is based on the baseline hazard $\lambda_{0}(t)=1$.

     * •

Generate $C|X$ from an distribution with rate $2$.

     * •

Compute $Y=\min\\{T,C\\}$ and $\Delta=1(T\leq C)$.

  3. 3.

Generate the missing pattern of covariates based on CCMV:

$\displaystyle\dfrac{P(R=01|x,y,\delta)}{P(R=11|x,y,\delta)}=\dfrac{P(R=01|x_{2},y,\delta)}{P(R=11|x_{2},y,\delta)}=\exp(-0.5x_{2}+0.75y+0.5\delta),$
$\displaystyle\dfrac{P(R=10|x,y,\delta)}{P(R=11|x,y,\delta)}=\dfrac{P(R=10|x_{1},y,\delta)}{P(R=11|x_{1},y,\delta)}=\exp(-1+0.5x_{1}+y+\delta),$

and remove the covariates accordingly.

Roughly, we expect approximately $40\%$, $30\%$, and $30\%$ of the
observations to fall in the classes $R=01$, $R=10$, and $R=11$, respectively.
This data generating mechanism completely specifies the joint model, and we
consider three different estimators. The first estimator is the IPW estimator,
which is based off the weighted Cox model from Section 3.1, and the complete
odds are estimated using logistic regression under CCMV. The second estimator
is multiple imputation (MI) / regression adjustment estimator as discussed in
Section 3.2 with $M=50$ imputations. Notably, the complete odds and the
imputation distribution are not variation independent in this problem (and
will generally not be in most problems of interest), so we have to fit the
imputation distribution using knowledge of the complete odds. The third
estimator is called the transformed MLE estimator. The specific parametric
model assumed for $p_{Y,\Delta|X}(y,\delta|x)$ implies that we can recover the
terms

$\nu_{1},\quad\nu_{1}\exp(\beta_{1}),\quad\nu_{1}\exp(\beta_{2}),\quad\text{and}\quad\nu_{1}\exp(\beta_{1}+\beta_{2})$

via maximum likelihood through the EM algorithm (details provided in Appendix
E.1). Let $\widehat{\gamma}_{i}$ denote the estimator for the $i$th term above
($i=1,2,3,4$). Two naive estimators, which we term the transformed MLEs, are
as follows:

$\displaystyle\widehat{\beta}_{1,\text{Transformed
MLE}}=\frac{1}{2}\log(\widehat{\gamma}_{2}/\widehat{\gamma}_{1})+\frac{1}{2}[\log(\widehat{\gamma}_{4}/\widehat{\gamma}_{1})-\log(\widehat{\gamma}_{3}/\widehat{\gamma}_{1})],$
$\displaystyle\widehat{\beta}_{2,\text{Transformed
MLE}}=\frac{1}{2}\log(\widehat{\gamma}_{3}/\widehat{\gamma}_{1})+\frac{1}{2}[\log(\widehat{\gamma}_{4}/\widehat{\gamma}_{1})-\log(\widehat{\gamma}_{2}/\widehat{\gamma}_{1})].$

We did not include the DR estimator in comparison because it has a very
complex form (see Lemma 6).

Figure 1: Estimation of $\beta$ in a simulated Cox model for varying sample
sizes.

For each of the sample sizes (350, 500, 1000, 2000), we generate 1000 data
sets and construct the three estimators. The means and the standard errors are
calculated across the 1000 replicates, and we show the 95% confidence
interval. The results are summarized in Figure 1. We first observe that all
three estimators are consistent and have standard errors that appear to
roughly decline at the appropriate $\sqrt{n}$-rate. The IPW estimator
(equivalent to the reweighted Cox model) has the highest variance, which is to
be expected because it uses the least information and inverse probabilities
are known to introduce instability. The multiple imputation / regression
adjustment (MI/RA) estimator uses slightly more information than the IPW
estimator by assuming a parametric model for $p_{Y,\Delta|X}(y,\delta|x)$. As
the imputation distribution and the complete odds are typically not variation
independent, we are forced to fit them together. We see that the RA estimator
appears to have similar variance to the transformed MLE, which is to be
expected for large enough $M$ (number of imputations).

###  4.2 Binary treatment problem

We perform another simulation under the binary treatment setting from causal
inference to show consistency under a new setting. The generating process is
described in Appendix E.2. For each of the samples sizes (1000, 2000, 5000,
10000), we generate 1000 data sets and construct $6$ estimators. We consider
two methods for accounting for the missing covariates (IPW-$R$ and RA-$R$) and
three methods for estimating the ATE (IPW-$A$, RA-$A$ and DR-$A$). So, there
are 6 combinations of one method for missing covariates and one method for the
ATE. The left panel of Figure 2 displays the results of three estimators
formed by using an IPW approach for the missing covariates: (IPW-$R$,
IPW-$A$), (IPW-$R$, RA-$A$), and (IPW-$R$, DR-$A$). In this case, we use the
bootstrap percentile method for computing the confidence intervals. The right
panel of Figure 2 provides the results of three estimators using a RA
(imputation) approach for the missing covariates: (RA-$R$, IPW-$A$), (RA-$R$,
RA-$A$), and (RA-$R$, DR-$A$). We use $M=50$ imputations and report the 95%
confidence intervals. We simply construct the confidence intervals by
computing the standard errors across all 1000 generated data sets. Again, we
also report the means across the 1000 replicates.

Figure 2: Estimation of the ATE in a simulated causal inference problem for
varying sample sizes.

In the left panel of Figure 2, we observe some surprising results. First, in
finite samples, the (IPW-$R$, IPW-$A$) estimator can lead to unstable
estimates for some data sets due to the fact that the complete odds are
estimated using a density ratio. Although we expect the doubly-robust
estimator (IPW-$R$, DR-$A$) to have smaller variance than the (IPW-$R$,
IPW-$A$) estimator, we see that it has larger variance. This is because the
unstable part of IPW-$R$ appears in each of the 3 terms of the DR-$A$
estimator.

In the right panel of Figure 2, we see very similar performance across all 3
estimators. These 3 estimators dominate the other 3 estimators in the right
panel, which suggests that a multiple imputation approach can provide more
stable results in practice. We discuss further sensitivity analysis for this
simulation in Appendix D.2.

##  5 Data analysis

We apply our methods on the data from the National Alzheimer’s Coordinating
Center’s (NACC) Uniform Data Set111https://naccdata.org/. This data set
contains longitudinal data collected at Alzheimer’s Disease Research Centers,
funded by the National Institutes of Health and starting from the year 2005.
Participating individuals have varying cognitive statuses from no cognitive
impairment to mild impairment to severe impairment, and data is collected on a
yearly basis with certain covariates measured at entry into the study.

We consider individuals who display mild cognitive impairment (MCI) and follow
their trajectory until they either display evidence of dementia (mild,
moderate, severe) or are censored. These individuals are binned into 5 age
groups based on when they first display MCI: 60-64, 65-69, 70-74, 75-79, and
80-84. In each age group, we fit a Cox model, regressing on the binary
variables: sex (female or not), race (white or not), education (have graduate
degree or not), depression within last 2 years (yes or no), and anxiety (yes
or no). Tracking the change in the $\beta$ across the different age groups
allows us to observe the age effect. The first three covariates (sex, race,
and education) are generally measured for all individuals at the baseline.
There are a few individuals for which the education status is unknown, but
because these represent less than $1\%$ of the observations, we exclude this
individuals from the analysis for simplicity. The last two covariates
(depression and anxiety) are subject to missingness, yielding a total of
$2^{2}=4$ covariate missing patterns. There are about $16\%$ participants in
our data with complete observations.

We consider three different estimators: complete case (CC) analysis, IPW and
RA estimators. The CC analysis is a common approach used by practitioners. In
the CC analysis, we remove any individuals with any missing covariates and fit
a standard Cox model. For the IPW estimator, we use the CCMV assumption using
logistic models for the complete odds. Since the variables are binary, the
imputation distribution in the RA estimator is the usual nonparametric model.
Point estimates and the 95% confidence intervals (constructed using the
nonparametric bootstrap) are provided in Figure 3.

Figure 3: Point estimates and the $95\%$ confidence intervals for each $\beta$
parameter. The total sample size and the number of complete cases is included
on the $x$-axis for each age group.

First, we observe that the results from a CC analysis are generally not
significant. For example, in the plot for $\beta_{3}$ (education), the CC
confidence intervals all strongly cover $0$. On the other hand, across the
plots, the RA approach generally yields narrower confidence intervals while
the IPW approach sometimes yields narrower confidence intervals; this is due
to the instability introduced by the inverse probability weights.
Specifically, for $\beta_{3}$ (education), the RA provides clear significant
results for several of the age groups. These plots suggest that often, when
applying a CC analysis, we suffer from such a loss in sample size that it
becomes difficult to draw a significant conclusion.

Moreover, it is known that higher education provides a protection against
cognitive impairment (known as the cognitive reserve; Meng and D’Arcy 2012;
Thow et al. 2018). The point estimates for $\beta_{3}$ from all three
estimation procedures are negative, which is what we would expect. As the RA
approach allows one to attain significant negative point estimates, it
highlights the potential benefits of accounting for the missing data. We
discuss further details for the analysis in Appendix F.

###  5.1 Sensitivity analysis

We conduct a sensitivity analysis to investigate the robustness of our
conclusion under CCMV. We perturb the complete odds using an exponential
tilting parameter $\rho$ such that

$O_{r,a}(X,W_{a})=Q_{r,a}(X_{r},W_{a})e^{\rho_{\bar{r}}X_{\bar{r}}}$

and every element in the vector $\rho_{\bar{r}}$ takes the same value $\rho$.
When $\rho=0$, there is no perturbation, but as $|\rho|\to\infty$, the scale
of the perturbation increases. For the imputation model, we use a sensitivity
analysis parameter $\zeta\in[0,1]$. When $\zeta\to 1$, the imputation
probability moves towards $0$, and when $\zeta\to 0$, the imputation
probability moves towards $1$. If $\zeta=1/2$, then the imputation probability
remains unchanged (i.e., the case of CCMV). We provide additional discussion
of these sensitivity parameters in Appendix D. The results are displayed in
Figures 4 and 5. The confidence intervals are shaded in gray. There is some
slight wiggles in the upper and lower bands because of the Monte Carlo errors
from the bootstrap (we construct each confidence interval using the
nonparametric bootstrap).

Figure 4: Point estimates and the $95\%$ confidence intervals for the IPW
estimator ($\beta_{3}$ and $\beta_{4}$) as we vary the sensitivity parameter
$\rho$ for the age groups $75-79$ and $80-84$.

Figure 5: Point estimates and the $95\%$ confidence intervals for the RA
estimator ($\beta_{3}$ and $\beta_{4}$) as we vary the sensitivity parameter
$\zeta$ for the age groups $75-79$ and $80-84$.

In Figure 4, we display the result of the IPW estimator. The IPW estimator is
robust when we perturb the missing data assumption in the sense that only the
education effect has a stronger dependency with the sensitivity parameter.
Across all $\rho\in[-1,1]$, the education effect is still negative, which
suggests that the cognitive reserve we obtain from the IPW estimator is
resilient.

In Figure 5, we observe that the RA estimator for $\beta_{3}$ (education
effect) is robust to changes in the sensitivity parameter, which again
reinforces our findings on the cognitive reserve.

On the other hand, we see that the confidence interval width of $\beta_{4}$
(depression) depends on the sensitivity parameter $\zeta$. While the sign of
the effect remains the same and the effect (value of $\beta_{4}$) is also
stable, the width of confidence intervals do change drastically. So, we have
to be careful when making claims regarding $\beta_{4}$.

##  6 Discussion

In this paper, we introduced a process for constructing a multistage estimator
for problems with missing not at random covariates and incomplete outcomes. We
discussed our method in the context of Cox regression from survival analysis
and a missing response problem (in Appendix A). Although we did not show the
specific details, the missing response problem provides a bridge to the binary
treatment problem from causal inference and a more general longitudinal study
problem with dropout. So, the arguments we provide for the missing response
problem naturally extend to these two problems.

Model Complexity (RA/Imputation)Computation (IPW)Efficiency (DR)

Figure 6: The tradeoff between model complexity, computation, and efficiency
in the CCMV setting.

Figure 6 shows an interesting tradeoff between model complexity, computation,
and efficiency in the special case of the CCMV assumption for the covariates.
The RA estimator is the easiest to model since we only need to place a model
on $p(x|1_{d},w_{a},a)$ for all $a$. But this comes at a cost of computation
since we generally need to perform imputation to fill in the missing data.
Also, RA does not achieve the (semiparametric) efficiency bound. The IPW
estimator is easy to compute because it is a simple reweighting of the cases
with completely observed covariates. But we need to model every
$Q_{r,a}(X_{r},W_{a})$, which would require many models. Also, IPW does not
achieve the (semiparametric) efficiency bound and may be unstable if the
weights are too large. The DR estimator achieves the semiparametric efficiency
bound by utilizing components from both the regression adjustment and inverse
probability weighted estimator. Due to the augmenting term, imputation must be
performed even for the complete cases, so the computational cost is high and
we have to place models for every $Q_{r,a}(X_{r},W_{a})$. Thus, efficiency
comes at the cost of complexity in modeling and high computation.

There are several possible areas for future work. One direction is to study
the problem of model compatibility issues between complete odds and regression
functions. When the assumption on missing covariates is not CCMV, the complete
odds and regression function may again not be variationally independent. Thus,
extra caution has to be taken when placing models on them. Another direction
is to study the use of nonparametric models in complete odds and regression
functions. This would lead to a more flexible modeling strategy and avoid the
problem of model compatibility, but overcoming the slow convergence rate will
be a challenge.

## References

  * Beesley and Taylor (2020) L. Beesley and J. Taylor.  A stacked approach for chained equations multiple imputation incorporating the substantive model.  _Biometrics_ , 09 2020.  doi: 10.1111/biom.13372.
  * Beesley et al. (2016) L. J. Beesley, J. W. Bartlett, G. T. Wolf, and J. M. G. Taylor.  Multiple imputation of missing covariates for the cox proportional hazards cure model.  _Statistics in Medicine_ , 35(26):4701–4717, 2016\.  doi: https://doi.org/10.1002/sim.7048.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.7048.
  * Berg and Yu (2021) E. Berg and C. Yu.  Estimation for nonignorable missing response or covariate using semi-parametric quantile regression imputation and a parametric response probability model.  _Statistica Sinica_ , 2021.  doi: 10.5705/ss.202020.0053.  URL https://par.nsf.gov/biblio/10297276.
  * Chen (2020) Y.-C. Chen.  Pattern graphs: a graphical approach to nonmonotone missing data.  _arXiv preprint arXiv:2004.00744_ , 2020.
  * Chen and Sadinle (2019) Y.-C. Chen and M. Sadinle.  Nonparametric pattern-mixture models for inference with missing data.  _arXiv preprint arXiv:1904.11085_ , 2019.
  * Crowe et al. (2010) B. J. Crowe, I. A. Lipkovich, and O. Wang.  Comparison of several imputation methods for missing baseline data in propensity scores analysis of binary outcome.  _Pharmaceutical statistics_ , 9(4):269–279, 2010\.
  * Kim and Yu (2011) J. K. Kim and C. L. Yu.  A semiparametric estimation of mean functionals with nonignorable missing data.  _Journal of the American Statistical Association_ , 106(493):157–165, 2011.
  * Kitching (2015) D. Kitching.  Depression in dementia.  _Australian Prescriber_ , 38:209–211, 12 2015.  doi: 10.18773/austprescr.2015.071.
  * Little (1993) R. J. Little.  Pattern-mixture models for multivariate incomplete data.  _Journal of the American Statistical Association_ , 88(421):125–134, 1993.
  * Little and Rubin (2002) R. J. A. Little and D. B. Rubin.  _Statistical Analysis with Missing Data_.  Wiley, Hoboken, New Jersey, 2nd edition, 2002.
  * Malinsky et al. (2020) D. Malinsky, I. Shpitser, and E. J. T. Tchetgen.  Semiparametric inference for non-monotone missing-not-at-random data: the no self-censoring model, 2020.
  * Mayer et al. (2020) I. Mayer, E. Sverdrup, T. Gauss, J.-D. Moyer, S. Wager, and J. Josse.  Doubly robust treatment effect estimation with missing attributes.  _The Annals of Applied Statistics_ , 14(3):1409 – 1431, 2020.  doi: 10.1214/20-AOAS1356.  URL https://doi.org/10.1214/20-AOAS1356.
  * Mayer et al. (2021) I. Mayer, J. Josse, and T. Group.  Transporting treatment effects with incomplete attributes, 2021.
  * Meng and D’Arcy (2012) X. Meng and C. D’Arcy.  Education and dementia in the context of the cognitive reserve hypothesis: A systematic review with meta-analyses and qualitative analyses.  _PLoS ONE_ , 7, 2012.
  * Meng et al. (1994) X.-L. Meng et al.  Multiple-imputation inferences with uncongenial sources of input.  _Statistical Science_ , 9(4):538–558, 1994.
  * Mohan and Pearl (2014) K. Mohan and J. Pearl.  Graphical models for recovering probabilistic and causal queries from missing data.  In _Advances in Neural Information Processing Systems_ , pages 1520–1528, 2014.
  * Mohan and Pearl (2021) K. Mohan and J. Pearl.  Graphical models for processing missing data.  _Journal of the American Statistical Association_ , 116(534):1023–1037, 2021.  doi: 10.1080/01621459.2021.1874961.  URL https://doi.org/10.1080/01621459.2021.1874961.
  * Mohan et al. (2013) K. Mohan, J. Pearl, and J. Tian.  Graphical models for inference with missing data.  In _Advances in neural information processing systems_ , pages 1277–1285, 2013.
  * Molenberghs et al. (1998) G. Molenberghs, B. Michiels, M. G. Kenward, and P. J. Diggle.  Monotone missing data and pattern-mixture models.  _Statistica Neerlandica_ , 52(2):153–161, 1998\.
  * Nabi et al. (2020) R. Nabi, R. Bhattacharya, and I. Shpitser.  Full law identification in graphical models of missing data: Completeness results.  _arXiv preprint arXiv:2004.04872_ , 2020.
  * Paik (1997) M. Paik.  Multiple imputation for the cox proportional hazards model with missing covariates.  _Lifetime data analysis_ , 3:289–98, 02 1997.  doi: 10.1023/A:1009657116403.
  * Qi et al. (2005) L. Qi, C. Wang, and R. L. Prentice.  Weighted estimators for proportional hazards regression with missing covariates.  _Journal of the American Statistical Association_ , 100(472):1250–1263, 2005.
  * Qu and Lipkovich (2009) Y. Qu and I. Lipkovich.  Propensity score estimation with missing values using a multiple imputation missingness pattern (mimp) approach.  _Statistics in medicine_ , 28(9):1402–1414, 2009\.
  * Robins and Gill (1997) J. M. Robins and R. D. Gill.  Non-response models for the analysis of non-monotone ignorable missing data.  _Statistics in medicine_ , 16(1):39–56, 1997\.
  * Robins and Wang (2000) J. M. Robins and N. Wang.  Inference for imputation estimators.  _Biometrika_ , 87(1):113–124, 2000.  ISSN 00063444.  URL http://www.jstor.org/stable/2673565.
  * Robins et al. (1994) J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American Statistical Association_ , 89(427):846–866, 1994.  ISSN 01621459.
  * Robins et al. (2000) J. M. Robins, A. Rotnitzky, and D. O. Scharfstein.  Sensitivity analysis for selection bias and unmeasured confounding in missing data and causal inference models.  In _Statistical models in epidemiology, the environment, and clinical trials_ , pages 1–94. Springer, 2000.
  * Rosenbaum and Rubin (1984) P. R. Rosenbaum and D. B. Rubin.  Reducing bias in observational studies using subclassification on the propensity score.  _Journal of the American statistical Association_ , 79(387):516–524, 1984.
  * Rosenbaum and Rubin (1985) P. R. Rosenbaum and D. B. Rubin.  Constructing a control group using multivariate matched sampling methods that incorporate the propensity score.  _The American Statistician_ , 39(1):33–38, 1985\.
  * Rubin (1974) D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688, 1974.
  * Rubin (1978) D. B. Rubin.  Bayesian inference for causal effects: The role of randomization.  _The Annals of Statistics_ , 6(1):34–58, 1978\.
  * Sadinle and Reiter (2017) M. Sadinle and J. P. Reiter.  Itemwise conditionally independent nonresponse modelling for incomplete multivariate data.  _Biometrika_ , 104(1):207–220, 2017.
  * Shao and Wang (2016) J. Shao and L. Wang.  Semiparametric inverse propensity weighting for nonignorable missing data.  _Biometrika_ , 103(1):175–187, 2016.
  * Shpitser (2016) I. Shpitser.  Consistent estimation of functions of data missing non-monotonically and not at random.  In _Advances in Neural Information Processing Systems_ , pages 3144–3152, 2016.
  * Splawa-Neyman et al. (1990) J. Splawa-Neyman, D. M. Dabrowska, and T. P. Speed.  On the application of probability theory to agricultural experiments. essay on principles. section 9.  _Statistical Science_ , 5(4):465–472, 1990.  ISSN 08834237.
  * Sugiyama et al. (2012) M. Sugiyama, T. Suzuki, and T. Kanamori.  _Density Ratio Estimation in Machine Learning_.  Cambridge University Press, 2012.  doi: 10.1017/CBO9781139035613.
  * Tchetgen Tchetgen et al. (2018) E. J. Tchetgen Tchetgen, L. Wang, and B. Sun.  Discrete choice models for nonmonotone nonignorable missing data: Identification and inference.  _Statistica Sinica_ , 28(4):2069–2088, 2018.
  * Thow et al. (2018) M. E. Thow, M. J. Summers, N. L. Saunders, J. J. Summers, K. Ritchie, and J. C. Vickers.  Further education improves cognitive reserve and triggers improvement in selective cognitive functions in older adults: The tasmanian healthy brain project.  _Alzheimer’s & dementia : diagnosis, assessment & disease monitoring_, 10(1):22–30, 2018.  ISSN 2352-8729.
  * Tsiatis (2007) A. Tsiatis.  _Semiparametric theory and missing data_.  Springer Science & Business Media, 2007.
  * van der Vaart (1998) A. W. van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 1998.
  * van der Vaart and Wellner (1996) A. W. van der Vaart and J. A. Wellner.  _Weak convergence_.  Springer, 1996.
  * Wang and Chen (2001) C. Y. Wang and H. Y. Chen.  Augmented inverse probability weighted estimator for cox missing covariate regression.  _Biometrics_ , 57(2):414–419, 2001.  ISSN 0006341X, 15410420.
  * Wang and Kim (2021) H. Wang and J. K. Kim.  Propensity score estimation using density ratio model under item nonresponse, 2021.
  * Wang and Robins (1998) N. Wang and J. M. Robins.  Large-sample theory for parametric multiple imputation procedures.  _Biometrika_ , 85(4):935–948, 1998.  ISSN 00063444.
  * White and Royston (2009) I. R. White and P. Royston.  Imputing missing covariate values for the cox model.  _Statistics in Medicine_ , 28(15):1982–1998, 2009\.  doi: https://doi.org/10.1002/sim.3618.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.3618.
  * Zhao et al. (2017) P. Zhao, N. Tang, A. Qu, and D. Jiang.  Semiparametric estimating equations inference with nonignorable missing data.  _Statistica Sinica_ , pages 89–113, 2017.

## Acknowledgements.

DS is supported by NSF DGE-2140004. YC is supported by NSF DMS-1952781, NSF
DMS-2112907, and NIH U24-AG072122.

Supplementary Materials: Multistage Estimators for Missing Covariates and
Incomplete Outcomes

Appendices

The appendix has several sections to provide additional theoretical results,
proofs, details on simulations and data analyses, and further generalizations
of work in the main body. In Appendix A, we construct a multistage estimator
for the missing response problem. The assumptions for the asymptotic theory
are provided in Appendix B. Since we utilize CCMV as the main missing data
assumption in the main body, we generalize this assumption in Appendix C and
show how to attain identifiable complete odds and imputation distributions.
Further details on the sensitivity analyses and simulations are included in
Appendices D and E, respectively. We include more data analyses in Appendix F.
The proofs for any theoretical results are included in Appendix G.

The table of contents is provided below.

  1. 1.

Appendix A: Missing response problem.

  2. 2.

Appendix B: Technical assumptions for the asymptotic theory.

  3. 3.

Appendix C: Beyond CCMV.

  4. 4.

Appendix D: Procedure of sensitivity analysis.

  5. 5.

Appendix E: Details of simulations.

  6. 6.

Appendix F: Additional analysis of the NACC data.

  7. 7.

Appendix G: Proofs of 1. identification theories (Appendix G.1), 2. efficiency
theories (Appendix G.2), 3. asymptotic theories (Appendix G.3)

##  Appendix A A multistage estimator for the missing response problem

We now turn our attention to the missing response problem and show how to
construct a multistage estimator in this setting. Note that this procedure can
be easily generalized to the binary treatment effect in causal inference. The
parameter of interest is the marginal mean of the outcome variable
$\mathbb{E}[Y]$. Let $\pi_{A}(X)=P(A=1|X)$ and $m_{1}(X)$ denote the
propensity score and regression function for $Y$, respectively. We start with
the classic DR estimator assuming fully-observed covariates for the missing
response problem. The efficient influence function takes the form

$\frac{1(A=1)Y}{\pi_{A}(X)}+m_{1}(X)-\frac{1(A=1)m_{1}(X)}{\pi_{A}(X)}.$

The expectation of the above quantity equals $\mathbb{E}[Y]$ as long as at
least one of $\pi_{A}(X)$ and $m_{1}(X)$ is correctly specified. Using
empirical average of the above quantity, we obtain the DR estimator under the
missing response problem.

###  A.1 Multiply-robust property

In this section, we highlight the multiply-robust property of the multistage
estimator by demonstrating the construction of one such multiply-robust
estimator in the missing response problem.

For this problem, we also define the following regression functions

$\displaystyle M_{m,r,0}(X_{r})=\mathbb{E}[m_{1}(X)|X_{r},R=r,A=0],$ (15)
$\displaystyle M_{m,r,1}(X_{r},Y)=\mathbb{E}[m_{1}(X)|X_{r},R=r,Y,A=1],$ (16)
$\displaystyle
M_{1/\pi_{A},r,1}(X_{r},Y)=\mathbb{E}\left[\frac{1}{\pi_{A}(X)}\biggr{|}X_{r},R=r,Y,A=1\right],$
(17) $\displaystyle
M_{m/\pi_{A},r,1}(X_{r},Y)=\mathbb{E}\left[\frac{m_{1}(X)}{\pi_{A}(X)}\biggr{|}X_{r},R=r,Y,A=1\right].$
(18)

######  Definition 4 (Semiparametric submodels).

Define the following semiparametric submodels.

  * •

$\mathcal{M}_{R}(r)$: The semiparametric model such that the parametric model
$p(x_{\bar{r}}|x_{r},r,w_{a},a;\eta_{r,a})$ is correctly specified.

  * •

$\mathcal{M}_{O}(r)$: The semiparametric model such that the parametric model
$O_{r,a}(x,w_{a};\alpha_{r,a})$ is correctly specified.

It follows from the previous definition that under model $\mathcal{M}_{R}(r)$,
the regression functions (15)-(18) are correctly specified. We specify the
semiparametric model in this manner to avoid the model congeniality issue
described in Section 3.2. Note that one could instead define
$\mathcal{M}_{R}(r)$ to be the semiparametric model under which these
regression functions are correctly specified.

######  Theorem 8 (Double robustness for missing response under CCMV for given
$r$).

Let $\theta_{r}=\int y\ P(dx,r,dy,da)$ be the parameter of interest. Consider
the uncentered influence function

$\displaystyle q_{r}(X,R,Y,A)=m_{1}(X)1(A=0)1(R=1_{d})Q_{r,0}(X_{r})$
$\displaystyle\quad+m_{1}(X)1(A=1)1(R=1_{d})Q_{r,1}(X_{r},Y)$
$\displaystyle\quad+M_{m,r,0}(X_{r})1(A=0)[1(R=r)-1(R=1_{d})Q_{r,0}(X_{r})]$
$\displaystyle\quad+M_{m,r,1}(X_{r},Y)1(A=1)[1(R=r)-1(R=1_{d})Q_{r,1}(X_{r},Y)]$
$\displaystyle\quad+\frac{1(R=1_{d})1(A=1)Y}{\pi_{A}(X)}Q_{r,1}(X_{r},Y)$
$\displaystyle\quad+M_{1/\pi_{A},r,1}(X_{r},Y)Y1(A=1)[1(R=r)-1(R=1_{d})Q_{r,1}(X_{r},Y)]$
$\displaystyle\quad-\frac{1(R=1_{d})1(A=1)m_{1}(X)}{\pi_{A}(X)}Q_{r,1}(X_{r},Y)$
$\displaystyle\quad-
M_{m/\pi_{A},r,1}(X_{r},Y)1(A=1)[1(R=r)-1(R=1_{d})Q_{r,1}(X_{r},Y)].$

Under the model $\mathcal{M}_{R}(r)\cup\mathcal{M}_{O}(r)$ (i.e., either the
complete odds or the regression functions are correct),
$\mathbb{E}[q_{r}(X,R,Y,A)]=\theta_{r}$.

Theorem 8 shows that the estimator for $\theta_{r}$ based on the provided
influence function is doubly robust. One key takeaway is that for a fairly
simple problem, the estimation procedure for the doubly robust estimator is
rather complex. In particular, there are 12 terms in total. With this in mind,
we expect that this doubly robust estimator to be more of a theoretical
construct than one that will generally be implemented in practice.

We further develop the theory by considering a doubly robust estimator for
$\theta_{r}$ for each $r$. Merging this together, we are able to estimate

$\theta=\int y\ P(dx,dr,dy,da)=\sum_{r}\theta_{r}.$

We see that in the following corollary, this yields a multiply-robust
estimator for $\theta=\mathbb{E}[Y]$.

######  Corollary 9 (Multiply-robustness for missing response under CCMV).

Define the parameter of interest

$\theta=\sum_{r}\theta_{r}=\sum_{r}\int y\ P(dx,r,dy,da)=\int y\
P(dx,dr,dy,da).$

The uncentered influence function $\sum_{r}q_{r}(X,R,Y,A)$ yields a consistent
estimator of $\theta=\sum_{r}\theta_{r}$ under the model
$\bigcap_{r}(\mathcal{M}_{R}(r)\cup\mathcal{M}_{O}(r))$.

Corollary 9 provides a strong result in that the estimator is
$2^{|\mathcal{R}|-1}$-robust, where $|\mathcal{R}|$ is the number of missing
covariate patterns. As shown in Theorem 8, for an estimator of $\theta_{r}$ to
be consistent, it is sufficient for at least one of $O_{r,a}(X,W_{a})$ or the
distribution $p(x_{\bar{r}}|x_{r},w_{a},a)$ to be correctly specified. In
other words, we have two chances to obtain a correct estimator for each
missing covariate pattern $r\neq 1_{d}$.

Since the missing response and binary treatment problems are related and well-
studied problems, we omit the construction of our multistage estimator for the
binary treatment problem. The included arguments can easily be extended to the
binary treatment problem and the problem of longitudinal data with monotone
missingness in the response.

###  A.2 Semiparametric efficiency

We now show that the multiply-robust estimator developed in Theorem 8 and
Corollary 9 is the semiparametric efficient estimator that attains the minimum
possible variance in the model
$\bigcap_{r}(\mathcal{M}_{R}(r)\cap\mathcal{M}_{O}(r))$, where all the
regression functions and complete odds are correctly specified.

######  Theorem 10 (EIF for missing response under CCMV for $X$).

For the mean statistical functional

$\psi(P)=\int y\ P(dx,dr,dy,da),$

the multiply-robust influence function in Corollary 9 is the efficient
influence function (EIF), and it achieves the semiparametric efficiency bound
in the model $\bigcap_{r}(\mathcal{M}_{R}(r)\cap\mathcal{M}_{O}(r))$.

This theorem further shows that the multiply-robust estimator is also full
semiparametric efficient.

##  Appendix B Asymptotic theory

In this section, we discuss the asymptotic theory for our multistage
estimators. We direct our attention to estimating the population quantity

$\theta=\mathbb{E}[h(X,W,A)]=\sum_{a\in\mathcal{A}}\mathbb{E}[f_{a}(X,W_{a})1(A=a)],$

where is an ODF as described in Definition 1. We will show that the estimation
methods that we have previously outlined lead to consistent, asymptotically
normal estimators under appropriate conditions. This theory will be useful for
constructing confidence intervals for data analysis and demonstrating
consistency. If $h$ is an estimating equation, then the asymptotic normality
of the solution to the finite-sample version of the estimating equation
follows under standard regularity conditions for $Z$-estimators.

###  B.1 Assumptions for estimating $f_{a}$

Fix an observed outcome pattern $a$. These assumptions hold for estimating the
function $f_{a}$.

  1. (A1)

There exists $\tau_{a}^{*}$ in the interior of the parameter space $T_{a}$
such that

$f_{a}(x,w_{a};\tau_{a}^{*})=f_{a}(x,w_{a}).$

We assume that there exists $\widehat{\tau}_{a}$, which is a regular
asymptotically linear estimator such that

$\widehat{\tau}_{a}-\tau_{a}^{*}=\frac{1}{n}\sum_{i=1}^{n}\phi_{r,a}(X_{i},R_{i},W_{i,A_{i}},A_{i})+o_{p}(1/\sqrt{n}),$

and the following convergence in distribution holds

$\widehat{\tau}_{a}-\tau_{a}^{*}\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma_{f_{a}}^{2})$

for some $\sigma_{f_{a}}^{2}>0$.

  2. (A2)

Assume that $f_{a}$ is uniformly bounded for all values of its parameter such
that

$\sup_{x,w_{a},\tau_{a}}|f_{a}(x,w_{a};\tau_{a})|<M$

for some constant $M>0$.

  3. (A3)

Let $\int(f_{a}(x,w_{a};\widehat{\tau}_{a})-f_{a}(x,w_{a};\tau_{a}^{*}))^{2}\
P(dx,dr,dw_{a},da)=o_{p}(1)$.

  4. (A4)

Let $f_{a}(x,r,w_{a};\tau_{a})$ have a nonzero first derivative with respect
to $\tau_{a}$, and assume that

$\int|\nabla_{\tau_{a}}f_{a}(x,w_{a};\tau_{a})|\ P(dx,dr,dw_{a},da)$

is finite in all of its entries.

The assumption (A1) says that the $f_{a}$ functions belong to a parametric
class and can be estimated at $\sqrt{n}$-consistency. We elect to use a
parametric assumption to make the analysis easier. Assumption (A2) is a
bounded assumption, which is reasonable if the input variables have a compact
support. Assumption (A3) describes the $L_{2}(P)$ convergence of the estimated
$f_{a}$ functions. Finally, Assumption (A4) is a technical assumption to
ensure that asymptotic normality can be achieved and allows for application of
the Dominated Convergence Theorem.

###  B.2 The IPW estimator

####  B.2.1 Assumptions

We now introduce the following 4 assumptions for the IPW theorem below. Each
assumption holds for each pair of fixed covariate pattern $r$ and observed-
outcome pattern $a$.

  1. (B1)

For all missing covariate patterns $r$ and observed outcome patterns $a$,
there exist constants $C_{1}$ and $C_{2}$ such that

$0<C_{1}\leq O_{r,a}(x,w_{a};\alpha_{r,a})\leq C_{2}<1.$

  2. (B2)

There exists $\alpha_{r,a}^{*}$ in the interior of the parameter space
$A_{r,a}$ such that

$O_{r,a}(x,w_{a};\alpha_{r,a}^{*})=\frac{P(R=r|x,a,w_{a})}{P(R=1_{d}|x,a,w_{a})}.$

We assume that there exists $\widehat{\alpha}_{r,a}$, which is a regular
asymptotically linear estimator

$\widehat{\alpha}_{r,a}-\alpha_{r,a}^{*}\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma_{\text{odds},r,a}^{2})$

for some $\sigma_{\text{odds},r,a}^{2}>0$.

  3. (B3)

The following class of functions is Donsker:

$\\{O_{r,a}(x,w_{a};\alpha_{r,a}:\alpha_{r,a}\in A_{r,a}\\}.$

  4. (B4)

Assume that $O_{r,a}(x,w_{a};\alpha_{r,a})$ have a nonzero first derivative
with respect to $\alpha_{r,a}$. Also, let

$\int|\nabla_{\alpha_{r,a}}O_{r,a}(x,w_{a};\alpha_{r,a})|\ P(dx,dr,dw_{a},da)$

be finite in all of its entries, and assume that

$\int(O_{r,a}(x,w_{a};\widehat{\alpha}_{r,a})-O_{r,a}(x,w_{a};\alpha_{r,a}^{*}))^{2}\
P(dx,dr,dw_{a},da)=o_{p}(1).$

Assumption (B1) is to ensure that the complete odds are sufficiently bounded
away from 0, which is true if the variables belong to a compact support.
Assumption (B2) says that the complete odds belong to a parametric model, and
there exists an estimator that is $\sqrt{n}$-consistent. This is reasonable
for most common parametric models under consideration. Assumption (B3) is a
Donsker condition, which is also reasonable if the parametric model is
sufficiently smooth. Assumption (B4) ensures that the estimated regression
function $O_{r,a}(x,w_{a};\widehat{\alpha}_{r,a})$ converges in the $L_{2}(P)$
norm to the true regression function $O_{r,a}(x,w_{a};\alpha_{r,a}^{*})$. The
finite first derivative ensures that we can apply the Dominated Convergence
Theorem to simplify the expression in the proof.

####  B.2.2 Asymptotic theory

We continue our goal of estimating $\theta$ by using an inverse probability
weighting approach. First, using Lemma 1, the parameter of interest writes as

$\displaystyle\theta$
$\displaystyle=\sum_{a}\sum_{r}\mathbb{E}[f_{a}(X,W_{a})O_{r,a}(X,W_{a})1(R=1_{d})1(A=a)],$

where $O_{r,a}(X,W_{a})=P(R=r|X,W_{a},A=a)/P(R=1_{d}|X,W_{a},A=a)$ is a
complete odds.

######  Theorem 11 (Inverse probability weighting).

Suppose the assumptions (A1)-(A4) and (B1)-(B4) hold. Then, the estimator

$\displaystyle\widehat{\theta}_{\text{IPW}}$
$\displaystyle=\sum_{r}\sum_{a}\widehat{\theta}_{IPW,r,a}$
$\displaystyle=\sum_{r}\sum_{a}\left(\frac{1}{n}\sum_{i=1}^{n}f_{a}(X_{i},W_{i,a})O_{r,a}(X_{i},W_{i,a})1(R_{i}=1_{d})1(A_{i}=a)\right)$

is consistent, and for some $\sigma^{2}_{\text{IPW}}>0$, we have the following
convergence in distribution

$\sqrt{n}(\widehat{\theta}_{\text{IPW}}-\theta)\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma^{2}_{\text{IPW}}).$

###  B.3 The RA estimator

####  B.3.1 Assumptions

We introduce the following 3 assumptions for the RA estimator.

  1. (C1)

For all covariate patterns $r$ and observed-outcome patterns $a$, there exists
$\eta_{r,a}^{*}$ in the interior of the parameter space $\mathcal{H}_{r,a}$
such that

$P_{ex}(x_{\bar{r}}|x_{r},r,w_{a},a)=P_{ex}(x_{\bar{r}}|x_{r},r,w_{a},a;\eta_{r,a}^{*}).$

Let $\widehat{P}_{ex}(x_{\bar{r}}|x_{r},r,w_{a},a)\equiv
P_{ex}(x_{\bar{r}}|x_{r},r,w_{a},a;\widehat{\eta}_{r,a})$. We assume that
$\widehat{\eta}_{r,a}$ is a regular asymptotically linear (RAL) estimator
satisfying

$\widehat{\eta}_{r,a}-\eta_{r,a}^{*}=\frac{1}{n}\sum_{i=1}^{n}\psi_{r,a}(X_{i},R_{i},W_{i,A_{i}},A_{i})+o_{p}(1/\sqrt{n}).$

Moreover,
$\sqrt{n}(\widehat{\eta}_{r,a}-\eta_{r,a}^{*})\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma_{\text{RA},r,a}^{2})$ for some
$\sigma_{\text{RA},r,a}^{2}>0$.

  2. (C2)

For a fixed covariate pattern $r$ and observed-outcome pattern $a$, define

$\displaystyle m_{r,a}(x_{r},w_{a};\eta_{r,a})$
$\displaystyle=\mathbb{E}[f_{a}(X,W_{a})|X_{r}=x_{r},R=r,W_{a}=w_{a},A=a]$
$\displaystyle=\int f_{a}(x,w_{a})\
P_{ex}(dx_{\bar{r}}|x_{r},r,w_{a},a;\eta_{r,a}).$

The following class of regression functions is Donsker:

$\\{m_{r,a}(x_{r},w_{a};\eta_{r,a}):\eta_{r,a}\in\mathcal{H}_{r,a}\\}.$

  3. (C3)

For all covariate patterns $r$ and observed-outcome patterns $a$, assume that
$m_{r,a}(x_{r},w_{a};\eta_{r,a})$ has a nonzero first derivative with respect
to $\eta_{r,a}$. Also, let

$\int|\nabla_{\eta_{r,a}}m_{r,a}(x_{r},w_{a};\eta_{r,a})|\
P_{obs}(dx_{r},r,dw_{a},a)$

be finite in all of its entries, and assume that

$\int(m_{r,a}(x_{r},w_{a};\widehat{\eta}_{r,a})-m_{r,a}(x_{r},w_{a};\eta_{r,a}^{*}))^{2}\
P_{obs}(dx_{r},r,dw_{a},a)=o_{p}(1).$

These are fairly mild assumptions. Assumption (C1) says that the extrapolation
distribution belongs to a parametric family and that the estimator
$\widehat{\eta}_{r,a}$ admits an asymptotically linear expansion. Provided the
parametric model holds, the $\sqrt{n}$-consistency of $\widehat{\beta}_{r,a}$
is a reasonable assumption since $\eta^{*}_{r,a}$ is often estimated using
$M$-estimators. Assumptions (C2) and (C3) are similar to assumptions (B3) and
(B4) made in the IPW estimation theorem.

###  B.4 Asymptotic theory

We can estimate $\theta$ using a regression adjustment approach. Using Lemma
5, the parameter writes as

$\displaystyle\theta$ $\displaystyle=\mathbb{E}[h(X,W,A)]$
$\displaystyle=\sum_{a}\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$
$\displaystyle=\sum_{a}\sum_{r}\mathbb{E}[m_{r,a}(X_{r},W_{a})1(R=r)1(A=a)]$
$\displaystyle=\sum_{a}\sum_{r}\theta_{\text{RA},r,a}.$

Thus, the main strategy will be to show asymptotic results for every fixed $r$
and $a$ pair.

######  Theorem 12 (Regression adjustment).

Suppose the assumptions (A1)-(A4) and (C1)-(C3) hold. Then, the estimator

$\displaystyle\widehat{\theta}_{\text{RA}}$
$\displaystyle=\sum_{r}\sum_{a}\widehat{\theta}_{\text{RA},r,a}$
$\displaystyle=\sum_{r}\sum_{a}\left(\frac{1}{n}\sum_{i=1}^{n}m_{r,a}(X_{r},W_{a};\widehat{\eta}_{r,a})1(R_{i}=r)1(A_{i}=a)\right)$

is consistent, and for some $\sigma^{2}_{\text{RA}}$, we have the following
convergence in distribution

$\sqrt{n}(\widehat{\theta}_{\text{RA}}-\theta)\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma^{2}_{\text{RA}}).$

Theorem 12 establishes that a population quantity expressible in the form (B)
can be estimated via a regression adjustment approach. The asymptotic
normality of the provided estimator can be used to construct a confidence
interval for inference. Although it is possible to construct a sandwich
estimator, we again recommend using the bootstrap to estimate the variance
because writing a closed form expression may not be so elegant.

##  Appendix C Beyond CCMV

We extend the work in the main body of the paper to see the analysis under an
identifying restriction that is not CCMV. We consider the case with 2
covariate variables $(X_{1},X_{2})$ subject to missingness for a total of
$2^{2}=4$ possible patterns. For simplicity, let $v=(w_{a},a).$ The CCMV
assumption (Definition 6) is equivalent to:

$\displaystyle p(x_{1}|x_{2},R=01,v)=p(x_{1}|x_{2},R=11,v),$ (19)
$\displaystyle p(x_{2}|x_{1},R=10,v)=p(x_{2}|x_{1},R=11,v),$ $\displaystyle
p(x_{1},x_{2}|R=00,v)=p(x_{1},x_{2}|R=11,v).$

The last equality can be relaxed by other assumptions. Here we consider the
following identifying restriction:

$\displaystyle p(x_{1}|x_{2},R=01,v)=p(x_{1}|x_{2},R=11,v),$ (20)
$\displaystyle p(x_{2}|x_{1},R=10,v)=p(x_{2}|x_{1},R=11,v),$ $\displaystyle
p(x_{1},x_{2}|R=00,v)=\kappa_{1}\cdot p(x_{1},x_{2}|R=01,v)+\kappa_{2}\cdot
p(x_{1},x_{2}|R=10,v)+\kappa_{3}\cdot p(x_{1},x_{2}|R=11,v)$

for nonnegative constants $\kappa_{1},\kappa_{2}$, and $\kappa_{3}$ satisfying
$\kappa_{1}+\kappa_{2}+\kappa_{3}=1$. When $\kappa_{3}=1$, equation (20)
reduces to the CCMV in equation (19). If there exist variables that are always
observed, then the above identifying restrictions can be generalized to allow
$\kappa_{1}$, $\kappa_{2}$, and $\kappa_{3}$ to be real-valued functions of
the observed variables. Under suitable choice of
$\kappa_{1},\kappa_{2},\kappa_{3}$ (as a function of observed variables), the
last restriction in equation (20) becomes

$\displaystyle p(x_{1},x_{2}|R=00,v)=p(x_{1},x_{2}|R\in\\{01,10,11\\},v),$

which shows connections to the available-case missing value assumption in
(Molenberghs et al., 1998).

###  C.1 IPW estimators

We now show that the (complete) odds admits a nice identifiable
representation. Let $V=(W_{A},A)$. The two odds

$\frac{P(R=10|X,V)}{P(R=11|X,V)}\quad\text{and}\quad\frac{P(R=01|X,V)}{P(R=11|X,V)}$

are identifiable via the same formulas under CCMV. We now decompose the odds
$P(R=00|X,V)/P(R\in\\{01,10,11\\}|X,V)$ via

$\displaystyle\kappa_{1}\cdot\frac{p(x_{1},x_{2}|R=01,v)P(R=00|v)}{p(x_{1},x_{2},R\in\\{01,10,11\\}|v)}+\kappa_{2}\cdot\frac{p(x_{1},x_{2}|R=10,v)P(R=00|v)}{p(x_{1},x_{2},R\in\\{01,10,11\\}|v)}$
$\displaystyle\quad+\kappa_{3}\cdot\frac{p(x_{1},x_{2}|R=11,v)P(R=00|v)}{p(x_{1},x_{2},R\in\\{01,10,11\\}|v)}.$

The first term can be decomposed by

$\displaystyle\kappa_{1}\cdot\frac{p(x_{1},x_{2}|R=01,v)P(R=00|v)}{p(x_{1},x_{2},R\in\\{01,10,11\\}|v)}$
$\displaystyle=\kappa_{1}\cdot\frac{p(x_{1}|x_{2},R=01,v)p(x_{2}|R=01,v)P(R=00|v)}{p(x_{2}|x_{1},R=11,v)p(x_{1},R=10|v)+p(x_{1}|x_{2},R=11|v)p(x_{2},R=01|v)+p(x_{1},x_{2},R=11|v)}$
$\displaystyle=\kappa_{1}\cdot\frac{p(x_{2}|R=01,v)P(R=00|v)}{\dfrac{p(x_{2}|x_{1},R=11,v)p(x_{1},R=10|v)}{p(x_{1}|x_{2},R=01,v)}+p(x_{2},R=01|v)+p(x_{2},R=11|v)}$
$\displaystyle=\kappa_{1}\cdot\frac{p(x_{2}|R=01,v)P(R=00|v)}{\dfrac{p(x_{1},x_{2},R=11|v)p(x_{1},R=10|v)p(x_{2},R=01|v)}{p(x_{1},x_{2},R=01|v)p(x_{1},R=11|v)}+p(x_{2},R=01|v)+p(x_{2},R=11|v)}$
$\displaystyle=\kappa_{1}\cdot\frac{p(x_{2}|R=01,v)P(R=00|v)}{\dfrac{P(R=11|x_{2},v)}{P(R=01|x_{2},v)}\cdot\dfrac{P(R=10|x_{1},v)}{P(R=11|x_{1},v)}\cdot
p(x_{2},R=01|v)+p(x_{2},R=01|v)+p(x_{2},R=11|v)}$
$\displaystyle=\kappa_{1}\cdot\frac{P(R=00|v)/P(R=01|v)}{\dfrac{P(R=11|x_{2},v)}{P(R=01|x_{2},v)}\cdot\dfrac{P(R=10|x_{1},v)}{P(R=11|x_{1},v)}+\dfrac{P(R=11|x_{2},v)}{P(R=01|x_{2},v)}+1},$

which is identifiable.

A similar argument for the second term yields

$\kappa_{2}\cdot\frac{P(R=00|v)/P(R=10|v)}{\dfrac{P(R=11|x_{1},v)}{P(R=10|x_{1},v)}\cdot\dfrac{P(R=01|x_{2},v)}{P(R=11|x_{2},v)}+\dfrac{P(R=11|x_{1},v)}{P(R=10|x_{1},v)}+1}.$

The third term yields

$\displaystyle\kappa_{3}\cdot\frac{p(x_{1},x_{2}|R=11,v)P(R=00|v)}{p(x_{1},x_{2},R\in\\{01,10,11\\}|v)}$
$\displaystyle=\kappa_{3}\cdot\frac{P(R=00|v)/P(R=11|v)}{\dfrac{P(R=10|x_{1},v)}{P(R=11|x_{1},v)}+\dfrac{P(R=01|x_{2},v)}{P(R=11|x_{2},v)}+1},$

which is identifiable.

Thus, this implies that $P(R=00|x,v)/P(R=11|x,v)$ is identifiable. We can
model the odds $P(R=r|x_{r},v)/P(R=11|x_{r},v)$, and then, plug them into the
formulas above to obtain an identifiable expression for
$P(R=00|x,v)/P(R=11|x,v)$.

###  C.2 RA estimator (imputation)

First, as the identifying restriction is expressed in terms of the
extrapolation distribution, the regression adjustment estimator is fairly
straightforward. Since the extrapolation distribution $p(x_{1},x_{2}|R=00)$ is
a mixture of three identifiable distributions, we can again utilize a Monte
Carlo approximation to construct the regression adjustment estimator. This
approximation is essentially the same as under CCMV except we sample from this
mixture distribution for pattern $R=00$.

##  Appendix D Sensitivity analysis

We present details on the sensitivity analyses utilized in the paper.

###  D.1 Sensitivity analysis for the odds under logistic regression

We have showed with Proposition 3 that under CCMV, the odds satisfies

$O_{r,a}(X,W_{a}):=\frac{P(R=r|X,W_{a},a)}{P(R=1_{d}|X,W_{a},a)}\stackrel{{\scriptstyle\text{CCMV}}}{{=}}\frac{P(R=r|X_{r},W_{a},a)}{P(R=1_{d}|X_{r},W_{a},a)}=:Q_{r,a}(X_{r},W_{a}).$

Thus, the complete odds depend only on the variables observed under pattern
$R=r$. We can introduce a sensitivity parameter to perturb this CCMV
assumption using an exponential tilting method with logistic regression (Kim
and Yu, 2011; Shao and Wang, 2016; Zhao et al., 2017). This results in

$O_{r,a}(X,W_{a})=Q_{r,a}(X_{r},W_{a})e^{\rho_{\bar{r}}X_{\bar{r}}},$ (21)

where $\rho_{\bar{r}}$ is a sensitivity parameter that consists of a vector of
real numbers (one for each missing covariate under pattern $r$). In practice,
we often let every element of $\rho_{\bar{r}}$ to be the same for simplicity.

We now consider the special case of logistic regression. Under CCMV with an
assumed logistic regression model, the complete odds satisfies

$\log O_{r,a}(X,W_{a};\alpha_{r,a})=\log
Q_{r,a}(X_{r},W_{a};\alpha_{r,a})=\alpha_{r,a}^{\top}X_{r}.$

Under the sensitivity model (equation (21)), this becomes

$\log
O_{r,a}(X,W_{a};\alpha_{r,a},\rho_{\bar{r}})=\alpha_{r,a}^{\top}X_{r}+\rho_{\bar{r}}^{\top}X_{\bar{r}}.$

Note that the first term on the RHS is the traditional term obtained by
fitting a logistic regression model, and the second term is the tilting factor
to push the missing data assumption away from CCMV.

Varying the elements of $\rho_{\bar{r}}$ will modify the missing data
assumption and change the weights. This allows one to determine the robustness
of the estimator to changes in this sensitivity parameter. This sensitivity
parameter has a nice interpretation because it is the coefficient in a linear
model, and the range can be specified based on one’s belief about the relation
of impact of the missing variables to the impact of the observed variables.
For example, if the coefficient of the $X_{1}$ variable is $1$ and if $X_{2}$
is subject to missingness, then the sensitivity parameter corresponding to
$X_{2}$ can be chosen to be within the interval $[-1,1]$ by assuming that
$X_{2}$ has at most the same impact in magnitude as $X_{1}$.

###  D.2 Sensitivity analysis for binary imputation via rejection sampling

We also present an approach for sensitivity analysis under imputation model.
We first consider the case where the variable to be imputed $X$ is binary.
Suppose that under an identifying restriction (such as the CCMV), we will
impute the variable $X$ by drawing $X\sim\text{Binomial}(p)$, where $p$ is the
probability parameter from a particular imputation model.

To perturb the imputation model, we consider the following procedure:

  * •

Draw $X\sim\text{Bernoulli}(p)$.

  * •

Consider a sensitivity parameter $\zeta\in[0,1]$. We iterate the following
steps until a draw is accepted:

    * –

If $X=1$, with probability $\zeta$, redraw $X\sim\text{Bernoulli}(p)$.
Otherwise, accept the draw.

    * –

If $X=0$, with probability $1-\zeta$, redraw $X\sim\text{Bernoulli}(p)$.
Otherwise, accept the draw.

We denote this new random variable by $Y$. We can see when $\zeta=1$, we force
$Y$ to be deterministically $0$, and when $\zeta=0$, we force $Y$ to be
deterministically $1$. When $\zeta=1/2$, we resample exactly half of the time,
so we obtain $Y=X$. Namely, when $\zeta=1/2$, we are under the unperturbed
model. Every other value of $\zeta$ interpolates between these extremes. Thus,
the parameter $\zeta$ can be viewed as a preference toward $X=0$.

It turns out that the above procedure admits a simpler expression. Given that
the final output is a binary number, the accepted value $Y$ follows from a
Bernoulli random variable. Since $Y$ is binary, $\mathbb{E}[Y]$ is the
probability of success. On the initial iteration, we draw $1$ and accept it
with probability $p(1-\zeta)$. Otherwise, we redraw with probability
$p\zeta+(1-p)(1-\zeta)$ and return to the same state. This yields the
following equation

$\mathbb{E}[Y]=p(1-\zeta)+p\zeta\mathbb{E}[Y]+(1-p)(1-\zeta)\mathbb{E}[Y],$

so we have
$\mathbb{E}[Y]=\dfrac{p(1-\zeta)}{1-p\zeta-(1-p)(1-\zeta)}=\dfrac{p(1-\zeta)}{p+\zeta-2p\zeta}$.
Thus, we actually do not need to perform the rejection sampling procedure.
Instead, we can simply draw

$Y\sim\text{Bernoulli}\left(\dfrac{p(1-\zeta)}{p+\zeta-2p\zeta}\right).$

###  D.3 Sensitivity analysis for continuous imputation via Gaussian tilting

Now, we discuss the case of sensitivity analysis on imputing a continuous
random variable. Suppose that without perturbations, we impute a continuous
random variable $X$ from a PDF $p(x)$. We propose the following modified
approach:

  * •

Draw $X\sim p(x)$.

  * •

Consider a sensitivity parameter $\xi\in[0,\infty)$. First, we sample $X\sim
p(x)$. We iterate the following steps until a draw is accepted:

    * –

With probability $\exp(-\xi X^{2})$, accept $X\sim p(x)$.

    * –

Otherwise, redraw $X$.

First, note that as $\xi\to\infty$, we almost always resample, and as $\xi\to
0$, we almost never resample. We denote this new random variable by $Y$. The
distribution of $Y$ is given by

$P(Y\leq y)=P(X\leq y\ |\ \text{accept }X)=P(X\leq y\ |\ U\leq\exp(-\xi
X^{2}))$

for a random variable $U\sim\text{Uniform}(0,1)$ independent of $X$. We see
that

$\displaystyle P(Y\leq y)$ $\displaystyle=P(X\leq y\ |\ U\leq\exp(-\xi
X^{2}))$ $\displaystyle\propto P(X\leq y,U\leq\exp(-\xi X^{2}))$
$\displaystyle=P(U\leq\exp(-\xi X^{2})\ |\ X\leq y)\cdot P(X\leq y)$
$\displaystyle=\int_{-\infty}^{y}\frac{P(U\leq\exp(-\xi x^{2})\ |\ x)}{P(X\leq
y)}\cdot p(x)\ dx\cdot P(X\leq y)$
$\displaystyle=\int_{-\infty}^{y}P(U\leq\exp(-\xi x^{2})\ |\ x)\cdot p(x)\ dx$
$\displaystyle=\int_{-\infty}^{y}\exp(-\xi x^{2})p(x)\ dx.$

Therefore, we see that this resampling method is equivalent to multiplying the
original PDF $p(x)$ by $\exp(-\xi x^{2})$ and renormalizing. Since the
multiplication factor $\exp(-\xi x^{2})$ is an exponential, this action can be
considered a form of Gaussian tilting.

$p(x)$ is Gaussian. Suppose that the original imputation distribution is
Gaussian such that $X\sim N(\mu,\sigma^{2})$. A simple completing the square
argument shows that

$\displaystyle-\frac{1}{2\sigma^{2}}(x-\mu)^{2}-\xi
x^{2}=-\left(\frac{1}{2\sigma^{2}}+\xi\right)\left(x-\frac{\mu}{2\sigma^{2}\xi+1}\right)^{2}+\left(\frac{\mu^{2}}{2\sigma^{2}+4\sigma^{4}\xi}-\frac{\mu^{2}}{2\sigma^{2}}\right).$

Therefore, Gaussian tilting by $\exp(-\xi x^{2})$ results in a random variable

$Y\sim
N\left(\frac{\mu}{2\sigma^{2}\xi+1},\frac{\sigma^{2}}{2\sigma^{2}\xi+1}\right).$

As $\xi\to\infty$, the distribution for $Y$ degenerates to a point mass of $1$
at $0$. As $\xi\to 0$, we approach the original distribution $p(x)$, which is
what we expect as we almost never resample.

$p(x)$ is Exponential. Suppose $p(x)$ is Exponential$(\lambda)$ such that
$p(x)=\lambda\exp(-\lambda x)$. Then, Gaussian tilting by $\exp(-\xi x^{2})$
results in a truncated Gaussian that is formed by taking the PDF over the
positive support of a Gaussian with mean $-\lambda/(2\xi)$ and variance
$1/\xi$. As $\xi\to\infty$, again, the distribution for $Y$ degenerates to a
point mass of $1$ at $0$.

##  Appendix E Details on survival analysis simulation and additional
simulation studies

###  E.1 Details on Cox model simulation

Here we provide details on the Cox model simulation. Note that the data
generating process is described in the main paper.

####  E.1.1 Estimation of the imputation model

In our simulated Cox model problem, there are two covariates
$X=(X_{1},X_{2})\in\\{0,1\\}^{2}$. In the complete case, the observed-data
distribution can be decomposed into the following terms:

$\displaystyle p(x,R=11,y,\delta)$
$\displaystyle=\underbrace{p(x,R=11)}_{\text{Empirical}}\cdot\frac{\overbrace{P(R=11|x,y,\delta)}^{\text{Propensity
score}}}{\underbrace{P(R=11|x)}_{\text{
Integration}}}\cdot\underbrace{p(y,\delta|x)}_{\text{EM}}.$

Each of the above four components in the RHS is estimated as follows:

  * •

Empirical. This component can be directly estimated from the data using
empirical proportion.

  * •

Propensity score. This is estimated by combining estimators of the complete
odds and use equation (5).

  * •

Integration. The term $P(R=11|x)$ can be rewritten as follows

$\displaystyle P(R=11|x)$
$\displaystyle=\sum_{\delta=0}^{1}\int_{0}^{\infty}P(R=11|x,y,\delta)p(y,\delta|x)\
dy.$

In practice, this integral does not have an analytic form, so we use numerical
integration in R to circumvent the problem.

  * •

EM. The estimation of this component is more involved. We use an EM-like
method to estimate it. We assume a parametric model
$p_{Y,\Delta|X}(y,\delta|x)=[\nu_{1}\exp(x^{\top}\beta)]^{\delta}\nu_{2}^{1-\delta}e^{-y(\nu_{1}\exp(x^{\top}\beta)+\nu_{2})}$;
this model is derived from our design of the simulation. Fitting
$p(y,\delta|x)$ is a simple maximum likelihood procedure when the covariates
are fully observed. Consider a specific $\widetilde{x}$. The model becomes

$\displaystyle p_{Y,\Delta|X}(y,\delta|\widetilde{x})$
$\displaystyle=[\nu_{1}\exp(\widetilde{x}^{\top}\beta)]^{\delta}\nu_{2}^{1-\delta}e^{-y(\nu_{1}\exp(\widetilde{x}^{\top}\beta)+\nu_{2})}$
$\displaystyle=\gamma_{\widetilde{x}}^{\delta}\nu_{2}^{1-\delta}e^{-y(\gamma_{\widetilde{x}}+\nu_{2})}$

for $\gamma_{x}=\nu_{1}\exp(x^{\top}\beta)$. The ML estimators are

$\widehat{\gamma}_{\widetilde{x}}=\frac{\sum_{i=1}^{n}\delta_{i}1(x_{i}=\widetilde{x})}{\sum_{i=1}^{n}y_{i}1(x_{i}=\widetilde{x})}\quad\text{and}\quad\widehat{\nu}_{2}=\frac{\sum_{i=1}^{n}(1(x_{i}=\widetilde{x})-\delta_{i}1(x_{i}=\widetilde{x}))}{\sum_{i=1}^{n}y_{i}1(x_{i}=\widetilde{x})}.$

In practice, we do not know the covariates for all of the individuals. So, we
fit the models using an EM algorithm. We replace the indicators
$1(x_{i}=\widetilde{x})$ with weights $w_{i}$.

Suppose $\widetilde{x}=(1,1)$. The weights are as follows:

$w_{i}=\begin{cases}1&\text{ if $x_{1}=x_{2}=1$}\\\ 0&\text{ if $x_{1}=0$ or
$x_{2}=0$}\\\ \widehat{P}(X_{2}=1|X_{1}=1,R=11,y,\delta)&\text{ if $x_{1}=1$
and $x_{2}$ missing}\\\ \widehat{P}(X_{1}=1|X_{2}=1,R=11,y,\delta)&\text{ if
$x_{2}=1$ and $x_{1}$ missing}\end{cases}.$

The weights depend on the fitted parameters. For the EM algorithm, we iterate
between updating $w_{i}$ (E-step) and computing the MLEs (M-step).

With the above procedure, we are able to estimate the distribution
$p(x,R=11,y,\delta)$, which describes all imputation models under the CCMV
assumption.

####  E.1.2 Sensitivity analysis

We include the following sensitivity analysis as part of our simulations. We
fix a sample size of $n=1000$ with $M=50$ imputations. For the IPW approach,
we vary the exponential tilting sensitivity parameter $\rho$ in the interval
$[-3,3]$. For the imputation/RA approach, we change the binary sensitivity
analysis parameter $\zeta$ in the interval $[0,1]$. For each corresponding
value of the sensitivity parameter, we compute the IPW and RA estimators,
respectively. The true value of the parameter is indicated with the red line.

Figure 7: Sensitivity analysis for the Cox model simulation with both the IPW
and RA estimators.

First, we observe that the IPW estimator appears to be fairly stable as we
vary the sensitivity parameter, and the size of the confidence interval does
not change drastically. Therefore, we conclude that the $\beta_{1}$ parameter
is significant even as we perturb the estimator. On the other hand, the point
estimates of the RA estimator do appear to change slightly as we vary the
imputation model. The $\beta_{2}$ parameter is large enough that any of these
perturbations do not result in a change in conclusion. However, the
$\beta_{1}$ parameter is relatively close in magnitude to $0$, so for extreme
values of the $\zeta$ parameter, the results are not significant. This
suggests that the choice of identifying restriction is critical, and a
sensitivity analysis is crucial for a sound data analysis. We note that the
extreme values of the $\zeta$ parameter generally imply near almost sure
deterministic imputation for both covariates, which generally will be unlikely
to be true.

###  E.2 Details on binary treatment simulation

####  E.2.1 Data generating process

We generate the data as follows:

  1. 1.

We sample $(R,Y,A)$ as discrete random variables as follows:

     * •

$P(R=11,Y=0,A=0)=1/16$

     * •

$P(R=11,Y=0,A=1)=1/12$

     * •

$P(R=11,Y=1,A=0)=1/16$

     * •

$P(R=11,Y=1,A=1)=1/24$

     * •

$P(R=10,Y=0,A=0)=1/12$

     * •

$P(R=10,Y=0,A=1)=1/8$

     * •

$P(R=10,Y=1,A=0)=1/24$

     * •

$P(R=10,Y=1,A=1)=1/12$

     * •

$P(R=01,Y=0,A=0)=1/16$

     * •

$P(R=01,Y=0,A=1)=1/24$

     * •

$P(R=01,Y=1,A=0)=1/16$

     * •

$P(R=01,Y=1,A=1)=1/12$

     * •

$P(R=00,Y=0,A=0)=1/36$

     * •

$P(R=00,Y=0,A=1)=1/24$

     * •

$P(R=00,Y=1,A=0)=1/24$

     * •

$P(R=00,Y=1,A=1)=1/18$

  2. 2.

We sample $(X_{1},X_{2})|R,Y,A$ via the following distributions:

     * •

$X_{1},X_{2}|R=11,Y=0,A=0\sim N\left(\begin{bmatrix}3\\\
4\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=11,Y=0,A=1\sim N\left(\begin{bmatrix}4\\\
4\end{bmatrix},\begin{bmatrix}0.5&0.2\\\ 0.2&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=11,Y=1,A=0\sim N\left(\begin{bmatrix}3\\\
2\end{bmatrix},\begin{bmatrix}0.4&0.1\\\ 0.1&0.4\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=11,Y=1,A=1\sim N\left(\begin{bmatrix}2\\\
2\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=10,Y=0,A=0\sim N\left(\begin{bmatrix}2\\\
19/5\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=10,Y=0,A=1\sim N\left(\begin{bmatrix}2\\\
16/5\end{bmatrix},\begin{bmatrix}0.5&0.2\\\ 0.2&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=10,Y=1,A=0\sim N\left(\begin{bmatrix}1\\\
3/2\end{bmatrix},\begin{bmatrix}0.4&0.1\\\ 0.1&0.4\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=10,Y=1,A=1\sim N\left(\begin{bmatrix}3\\\
11/5\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=01,Y=0,A=0\sim N\left(\begin{bmatrix}13/5\\\
2\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=01,Y=0,A=1\sim N\left(\begin{bmatrix}14/5\\\
1\end{bmatrix},\begin{bmatrix}0.5&0.2\\\ 0.2&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=01,Y=1,A=0\sim N\left(\begin{bmatrix}3.125\\\
2.5\end{bmatrix},\begin{bmatrix}0.4&0.1\\\ 0.1&0.4\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=01,Y=1,A=1\sim N\left(\begin{bmatrix}1.9\\\
1.5\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=00,Y=0,A=0\sim N\left(\begin{bmatrix}3\\\
4\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=00,Y=0,A=1\sim N\left(\begin{bmatrix}4\\\
4\end{bmatrix},\begin{bmatrix}0.5&0.2\\\ 0.2&0.5\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=00,Y=1,A=0\sim N\left(\begin{bmatrix}3\\\
2\end{bmatrix},\begin{bmatrix}0.4&0.1\\\ 0.1&0.4\end{bmatrix}\right)$

     * •

$X_{1},X_{2}|R=00,Y=1,A=1\sim N\left(\begin{bmatrix}2\\\
2\end{bmatrix},\begin{bmatrix}0.5&0.1\\\ 0.1&0.5\end{bmatrix}\right)$

The generative model for $(X_{1},X_{2})|R,Y,A$ is constructed in such a way to
ensure that the CCMV assumption holds.

####  E.2.2 Details on the estimation

Note that in the binary treatment problem, we have two estimators; one
estimator to account for the missing covariates (denoted as IPW-R or RA-R) and
the other estimator for the ATE (denoted as IPW-A, RA-A, DR-A).

Final estimators. For ease of simulation, we can consider the following 6
estimators for $\mathbb{E}[Y(1)]$. Note that here we describe the linear
functional form of the estimator. We take the empirical average of them to
form the final estimate.

  1. 1.

_IPW-R, IPW-A._

$\frac{1(R=11)1(A=1)Y}{P(A=1|X)}\sum_{r}Q_{r,1}(X_{r},Y)$

  2. 2.

_IPW-R, RA-A._

$\frac{1(R=11)m_{1}(X)}{P(A=1|X)}\sum_{r}\sum_{a}1(A=a)Q_{r,a}(X_{r},Y)$

  3. 3.

_IPW-R, DR-A._

$\displaystyle\frac{1(R=11)1(A=1)(Y-m_{1}(X))}{P(A=1|X)}\sum_{r}Q_{r,1}(X_{r},Y)$
$\displaystyle\quad+\frac{1(R=11)m_{1}(X)}{P(A=1|X)}\sum_{r}\sum_{a}1(A=a)Q_{r,a}(X_{r},Y)$

  4. 4.

_RA-R, IPW-A._

$1(A=1)Y\sum_{r}1(R=r)M_{1/\pi_{A},r,1}(X_{r},Y)$

  5. 5.

_RA-R, RA-A._

$\sum_{r}\sum_{a}1(R=r)1(A=a)M_{m,r,a}(X_{r},Y)$

  6. 6.

_RA-R, DR-A._

$\displaystyle 1(A=1)Y\sum_{r}1(R=r)M_{m,r,1}(X_{r},Y)$
$\displaystyle\quad+\sum_{r}\sum_{a}1(R=r)1(A=a)M_{m,r,a}(X_{r},Y)$
$\displaystyle\quad-1(A=1)\sum_{r}1(R=r)M_{m/\pi_{A},r,1}(X_{r},Y)$

As a reminder, the notation $m_{1}$, $M_{m,r,0}$, $M_{m,r,1}$,
$M_{1/\pi_{A},r,1}$, and $M_{m/\pi_{A},r,1}$ refer to the following regression
functions:

$\displaystyle m_{1}(X)=\mathbb{E}[Y|A=1,X],$ $\displaystyle
M_{m,r,0}(X_{r})=\mathbb{E}[m_{1}(X)|X_{r},R=r,A=0],$ $\displaystyle
M_{m,r,1}(X_{r},Y)=\mathbb{E}[m_{1}(X)|X_{r},R=r,Y,A=1],$ $\displaystyle
M_{1/\pi_{A},r,1}(X_{r},Y)=\mathbb{E}\left[\frac{1}{\pi_{A}(X)}\biggr{|}X_{r},R=r,Y,A=1\right],$
$\displaystyle
M_{m/\pi_{A},r,1}(X_{r},Y)=\mathbb{E}\left[\frac{m_{1}(X)}{\pi_{A}(X)}\biggr{|}X_{r},R=r,Y,A=1\right].$

Similar estimators can be constructed for $\mathbb{E}[Y(0)]$. Merging
estimators for $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$ together, we can
obtain $6$ estimators similar to the ones above for the ATE. We choose not to
construct the multiply-robust estimator for simplicity. As seen in Theorem 8,
there will be numerous terms, so the multiply-robust estimator will generally
be considered a theoretical construct.

Here we describe how to estimate the components in the above estimators.

IPW-R component. The above construction implies that the complete odds has the
following form:

$\displaystyle Q_{r,a}(X_{r},Y)$
$\displaystyle=\frac{P(R=r|X_{r},Y,A)}{P(R=11|X_{r},Y,A)}$
$\displaystyle=\frac{p(x_{r}|r,y,a)p(r,y,a)}{p(x_{r}|11,y,a)p(11,y,a)}.$

In the special case where covariances are assumed to be the same for a fixed
$y$ and $a$, one may simplify the complete odds accordingly

$\displaystyle Q_{10,a}(x_{1},y)$
$\displaystyle=\frac{p(x_{1}|10,y,a)p(r,y,a)}{p(x_{1}|11,y,a)p(11,y,a)}$
$\displaystyle=\frac{\exp(-(x_{1}-\mu_{1,10})^{2}/(2\sigma^{2}))\cdot
p(10,y,a)}{\exp(-(x_{1}-\mu_{1,11})^{2}/(2\sigma^{2}))\cdot p(11,y,a)}$
$\displaystyle=\exp\left(\frac{\mu_{1,10}-\mu_{1,11}}{\sigma^{2}}x_{1}+\frac{\mu_{1,11}^{2}-\mu_{1,10}^{2}}{2\sigma^{2}}\right)\cdot\frac{p(10,y,a)}{p(11,y,a)},$
$\displaystyle Q_{01,a}(x_{2},y)$
$\displaystyle=\frac{p(x_{2}|01,y,a)p(r,y,a)}{p(x_{2}|11,y,a)p(11,y,a)}$
$\displaystyle=\frac{\exp(-(x_{2}-\mu_{2,01})^{2}/(2\sigma^{2}))\cdot
p(01,y,a)}{\exp(-(x_{2}-\mu_{2,11})^{2}/(2\sigma^{2}))\cdot p(11,y,a)}$
$\displaystyle=\exp\left(\frac{\mu_{2,01}-\mu_{2,11}}{\sigma^{2}}x_{1}+\frac{\mu_{2,11}^{2}-\mu_{2,01}^{2}}{2\sigma^{2}}\right)\cdot\frac{p(01,y,a)}{p(11,y,a)},$
$\displaystyle Q_{00,a}(y)$ $\displaystyle=\frac{p(00,y,a)}{p(11,y,a)}.$

The probabilities $p(r,y,a)$ can be estimated by the empirical frequencies,
and the means and variances are estimated using the sample versions. Finally,
the propensity score of $R=11$ is

$\pi_{R}(x,y,a)=P(R=11|x,y,a)=\frac{1}{1+Q_{10,a}(x_{1},y)+Q_{01,a}(x_{1},y)+Q_{00,a}(y)}.$

RA-R component. We estimate the regression functions $M_{m,r,0}$, $M_{m,r,1}$,
$M_{1/\pi_{A},r,1}$, and $M_{m/\pi_{A},r,1}$ using a Monte Carlo approach,
described in Section 3.2. Specifically, in the simulation procedure, we make
the assumption that the distribution $p(x_{1},x_{2}|r,y,a)$ is a multivariate
Gaussian and that the missing mechanism for the covariates satisfies the CCMV
assumption. This implies that the conditional distributions
$p(x_{1}|x_{2},r,y,a)$ and $p(x_{2}|x_{1},r,y,a)$ are conditional Gaussians.
Therefore, we first can estimate the distributions $p(x_{1}|R=10,r,y,a)$,
$p(x_{2}|R=01,r,y,a)$, and $p(x_{1},x_{2}|R=11,r,y,a)$ with a maximum
likelihood procedure directly (plug in the sample means and
variances/covariances as parameters). Then, under the CCMV assumption, we can
construct the estimates of the distribution $p(x_{1},x_{2}|r,y,a)$ in the
following way:

$\displaystyle\widehat{p}(x_{1},x_{2}|R=00,y,a)=\widehat{p}(x_{1},x_{2}|R=11,y,a),$
$\displaystyle\widehat{p}(x_{1},x_{2}|R=01,y,a)=\widehat{p}(x_{1}|x_{2},R=11,y,a)\widehat{p}(x_{2}|R=01,y,a),$
$\displaystyle\widehat{p}(x_{1},x_{2}|R=10,y,a)=\widehat{p}(x_{2}|x_{1},R=11,y,a)\widehat{p}(x_{1}|R=10,y,a).$

This further implies that we can use the estimated distributions

$\displaystyle\widehat{p}(x_{1},x_{2}|R=00,y,a)=\widehat{p}(x_{1},x_{2}|R=11,y,a),$
$\displaystyle\widehat{p}(x_{1}|x_{2},R=01,y,a)=\widehat{p}(x_{1}|x_{2},R=11,y,a),$
$\displaystyle\widehat{p}(x_{2}|x_{1},R=10,y,a)=\widehat{p}(x_{2}|x_{1},R=11,y,a)$

to complete our data. We use the above extrapolation densities as an estimated
imputation model. We then sample the missing covariates $M=50$ times to
generate multiple complete data sets
$\widetilde{D}^{(j)}=\\{(\widetilde{X}_{i}^{(j)},R_{i},Y_{i},A_{i})\\}_{i=1}^{n}$,
where
$\widetilde{X}_{i}^{(j)}=\left(X_{i,R_{i}},\widetilde{X}_{i,\bar{R}_{i}}^{(j)}\right)$
refers to the $j$th imputed missing covariates appended to the observed
covariates corresponding to the $i$th observation.

We pool all $\widetilde{D}^{(j)}$ data sets together to form one large data
set $\check{D}$. This enables us to take the empirical average of the RA-$R$
estimators directly using all of the completed data $\check{D}$. For example,
for the (RA-$R$, IPW-$A$) estimator, we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[1(A_{i}=1)Y_{i}\sum_{r}1(R_{i}=r)\widehat{M}_{1/\pi_{A},r,1}(X_{R_{i},i},Y_{i})\right]$
$\displaystyle\quad=\frac{1}{n}\sum_{i=1}^{n}\left[1(A_{i}=1)Y_{i}\cdot\frac{1}{M}\sum_{j=1}^{M}\frac{1}{\widehat{\pi}_{A}(\widetilde{X}_{i}^{(j)})}\right]$
$\displaystyle\quad=\frac{1}{Mn}\sum_{i=1}^{n}\sum_{j=1}^{M}\left[1(A_{i}=1)Y_{i}\cdot\frac{1}{\widehat{\pi}_{A}(\widetilde{X}_{i}^{(j)})}\right].$

Similar logic implies that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{r}\sum_{a}1(R_{i}=r)1(A_{i}=a)\widehat{M}_{m,r,a}(X_{R_{i},i},Y_{i})\right]$
$\displaystyle\quad=\frac{1}{n}\sum_{i=1}^{n}\left[\sum_{a}1(A_{i}=a)\frac{1}{M}\sum_{j=1}^{M}\widehat{m}_{1}(\widetilde{X}_{i}^{(j)})\right]$
$\displaystyle\quad=\frac{1}{Mn}\sum_{i=1}^{n}\sum_{j=1}^{M}\widehat{m}_{1}(\widetilde{X}_{i}^{(j)}),$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left[1(A_{i}=1)\sum_{r}1(R_{i}=r)\widehat{M}_{m/\pi_{A},r,a}(X_{R_{i},i},Y_{i})\right]$
$\displaystyle\quad=\frac{1}{n}\sum_{i=1}^{n}\left[1(A_{i}=1)\frac{1}{M}\sum_{j=1}^{M}\frac{\widehat{m}_{1}(\widetilde{X}_{i}^{(j)})}{\widehat{\pi}_{A}(\widetilde{X}_{i}^{(j)})}\right]$
$\displaystyle\quad=\frac{1}{Mn}\sum_{i=1}^{n}\sum_{j=1}^{M}1(A_{i}=1)\frac{\widehat{m}_{1}(\widetilde{X}_{i}^{(j)})}{\widehat{\pi}_{A}(\widetilde{X}_{i}^{(j)})}.$

With these equations, we can construct the (RA-$R$, RA-$A$) and (RA-$R$,
DR-$A$) estimators.

IPW-A component. The propensity score can be calculated via

$\displaystyle\pi_{A}(x)=P(A=1|x_{1},x_{2})$
$\displaystyle=\frac{\sum_{r}\sum_{y=0}^{1}p(x_{1},x_{2}|r,y,A=1)p(r,y,A=1)}{\sum_{r}\sum_{y=0}^{1}\sum_{a=0}^{1}p(x_{1},x_{2}|r,y,a)p(r,y,a)}.$

We describe how to estimate the distribution $\widehat{p}(x_{1},x_{2}|r,y,a)$
in the previous discussion on the RA-R component.

RA-A component. Since $Y$ is binary, the regression function can be calculated
via

$\displaystyle m_{a}(x)=\mathbb{E}[Y|A=a,x_{1},x_{2}]$
$\displaystyle=P(Y=1|A=a,x_{1},x_{2})$
$\displaystyle=\frac{\sum_{r}p(x_{1},x_{2}|r,Y=1,A=a)p(r,Y=1,A=a)}{\sum_{r}\sum_{y=0}^{1}p(x_{1},x_{2}|r,y,a)p(r,y,a)}.$

We describe how to estimate the distribution $\widehat{p}(x_{1},x_{2}|r,y,a)$
in the previous discussion on the RA-R component.

####  E.2.3 Sensitivity analysis

We perform a sensitivity analysis for this simulation by varying the
sensitivity parameters $\rho\in[-2,0.5]$ and $\xi\in[0,4]$ based on Appendices
D.1 and D.3. For each sensitivity parameter, we generate 1000 data sets with
$n=1000$ and perturb the 6 estimators. We use $M=50$ imputations for the RA
estimator. In keeping with what was done in the primary analysis, the
confidence intervals for the IPW-$R$ estimators are constructed by pooling the
results across all 1000 iterations and then taking the 0.025% and 97.5%
quantiles. This was done to reduce the effect of outliers (created by unstable
inverse probability weights) on the standard error estimate. The confidence
intervals for the RA-$R$ estimators are generated by calculating the standard
error across all 1000 iterations. The true ATE is plotted in red.

Figure 8: Sensitivity analysis for the Cox model simulation with both the IPW
and RA estimators.

First, we observe that when $\xi$ reaches $4$, we are very close to
deterministic imputation, so often in most practical applications, $\xi$ does
not need to belong to a large interval. The estimator where we use RA in both
stages (RA-$R$, RA-$A$) performs the best and appears fairly robust to changes
in $\xi$. The (RA-$R$, IPW-$A$) estimator also appears to have minimal bias as
we vary $\xi$, but the variance increases. The IPW estimators are fairly
sensitive to changes in $\rho$, and in particular, the variance can severely
increase for moderate changes in $\rho$. This sensitivity analysis suggests
that RA tends to provide the most stable estimates when the assumptions are
perturbed. Intuitively, this makes sense because a perturbation in the missing
data assumption changes the complete odds, and therefore, can drastically
inflate the inverse probability weights in some scenarios.

##  Appendix F Additional analysis of the NACC data

We performed sensitivity analysis under the IPW and RA settings by perturbing
both the IPW and RA estimators with sensitivity parameters $\rho\in[-1,1]$ and
$\zeta\in[0,1]$ (since the covariates in the NACC data are binary). The
details for sensitivity analysis are outlined Appendices D.1 and D.2.

We include the results of all sensitivity analysis of the depression and
anxiety in this section. Figures 9 and 10 show that for three age groups
($60-64$, $65-69$, and $70-74$), the point estimates given by the IPW and RA
estimators are generally fairly stable even as we vary the sensitivity
parameters. Moreover, the confidence intervals for the age groups $60-64$ and
$70-74$ for both the IPW and RA estimators either exhibit significance or have
endpoints close to $0$ across much of the sensitivity parameters’ ranges. This
preliminary analysis suggests that at least for some age groups, recent
depression is associated with shorter time to dementia. There is evidence in
the literature to support this claim as well (Kitching, 2015).

Figure 9: Sensitivity analysis results for the IPW estimator for $\beta_{4}$.

Figure 10: Sensitivity analysis results for the RA estimator for $\beta_{4}$.

Figures 11 and 12 display the results for anxiety. We see across all age
groups, as we vary the respective sensitive parameters $\rho$ and $\zeta$, the
point estimates for $\beta_{5}$ are all close to $0$ with all of the
confidence intervals covering $0$. This suggests some evidence that the self-
reported anxiety is not associated with shorter time to dementia.

Figure 11: Sensitivity analysis results for the IPW estimator for $\beta_{5}$.

Figure 12: Sensitivity analysis results for the RA estimator for $\beta_{5}$.

##  Appendix G Proofs

The proofs of the results in the paper are provided here. We separate the
proofs into three sections by topic: 1. estimator construction and
identification, 2. multiple robustness and efficiency theory, and 3.
asymptotic theory.

###  G.1 Estimator construction and identification

This section contains proofs of Lemma 1, Proposition 2, Proposition 3, Theorem
4, Lemma 5.

Proof of Lemma 1. Let the observed outcome pattern $a$ be fixed. We have

$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$ $\displaystyle=\int
f_{a}(x,w_{a})1(a^{\prime}=a)1(r^{\prime}=r)\
p(x,r^{\prime},w_{a},a^{\prime})\ dxdr^{\prime}dw_{a}da^{\prime}$
$\displaystyle=\int f_{a}(x,w_{a})1(a^{\prime}=a)\ p(x,r,w_{a},a^{\prime})\
dxdw_{a}da^{\prime}$ $\displaystyle=\int f_{a}(x,w_{a})1(a^{\prime}=a)\
\frac{p(x,r,w_{a},a^{\prime})}{p(x,1_{d},w_{a},a^{\prime})}\
p(x,1_{d},w_{a},a^{\prime})\ dxdw_{a}da^{\prime}$ $\displaystyle=\int
f_{a}(x,w_{a})1(a^{\prime}=a)\frac{P(R=r|x,w_{a},a^{\prime})}{P(R=1_{d}|x,w_{a},a^{\prime})}$
$\displaystyle\qquad\times 1(s=1_{d})\ p(x,s,w_{a},a^{\prime})\
dxdsdw_{a}da^{\prime}$
$\displaystyle=\mathbb{E}[f_{a}(X,W_{a})1(A=a)O_{r,a}(X,W_{a})1(R=1_{d})],$

as desired.

Finally, note that summing over all patterns $r$ yields the second
decomposition

$\mathbb{E}[f_{a}(X,W_{a})1(A=a)]=\sum_{r}\mathbb{E}[f_{a}(X,W_{a})1(A=a)O_{r,a}(X,W_{a})1(R=1_{d})].$

$\square$

Proof of Proposition 2. Fix an observed outcome pattern $a$. We start with the
forward direction. Suppose the complete odds $O_{r,a}(X,W_{a})$ is
identifiable for all $r$. Then,

$\sum_{r}O_{r,a}(X,W_{a})=\sum_{r}\frac{P(R=r|X,W_{a},A=a)}{P(R=1_{d}|X,W_{a},A=a)}=\frac{1}{P(R=1_{d}|X,W_{a},A=a)}$

is identifiable. The desired selection probabilities can be recovered because

$P(R=r|X,W_{a},A=a)=P(R=1_{d}|X,W_{a},A=a)\cdot O_{r,a}(X,W_{a})$

for any covariate missingness pattern $r$.

For the converse, an identifiable set of selection probabilities
$P(R=r|X,W_{a},A=a)$ for all $r$ immediately implies that the complete odds
can be formed and is identifiable since the complete odds is the exactly the
ratio of two selection probabilities. $\square$

Proof of Proposition 3. This result is best seen by showing that the CCMV
formulation in (6) is equivalent to the CCMV formulation in (7), so we will
proceed with a series of equivalences. We have

$\displaystyle
p(x_{\bar{r}}|x_{r},r,w_{a},a)=p(x_{\bar{r}}|x_{r},R=1_{d},w_{a},a)$
$\displaystyle\quad\iff\frac{p(x,r,w_{a},a)}{p(x_{r},r,w_{a},a)}=\frac{p(x,1_{d},w_{a},a)}{p(x_{r},1_{d},w_{a},a)}$
$\displaystyle\quad\iff\frac{p(x,r,w_{a},a)}{p(x,1_{d},w_{a},a)}=\frac{p(x_{r},r,w_{a},)}{p(x_{r},1_{d},w_{a},a)}$
$\displaystyle\quad\iff\frac{P(R=r|X,W_{a},a)}{P(R=1_{d}|X,W_{a},a)}=\frac{P(R=r|X_{r},W_{a},a)}{P(R=1_{d}|X_{r},W_{a},a)}$
$\displaystyle\quad\iff O_{r,a}(X,W_{a})=Q_{r,a}(X_{r},W_{a}).$

Thus, the CCMV missing data assumption on the covariates is equivalent through
the pattern mixture model and complete odds formulations. $\square$

Proof of Theorem 4. We wish to show that under the CCMV assumption, the full-
covariate distribution $p(x,r,w_{a},a)$ is nonparametrically identified.
First, we factorize the target distribution using a pattern mixture model
approach

$p(x,r,w_{a},a)=p(x_{\bar{r}}|x_{r},r,w_{a},a)p(x_{r},w_{a},a).$

Note that the distribution $p(x_{r},r,w_{a},a)$ is identifiable from the
observed data, so the only part we need to handle is the extrapolation density
$p(x_{\bar{r}}|x_{r},r,w_{a},a)$. By the CCMV assumption, we have the equality

$p(x_{\bar{r}}|x_{r},r,w_{a},a)\stackrel{{\scriptstyle\text{CCMV}}}{{=}}p(x_{\bar{r}}|x_{r},R=1_{d},w_{a},a).$

Since the distribution $p(x,R=1_{d},w_{a},a)$ is identifiable from the
observed data (consider only the observations with fully-observed covariates
and observed outcome pattern $A=a$), the conditional distribution
$p(x_{\bar{r}}|x_{r},R=1_{d},w_{a},a)$ is identifiable. This means that the
extrapolation density is equated to an identifiable quantity.

The distribution $p(x,r,w_{a},a)$ derived from the CCMV identifying
restriction does not conflict with the observed data because after
marginalizing out the unobserved variable $x_{\bar{r}}$, we obtain
$p(x_{r},w_{a},a)$, so we have nonparametric identification.

$\square$

Proof of Lemma 5. The expectation rewrites as

$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$ $\displaystyle=\int
f_{a}(x,w_{a})1(a^{\prime}=a)1(r^{\prime}=r)\ p(x,r,w_{a},a^{\prime})\
dxdr^{\prime}dw_{a}da^{\prime}$ $\displaystyle=\int f_{a}(x,w_{a})\
p(x,r,w_{a},a)\ dxdw_{a}$ $\displaystyle=\int
f_{a}(x,w_{a})p(x_{\bar{r}}|x_{r},r,w_{a},a)\ dx_{\bar{r}}\
p(x_{r},r,w_{a},a)\ dx_{r}dw_{a}.$

Recall that $m_{r,a}(x_{r},w_{a})=\int f_{a}(x,w_{a})\
p(x_{\bar{r}}|x_{r},r,w_{a},a)\ dx_{\bar{r}}$ is the regression function.
Substituting it into the previous expression yields

$\displaystyle\int m_{r,a}(x_{r},w_{a})\ p(x_{r},r,w_{a},a)\ dx_{r}dw_{a}$
$\displaystyle=\mathbb{E}[m_{r,a}(x_{r},w_{a})1(R=r)1(A=a)].$

Finally, summing over all patterns $r$ yields the second result

$\displaystyle\mathbb{E}[f_{a}(X,W_{a})1(A=a)]$
$\displaystyle=\sum_{r}\mathbb{E}[m_{r,a}(x_{r},w_{a})1(R=r)1(A=a)].$

$\square$

###  G.2 Multiple robustness and efficiency theory

This section contains proofs of Lemma 6, Theorem 7, Theorem 8, Corollary 9,
Theorem 10.

Proof of Lemma 6. Fix the observed outcome pattern $a$, and suppose the CCMV
assumption holds for $X$. We will show that
$\mathbb{E}[g_{r}(X,W_{a})]=\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$ provided
either $O_{r,a}(X,W_{a})$ or $m_{r,a}(X_{r},W_{a})$ is correctly specified. By
Lemmas 1 and 5, the terms (12) and (13) are known to have expected value equal
to $\mathbb{E}[f_{a}(X,W_{a})1(A=a)1(R=r)]$ when the complete odds
$O_{r,a}(X,W_{a})$ and the regression functions $m_{r,a}(X_{r},W_{a})$ are
correctly specified, respectively.

Therefore, we focus our attention on term (14). There are two cases.

Case 1. In the first case, we will assume that the complete odds
$O_{r,a}(X,W_{a})$ is correctly specified in that it equals
$P(R=r|X,W_{a},A=a)/P(R=1_{d}|X,W_{a},A=a)$. Since term (12) has the desired
expectation when the complete odds are correctly specified, it suffices to
show that term (13) minus term (14) has mean 0. The difference of the terms is

$\displaystyle
m_{r,a}(X_{r},W_{a})1(R=r)1(A=a)-m_{r,a}(X_{r},W_{a})1(A=a)1(R=1_{d})O_{r,a}(X,W_{a}).$

Taking the expectation yields

$\displaystyle\int m_{r,a}(x_{r},w_{a})p(x_{r},r,w_{a},a)\ dx_{r}dw_{a}$
$\displaystyle\quad-\int
m_{r,a}(x_{r},w_{a})1(a^{\prime}=a)1(s=1_{d})O_{r,a}(x,w_{a})\
p(x,s,w_{a},a^{\prime})\ dxdsdw_{a}da^{\prime}$ $\displaystyle=\int
m_{r,a}(x_{r},w_{a})p(x_{r},r,w_{a},a)\ dx_{r}dw_{a}$
$\displaystyle\qquad-\int m_{r,a}(x_{r},w_{a})O_{r,a}(x,w_{a})\
p(x,1_{d},w_{a},a)\ dxdw_{a}.$

Since $O_{r,a}(X,W_{a})$ is assumed to properly specified, this further