# Doubly Flexible Estimation under Label Shift

Seong-ho Lee Yanyuan Ma Jiwei Zhao Penn State UniversityPenn State
UniversityUniversity of Wisconsin-Madison; e-mail<EMAIL_ADDRESS>
###### Abstract

In studies ranging from clinical medicine to policy research, complete data
are usually available from a population ${\mathcal{P}}$, but the quantity of
interest is often sought for a related but different population
${\mathcal{Q}}$ which only has partial data. In this paper, we consider the
setting that both outcome $Y$ and covariate ${\bf X}$ are available from
${\mathcal{P}}$ whereas only ${\bf X}$ is available from ${\mathcal{Q}}$,
under the so-called label shift assumption, i.e., the conditional distribution
of ${\bf X}$ given $Y$ remains the same across the two populations. To
estimate the parameter of interest in population ${\mathcal{Q}}$ via
leveraging the information from population ${\mathcal{P}}$, the following
three ingredients are essential: (a) the common conditional distribution of
${\bf X}$ given $Y$, (b) the regression model of $Y$ given ${\bf X}$ in
population ${\mathcal{P}}$, and (c) the density ratio of the outcome $Y$
between the two populations. We propose an estimation procedure that only
needs some standard nonparametric regression technique to approximate the
conditional expectations with respect to (a), while by no means needs an
estimate or model for (b) or (c); i.e., doubly flexible to the possible model
misspecifications of both (b) and (c). This is conceptually different from the
well-known doubly robust estimation in that, double robustness allows at most
one model to be misspecified whereas our proposal here can allow both (b) and
(c) to be misspecified. This is of particular interest in our setting because
estimating (c) is difficult, if not impossible, by virtue of the absence of
the $Y$-data in population ${\mathcal{Q}}$. Furthermore, even though the
estimation of (b) is sometimes off-the-shelf, it can face curse of
dimensionality or computational challenges. We develop the large sample theory
for the proposed estimator, and examine its finite-sample performance through
simulation studies as well as an application to the MIMIC-III database.

Doubly Flexible Estimation under Label Shift

Key Words: Distribution shift, label shift, doubly flexible, semiparametric
statistics, efficient influence function, model misspecification.

##  1 Introduction

In studies ranging from clinical medicine to policy research, there often
exist data and information from a population ${\mathcal{P}}$, while the
quantity of interest is defined on a particular target ${\mathcal{Q}}$,
relevant but different from ${\mathcal{P}}$. For instance, in a clinical trial
setting, physicians may be left interpreting evidence from a randomized
controlled trial consisting of patients who have demographics and
comorbidities that are quite different from those of their own patients
(population ${\mathcal{Q}}$). As another example, to build a predictive model
on pneumonia outbreak for the flu season (population ${\mathcal{Q}}$),
researchers might find a similar model during the non-flu season relevant and
useful. In these scenarios, there is a discrepancy between the distributions
of ${\mathcal{P}}$ and ${\mathcal{Q}}$, termed distribution shift throughout.
Distribution shift can also refer to the fact that the distribution of the
training sample is different from that of the testing sample, in the
evaluation of a learning algorithm.

In all of these situations, it is of vital interest to propose methods that
can appropriately leverage the information from ${\mathcal{P}}$ into the
statistical tasks for ${\mathcal{Q}}$. Our methodology will use the
information from both outcome (output, response, label) $Y$ and covariate
(input, predictor, feature) ${\bf X}$ in population ${\mathcal{P}}$ as well as
covariate ${\bf X}$ in population ${\mathcal{Q}}$. This setting is also named
unsupervised domain adaptation (Quinonero-Candela et al., 2008; Moreno-Torres
et al., 2012; Kouw & Loog, 2019).

Without any assumptions on the nature of shift, it is certainly impossible to
leverage information between two heterogeneous populations. Two major types of
distribution shifts have been defined in the literature. The first is called
covariate shift where the shift happens between the marginal distributions of
${\bf X}$ while the conditional distribution of $Y$ given ${\bf X}$ does not
change; i.e., $p_{{\bf X}}({\bf x})\neq q_{{\bf X}}({\bf x})$ and
$p_{Y\mid{\bf X}}(y,{\bf x})=q_{Y\mid{\bf X}}(y,{\bf x})$. The difference
between ${\mathcal{P}}$ and ${\mathcal{Q}}$ can be summarized as a density
ratio $q_{{\bf X}}({\bf x})/p_{{\bf X}}({\bf x})$, which is, fortunately,
estimable since covariate ${\bf X}$ is available from both populations.
Covariate shift aligns with the causal learning setting (Schölkopf et al.,
2012) where ${\bf X}$ is the cause and $Y$ is the effect. Covariate shift has
attracted a great deal of attention and has been investigated in many
literatures, such as Shimodaira (2000); Huang et al. (2006); Sugiyama et al.
(2008); Gretton et al. (2009); Sugiyama & Kawanabe (2012); Kpotufe & Martinet
(2021) and the references therein.

The second type, which is the focus of this paper, is named label shift,
because it assumes that the shift is induced by the marginal distributions of
$Y$ while the process generating ${\bf X}$ given $Y$ is identical in both
populations. Formally, it assumes

$\displaystyle p_{Y}(y)\neq q_{Y}(y),\mbox{ and }p_{{\bf X}\mid Y}({\bf
x},y)=q_{{\bf X}\mid Y}({\bf x},y)\equiv g({\bf x},y).$

Label shift is also called prior probability shift (Storkey, 2009; Tasche,
2017), target shift (Zhang et al., 2013; Nguyen et al., 2016), or class prior
change (Du Plessis & Sugiyama, 2014; Iyer et al., 2014). Label shift aligns
with the anticausal learning setting in which the outcome/label $Y$ causes the
covariate/feature ${\bf X}$; for example, diseases cause symptoms or objects
cause sensory observations. Consider the situation that one fits a model to
predict whether a patient has pneumonia based on observed symptoms, and that
this model predicts reliably when deployed in the clinic during the non-flu
season. When the flu season starts, there is a sudden surge of pneumonia cases
and the probability of patients developing pneumonia given that they show
symptoms rises, while the mechanism of showing symptoms of pneumonia is rather
stable. Label shift also exists in many computer vision applications, such as
predicting object locations and directions, and human poses; see Martinez et
al. (2017); Yang et al. (2018); Guo et al. (2020).

In the label shift framework, one fundamental problem (Garg et al., 2020) is
determining whether the shift has occurred and estimating the label
distribution $q_{Y}(y)$, or equivalently, assessing the density ratio
$q_{Y}(y)/p_{Y}(y)\equiv\rho(y)$. In contrast to estimating the density ratio
$q_{{\bf X}}({\bf x})/p_{{\bf X}}({\bf x})$ under covariate shift, estimating
$\rho(y)$ is a daunting task due to the absence of the $Y$ observations in
population ${\mathcal{Q}}$. Works in the label shift framework are mainly
limited to the classification problems in the machine learning literature.
Saerens et al. (2002) proposed a simple Expectation-Maximization (EM)
(Dempster et al., 1977) procedure, named maximum likelihood label shift
(MLLS), to estimate $q_{Y}(y)$ assuming access to a classifier that outputs
the true conditional probabilities of the population ${\mathcal{P}}$,
$p_{Y\mid{\bf X}}(y,{\bf x})$. Later on, Chan & Ng (2005) proposed an EM
algorithm that requires the estimation of $g({\bf x},y)$, which is
unfortunately difficult for high-dimensional ${\bf X}$ and moreover, it does
not apply to regression problems. Alternatively, Lipton et al. (2018) and
Azizzadenesheli et al. (2019) proposed moment-matching based estimators, named
black box shift learning (BBSL) and regularized learning under label shift
(RLLS), that make use of the invertible confusion matrix of a classifier
learned from population ${\mathcal{P}}$. The connection and comparison of
these two lines of research, either empirical or theoretical, remain unclear.
To our best knowledge, neither BBSL nor RLLS has been benchmarked against EM.
Alexandari et al. (2020) showed that, in combination with a calibration named
bias-corrected temperature scaling, MLLS outperforms BBSL and RLLS
empirically; whereas MLLS underperforms BBSL when applied naively. Under label
shift, Maity et al. (2020) also studied the minimax rate of convergence for
nonparametric classification.

For continuous $Y$ in regression problems, estimating $q_{Y}(y)$ becomes the
problem of estimating a function instead of a finite number of parameters. Not
surprisingly, its literature is quite scarce. Zhang et al. (2013) proposed a
nonparametric method to estimate the density ratio by kernel mean matching of
distributions. However, this approach does not scale to large data as the
computational cost is quadratic in the sample size. Nguyen et al. (2016)
considered continuous label shift adaptation and studied an importance weight
estimator, but their approach relies on a parametric model for $p_{Y\mid{\bf
X}}(y,{\bf x})$ hence can only be applied in supervised learning.

In this paper, we take a completely distinct approach from all of the existing
literature. Different from the current majority, our methodology is devised to
accommodate both classification and regression. Whether the outcome $Y$ is
discrete or continuous is not essential in our proposal. We directly estimate
a characteristic of the population ${\mathcal{Q}}$. Specifically, we estimate
the parameter ${\bm{\theta}}$ such that $\mathrm{E}_{q}\\{{\bf U}({\bf
X},Y,{\bm{\theta}})\\}={\bf 0}$ where ${\bf U}(\cdot)$ is a user specified
function and $\mathrm{E}_{q}(\cdot)$ stands for the expectation with respect
to $q_{Y}(y)g({\bf x},y)$ or equivalently to $q_{Y\mid{\bf X}}(y,{\bf
x})q_{{\bf X}}({\bf x})$. This is a general framework, including estimating
the mean of $Y$ or the $t$-th quantile of $Y$ as special cases. According to
how the nuisance components are estimated, detailed in the next three
paragraphs, we propose various estimators for ${\bm{\theta}}$, and develop
large sample theory for these estimators to quantify the estimation
uncertainties and to conduct statistical inference.

To estimate ${\bm{\theta}}$, three nuisance components are involved. First and
foremost is the density ratio $\rho(y)$, which is almost infeasible to
estimate based on the observed data due to the lack of $Y$-observations in
population ${\mathcal{Q}}$. Our intention is to bypass the challenging task of
estimating $\rho(y)$. This turns out achievable through careful manipulation
of other components of the influence function. In fact, a unique feature of
our work is that, we do not need to estimate $\rho(y)$ throughout the
estimation procedure. Instead, only a working model, denoted as $\rho^{*}(y)$,
is needed.

The second one is $p_{Y\mid{\bf X}}(y,{\bf x})$, or some dependent quantities
such as $\mathrm{E}_{p}(\cdot\mid{\bf x})$. In contrast to $\rho(y)$,
estimating $\mathrm{E}_{p}(\cdot\mid{\bf x})$ is blessed with the observed
data in population $\cal P$. Indeed, we can use off-the-shelf machine learning
methods or nonparametric regression methods to obtain the corresponding
estimator $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$. Nonetheless, we can
also choose to give up estimating $\mathrm{E}_{p}(\cdot\mid{\bf x})$ even
though we can do it. This means that we can misspecify the conditional
distribution $p_{Y\mid{\bf X}}(y,{\bf x})$ while we also misspecify the
density ratio $\rho(y)$. We call such an estimator $\widehat{\bm{\theta}}$
doubly flexible—the working density ratio model $\rho^{*}(y)$ is flexible, so
is the working conditional distribution model $p_{Y\mid{\bf X}}^{\star}(y,{\bf
x})$. Note that our superscripts here are different: superscript ∗ stands for
the working model of the density ratio whereas ⋆ is for the conditional
distribution model. This doubly flexible property is much more favorable than
the classic “doubly robust” in the literature. The standard double robustness
means that one can misspecify either one of two models but not both, while
here, we can misspecify both models. As an alternative, if one chooses to
estimate $\mathrm{E}_{p}(\cdot\mid{\bf x})$, say,
$\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$, we name the corresponding
estimator $\widetilde{\bm{\theta}}$ singly flexible—only flexible in working
model $\rho^{*}(y)$.

The third nuisance is the conditional density function $g({\bf x},y)$, whose
estimation might be subject to the curse of dimensionality. Fortunately, in
our estimation procedure, $g({\bf x},y)$ only affects quantities of the form
$\mathrm{E}(\cdot\mid y)$, which are one dimensional regression problems hence
can be easily solved via the most basic nonparametric regression procedure
such as the Nadaraya-Watson estimation.

The remaining of the paper is structured as follows. In Section 2, we first
outline our strategy of how to incorporate samples from two heterogeneous
populations. The proposed doubly flexible estimator is presented in Section 3,
and the alternative singly flexible estimator is contained in Section 4. For
easier understanding and improved readability, we present both the methodology
and the theory for a special parameter $\theta=\mathrm{E}_{q}(Y)$ in the main
text, while defer the general results for ${\bm{\theta}}$ such that
$\mathrm{E}_{q}\\{{\bf U}({\bf X},Y,{\bm{\theta}})\\}={\bf 0}$ to the
Supplement. Section 5 contains empirical results for extensive simulation
studies. We present an application to the MIMIC-III database in Section 6. The
paper is concluded with discussions in Section 7. All the technical details
are also included in the Supplement.

##  2 Model Structure

We consider independent and identically distributed (iid) observations
$\\{Y_{i},{\bf X}_{i}\\},i=1,\ldots,n_{1}$ from population ${\mathcal{P}}$,
and iid observations ${\bf X}_{j},j=n_{1}+1,\ldots,n_{1}+n_{0}=n$ from
population ${\mathcal{Q}}$. To use the information in population
${\mathcal{P}}$ under label shift, we stack the two random samples together
and assemble a new data set of size $n$, which represents a random sample for
an imaginary population consisting of $100\pi\%$ population ${\mathcal{P}}$
members and $100(1-\pi)\%$ population ${\mathcal{Q}}$ members. Here we define
$\pi\equiv n_{1}/n$. Throughout our derivation, other than
$\mathrm{E}_{p}(\cdot)$ and $\mathrm{E}_{q}(\cdot)$, we also compute
$\mathrm{E}(\cdot)$ that is with respect to this imaginary population;
however, this imaginary population is only used as an intermediate tool to
leverage information from two heterogeneous populations under label shift. Our
final conclusion will only be made for the target population ${\mathcal{Q}}$.

For convenience, we introduce a binary indicator $R$ in this stacked random
sample, where $R=1$ means the subject is from population ${\mathcal{P}}$ and
$R=0$ population ${\mathcal{Q}}$. Thus, the likelihood of one observation from
the stacked random sample is

$\displaystyle\\{g({\bf x},y)p_{Y}(y)\\}^{r}\left\\{\int g({\bf
x},y)q_{Y}(y)\differential y\right\\}^{1-r}\pi^{r}(1-\pi)^{1-r}$ (1)
$\displaystyle=$ $\displaystyle\\{g({\bf x},y)p_{Y}(y)\\}^{r}\left\\{\int
g({\bf x},y)\rho(y)p_{Y}(y)\differential
y\right\\}^{1-r}\pi^{r}(1-\pi)^{1-r}.$ (2)

Although $g({\bf x},y)$ and $p_{Y}(y)$ can be identified from (1),
unfortunately $q_{Y}(y)$ may not. Below is a simple example illustrating the
possible nonidentifiability of $q_{Y}(y)$.

######  Example 1.

Consider a discrete $Y$ with three supporting values $0,1,2$ and a discrete
$X$ with two supporting values $0,1$. In both populations ${\mathcal{P}}$ and
${\mathcal{Q}}$, $g(x,y)$ is given as

$\displaystyle\hbox{pr}(X=0\mid Y=0)=1/5,\hbox{pr}(X=0\mid Y=1)=1/8,\mbox{ and
}\hbox{pr}(X=0\mid Y=2)=2/3.$

In population ${\mathcal{P}}$, the marginal distribution of $Y$, $p_{Y}(y)$,
is given as

$\displaystyle\hbox{pr}(Y=0)=5/16,\hbox{pr}(Y=1)=1/2,\mbox{ and
}\hbox{pr}(Y=2)=3/16.$

In population ${\mathcal{Q}}$, the marginal distribution of $Y$, $q_{Y}(y)$,
is given as

$\displaystyle\hbox{pr}(Y=0)=\frac{5(25-416t)}{336},\hbox{pr}(Y=1)=\frac{32t-1}{6},\mbox{
and }\hbox{pr}(Y=2)=\frac{89+96t}{112},$

where $t\in(1/32,25/336)$. Clearly this satisfies the label shift assumption,
and the marginal distribution of $X$ in population ${\mathcal{Q}}$ is
identifiable since $\hbox{pr}(X=0)=7/12$ is free of $t$; however, the marginal
distribution of $Y$ in population ${\mathcal{Q}}$, $q_{Y}(y)$, is not
identifiable.

The following result demonstrates that the completeness condition on
$p_{Y\mid{\bf X}}(y,{\bf x})$ would ensure the identifiability. Its proof is
contained in Section S.1 of the Supplement.

######  Lemma 1.

If the conditional pdf/pmf $p_{Y\mid{\bf X}}(y,{\bf x})$ of population
${\mathcal{P}}$ satisfies the completeness condition in the sense that, for
any function $h(Y)$ with finite mean, $\mathrm{E}_{p}\\{h(Y)\mid{\bf
X}\\}=\int h(y)p_{Y\mid{\bf X}}(y,{\bf x})\differential y=0$ implies $h(Y)=0$
almost surely, then all the unknown components in (2), i.e., $g({\bf x},y)$,
$p_{Y}(y)$ and $\rho(y)$, are identifiable. Subsequently, $q_{Y}(y)$ is also
identifiable.

The completeness condition in Lemma 1 is mild and has been widely assumed in
instrumental variables, measurement error models, and econometrics; see, e.g,
Newey & Powell (2003); d’Haultfoeuille (2011); Hu & Shiu (2018). Because the
condition is imposed on population ${\mathcal{P}}$ and we have random
observations $({\bf X}_{i},Y_{i})$’s from population ${\mathcal{P}}$, it can
be examined and verified in empirical studies. One can easily check that many
commonly-used distributions such as exponential families satisfy the
completeness condition. In particular, if the outcome $Y$ is discrete with
finitely many supporting values, Newey & Powell (2003) pointed out that the
completeness condition only means that the covariate ${\bf X}$ has a support
whose cardinality is no smaller than that of $Y$.

In this paper, we focus on estimating a characteristic of population
${\mathcal{Q}}$. For better clarity, we will present the results for
$\theta=\mathrm{E}_{q}(Y)=\int yg({\bf x},y)q_{Y}(y)\differential{\bf
x}\differential y$ in the main text, then generalize the results to
${\bm{\theta}}$ that satisfies $\mathrm{E}_{q}\\{{\bf U}(Y,{\bf
X},{\bm{\theta}})\\}={\bf 0}$ in the Supplement. The main challenge in
estimating $\theta$ is caused by the lack of knowledge and data on $q_{Y}(y)$,
or equivalently, $\rho(y)$. Nevertheless, we will construct an estimator that
bypasses the difficulty of assessing $\rho(y)$. We will show that we only need
a working model of $\rho(y)$, denoted as $\rho^{*}(y)$, that can be flexible.
Furthermore, we find that our procedure can also simultaneously avoid
estimating $p_{Y\mid{\bf X}}(y,{\bf x})$, in that we can insert a possibly
misspecified working model $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$. Thus, our
procedure is flexible with respect to both $\rho(y)$ and $p_{Y\mid{\bf
X}}(y,{\bf x})$—doubly flexible. This is a property different from the classic
“double robustness” which means that one can only misspecify one of two models
but not both. In contrast, here, we can misspecify both.

##  3 Proposed Doubly Flexible Estimation for $\theta=\mathrm{E}_{q}(Y)$

If the density ratio function $\rho(y)$ were known, an intuitive estimator of
$\theta=\mathrm{E}_{q}(Y)$ can be created by noticing the relation
$\theta=\mathrm{E}_{q}(Y)=\mathrm{E}_{p}\\{\rho(Y)Y\\}=\mathrm{E}\\{R\rho(Y)Y\\}/\pi$;
that is,

$\displaystyle\widecheck{\theta}=\frac{1}{n}\sum_{i=1}^{n}\frac{r_{i}}{\pi}\rho(y_{i})y_{i}.$
(3)

We call this estimator shift-dependent since it requires the correct
specification of $\rho(y)$. Clearly, if a working model $\rho^{*}(y)$ is
adopted, the corresponding estimator $\widecheck{\theta}^{*}$ is likely
biased.

###  3.1 General Approach to Estimating $\theta$

The creation of an estimator that is not solely shift-dependent is possible.
To motivate our proposed estimator, we first make some simple observations via
balancing the samples from populations ${\mathcal{P}}$ and ${\mathcal{Q}}$.
Recognizing the relation between $\mathrm{E}_{p}(\cdot)$,
$\mathrm{E}_{q}(\cdot)$ and $\mathrm{E}(\cdot)$, the balancing of $Y$ is

$\displaystyle\mathrm{E}\left\\{\frac{R}{\pi}\rho(Y)Y\right\\}=\mathrm{E}_{p}\\{\rho(Y)Y\\}=\mathrm{E}_{q}(Y)=\mathrm{E}\left(\frac{1-R}{1-\pi}\theta\right).$

Further, replacing the variable $Y$ above by an arbitrary function of ${\bf
X}$, we obtain another balancing function

$\displaystyle\mathrm{E}\left\\{\frac{R}{\pi}\rho(Y)b({\bf
X})\right\\}=\mathrm{E}\left\\{\frac{1-R}{1-\pi}b({\bf X})\right\\}.$

Certainly, we also have

$\displaystyle\mathrm{E}\left(\frac{R}{\pi}c\right)=\mathrm{E}\left(\frac{1-R}{1-\pi}c\right)$

for any constant $c$. Combining the above three, we can obtain a family of
mean zero functions

$\displaystyle\frac{r}{\pi}\\{\rho(y)y-b({\bf
x})\rho(y)+c\\}+\frac{1-r}{1-\pi}\\{b({\bf x})-\theta-c\\}:\forall b({\bf
x}),\forall c.$ (4)

Note that the model in (2) contains three unknown functions $p_{Y}(y)$,
$g({\bf x},y)$ and $\rho(y)$. For this model, in Section S.2 of the
Supplement, we establish that

$\displaystyle{\cal F}\equiv\left[\frac{r}{\pi}\\{\rho(y)y-b({\bf
x})\rho(y)+c\\}+\frac{1-r}{1-\pi}\\{b({\bf x})-\theta-c\\}:\mathrm{E}\\{b({\bf
X})\mid y\\}=y,\forall c\right]$

is the family of all influence functions (Bickel et al., 1993; Tsiatis, 2006)
for estimating $\theta$. According to the definition of the influence
function, ${\cal F}$ is sufficiently comprehensive since it can generate any
regular asymptotically linear estimator of $\theta$. The requirement
$\mathrm{E}\\{b({\bf X})\mid y\\}=y$ in the definition of ${\cal F}$ is
pivotal. Different from the mean zero function in (4), which critically relies
on the correct specification of $\rho(y)$, the element in ${\cal F}$ preserves
its zero mean even if $\rho(y)$ is misspecified as long as an appropriate
$b({\bf x})$ is chosen so that $\mathrm{E}\\{b({\bf X})\mid y\\}=y$. To
further discover a wise choice of such a $b({\bf x})$, we first derive a
special element in ${\cal F}$, the efficient influence function $\phi_{\rm
eff}({\bf x},r,ry)$, that corresponds to the semiparametric efficiency bound
and that provides guidance on constructing flexible estimators for $\theta$.

######  Proposition 1.

The efficient influence function $\phi_{\rm eff}({\bf x},r,ry)$ for $\theta$
is

$\displaystyle\phi_{\rm eff}({\bf x},r,ry)$ $\displaystyle=$
$\displaystyle\frac{r}{\pi}\rho(y)\left[y-\frac{\mathrm{E}_{p}\\{a(Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho(Y)\mid{\bf x}\\}}\right]$
$\displaystyle+\frac{1-r}{1-\pi}\left[\frac{\mathrm{E}_{p}\\{a(Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho(Y)\mid{\bf x}\\}}-\theta\right],$

where $a(y)$ satisfies

$\displaystyle\mathrm{E}\left[\frac{\mathrm{E}_{p}\\{a(Y)\rho(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho(Y)\mid{\bf X}\\}}\mid y\right]=y.$ (5)

The detailed derivation of the efficient influence function in Proposition 1
is provided in Section S.3 of the Supplement. Clearly, the unique $b({\bf x})$
that leads to the efficient influence function is

$\displaystyle b({\bf x})\equiv\frac{\mathrm{E}_{p}\\{a(Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho(Y)\mid{\bf
x}\\}}=\frac{\mathrm{E}_{q}\\{a(Y)\mid{\bf
X}\\}}{\mathrm{E}_{q}\\{\rho(Y)\mid{\bf X}\\}+\pi/(1-\pi)}.$

In principle, if both $\rho(y)$ and $b({\bf x})$ were known, we can estimate
$\theta$ by solving the estimating equation
$\sum_{i=1}^{n}\phi_{\rm eff}({\bf x}_{i},r_{i},r_{i}y_{i})=0$, which leads to

$\displaystyle\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.62889pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.87778pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\theta\hss$\crcr}}}\limits}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho(y_{i})\\{y_{i}-b({\bf
x}_{i})\\}+\frac{1-r_{i}}{1-\pi}b({\bf x}_{i})\right].$ (14)

However, the estimator $\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.62889pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.87778pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\theta\hss$\crcr}}}\limits}$
is impractical because of the following three obstacles. First, as we pointed
out, $\rho(y)$ is almost infeasible to estimate based on the observed data.
Second, $\mathrm{E}_{p}(\cdot\mid{\bf x})$ is unknown and needs to be
estimated. Though various off-the-shelf machine learning or nonparametric
regression methods are available, when the dimension of ${\bf x}$ is high,
their performances are not always satisfactory and their computation can be
expensive. The third obstacle lies in solving $a(y)$ from the integral
equation (5), which requires $g({\bf x},y)$ to evaluate its left hand side.
Estimating conditional density $g({\bf x},y)$ could be even more difficult
than estimating $\mathrm{E}_{p}(\cdot\mid{\bf x})$, due to the curse of
dimensionality.

Our proposed estimator will bypass the challenging task of estimating
$\rho(y)$. Throughout the estimation procedure, only a working model
$\rho^{*}(y)$ is needed, which can be arbitrarily misspecified hence is
flexible. This turns out achievable through careful manipulation of other
components of the efficient influence function. Our proposed estimator can
also avoid estimating $\mathrm{E}_{p}(\cdot\mid{\bf x})$ even though we can do
it if we decide to. This means that we can misspecify the conditional density
model $p_{Y\mid{\bf X}}(y,{\bf x})$, encoded as $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x})$, while we also misspecify the density ratio $\rho(y)$.
We call such an estimation procedure doubly flexible. To overcome the third
obstacle, we recognize that $g({\bf x},y)$ only affects quantities of the form
$\mathrm{E}(\cdot\mid y)$, which are one dimensional regression problems hence
can be easily solved via the most basic nonparametric regression procedure
such as the Nadaraya-Watson estimator.

In a nutshell, a unique feature of our work is the tolerance of both
$\rho^{*}(y)$ and $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$, which can be
simultaneously misspecified. We thus name the procedure doubly flexible.

###  3.2 Proposed Estimator $\widehat{\theta}$: Doubly Flexible in
$\rho^{*}(y)$ and $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$

Interestingly and critically, we discover that, even when both $\rho^{*}(y)$
and $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ are misspecified, the corresponding
estimator following the implementation of
$\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.75555pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.62889pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\theta\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{1.87778pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\theta\hss$\crcr}}}\limits}$
in (14) is still consistent for $\theta$. We summarize this result in
Proposition 2 and give its proof in Section S.4 of the Supplement. Below, we
use superscripts ∗ and ⋆ to indicate that the corresponding quantities are
calculated based on the working models $\rho^{*}(y)$ and $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x})$ respectively.

######  Proposition 2.

Define

$\displaystyle\widehat{\theta}_{t}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{y_{i}-b^{*\star}({\bf
x}_{i})\\}+\frac{1-r_{i}}{1-\pi}b^{*\star}({\bf x}_{i})\right],$

where

$\displaystyle b^{*\star}({\bf
x})\equiv\frac{\mathrm{E}_{p}^{\star}\\{a^{*\star}(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf x}\\}},$

and $a^{*\star}(y)$ is a solution to

$\displaystyle\mathrm{E}\left[\frac{\mathrm{E}_{p}^{\star}\\{a^{*\star}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid
y\right]=y.$ (15)

Then $\widehat{\theta}_{t}$ is a consistent estimator of $\theta$.

In Proposition 2, the subscript t in $\widehat{\theta}_{t}$ indicates the
conditional density $g({\bf x},y)$ in (15) is the truth. In reality, note that
$g({\bf x},y)$ only involves in the evaluation of the conditional expectation
$\mathrm{E}(\cdot\mid y)$ on the left hand side of (15). This is a one
dimensional regression problem and can be easily estimated by many basic
nonparametric regression procedures such as the Nadaraya-Watson estimator.
Specifically, we approximate the integral equation (15) by

$\displaystyle y$ $\displaystyle=$
$\displaystyle\widehat{\mathrm{E}}\left[\frac{\mathrm{E}_{p}^{\star}\\{a^{*\star}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid
y\right]$ (16) $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{\mathrm{E}_{p}^{\star}\\{a^{*\star}(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})}$
$\displaystyle=$ $\displaystyle\int
a^{*\star}(t)\rho^{*}(t)\sum_{i=1}^{n}\frac{p_{Y\mid{\bf X}}^{\star}(t,{\bf
x}_{i})}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})}\differential
t,$

where $K_{h}(\cdot)\equiv K(\cdot/h)/h$, $K(\cdot)$ is a kernel function and
$h$ is a bandwidth, with conditions imposed later in our theoretical
investigation. (16) is a Fredholm integral equation of the first type, which
is ill-posed. Numerical methods to provide stable and reliable solutions have
been well studied in the literature (Hansen, 1992). In our numerical
implementations in Sections 5 and 6, we use Landweber’s iterative method
(Landweber, 1951) that is well-known to produce a convergent solution. We
provide those technical details in Section S.5 of the Supplement.

We summarize the complete estimation procedure in Algorithm 1.

Input: data from population ${\mathcal{P}}$: $(y_{i},{\bf x}_{i},r_{i}=1)$,
$i=1,\ldots,n_{1}$, data from population ${\mathcal{Q}}$: $({\bf
x}_{j},r_{j}=0)$, $j=n_{1}+1,\ldots,n$, and value $\pi=n_{1}/n$.

do

(a) adopt a working model for $\rho(y)$, denoted as $\rho^{*}(y)$;

(b) adopt a working model for $p_{Y\mid{\bf X}}(y,{\bf x})$, denoted as
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ or $p_{Y\mid{\bf X}}^{\star}(y,{\bf
x},\widehat{\bm{\zeta}})$;

(c) compute $w_{i}=[\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}]^{-1}$ for $i=1,\dots,n$;

(d) obtain $\widehat{a}^{*\star}(\cdot)$ by solving the integral equation
(16);

(e) compute $\widehat{b}^{*\star}({\bf
x}_{i})=w_{i}\mathrm{E}_{p}^{\star}\\{\widehat{a}^{*\star}(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}$ for $i=1,\dots,n$;

(f) obtain $\widehat{\theta}$ as
$\displaystyle\widehat{\theta}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{y_{i}-\widehat{b}^{*\star}({\bf
x}_{i})\\}+\frac{1-r_{i}}{1-\pi}\widehat{b}^{*\star}({\bf x}_{i})\right].$
(17)

Output: $\widehat{\theta}$.

Algorithm 1 Proposed Estimator $\widehat{\theta}$: Doubly Flexible in
$\rho^{*}(y)$ and $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$

######  Remark 1.

In step (b) of Algorithm 1, one may adopt a completely specified $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x})$ or a partially specified model $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x},{\bm{\zeta}})$ with an unknown parameter ${\bm{\zeta}}$.
If the latter case, a natural strategy is to estimate ${\bm{\zeta}}$ first
based on the observed samples from $\cal P$ via, say MLE, to obtain
$\widehat{\bm{\zeta}}$, then use $p_{Y\mid{\bf X}}^{\star}(y,{\bf
x},\widehat{\bm{\zeta}})$ to replace the completely fixed $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x})$. In fact, we will show that the action of estimating
${\bm{\zeta}}$ has no consequence in terms of estimating $\theta$. This is an
important discovery, because this means one can always include a reasonably
flexible model $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},{\bm{\zeta}})$ so that it
has a good chance of approximating the true $p_{Y\mid{\bf X}}(y,{\bf x})$. If
$p_{Y\mid{\bf X}}(y,{\bf x})=p_{Y\mid{\bf X}}^{\star}(y,{\bf
x},{\bm{\zeta}}_{0})$ for certain ${\bm{\zeta}}_{0}$, then even though the
additional parameter ${\bm{\zeta}}$ causes extra work, the reward is that
$\theta$ can be estimated as efficiently as if we knew $p_{Y\mid{\bf
X}}(y,{\bf x})$ completely. In all the subsequent steps, we replace
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ by $p_{Y\mid{\bf X}}^{\star}(y,{\bf
x},\widehat{\bm{\zeta}})$ for its generality, bearing in mind that
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ degenerates to
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ when the parameter ${\bm{\zeta}}$
vanishes.

We now study the theoretical properties of $\widehat{\theta}$ defined in (17).
The main technical challenge is quantifying the gap between the solutions for
the integral equations (15) and (16), encoded as $a^{*\star}(y)$ and
$\widehat{a}^{*\star}(y)$ respectively. To facilitate the derivation, we
define the linear operator

$\displaystyle{\cal L}^{*\star}(a)(y)$ $\displaystyle\equiv$ $\displaystyle
p_{Y}(y)\mathrm{E}\left[\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid
y\right]=\int a(t)u^{*\star}(t,y)\differential t,\mbox{ where}$ $\displaystyle
u^{*\star}(t,y)$ $\displaystyle\equiv$ $\displaystyle
p_{Y}(y)\int\frac{\rho^{*}(t)p_{Y\mid{\bf X}}^{\star}(t,{\bf
x},{\bm{\zeta}})}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf x}\\}}g({\bf
x},y)\differential{\bf x}.$

Apparently, $a^{*\star}(y)$ satisfies

$\displaystyle{\cal L}^{*\star}(a^{*\star})(y)=v(y),\mbox{ where }v(y)\equiv
p_{Y}(y)y.$

Similarly, $\widehat{a}^{*\star}(y)$ satisfies $\widehat{\cal
L}^{*\star}(\widehat{a}^{*\star})(y)=\widehat{v}(y)$, where

$\displaystyle\widehat{\cal L}^{*\star}(a)(y)$ $\displaystyle\equiv$
$\displaystyle
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}_{i},\widehat{\bm{\zeta}}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i},\widehat{\bm{\zeta}}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}_{i},\widehat{\bm{\zeta}}\\}},\mbox{ and}$ $\displaystyle\widehat{v}(y)$
$\displaystyle\equiv$ $\displaystyle n_{1}^{-1}\sum_{i=1}^{n}v_{i,h}(y)\equiv
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})y.$

We first establish in Lemma 2 that given regularity conditions (A1)-(A4), the
linear operator ${\cal L}^{*\star}$, as well as its inverse, is well behaved.

  1. (A1)

The working model $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ or $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ satisfies the completeness
condition stated in Lemma 1.

  2. (A2)

$\rho^{*}(y)>\delta$ for all $y$ on the support of $p_{Y}(y)$ where $\delta$
is a positive constant, and $\rho^{*}(y)$ is twice differentiable and its
derivative is bounded.

  3. (A3)

The function $u^{*\star}(t,y)$ is bounded and has bounded derivatives with
respect to $t$ and $y$ on its support. The function $a^{*\star}(y)$ in (15) is
bounded.

  4. (A4)

The support sets of $g({\bf x},y),p_{Y}(y),\rho^{*}(y)$ are compact.

######  Lemma 2.

Let $\|a\|_{\infty}\equiv\sup_{y}|a(y)|$. Under Conditions (A1)-(A4), the
linear operator ${\cal L}^{*\star}:L^{\infty}(R)\to L^{\infty}(R)$ is
invertible. In addition, there exist positive finite constants $c_{1},c_{2}$
such that for all $a(y)\in L^{\infty}(R)$,

  1. (i)

$c_{1}\|a\|_{\infty}\leq\|{\cal L}^{*\star}(a)\|_{\infty}\leq
c_{2}\|a\|_{\infty}$,

  2. (ii)

$\|{\cal L}^{*\star-1}(a)\|_{\infty}\leq c_{1}^{-1}\|a\|_{\infty}$.

The proof of Lemma 2 is in Section S.6 of the Supplement. To analyze the
asymptotic normality of the estimator $\widehat{\theta}$, we add two more
regularity conditions on the kernel function and the bandwidth $h$.

  1. (A5)

The kernel function $K(\cdot)\geq 0$ is symmetric, bounded, and twice
differentiable with bounded first derivative. It has support on $(-1,1)$ and
satisfies $\int_{-1}^{1}K(t)\differential t=1$.

  2. (A6)

The bandwidth $h$ satisfies $n_{1}(\log n_{1})^{-4}h^{2}\to\infty$ and
$n^{2}n_{1}^{-1}h^{4}\to 0$.

Condition (A5) is standard for kernel functions. Note that we only need a one-
dimensional kernel function $K(\cdot)$ in our estimation procedure. Condition
(A6) specifies the requirement of the bandwidth $h$ associated with kernel
function $K(\cdot)$. In general we need both $h^{-1}=o\\{n_{1}^{1/2}(\log
n_{1})^{-2}\\}$ and $h=o(n^{-1/2}n_{1}^{1/4})$. If $\pi$ is further assumed to
be bounded away from zero, the second requirement becomes $h=o(n_{1}^{-1/4})$
and one can simply choose $h=n_{1}^{-1/3}$ to meet both requirements. We are
now ready to present the asymptotic normality of the estimator
$\widehat{\theta}$ below. Its proof is contained in Section S.7 of the
Supplement.

######  Theorem 1.

Assume $\widehat{\bm{\zeta}}$ satisfies
$\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$ and
$\mathrm{E}_{p}^{\star}\\{\|{\bf S}_{\bm{\zeta}}^{\star}(Y,{\bf
x},{\bm{\zeta}})\|_{2}\mid{\bf x}\\}$ is bounded, where ${\bf
S}_{\bm{\zeta}}^{\star}(y,{\bf x},{\bm{\zeta}})\equiv\partial\log p_{Y\mid{\bf
X}}^{\star}(y,{\bf x},{\bm{\zeta}})/\partial{\bm{\zeta}}$. For any choice of
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x},{\bm{\zeta}})$ and $\rho^{*}(y)$, under
Conditions (A1)-(A6),

$\displaystyle\sqrt{n_{1}}(\widehat{\theta}-\theta)\to
N(0,\sigma_{\theta}^{2})$

in distribution as $n_{1}\to\infty$, where

$\displaystyle\sigma_{\theta}^{2}=\hbox{var}\left(\sqrt{\pi}\phi_{\rm
eff}^{*\star}({\bf
X},R,RY)+\frac{R}{\sqrt{\pi}}\left[\frac{\mathrm{E}_{p}^{\star}\\{a^{*\star}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
X}\\}}-Y\right]\\{\rho^{*}(Y)-\rho(Y)\\}\right).$

In Theorem 1, the only requirement on $\widehat{\bm{\zeta}}$ is
$\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$. Thus, the
asymptotic variance of $\sqrt{n_{1}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})$ does
not affect the result in Theorem 1 as long as $\widehat{\bm{\zeta}}$ is
$\sqrt{n_{1}}$-consistent for ${\bm{\zeta}}$. This is easily achievable by
constructing a standard MLE or moment based estimator for ${\bm{\zeta}}$ in
the regression model of $Y$ given ${\bf X}$, $p_{Y\mid{\bf X}}^{\star}(y,{\bf
x},{\bm{\zeta}})$, based on the $n_{1}$ observations from population $\cal P$.

In addition, it is clear that the estimator $\widehat{\theta}$ is
$\sqrt{n_{1}}$-consistent, even if $n$ goes to infinity much faster than
$n_{1}$ does. The intuition is that when we only have $n_{1}$ complete
observations in this problem, although a much larger $n_{0}$ can help us
better understand the label shift mechanism, it cannot improve the convergence
rate of $\widehat{\theta}$.

Finally, Theorem 1 also indicates that, when $\widehat{\bm{\zeta}}$ is a
$\sqrt{n_{1}}$-consistent estimator for ${\bm{\zeta}}$ such that $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x},{\bm{\zeta}})=p_{Y\mid{\bf X}}(y,{\bf x})$, and
$\rho^{*}(y)$ is correctly specified as $\rho(y)$, the corresponding estimator
$\widehat{\theta}$ achieves the semiparametric efficiency bound and is the
efficient estimator. We state this result formally as Corollary 1. Since it is
a special case of Theorem 1, its proof is omitted.

######  Corollary 1.

Assume $\widehat{\bm{\zeta}}$ satisfies
$\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$ and
$\mathrm{E}_{p}^{\star}\\{\|{\bf S}_{\bm{\zeta}}^{\star}(Y,{\bf
x},{\bm{\zeta}})\|_{2}\mid{\bf x}\\}$ is bounded. If $p_{Y\mid{\bf
X}}^{\star}(y,{\bf x},{\bm{\zeta}})=p_{Y\mid{\bf X}}(y,{\bf x})$ and
$\rho^{*}(y)=\rho(y)$, under Conditions (A1)-(A6),

$\displaystyle\sqrt{n_{1}}(\widehat{\theta}_{\rm eff}-\theta)\to
N[0,\hbox{var}\\{\sqrt{\pi}\phi_{\rm eff}({\bf X},R,RY)\\}]$

in distribution as $n_{1}\to\infty$.

##  4 Alternative Estimator $\widetilde{\theta}$: Singly Flexible in
$\rho^{*}(y)$

Because the assessment of $\mathrm{E}_{p}(\cdot\mid{\bf x})$ only relies on
the observed data, instead of adopting an arbitrary known model
$\mathrm{E}_{p}^{\star}(\cdot\mid{\bf x})$ or parametric model
$\mathrm{E}_{p}^{\star}(\cdot\mid{\bf x},{\bm{\zeta}})$, one might be willing
to estimate $\mathrm{E}_{p}(\cdot\mid{\bf x})$ in a model free fashion and
replace $\mathrm{E}_{p}^{\star}(\cdot\mid{\bf x})$ in the estimation procedure
presented in Section 3.2 by a well-behaved estimator
$\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$. Here we consider a general
estimator $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ which has convergence
rate faster than $n_{1}^{-1/4}$. This rate is achievable for many
nonparametric regression or machine learning algorithms (Chernozhukov et al.,
2018), see for example, Chen & White (1999) for a class of neural network
models, Wager & Walther (2015) for a class of regression trees and random
forests, and Bickel et al. (2009); Bühlmann & Van De Geer (2011); Belloni &
Chernozhukov (2011, 2013) for a variety of sparse models. Meanwhile, we still
do not aim to estimate $\rho(y)$ since we do not have the $Y$-data in
population ${\mathcal{Q}}$. We denote the corresponding estimator
$\widetilde{\theta}$ and call it singly flexible because of its flexibility in
using a working model $\rho^{*}(y)$.

The idea behind the estimator $\widetilde{\theta}$ is similar to
$\widehat{\theta}$, therefore we only emphasize the difference from Section
3.2. Similar to (15), we define $a^{*}(y)$ as the solution of

$\displaystyle\mathrm{E}\left[\frac{\mathrm{E}_{p}\\{a^{*}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid y\right]=y.$
(18)

Equivalently, $a^{*}(y)$ satisfies ${\cal L}^{*}(a^{*})(y)=v(y)$, where

$\displaystyle{\cal L}^{*}(a)(y)$ $\displaystyle\equiv$ $\displaystyle
p_{Y}(y)\mathrm{E}\left[\frac{\mathrm{E}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid y\right]=\int
a(t)u^{*}(t,y)\differential t,\mbox{ and}$ $\displaystyle u^{*}(t,y)$
$\displaystyle\equiv$ $\displaystyle p_{Y}(y)\int\frac{\rho^{*}(t)p_{Y\mid{\bf
X}}(t,{\bf x})}{\mathrm{E}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf x}\\}}g({\bf
x},y)\differential{\bf x}.$

Using the estimator $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$, we
approximate the integral equation (18) as

$\displaystyle y$ $\displaystyle=$
$\displaystyle\widehat{\mathrm{E}}\left[\frac{\widehat{\mathrm{E}}_{p}\\{a^{*}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf X}\\}}\mid
y\right]$ (19) $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{\widehat{\mathrm{E}}_{p}\\{a^{*}(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}}{\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})},$

and we write $\widehat{a}^{*}(y)$ as the solution to $\widehat{\cal
L}^{*}(\widehat{a}^{*})(y)=\widehat{v}(y)$, where

$\displaystyle\widehat{\cal L}^{*}(a)(y)\equiv
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\frac{\widehat{\mathrm{E}}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}}{\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}.$

We summarize the algorithm for computing the estimator $\widetilde{\theta}$
below.

Input: data from population ${\mathcal{P}}$: $(y_{i},{\bf x}_{i},r_{i}=1)$,
$i=1,\ldots,n_{1}$, data from population ${\mathcal{Q}}$: $({\bf
x}_{j},r_{j}=0)$, $j=n_{1}+1,\ldots,n$, and value $\pi=n_{1}/n$.

do

(a) adopt a working model for $\rho(y)$, denoted as $\rho^{*}(y)$;

(b) adopt a nonparametric or machine learning algorithm for estimating
$\mathrm{E}_{p}(\cdot\mid{\bf x})$, denoted as
$\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$;

(c) compute $\widehat{w}_{i}=[\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}]^{-1}$ for $i=1,\dots,n$;

(d) obtain $\widehat{a}^{*}(\cdot)$ by solving the integral equation (19);

(e) compute $\widehat{b}^{*}({\bf
x}_{i})=\widehat{w}_{i}\widehat{\mathrm{E}}_{p}\\{\widehat{a}^{*}(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}$ for $i=1,\dots,n$;

(f) obtain $\widetilde{\theta}$ as
$\displaystyle\widetilde{\theta}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{y_{i}-\widehat{b}^{*}({\bf
x}_{i})\\}+\frac{1-r_{i}}{1-\pi}\widehat{b}^{*}({\bf x}_{i})\right].$ (20)

Output: $\widetilde{\theta}$.

Algorithm 2 Alternative Estimator $\widetilde{\theta}$: Single Flexible in
$\rho^{*}(y)$

To develop the asymptotic normality of the estimator $\widetilde{\theta}$,
instead of Condition (A3), we need

  1. (A7)

The function $u^{*}(t,y)$ is bounded and has bounded derivatives with respect
to $t$ and $y$ on its support. The function $a^{*}(y)$ in (18) is bounded.

We present Theorem 2, with its proof contained in Section S.8 of the
Supplement.

######  Theorem 2.

Assume $\widehat{\mathrm{E}}_{p}$ satisfies
$|\widehat{\mathrm{E}}_{p}\\{a(Y)\mid{\bf x}\\}-\mathrm{E}_{p}\\{a(Y)\mid{\bf
x}\\}|=o_{p}(n_{1}^{-1/4})$ for any bounded function $a(y)$. For any choice of
$\rho^{*}(y)$, under Conditions (A2), (A4)-(A7),

$\displaystyle\sqrt{n_{1}}(\widetilde{\theta}-\theta)\to
N(0,\sigma_{\theta}^{2})$

in distribution as $n_{1}\to\infty$, where

$\displaystyle\sigma_{\theta}^{2}=\hbox{var}\left(\sqrt{\pi}\phi_{\rm
eff}^{*}({\bf
X},R,RY)+\frac{R}{\sqrt{\pi}}\left[\frac{\mathrm{E}_{p}\\{a^{*}(Y)\rho^{*}(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)\mid{\bf
X}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf
X}\\}}-Y\right]\\{\rho^{*}(Y)-\rho(Y)\\}\right).$

It is direct from Theorem 2 that when the posited model $\rho^{*}(y)$ is
correctly specified, the estimator $\widetilde{\theta}$ becomes the efficient
estimator for $\theta$. We point out this consequence as Corollary 2 below.

######  Corollary 2.

Assume $\widehat{\mathrm{E}}_{p}$ satisfies
$|\widehat{\mathrm{E}}_{p}\\{a(Y)\mid{\bf x}\\}-\mathrm{E}_{p}\\{a(Y)\mid{\bf
x}\\}|=o_{p}(n_{1}^{-1/4})$ for any bounded function $a(y)$. If
$\rho^{*}(y)=\rho(y)$, under Conditions (A2), (A4)-(A7),

$\displaystyle\sqrt{n_{1}}(\widetilde{\theta}_{\rm eff}-\theta)\to
N[0,\hbox{var}\\{\sqrt{\pi}\phi_{\rm eff}({\bf X},R,RY)\\}]$

in distribution as $n_{1}\to\infty$.

Last but not least, Sections 3 and 4 here only present the results for
estimating $\theta=\mathrm{E}_{q}(Y)$. The whole story can be extended to a
general parameter ${\bm{\theta}}$ such that $\mathrm{E}_{q}\\{{\bf U}({\bf
X},Y,{\bm{\theta}})\\}={\bf 0}$, and the results are stated in Sections S.9
and S.10 of the Supplement. In our numerical studies in Sections 5 and 6
below, we analyze both $\mathrm{E}_{q}(Y)$ and the $t$-th quantile of
population ${\mathcal{Q}}$, defined as

$\displaystyle\tau_{q,t}(Y)=\inf\left[y:\mathrm{E}_{q}\\{I(Y\leq y)\\}\geq
t\right],$

where $0<t<1$. This corresponds to
$\mathrm{E}_{q}\left[\eta_{t}\left\\{Y-\tau_{q,t}(Y)\right\\}\right]=0$ where
$\eta_{t}(r)=t-I(r<0)$.

##  5 Simulation Studies

We conduct simulation studies to assess the finite sample performance of our
proposed methods. We report the results for the mean $\mathrm{E}_{q}(Y)$ and
the median $\tau_{q,0.5}(Y)$ of the outcome $Y$ in population ${\mathcal{Q}}$.

We first generate a binary indicator $R_{i},i=1,\dots,n$ from the Bernoulli
distribution with probability $0.5$, and record $n_{1}=\sum_{i=1}^{n}r_{i}$,
$\pi=n_{1}/n$. Then we generate $Y_{i}$ from $N(0,1)$ if $R_{i}=1$ and from
$N(1,1)$ if $R_{i}=0$. The ${\bf X}\mid Y$ distribution is generated from a
3-dimensional normal with mean $(-0.5,0.5,1)^{\rm T}Y_{i}$ and covariance
${\bf I}$, the identity matrix. This implies, $\mathrm{E}_{q}(Y)=1$,
$\tau_{q,0.5}=1$ and the true density ratio model

$\displaystyle\rho(y)=\exp(-0.5+y).$

One can derive that $p_{Y\mid{\bf X}}(y,{\bf x},{\bm{\zeta}})$ follows normal
with mean $(1,{\bf x}^{\rm T}){\bm{\beta}}$ where
${\bm{\beta}}=(0,-0.2,0.2,0.4)^{\rm T}$ and variance $\sigma^{2}=0.4$. Here we
denote ${\bm{\zeta}}\equiv({\bm{\beta}}^{\rm T},\sigma^{2})^{\rm T}$.

We use the following misspecified working models. We define

$\displaystyle\rho^{*}(y)\equiv c^{*}\exp(-0.7+1.2y),$

where $c^{*}\equiv\pi/\\{n^{-1}\sum_{i=1}^{n}r_{i}\exp(-0.7+1.2y_{i})\\}$ in
order to satisify $\mathrm{E}\\{R\rho^{*}(Y)\\}=\pi$. For the working model
$p_{Y\mid{\bf X}}^{\star}$, we define ${\bf
x}^{\star}\equiv[x_{1},\exp(x_{2}/2),x_{3}/\\{1+\exp(x_{2})\\}+10]^{\rm T}$
and define $p_{Y\mid{\bf X}}^{\star}(y,{\bf x}^{\star},{\bm{\zeta}}^{\star})$
as the normal distribution with mean $(1,{\bf x}^{\star\rm
T}){\bm{\beta}}^{\star}$ where
${\bm{\beta}}^{\star}=(-7.000,-0.223,0.363,0.664)^{\rm T}$ and variance
$\sigma^{\star 2}=0.449$. The parameter
${\bm{\zeta}}^{\star}\equiv({\bm{\beta}}^{\star\rm T},\sigma^{\star 2})^{\rm
T}$ is obtained by minimizing the Kullback-Leibler distance
$D_{kl}(p_{Y\mid{\bf X}}\|p_{Y\mid{\bf X}}^{\star})$.

We implement the following seven estimators:

  * (1)

`shift-dependent`∗: $\widecheck{\theta}$ in (3) with $\rho^{*}(y)$;

  * (2)

`doubly-flexible`∗⋆: $\widehat{\theta}$ in (17) with $\rho^{*}(y)$ and
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x},{\bm{\zeta}}^{\star})$, theoretically
analyzed in Theorem 1;

  * (3)

`singly-flexible`∗: $\widetilde{\theta}$ in (20) with $\rho^{*}(y)$,
theoretically analyzed in Theorem 2;

  * (4)

`shift-dependent`0: $\widecheck{\theta}$ in (3) with correct $\rho(y)$;

  * (5)

`doubly-flexible`0: $\widehat{\theta}_{\rm eff}$ with correct $\rho(y)$ and
$p_{Y\mid{\bf X}}(y,{\bf x},{\bm{\zeta}})$, theoretically analyzed in
Corollary 1;

  * (6)

`singly-flexible`0: $\widetilde{\theta}_{\rm eff}$ with correct $\rho(y)$,
theoretically analyzed in Corollary 2;

  * (7)

`oracle`: the $\sqrt{n_{0}}$-consistent estimator
$\frac{1}{n}\sum_{i=1}^{n}\frac{1-r_{i}}{1-\pi}y_{i}$.

Note that the last four estimators (shown as “gray” in Figures 1 and 2) are
unrealistic since they either use the unknown models $\rho(y)$ and
$p_{Y\mid{\bf X}}(y,{\bf x},{\bm{\zeta}})$ or the $Y$-data in population
${\mathcal{Q}}$.

In implementing estimators `doubly-flexible`∗⋆, `singly-flexible`∗, `doubly-
flexible`0 and `singly-flexible`0, we solve the integral equations (16) and
(19) using the Nadaraya-Watson estimator for $\mathrm{E}(\cdot\mid y)$ with
Gaussian kernel and bandwidth $h=n_{1}^{-1/3}$ that is discussed in Condition
(A6). Numerically, the integrations are approximated by the Gauss-Legendre
quadrature with 50 points on the interval $[-5,5]$ and the integral equations
are evaluated at $y_{i},i=1,\dots,n_{1}$. In addition, for estimators `singly-
flexible`∗ and `singly-flexible`0, we estimate $\mathrm{E}_{p}(\cdot\mid{\bf
x})$ using the Nadaraya-Watson estimator based on the product Gaussian kernel
with bandwidth $2.5n_{1}^{-1/7}$, where the order comes from the optimal
bandwidth $n_{1}^{-1/(4+d)}$ with $d$ the dimensionality of covariate ${\bf
X}$. See Section S.5 of the Supplement for technical details on the numerical
implementation.

Based on 1000 simulation replicates, Figures 1 and 2 illustrate the boxplots
of the estimates for the mean and the median, respectively. With the
misspecified working model $\rho^{*}(y)$, the estimator `shift-dependent`∗ is
biased; in contrast, the proposed estimators `doubly-flexible`∗⋆ and `singly-
flexible`∗ are both unbiased. When the correct model $\rho(y)$ is used, not
surprisingly, all of the estimators `shift-dependent`0, `doubly-flexible`0 and
`singly-flexible`0 are unbiased. It is also clear that the two proposed
flexible estimators are always more efficient than the shift-dependent
estimator no matter the correct model $\rho(y)$ is used or not.

To further demonstrate the efficiency comparison and the inference results, in
Tables 1 and 2, we report the mean squared error (MSE), the empirical bias
(Bias), the empirical standard error (SE), the average of estimated standard
error ($\widehat{\rm SE}$), and the empirical coverage at $95\%$ confidence
level (CI), for each of the estimators. The estimator `shift-dependent`∗ has
an incorrect coverage (over-coverage in Table 1 and under-coverage in Table 2)
because of its severe bias. This issue is not mitigated at all in Table 1 or
becomes even worse in Table 2 when we increase the size of the stacked random
sample from 500 to 1000. On the contrary, the estimators `doubly-flexible`∗⋆
and `singly-flexible`∗ are correctly covered. Though there is no theoretical
justification, `doubly-flexible`∗⋆ is slightly less efficient than `singly-
flexible`∗ in this setting. This indicates that the effort of correctly
estimating $\mathrm{E}_{p}(\cdot\mid{\bf x})$ pays off in the sense of
improving the estimation efficiency. With the correct $\rho(y)$ model used,
each of the estimators `doubly-flexible`0 and `singly-flexible`0 is more
efficient than its counterpart `doubly-flexible`∗⋆ and `singly-flexible`∗. The
$\sqrt{n_{0}}$-consistent estimator `oracle` is the most efficient one in this
simulation setting.

Figure 1: Boxplots of estimates for mean in the simulation study. Dashed line:
the true estimand.

(a) mean, $n=500$

(b) mean, $n=1000$

Table 1: Summary of mean estimation results in the simulation study. $n$ | Estimator | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x})$ | MSE | Bias | SE | $\widehat{\rm SE}$ | CI
---|---|---|---|---|---|---|---|---
500 | shift-dependent∗ | $\rho^{*}(y)$ | - | .0699 | .1840 | .1899 | .2791 | 1.000
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0173 | -.0120 | .1311 | .1287 | .943
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0162 | -.0201 | .1258 | .1230 | .941
shift-dependent0 | $\rho(y)$ | - | .0533 | .0049 | .2309 | .2119 | .899
doubly-flexible0 | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0153 | -.0231 | .1214 | .1212 | .941
singly-flexible0 | $\rho(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0138 | -.0221 | .1155 | .1169 | .939
oracle | - | - | .0040 | .0006 | .0636 | .0633 | .943
1000 | shift-dependent∗ | $\rho^{*}(y)$ | - | .0561 | .1906 | .1406 | .2077 | .999
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0085 | .0013 | .0922 | .0912 | .955
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0081 | -.0031 | .0899 | .0850 | .952
shift-dependent0 | $\rho(y)$ | - | .0275 | .0024 | .1660 | .1533 | .927
doubly-flexible0 | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0075 | -.0125 | .0856 | .0861 | .958
singly-flexible0 | $\rho(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0069 | -.0094 | .0824 | .0827 | .955
oracle | - | - | .0020 | .0008 | .0451 | .0447 | .945
Figure 2: Boxplots of estimates for median in the simulation study. Dashed
line: the true estimand.

(a) median, $n=500$

(b) median, $n=1000$

Table 2: Summary of median estimation results in the simulation study. $n$ | Estimator | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x})$ | MSE | Bias | SE | $\widehat{\rm SE}$ | CI
---|---|---|---|---|---|---|---|---
500 | shift-dependent∗ | $\rho^{*}(y)$ | - | .0695 | .2025 | .1687 | .1547 | .802
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0368 | -.0072 | .1918 | .1764 | .940
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0211 | .0024 | .1453 | .1390 | .941
shift-dependent0 | $\rho(y)$ | - | .0178 | .0011 | .1336 | .1286 | .947
doubly-flexible0 | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0093 | -.0033 | .0964 | .0950 | .951
singly-flexible0 | $\rho(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0071 | -.0162 | .0827 | .0881 | .956
oracle | - | - | .0064 | -.0020 | .0799 | .0821 | .946
1000 | shift-dependent∗ | $\rho^{*}(y)$ | - | .0554 | .2018 | .1210 | .1104 | .560
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0229 | -.0013 | .1515 | .1444 | .943
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0128 | -.0023 | .1130 | .1124 | .954
shift-dependent0 | $\rho(y)$ | - | .0091 | .0005 | .0955 | .0915 | .934
doubly-flexible0 | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x},\widehat{\bm{\zeta}})$ | .0047 | .0000 | .0688 | .0669 | .954
singly-flexible0 | $\rho(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | .0036 | -.0098 | .0594 | .0625 | .948
oracle | - | - | .0031 | .0004 | .0558 | .0578 | .959

##  6 Data Application

We now illustrate the numerical performance of our proposed method through
analyzing the Medical Information Mart for Intensive Care III (MIMIC-III), an
openly available electronic health records database, developed by the MIT Lab
for Computational Physiology (Johnson et al., 2016). It comprises deidentified
health-related records including demographics, vital signs, laboratory test
and medications, for 46,520 patients who admitted to the intensive care unit
of the Beth Israel Deaconess Medical Center between 2001 and 2012.

The outcome of interest $Y$ in our analysis is the sequential organ failure
assessment (SOFA) score (Singer et al., 2016), used to track a patient’s
status during the stay in an intensive care unit to determine the extent of a
patient’s organ function or rate of failure. The score is based on six
different sub-scores, with one of each for the respiratory, cardiovascular,
hepatic, coagulation, renal and neurological systems. The SOFA score ranges
from 0 (best) to 24 (worst). We include 16 covariates from either chart events
(6 variables, diastolic blood pressure, systolic blood pressure, blood
glucose, respiratory rate per minute, and two measures from body temperature)
or laboratory tests (10 variables, peripheral caillary oxygen saturation, two
measures from each of hematocrit level, platelets count and red blood cell
count, and three measures from blood urea nitrogen). We choose these
covariates through assessing whether the absolute correlation with the outcome
$Y$ is greater than 0.2 and whether the missing rate is less than 1%. In our
analysis, we only include the first admission if the patient was admitted to
the intensive care unit more than once. We also exclude patients whose outcome
$Y$ is greater than or equal to 20, whose age is greater than or equal to 65,
and who has missing values in any of the covariates. This results in a total
of $n=16,691$ records.

In our analysis, we define the population ${\mathcal{P}}$ as patients with
private, goverment, and self-pay insurances ($R=1,n_{1}=11,695$), and
population ${\mathcal{Q}}$ as patients whose insurance type is either Medicaid
or Medicare ($R=0,n_{0}=4,996$). The label shift assumption that the
conditional distribution of ${\bf X}$ given $Y$ remains the same can be tested
via the conditional independence of ${\bf X}$ and $R$ given $Y$. In our
analysis, we test the conditional independence between $R$ and each of
covariates by the invariant environment prediction test (Heinze-Deml et al.,
2018) in R package `CondIndTests`, and the p-values range from 0.460 to 0.628.
This indicates the label shift assumption is indeed sensible in our analysis.
We first compute the sample mean (3.7409) and sample $t$-th quantiles
(1,1,3,5,8 for $t=(10,25,50,75,90)\%$) of SOFA scores among patients whose
insurance type is either Medicaid or Medicare. We regard these estimates as
`oracle` in order to compare with our proposed methods.

To identify a reasonable working model $\rho^{*}(y)$, we model the data
$Y+0.001$ from population ${\mathcal{P}}$ as a parametric gamma distribution
$f(y,\alpha,\beta)=\Gamma(\alpha)^{-1}\beta^{-\alpha}y^{\alpha-1}\exp(-y/\beta)$
with $\Gamma(\cdot)$ the $\Gamma$-function, the shape parameter $\alpha>0$ and
the scale parameter $\beta>0$. We estimate the unknown parameters $\alpha$ and
$\beta$ as $\widehat{\alpha}$ and $\widehat{\beta}$ using the MLE. For the
$Y$-data in population ${\mathcal{Q}}$, we assume $Y+0.001$ follows a similar
gamma distribution with shape parameter $\widehat{\alpha}+1$ and scale
parameter $\widehat{\beta}$. Hence, we use the working model

$\displaystyle\rho^{*}(y)=\frac{f(y+0.001,\widehat{\alpha}+1,\widehat{\beta})}{f(y+0.001,\widehat{\alpha},\widehat{\beta})}.$

To implement the estimator `doubly-flexible`∗⋆, we impose a parametric model
$p_{Y\mid{\bf X}}^{\star}(y,{\bf x},{\bm{\zeta}})$ by regressing $Y+0.001$ on
${\bf X}$ as a generalized linear model with gamma distribution, and estimate
${\bm{\zeta}}$ using the MLE. To implement the estimator `singly-flexible`∗,
we identify the first three principal components from the 16 covariates, and
then estimate $\mathrm{E}_{p}(\cdot\mid{\bf x})$ as a function of those three
principal components using the Nadaraya-Watson estimator based on the product
Gaussian kernel with bandwidth $0.5n_{1}^{-1/7}$. To solve the corresponding
integral equations, similar to Section 5, we approximate $\mathrm{E}(\cdot\mid
y)$ by its Nadaraya-Watson estimator with the Gaussian kernel and bandwidth
$h=n_{1}^{-1/3}$. In addition, for numerical implementation, the integration
is approximated at 50 equally-spaced points on the interval $[0,19]$, and the
integral equations are evaluated at each supporting point of
$\\{y_{i}:i=1,\dots,n_{1}\\}$. See Section S.5 of the Supplement for technical
details on the implementation.

The results are summarized in Table 3 for the mean and in Table 4 for the
quantiles. The estimator `shift-dependent`∗ that relies on a misspecified
model $\rho^{*}(y)$ severely over-estimates the quantities compared to the
`oracle` estimate, in most of the scenarios including estimating the mean,
25%, 75% and 90% quantiles. As a consequence, the `oracle` estimate cannot be
covered by the confidence intervals. In contrast, the proposed estimators
`doubly-flexible`∗⋆ and `singly-flexible`∗, although also rely on
$\rho^{*}(y)$, provide almost identical estimates as `oracle`. Accordingly,
the confidence intervals from the proposed methods all cover the `oracle`
estimate. In these scenarios, `singly-flexible`∗ is more efficient than
`doubly-flexible`∗⋆, which echoes our findings in Section 5.

When estimating 10% and 50% quantiles, we find that `shift-dependent`∗ gives
almost the same estimate as `oracle`. It might be plausible that the
difference between $\rho^{*}(y)$ and the true $\rho(y)$ is minor for
estimating these two quantities. Nevertheless, the proposed estimators
`doubly-flexible`∗⋆ and `singly-flexible`∗ are still more efficient than the
estimator `shift-dependent`∗.

Finally it is interesting to observe that, the estimator `oracle` is even less
efficient than the proposed estimators `doubly-flexible`∗⋆ and `singly-
flexible`∗, in estimating all of the quantiles. This is because `oracle` is
$\sqrt{n_{0}}$-consistent whereas the two proposed estimators are both
$\sqrt{n_{1}}$-consistent. In this application, $n_{1}$ is 2.34 times greater
than $n_{0}$, which might result in the situation that the `oracle` estimate
being less efficient. In Section 5, we also considered situations that $n_{1}$
is much larger than $n_{0}$ and similar phenomenon was observed as well, with
detailed results omitted.

##  7 Conclusion

In this paper, we estimate a characteristic of a target population
${\mathcal{Q}}$, via exploiting the data and information from a different but
relevant population ${\mathcal{P}}$, under the label shift assumption.
Different from most existing literatures, our proposal is devised to
accommodate both classification and regression problems. We primarily propose
the doubly flexible estimate, whose unique feature is to simultaneously allow
both models to be misspecified thus is flexible: the density ratio model
$\rho(y)$ that governs the label shift mechanism, and the conditional
distribution model $p_{Y\mid{\bf X}}(y,{\bf x})$ of population
${\mathcal{P}}$. While the estimation of the latter can be done via off-the-
shelf procedures sometimes, it often faces curse of dimensionality or
computational challenges. Further, estimating $\rho(y)$ is even more difficult
because the $Y$-data in population ${\mathcal{Q}}$ is not accessible in our
procedure.

Table 3: Mean estimation results in the data application. Estimator | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x})$ | Estimate | diff. with oracle | $\widehat{\rm SE}$ | CI
---|---|---|---|---|---|---
shift-dependent∗ | $\rho^{*}(y)$ | - | 4.0529 | .3120 | .0593 | [3.9367, 4.1691]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 3.7579 | .0170 | .0803 | [3.6005, 3.9153]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 3.7542 | .0133 | .0496 | [3.6570, 3.8514]
oracle | - | - | 3.7409 | - | .0405 | [3.6616, 3.8202]
Table 4: Quantile estimation results in the data application. $\tau$ | Estimator | $\rho(y)$ | $p_{Y\mid{\bf X}}(y,{\bf x})$ | Estimate | diff. with oracle | $\widehat{\rm SE}$ | CI
---|---|---|---|---|---|---|---
10% |  shift-dependent∗ | $\rho^{*}(y)$ | - | 1.0000 | .0000 | .0102 | [0.9801, 1.0199]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 0.9998 | -.0002 | .0095 | [0.9812, 1.0185]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 0.9998 | -.0002 | .0099 | [0.9803, 1.0192]
oracle | - | - | 1.0000 | - | .0218 | [0.9572, 1.0428]
25% |  shift-dependent∗ | $\rho^{*}(y)$ | - | 1.9995 | .9995 | .0249 | [1.9508, 2.0483]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 1.0002 | .0002 | .0276 | [0.9460, 1.0543]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 1.0004 | .0004 | .0196 | [0.9620, 1.0389]
oracle | - | - | 1.0000 | - | .0315 | [0.9383, 1.0617]
50% |  shift-dependent∗ | $\rho^{*}(y)$ | - | 3.0002 | .0002 | .0408 | [2.9203, 3.0801]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 3.0001 | .0001 | .0333 | [2.9348, 3.0654]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 3.0005 | .0005 | .0334 | [2.9350, 3.0659]
oracle | - | - | 3.0000 | - | .0550 | [2.8923, 3.1077]
75% |  shift-dependent∗ | $\rho^{*}(y)$ | - | 5.9999 | .9999 | .0742 | [5.8544, 6.1454]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 5.0001 | .0001 | .0759 | [4.8512, 5.1489]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 5.0001 | .0001 | .0381 | [4.9254, 5.0748]
oracle | - | - | 5.0000 | - | .0591 | [4.8842, 5.1158]
90% |  shift-dependent∗ | $\rho^{*}(y)$ | - | 8.9999 | .9999 | .1380 | [8.7295, 9.2703]
doubly-flexible∗⋆ | $\rho^{*}(y)$ | $p_{Y\mid{\bf X}}^{\star}(y,{\bf x},\widehat{\bm{\zeta}})$ | 8.0004 | .0004 | .1003 | [7.8039, 8.1969]
singly-flexible∗ | $\rho^{*}(y)$ | $\widehat{\mathrm{E}}_{p}(\cdot\mid{\bf x})$ | 8.0004 | .0004 | .0638 | [7.8754, 8.1253]
oracle | - | - | 8.0000 | - | .1093 | [7.7858, 8.2142]

## References

  * (1)
  * Alexandari et al. (2020) Alexandari, A. M., Kundaje, A. & Shrikumar, A. (2020), Maximum likelihood with bias-corrected calibration is hard-to-beat at label shift adaptation, in ‘Proceedings of the 37th International Conference on Machine Learning’.
  * Azizzadenesheli et al. (2019) Azizzadenesheli, K., Liu, A., Yang, F. & Anandkumar, A. (2019), ‘Regularized learning for domain adaptation under label shifts’, arXiv preprint arXiv:1903.09734 .
  * Belloni & Chernozhukov (2011) Belloni, A. & Chernozhukov, V. (2011), ‘L1-penalized quantile regression in high-dimensional sparse models’, The Annals of Statistics 39(1), 82–130.
  * Belloni & Chernozhukov (2013) Belloni, A. & Chernozhukov, V. (2013), ‘Least squares after model selection in high-dimensional sparse models’, Bernoulli 19(2), 521–547.
  * Bickel et al. (1993) Bickel, P. J., Klaassen, J., Ritov, Y. & Wellner, J. A. (1993), Efficient and Adaptive Estimation for Semiparametric Models, Johns Hopkins University Press Baltimore.
  * Bickel et al. (2009) Bickel, P. J., Ritov, Y. & Tsybakov, A. B. (2009), ‘Simultaneous analysis of lasso and dantzig selector’, The Annals of Statistics 37(4), 1705–1732.
  * Bühlmann & Van De Geer (2011) Bühlmann, P. & Van De Geer, S. (2011), Statistics for high-dimensional data: methods, theory and applications, Springer Science & Business Media.
  * Chan & Ng (2005) Chan, Y. S. & Ng, H. T. (2005), Word sense disambiguation with distribution estimation., in ‘IJCAI’, Vol. 5, Citeseer, pp. 1010–5.
  * Chen & White (1999) Chen, X. & White, H. (1999), ‘Improved rates and asymptotic normality for nonparametric neural network estimators’, IEEE Transactions on Information Theory 45(2), 682–691.
  * Chernozhukov et al. (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. & Robins, J. (2018), ‘Double/debiased machine learning for treatment and structural parameters: Double/debiased machine learning’, The Econometrics Journal 21(1).
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M. & Rubin, D. B. (1977), ‘Maximum likelihood from incomplete data via the em algorithm’, Journal of the Royal Statistical Society: Series B (Methodological) 39(1), 1–22.
  * d’Haultfoeuille (2011) d’Haultfoeuille, X. (2011), ‘On the completeness condition in nonparametric instrumental problems’, Econometric Theory 27(3), 460–471.
  * Du Plessis & Sugiyama (2014) Du Plessis, M. C. & Sugiyama, M. (2014), ‘Semi-supervised learning of class balance under class-prior change by distribution matching’, Neural Networks 50, 110–119.
  * Garg et al. (2020) Garg, S., Wu, Y., Balakrishnan, S. & Lipton, Z. (2020), ‘A unified view of label shift estimation’, Advances in Neural Information Processing Systems 33, 3290–3300.
  * Gretton et al. (2009) Gretton, A., Smola, A., Huang, J., Schmittfull, M., Borgwardt, K. & Schölkopf, B. (2009), ‘Covariate shift by kernel mean matching’, Dataset Shift in Machine Learning 3(4), 5.
  * Guo et al. (2020) Guo, J., Gong, M., Liu, T., Zhang, K. & Tao, D. (2020), Ltf: A label transformation framework for correcting label shift, in ‘International Conference on Machine Learning’, PMLR, pp. 3843–3853.
  * Hansen (1992) Hansen, P. C. (1992), ‘Numerical tools for analysis and solution of fredholm integral equations of the first kind’, Inverse problems 8(6), 849.
  * Heinze-Deml et al. (2018) Heinze-Deml, C., Peters, J. & Meinshausen, N. (2018), ‘Invariant causal prediction for nonlinear models’, Journal of Causal Inference 6(2).
  * Hu & Shiu (2018) Hu, Y. & Shiu, J.-L. (2018), ‘Nonparametric identification using instrumental variables: sufficient conditions for completeness’, Econometric Theory 34(3), 659–693.
  * Huang et al. (2006) Huang, J., Gretton, A., Borgwardt, K., Schölkopf, B. & Smola, A. (2006), ‘Correcting sample selection bias by unlabeled data’, Advances in neural information processing systems 19.
  * Iyer et al. (2014) Iyer, A., Nath, S. & Sarawagi, S. (2014), Maximum mean discrepancy for class ratio estimation: Convergence bounds and kernel selection, in ‘International Conference on Machine Learning’, PMLR, pp. 530–538.
  * Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Lehman, L.-w. H., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Anthony Celi, L. & Mark, R. G. (2016), ‘Mimic-iii, a freely accessible critical care database’, Scientific Data 3(1), 1–9.
  * Kouw & Loog (2019) Kouw, W. M. & Loog, M. (2019), ‘A review of domain adaptation without target labels’, IEEE Transactions on Pattern Analysis and Machine Intelligence 43(3), 766–785.
  * Kpotufe & Martinet (2021) Kpotufe, S. & Martinet, G. (2021), ‘Marginal singularity and the benefits of labels in covariate-shift’, The Annals of Statistics 49(6), 3299–3323.
  * Landweber (1951) Landweber, L. (1951), ‘An iteration formula for fredholm integral equations of the first kind’, American Journal of Mathematics 73(3), 615–624.
  * Lipton et al. (2018) Lipton, Z., Wang, Y.-X. & Smola, A. (2018), Detecting and correcting for label shift with black box predictors, in ‘International Conference on Machine Learning’, PMLR, pp. 3122–3130.
  * Maity et al. (2020) Maity, S., Sun, Y. & Banerjee, M. (2020), ‘Minimax optimal approaches to the label shift problem’, arXiv preprint arXiv:2003.10443 .
  * Martinez et al. (2017) Martinez, J., Hossain, R., Romero, J. & Little, J. J. (2017), A simple yet effective baseline for 3d human pose estimation, in ‘Proceedings of the IEEE International Conference on Computer Vision’, pp. 2640–2649.
  * Moreno-Torres et al. (2012) Moreno-Torres, J. G., Raeder, T., Alaiz-Rodríguez, R., Chawla, N. V. & Herrera, F. (2012), ‘A unifying view on dataset shift in classification’, Pattern Recognition 45(1), 521–530.
  * Newey & McFadden (1994) Newey, W. K. & McFadden, D. (1994), ‘Large sample estimation and hypothesis testing’, Handbook of econometrics 4, 2111–2245.
  * Newey & Powell (2003) Newey, W. K. & Powell, J. L. (2003), ‘Instrumental variable estimation of nonparametric models’, Econometrica 71(5), 1565–1578.
  * Nguyen et al. (2016) Nguyen, T. D., Christoffel, M. & Sugiyama, M. (2016), Continuous target shift adaptation in supervised learning, in ‘Asian Conference on Machine Learning’, PMLR, pp. 285–300.
  * Quinonero-Candela et al. (2008) Quinonero-Candela, J., Sugiyama, M., Schwaighofer, A. & Lawrence, N. D. (2008), Dataset Shift in Machine Learning, MIT Press.
  * Saerens et al. (2002) Saerens, M., Latinne, P. & Decaestecker, C. (2002), ‘Adjusting the outputs of a classifier to new a priori probabilities: a simple procedure’, Neural Computation 14(1), 21–41.
  * Schölkopf et al. (2012) Schölkopf, B., Janzing, D., Peters, J., Sgouritsa, E., Zhang, K. & Mooij, J. (2012), On causal and anticausal learning, in ‘29th International Conference on Machine Learning (ICML 2012)’, Omnipress, pp. 1255–1262.
  * Shimodaira (2000) Shimodaira, H. (2000), ‘Improving predictive inference under covariate shift by weighting the log-likelihood function’, Journal of Statistical Planning and Inference 90(2), 227–244.
  * Singer et al. (2016) Singer, M., Deutschman, C. S., Seymour, C. W., Shankar-Hari, M., Annane, D., Bauer, M., Bellomo, R., Bernard, G. R., Chiche, J.-D., Coopersmith, C. M. et al. (2016), ‘The third international consensus definitions for sepsis and septic shock (sepsis-3)’, JAMA 315(8), 801–810.
  * Storkey (2009) Storkey, A. (2009), ‘When training and test sets are different: characterizing learning transfer’, Dataset shift in machine learning 30, 3–28.
  * Sugiyama & Kawanabe (2012) Sugiyama, M. & Kawanabe, M. (2012), Machine learning in non-stationary environments: Introduction to covariate shift adaptation, MIT press.
  * Sugiyama et al. (2008) Sugiyama, M., Suzuki, T., Nakajima, S., Kashima, H., von Bünau, P. & Kawanabe, M. (2008), ‘Direct importance estimation for covariate shift adaptation’, Annals of the Institute of Statistical Mathematics 60(4), 699–746.
  * Tasche (2017) Tasche, D. (2017), ‘Fisher consistency for prior probability shift’, The Journal of Machine Learning Research 18(1), 3338–3369.
  * Tsiatis (2006) Tsiatis, A. A. (2006), ‘Semiparametric theory and missing data’.
  * Wager & Walther (2015) Wager, S. & Walther, G. (2015), ‘Adaptive concentration of regression trees, with application to random forests’, arXiv preprint arXiv:1503.06388 .
  * Yang et al. (2018) Yang, X., Sun, H., Sun, X., Yan, M., Guo, Z. & Fu, K. (2018), ‘Position detection and direction prediction for arbitrary-oriented ships via multitask rotation region convolutional neural network’, IEEE Access 6, 50839–50849.
  * Zhang et al. (2013) Zhang, K., Schölkopf, B., Muandet, K. & Wang, Z. (2013), Domain adaptation under target and conditional shift, in ‘International Conference on Machine Learning’, PMLR, pp. 819–827.

## SUPPLEMENTARY MATERIAL

###  S.1 Proof of Lemma 1

Suppose we have two different sets of models: $g({\bf x},y)$, $p_{Y}(y)$,
$\rho(y)$ and $\widetilde{g}({\bf x},y)$, $\widetilde{p}_{Y}(y)$,
$\widetilde{\rho}(y)$. From the likelihood (1), it is obvious that

$\displaystyle g({\bf x},y)p_{Y}(y)$ $\displaystyle=$
$\displaystyle\widetilde{g}({\bf x},y)\widetilde{p}_{Y}(y),$
$\displaystyle\int g({\bf x},y)p_{Y}(y)\rho(y)\differential y$
$\displaystyle=$ $\displaystyle\int\widetilde{g}({\bf
x},y)\widetilde{p}_{Y}(y)\widetilde{\rho}(y)\differential y.$

Integrate with respect to ${\bf x}$ on both sides of the first equation above,
it is apparent that $p_{Y}(y)=\widetilde{p}_{Y}(y)$ and $g({\bf
x},y)=\widetilde{g}({\bf x},y)$. Therefore, the second equation above becomes

$\displaystyle\int g({\bf x},y)p_{Y}(y)\rho(y)\differential y=\int g({\bf
x},y)p_{Y}(y)\widetilde{\rho}(y)\differential y$ $\displaystyle=$
$\displaystyle\int p_{Y\mid{\bf X}}(y\mid{\bf x})p_{{\bf X}}({\bf
x})\rho(y)\differential y=\int p_{Y\mid{\bf X}}(y\mid{\bf x})p_{{\bf X}}({\bf
x})\widetilde{\rho}(y)\differential y.$

Hence, $\int p_{Y\mid{\bf X}}(y\mid{\bf x})\rho(y)\differential y=\int
p_{Y\mid{\bf X}}(y\mid{\bf x})\widetilde{\rho}(y)\differential y$. Using the
completeness condition, clearly $\rho(y)=\widetilde{\rho}(y)$. ∎

###  S.2 Derivation of ${\cal F}$

We first state and prove a result regarding the tangant space.

######  Proposition S.1.

The tangent space of (1) is ${\cal T}\equiv{\cal T}_{\bm{\alpha}}\oplus({\cal
T}_{\bm{\beta}}+{\cal T}_{\bm{\gamma}})$, where

$\displaystyle{\cal T}_{\bm{\alpha}}$ $\displaystyle=$
$\displaystyle\left[r{\bf a}_{1}(y):\mathrm{E}_{p}\\{{\bf a}_{1}(Y)\\}={\bf
0}\right],$ $\displaystyle{\cal T}_{\bm{\beta}}$ $\displaystyle=$
$\displaystyle\left[r{\bf a}_{2}({\bf x},y)+(1-r)\mathrm{E}_{q}\\{{\bf
a}_{2}({\bf x},Y)\mid{\bf x}\\}:\mathrm{E}\\{{\bf a}_{2}({\bf X},y)\mid
y\\}={\bf 0}\right],$ $\displaystyle{\cal T}_{\bm{\gamma}}$ $\displaystyle=$
$\displaystyle\left[(1-r)\mathrm{E}_{q}\\{{\bf a}_{3}(Y)\mid{\bf
x}\\}:\mathrm{E}_{q}\\{{\bf a}_{3}(Y)\\}={\bf 0}\right].$

Proof. Consider a parametric submodel of (1),

$\displaystyle f_{{\bf X},R,RY}({\bf x},r,ry,{\bm{\delta}})$ $\displaystyle=$
$\displaystyle\pi^{r}(1-\pi)^{1-r}\\{g({\bf
x},y,{\bm{\beta}})p_{Y}(y,{\bm{\alpha}})\\}^{r}\left\\{\int g({\bf
x},y,{\bm{\beta}})q_{Y}(y,{\bm{\gamma}})dy\right\\}^{1-r},$ (S.1)

where ${\bm{\delta}}=({\bm{\alpha}}^{\rm T},{\bm{\beta}}^{\rm
T},{\bm{\gamma}}^{\rm T})^{\rm T}$. We can derive that the score function
associated with an arbitrary ${\bm{\delta}}$ is ${\bf
S}_{\bm{\delta}}\equiv({\bf S}_{\bm{\alpha}}^{\rm T},{\bf S}_{\bm{\beta}}^{\rm
T},{\bf S}_{\bm{\gamma}}^{\rm T})^{\rm T}$, where

$\displaystyle{\bf S}_{\bm{\alpha}}({\bf x},r,ry)$ $\displaystyle\equiv$
$\displaystyle r{\bf a}_{\bm{\alpha}}(y),$ $\displaystyle{\bf
S}_{\bm{\beta}}({\bf x},r,ry)$ $\displaystyle\equiv$ $\displaystyle r{\bf
a}_{\bm{\beta}}({\bf x},y)+(1-r)\mathrm{E}_{q}\\{{\bf a}_{\bm{\beta}}({\bf
x},Y)\mid{\bf x}\\},$ $\displaystyle{\bf S}_{\bm{\gamma}}({\bf x},r,ry)$
$\displaystyle\equiv$ $\displaystyle(1-r)\mathrm{E}_{q}\\{{\bf
a}_{\bm{\gamma}}(Y)\mid{\bf x}\\},$

$\mathrm{E}_{p}\\{{\bf a}_{\bm{\alpha}}(Y)\\}={\bf 0}$, $\mathrm{E}_{p}\\{{\bf
a}_{\bm{\beta}}({\bf X},y)|y\\}=\mathrm{E}_{q}\\{{\bf a}_{\bm{\beta}}({\bf
X},y)|y\\}=\mathrm{E}\\{{\bf a}_{\bm{\beta}}({\bf X},y)|y\\}={\bf 0}$, and
$\mathrm{E}_{q}\\{{\bf a}_{\bm{\gamma}}(Y)\\}={\bf 0}$. The above derivation
directly leads to Proposition S.1.∎

We are now ready to establish the result regarding ${\cal F}$, which we
explicitly write out as Proposition S.2.

######  Proposition S.2.

The set of the influence functions for $\theta$ is

$\displaystyle{\cal F}\equiv\left[\frac{r}{\pi}\\{y\rho(y)-b({\bf
x})\rho(y)+c\\}+\frac{1-r}{1-\pi}\\{b({\bf x})-\theta-c\\}:\mathrm{E}\\{b({\bf
X})\mid y\\}=y,\forall c\right].$ (S.2)

Proof. Note that $\theta=\mathrm{E}_{q}(Y)$, hence

$\displaystyle\frac{\partial\theta}{\partial{\bm{\alpha}}^{\rm T}}$
$\displaystyle=$ $\displaystyle{\bf 0}^{\rm T},$
$\displaystyle\frac{\partial\theta}{\partial{\bm{\beta}}^{\rm T}}$
$\displaystyle=$ $\displaystyle{\bf 0}^{\rm T},$
$\displaystyle\frac{\partial\theta}{\partial{\bm{\gamma}}^{\rm T}}$
$\displaystyle=$ $\displaystyle\mathrm{E}_{q}\\{Y{\bf a}_{\bm{\gamma}}^{\rm
T}(Y)\\}.$

Let $\phi({\bf x},r,ry)$ be

$\displaystyle\phi({\bf x},r,ry)\equiv\frac{r}{\pi}\phi_{1}({\bf
x},y)+\frac{1-r}{1-\pi}\phi_{2}({\bf x}).$

For $\phi({\bf x},r,ry)$ to be an influence function, it must satisfy

$\displaystyle\mathrm{E}(\phi)=\mathrm{E}_{p}\\{\phi_{1}({\bf
X},Y)\\}+\mathrm{E}_{q}\\{\phi_{2}({\bf X})\\}=0$ (S.3)

and $\mathrm{E}(\phi{\bf S}_{\bm{\delta}}^{\rm
T})=\partial\theta/\partial{\bm{\delta}}^{\rm T}$. $\mathrm{E}(\phi{\bf
S}_{\bm{\alpha}}^{\rm T})=\partial\theta/\partial{\bm{\alpha}}^{\rm T}={\bf
0}^{\rm T}$ is equivalent to

$\displaystyle\mathrm{E}\\{\phi_{1}({\bf X},y)\mid y\\}=c$ (S.4)

for some constant $c$. In addition, since

$\displaystyle\mathrm{E}(\phi{\bf S}_{\bm{\beta}}^{\rm T})$ $\displaystyle=$
$\displaystyle\mathrm{E}_{p}\\{\phi_{1}({\bf X},Y){\bf a}_{\bm{\beta}}^{\rm
T}({\bf X},Y)\\}+\mathrm{E}_{q}\\{\phi_{2}({\bf X}){\bf a}_{\bm{\beta}}^{\rm
T}({\bf X},Y)\\}$ $\displaystyle=$
$\displaystyle\mathrm{E}_{p}[\\{\phi_{1}({\bf X},Y)+\phi_{2}({\bf
X})\rho(Y)\\}{\bf a}_{\bm{\beta}}^{\rm T}({\bf X},Y)]$

and ${\bf a}_{\bm{\beta}}({\bf x},y)$ is an arbitrary function which satisfies
$\mathrm{E}\\{{\bf a}_{\bm{\beta}}({\bf X},y)|y\\}={\bf 0}$,
$\mathrm{E}(\phi{\bf S}_{\bm{\beta}}^{\rm
T})=\partial\theta/\partial{\bm{\beta}}^{\rm T}={\bf 0}^{\rm T}$ implies
$\phi_{1}({\bf x},y)+\phi_{2}({\bf x})\rho(y)=a(y)$ for some function $a(y)$.
Then (S.4) yields

$\displaystyle\phi_{1}({\bf x},y)=[\mathrm{E}\\{\phi_{2}({\bf X})\mid
y\\}-\phi_{2}({\bf x})]\rho(y)+c.$

Also, noting that $\mathrm{E}(\phi{\bf S}_{\bm{\gamma}}^{\rm
T})=\partial\theta/\partial{\bm{\gamma}}^{\rm T}=\mathrm{E}_{q}\\{Y{\bf
a}_{\bm{\gamma}}^{\rm T}(Y)\\}$ is equivalent to $\mathrm{E}\\{\phi_{2}({\bf
X})\mid y\\}=y+c^{*}$ for some constant $c^{*}$, we have from (S.3) and (S.4)
that

$\displaystyle\mathrm{E}\\{\phi_{2}({\bf X})\mid y\\}=y-\mathrm{E}_{q}(Y)-c.$

Therefore, defining $b({\bf x})\equiv\phi_{2}({\bf x})+\mathrm{E}_{q}(Y)+c$,
the summary description of the influence function is

$\displaystyle\phi({\bf x},r,ry)=\frac{r}{\pi}\\{y\rho(y)-b({\bf
x})\rho(y)+c\\}+\frac{1-r}{1-\pi}\\{b({\bf x})-\mathrm{E}_{q}(Y)-c\\},$

where $b({\bf x})$ satisfies $\mathrm{E}\\{b({\bf X})\mid y\\}=y$ and $c$ is a
constant. Hence, we get the result (S.2). ∎

###  S.3 Proof of Proposition 1

Note that

$\displaystyle\frac{\mathrm{E}_{p}\\{a({\bf x},Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho(Y)\mid{\bf x}\\}}=\frac{\int a({\bf
x},y)\rho(y)g({\bf x},y)p_{Y}(y)dy}{\int\rho(y)g({\bf
x},y)p_{Y}(y)dy}=\frac{\int a({\bf x},y)g({\bf x},y)q_{Y}(y)dy}{\int g({\bf
x},y)q_{Y}(y)dy}=\mathrm{E}_{q}\\{a({\bf x},Y)\mid{\bf x}\\},$

then $\phi_{\rm eff}({\bf x},r,ry)$ can be alternatively written as

$\displaystyle\phi_{\rm eff}({\bf x},r,ry)$ $\displaystyle=$
$\displaystyle\frac{r}{\pi}\rho(y)\left[y-\frac{\mathrm{E}_{q}\\{a(Y)\mid{\bf
x}\\}}{\mathrm{E}_{q}\\{\rho(Y)\mid{\bf
x}\\}+\pi/(1-\pi)}\right]+\frac{1-r}{1-\pi}\left[\frac{\mathrm{E}_{q}\\{a(Y)\mid{\bf
x}\\}}{\mathrm{E}_{q}\\{\rho(Y)\mid{\bf x}\\}+\pi/(1-\pi)}-\theta\right],$

where $a(y)$ satisfies

$\displaystyle\mathrm{E}\left[\frac{\mathrm{E}_{q}\\{a(Y)\mid{\bf
X}\\}}{\mathrm{E}_{q}\\{\rho(Y)\mid{\bf X}\\}+\pi/(1-\pi)}\mid y\right]=y.$

First, it is immediate that $\phi_{\rm eff}({\bf x},r,ry)$ is an influence
function for $\theta$, i.e., belongs to $\cal F$ given in (S.2) from letting

$\displaystyle b({\bf x})$ $\displaystyle\equiv$
$\displaystyle\frac{\mathrm{E}_{q}\\{a(Y)\mid{\bf
x}\\}}{\mathrm{E}_{q}\\{\rho(Y)\mid{\bf x}\\}+\pi/(1-\pi)},$ $\displaystyle c$
$\displaystyle\equiv$ $\displaystyle 0.$

Next, we show that $\phi_{\rm eff}({\bf x},r,ry)$ is in the tangent space
${\cal T}$ of (S.1). We decompose $\phi_{\rm eff}({\bf x},r,ry)$ into

$\displaystyle\phi_{\rm eff}({\bf x},r,ry)$ $\displaystyle=$
$\displaystyle\frac{r}{\pi}\left\\{y\rho(y)-b({\bf
x})\rho(y)\right\\}+\frac{1-r}{1-\pi}\left\\{b({\bf x})-\theta\right\\}$
$\displaystyle=$ $\displaystyle r\\{a_{1}(y)+a_{2}({\bf
x},y)\\}+(1-r)[\mathrm{E}_{q}\\{a_{2}({\bf x},Y)\mid{\bf
x}\\}+\mathrm{E}_{q}\\{a_{3}(Y)\mid{\bf x}\\}],$

where

$\displaystyle a_{1}(y)$ $\displaystyle\equiv$ $\displaystyle 0,$
$\displaystyle a_{2}({\bf x},y)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\pi}\\{y\rho(y)-b({\bf x})\rho(y)\\},$ $\displaystyle
a_{3}(y)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\pi}\\{a(y)-y\rho(y)\\}-\frac{1}{1-\pi}\theta.$

Then it is easy to show $ra_{1}(y)\in{\cal T}_{\bm{\alpha}}$, and $ra_{2}({\bf
x},y)+(1-r)\mathrm{E}_{q}\\{a_{2}({\bf x},Y)\mid{\bf x}\\}\in{\cal
T}_{\bm{\beta}}$ because $\mathrm{E}\\{b({\bf X})\mid y\\}=y$. Further,
$(1-r)\mathrm{E}_{q}\\{a_{3}(Y)\mid{\bf x}\\}\in{\cal T}_{\bm{\gamma}}$ since

$\displaystyle\mathrm{E}_{q}\\{a_{3}(Y)\\}$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\mathrm{E}_{q}[\mathrm{E}_{q}\\{a(Y)\mid{\bf
X}\\}]-\frac{1}{\pi}\mathrm{E}_{q}\\{Y\rho(Y)\\}-\frac{1}{1-\pi}\theta$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\mathrm{E}_{q}[b({\bf
X})\mathrm{E}_{q}\\{\rho(Y)\mid{\bf
X}\\}]+\frac{1}{1-\pi}\mathrm{E}_{q}\\{b({\bf
X})\\}-\frac{1}{\pi}\mathrm{E}_{q}\\{Y\rho(Y)\\}-\frac{1}{1-\pi}\theta$
$\displaystyle=$ $\displaystyle 0.$

Hence $\phi_{\rm eff}({\bf x},r,ry)$ belongs to ${\cal T}$, which proves the
result. ∎

###  S.4 Proof of Proposition 2

The definition of $b^{*\star}({\bf x})$ immediately leads to
$\mathrm{E}\\{b^{*\star}({\bf X})\mid y\\}=y$. Therefore,

$\displaystyle\widehat{\theta}_{t}$ $\displaystyle\xrightarrow{P}$
$\displaystyle\mathrm{E}\left[\frac{R}{\pi}\rho^{*}(Y)\\{Y-b^{*\star}({\bf
X})\\}+\frac{1-R}{1-\pi}b^{*\star}({\bf X})\right]$ $\displaystyle=$
$\displaystyle\mathrm{E}_{p}[\rho^{*}(Y)\\{Y-b^{*\star}({\bf
X})\\}]+\mathrm{E}_{q}\\{b^{*\star}({\bf X})\\}$ $\displaystyle=$
$\displaystyle\mathrm{E}_{p}[\rho^{*}(Y)[Y-\mathrm{E}\\{b^{*\star}({\bf
X})\mid Y\\}])+\mathrm{E}_{q}[\mathrm{E}\\{b^{*\star}({\bf X})\mid Y\\}]$
$\displaystyle=$ $\displaystyle\mathrm{E}_{q}(Y).$

∎

###  S.5 Algorithm to Solve Equations (16) and (19) in Algorithms 1 and 2

We first illustrate how to solve (16). Let

$\displaystyle\phi^{*\star}(y,t)\equiv\rho^{*}(t)\sum_{i=1}^{n}\frac{p_{Y\mid{\bf
X}}^{\star}(t,{\bf x}_{i})}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})},$
(S.5)

then (16) is equivalently written as
$y=\int\phi^{*\star}(y,t)a^{*\star}(t)\differential t$, which is a Fredholm
integral equation of the first type. We can find the solution
$\widehat{a}^{*\star}$ through Landweber’s iterative method (Landweber 1951).
Choose a starting point $a_{0}$, then iterate the formula

$\displaystyle a_{k+1}(y)\leftarrow
a_{k}(y)+\int\phi^{*\star}(t,y)t\differential
t-\int\left\\{\int\phi^{*\star}(z,y)\phi^{*\star}(z,t)\differential
z\right\\}a_{k}(t)\differential t$

while $\int\\{a_{k+1}(t)-a_{k}(t)\\}^{2}\differential
t/\int\\{a_{k}(t)\\}^{2}\differential t$ is greater than a chosen tolerance.

For practical implementation, we can approximate the integration in (16) using
quadrature method. For a given quadrature rule, let the quadrature nodes be
${\bf t}\equiv(t_{1},\dots,t_{m})^{\rm T}$ and weights be ${\bf
w}\equiv(w_{1},\dots,w_{m})^{\rm T}$. Also, let $\widetilde{\bf
y}\equiv(\widetilde{y}_{1},\dots,\widetilde{y}_{l})^{\rm T},{\bf
a}^{*\star}\equiv\\{a^{*\star}(t_{1}),\dots,a^{*\star}(t_{m})\\}^{\rm T}$, and
${\bm{\Phi}}^{*\star}$ be an $l\times m$ matrix whose $(i,j)$ component is
$\phi^{*\star}(\widetilde{y}_{i},t_{j})$. Then (16) can be discretized as
$\widetilde{\bf y}={\bm{\Phi}}^{*\star}({\bf a}^{*\star}\cdot{\bf w})$ where
${\bf a}\cdot{\bf b}\equiv(a_{1}b_{1},\dots,a_{m}b_{m})^{\rm T}$ is the
Hadamard product. Also, the above iterative formula can be approximated as

$\displaystyle{\bf a}_{k+1}\leftarrow{\bf a}_{k}+l^{-1}{\bm{\Phi}}^{*\star\rm
T}\widetilde{\bf y}-l^{-1}{\bm{\Phi}}^{*\star\rm T}{\bm{\Phi}}^{*\star}({\bf
a}_{k}\cdot{\bf w}).$

We iterately use this formula with the starting point ${\bf a}_{0}$ as long
as$\|{\bf a}_{k+1}-{\bf a}_{k}\|_{2}^{2}/\|{\bf a}_{k}\|_{2}^{2}$ is greater
than a chosen tolerance. We summarize the algorithm in Algorithm S.1.

Input: a function $\phi^{*\star}(y,t)$ in (S.5) and a tolerance $\Delta$.

do

(a) adopt a set of evaluation points $\widetilde{\bf
y}=(\widetilde{y}_{1},\dots,\widetilde{y}_{l})^{\rm T}$;

(b) adopt a quadrature rule $({\bf t},{\bf w})$, where nodes ${\bf
t}=(t_{1},\dots,t_{m})^{\rm T}$ and weights ${\bf w}=(w_{1},\dots,w_{m})^{\rm
T}$;

(c) compute a matrix ${\bm{\Phi}}^{*\star}$ where
${\bm{\Phi}}^{*\star}_{(i,j)}=\phi^{*\star}(\widetilde{y}_{i},t_{j}),i=1,\dots,l,j=1,\dots,m$;

(d) Declare a starting point ${\bf a}_{0}=(a_{01},\dots,a_{0m})^{\rm T}$;

while  _$\|{\bf a}_{k+1}-{\bf a}_{k}\|_{2}^{2}/\|{\bf a}_{k}\|_{2}^{2}
>\Delta$_ do

${\bf a}_{k+1}\leftarrow{\bf a}_{k}+l^{-1}{\bm{\Phi}}^{*\star\rm
T}\widetilde{\bf y}-l^{-1}{\bm{\Phi}}^{*\star\rm T}{\bm{\Phi}}^{*\star}({\bf
a}_{k}\cdot{\bf w})$;

end while

Output: $\widehat{\bf a}^{*\star}\leftarrow{\bf a}_{k+1}$.

Algorithm S.1 Solving the integral equation (16)

We now describe how to solve (19). If $\widehat{\mathrm{E}}_{p}\\{a(Y)\mid{\bf
x}\\}=\int a(t)\widehat{p}_{Y\mid{\bf X}}(t,{\bf x}_{i})\differential t$ for
some $\widehat{p}_{Y\mid{\bf X}}$, then we can follow the above procedure to
obtain the solution while replacing $p_{Y\mid{\bf X}}^{\star}(t,{\bf x}_{i})$
in (S.5) by $\widehat{p}_{Y\mid{\bf X}}(t,{\bf x}_{i})$. Now, suppose
$\widehat{\mathrm{E}}_{p}\\{a(Y)\mid{\bf
x}\\}=\sum_{i=1}^{n_{1}}a(y_{i})w_{i}({\bf x})$ for some
$w_{i},i=1,\dots,n_{1}$. This form includes a general class of nonparametric
regression estimators, for instance, the Nadaraya-Watson estimator is of this
form with $w_{i}({\bf x})=r_{i}K_{h}({\bf x}-{\bf
x}_{i})/\sum_{j=1}^{n}r_{j}K_{h}({\bf x}-{\bf x}_{j})$. Then (19) is
equivalent to

$\displaystyle y$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\frac{\sum_{k=1}^{n_{1}}a^{*}(y_{k})\rho^{*}(y_{k})w_{k}({\bf
x}_{i})}{\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{n_{1}}\phi^{*}(y,y_{k},w_{k})a^{*}(y_{k}),$

where

$\displaystyle\phi^{*}(y,y_{k},w_{k})\equiv\sum_{i=1}^{n}\frac{\rho^{*}(y_{k})w_{k}({\bf
x}_{i})}{\widehat{\mathrm{E}}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\widehat{\mathrm{E}}_{p}\\{\rho^{*}(Y)\mid{\bf
x}_{i}\\}}\frac{r_{i}K_{h}(y-y_{i})}{\sum_{j=1}^{n}r_{j}K_{h}(y-y_{j})}.$
(S.6)

Now, let $\widetilde{\bf
y}\equiv(\widetilde{y}_{1},\dots,\widetilde{y}_{l})^{\rm T},{\bf
a}^{*}\equiv\\{a^{*}(y_{1}),\dots,a^{*}(y_{n_{1}})\\}^{\rm T}$, and
${\bm{\Phi}}^{*}$ be a $l\times n_{1}$ matrix whose $(i,j)$th component is
$\phi^{*}(\widetilde{y}_{i},y_{j},w_{j})$. Then we can find the approximate
solution $\widehat{\bf a}^{*}$ via iteratively using the formula

$\displaystyle{\bf a}_{k+1}\leftarrow{\bf a}_{k}+l^{-1}{\bm{\Phi}}^{*\rm
T}\widetilde{\bf y}-l^{-1}{\bm{\Phi}}^{*\rm T}{\bm{\Phi}}^{*}{\bf a}_{k}$

with a starting point ${\bf a}_{0}$, as long as $\|{\bf a}_{k+1}-{\bf
a}_{k}\|_{2}^{2}/\|{\bf a}_{k}\|_{2}^{2}$ is greater than a chosen tolerance.
We summarize the algorithm in Algorithm S.2.

Input: a function $\phi^{*}(y,y_{k},w_{k})$ in (S.6) and a tolerance $\Delta$.

do

(a) adopt a set of evaluation points $\widetilde{\bf
y}=(\widetilde{y}_{1},\dots,\widetilde{y}_{l})^{\rm T}$;

(b) compute a matrix ${\bm{\Phi}}^{*}$ where
${\bm{\Phi}}^{*}_{(i,j)}=\phi^{*}(\widetilde{y}_{i},y_{j},w_{j}),i=1,\dots,l,j=1,\dots,n_{1}$;

(c) Declare a starting point ${\bf a}_{0}=(a_{01},\dots,a_{0m})^{\rm T}$;

while  _$\|{\bf a}_{k+1}-{\bf a}_{k}\|_{2}^{2}/\|{\bf a}_{k}\|_{2}^{2}
>\Delta$_ do

${\bf a}_{k+1}\leftarrow{\bf a}_{k}+l^{-1}{\bm{\Phi}}^{*\rm T}\widetilde{\bf
y}-l^{-1}{\bm{\Phi}}^{*\rm T}{\bm{\Phi}}^{*}({\bf a}_{k}\cdot{\bf w})$;

end while

Output: $\widehat{\bf a}^{*}\leftarrow{\bf a}_{k+1}$.

Algorithm S.2 Solving the integral equation (19)

###  S.6 Proof of Lemma 2

From the definition of ${\cal L}^{*\star}(a)$, Conditions (A3), and (A4), it
is immediate that there exists a constant $0<c_{2}<\infty$ such that $\|{\cal
L}^{*\star}(a)\|_{\infty}\leq c_{2}\|a\|_{\infty}$. Now we show there is a
constant $0<c_{1}<\infty$ such that $c_{1}\|a\|_{\infty}\leq\|{\cal
L}^{*\star}(a)\|_{\infty}$. Note that if ${\cal L}^{*\star}$ is invertible, by
the bounded inverse theorem we have $\|{\cal L}^{*\star-1}(v)\|_{\infty}\leq
c_{1}^{-1}\|v\|_{\infty}$ for some constant $0<c_{1}<\infty$, i.e.,
$c_{1}\|a\|_{\infty}\leq\|{\cal L}^{*\star}(a)\|_{\infty}$. Hence it suffices
to show that ${\cal L}^{*\star}$ is invertible. We prove this by
contradiction. Suppose there are $a_{1}(y)$ and $a_{2}(y)$ such that ${\cal
L}^{*\star}(a_{1})(y)={\cal L}^{*\star}(a_{2})(y)=v(y)$ and $a_{1}\neq a_{2}$.
Then by Conditions (A1) and (A2), we have

$\displaystyle\mathrm{E}_{p}^{\star}\\{a_{1}(Y)\rho^{*}(Y)\mid{\bf
x}\\}\neq\mathrm{E}_{p}^{\star}\\{a_{2}(Y)\rho^{*}(Y)\mid{\bf x}\\}.$

Now, the efficient score calculated under the posited models is

$\displaystyle\phi_{\rm eff}^{*\star}({\bf x},r,ry)$ $\displaystyle=$
$\displaystyle\frac{r}{\pi}\rho^{*}(y)\left[y-\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf x}\\}}\right]$
$\displaystyle+\frac{1-r}{1-\pi}\left[\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}\\}}-\theta\right],$

where $a(y)$ satisfies ${\cal L}^{*\star}(a)(y)=v(y)$. Then letting $a=a_{1}$
and $a=a_{2}$ gives two distinct efficient scores, which contradicts the
uniqueness of the efficient score. Therefore, there is a unique solution
$a^{*\star}(y)$ for ${\cal L}^{*\star}(a^{*\star})(y)=v(y)$, hence ${\cal
L}^{*\star}$ is invertible. ∎

###  S.7 Proof of Theorem 1

We define

$\displaystyle b^{*\star}({\bf
x},a,{\bm{\zeta}})\equiv\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
x},{\bm{\zeta}}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x},{\bm{\zeta}}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x},{\bm{\zeta}}\\}},$

and analyze $b^{*\star}({\bf x},\widehat{a}^{*\star},\widehat{\bm{\zeta}})$.
First, under Conditions (A3)-(A6) and
$\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$, we have

$\displaystyle\widehat{\cal L}^{*\star}(a^{*\star})(y)$ $\displaystyle=$
$\displaystyle n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})b^{*\star}({\bf
x}_{i},a^{*\star},\widehat{\bm{\zeta}})$ (S.7) $\displaystyle=$ $\displaystyle
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\left\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})+\frac{\partial b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(n_{1}^{-1/2})\right\\}$
$\displaystyle=$ $\displaystyle n_{1}^{-1}\sum_{i=1}^{n}{\cal
L}_{i,h}^{*\star}(a^{*\star})(y)+\left[p_{Y}(y)\frac{\partial
E\\{b^{*\star}({\bf X},a^{*\star},{\bm{\zeta}})\mid
y\\}}{\partial{\bm{\zeta}}^{\rm
T}}+o_{p}(1)\right](\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle n_{1}^{-1}\sum_{i=1}^{n}{\cal
L}_{i,h}^{*\star}(a^{*\star})(y)+o_{p}(n_{1}^{-1/2})$

uniformly in $y$ since $E\\{b^{*\star}({\bf X},a^{*\star},{\bm{\zeta}})\mid
y\\}=y$. Here, we define

$\displaystyle{\cal L}_{i,h}^{*\star}(a)(y)\equiv
r_{i}K_{h}(y-y_{i})\frac{\mathrm{E}_{p}^{\star}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf x}_{i}\\}}.$

We also define $g_{i,h}(y)\equiv{\cal L}^{*\star-1}\\{v_{i,h}-{\cal
L}_{i,h}^{*\star}(a^{*\star})\\}(y)$. Note that
$\|E(V_{i,h})-v\|_{\infty}=O(h^{2})$ by Conditions (A4),(A5), then since
${\cal L}^{*\star-1}$ is a bounded linear operator by Lemma 2,

$\displaystyle\|E\\{{\cal L}^{*\star-1}(V_{i,h})\\}-{\cal
L}^{*\star-1}(v)\|_{\infty}$ $\displaystyle=$ $\displaystyle\|{\cal
L}^{*\star-1}\\{E(V_{i,h})\\}-{\cal L}^{*\star-1}(v)\|_{\infty}$
$\displaystyle=$ $\displaystyle\|{\cal
L}^{*\star-1}\\{E(V_{i,h})-v\\}\|_{\infty}$ $\displaystyle=$ $\displaystyle
O(h^{2}).$

Similarly, $\|E\\{({\cal L}^{*\star-1}{\cal
L}_{i,h}^{*\star})(a^{*\star})\\}-a^{*\star}\|_{\infty}=O(h^{2})$. Then using
${\cal L}^{*\star-1}(v)(y)=a^{*\star}(y)$, we get

$\displaystyle\|E(G_{i,h})\|_{\infty}$ $\displaystyle=$
$\displaystyle\|E\\{{\cal L}^{*\star-1}(V_{i,h})\\}-E\\{({\cal
L}^{*\star-1}{\cal L}_{i,h}^{*\star})(a^{*\star})\\}-\\{{\cal
L}^{*\star-1}(v)-a^{*\star}\\}\|_{\infty}$ (S.8) $\displaystyle\leq$
$\displaystyle\|E\\{{\cal L}^{*\star-1}(V_{i,h})\\}-{\cal
L}^{*\star-1}(v)\|_{\infty}+\|E\\{({\cal L}^{*\star-1}{\cal
L}_{i,h})(a^{*\star})\\}-a^{*\star}\|_{\infty}$ $\displaystyle=$
$\displaystyle O(h^{2}).$

Now, for any bounded function $a(y)$, $\|(\widehat{\cal L}^{*\star}-{\cal
L}^{*\star})(a)\|_{\infty}=O_{p}\left\\{(n_{1}h)^{-1/2}\log
n_{1}+h^{2}\right\\}=o_{p}(n_{1}^{-1/4})$ under Conditions (A3)-(A6) and the
assumption $\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$.
Similarly, $\|\widehat{v}-v\|_{\infty}=o_{p}(n_{1}^{-1/4})$. Then using that
${\cal L}^{*\star-1}$ is a bounded linear operator by Lemma 2,
$\widehat{a}^{*\star}(y)$ can be expressed as

$\displaystyle\widehat{a}^{*\star}(y)$ $\displaystyle=$
$\displaystyle\widehat{\cal L}^{*\star-1}(\widehat{v})(y)$ $\displaystyle=$
$\displaystyle\\{{\cal L}^{*\star}+(\widehat{\cal L}^{*\star}-{\cal
L}^{*\star})\\}^{-1}(\widehat{v})(y)$ $\displaystyle=$ $\displaystyle\\{{\cal
L}^{*\star-1}-{\cal L}^{*\star-1}(\widehat{\cal L}^{*\star}-{\cal
L}^{*\star}){\cal
L}^{*\star-1}\\}\\{v+(\widehat{v}-v)\\}(y)+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle a^{*\star}(y)+{\cal
L}^{*\star-1}(\widehat{v}-v)(y)-\\{{\cal L}^{*\star-1}(\widehat{\cal
L}^{*\star}-{\cal L}^{*\star})\\}(a^{*\star})(y)+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle a^{*\star}(y)+{\cal
L}^{*\star-1}\\{\widehat{v}-\widehat{\cal
L}^{*\star}(a^{*\star})\\}(y)+o_{p}(n_{1}^{-1/2})$ $\displaystyle=$
$\displaystyle
a^{*\star}(y)+n_{1}^{-1}\sum_{i=1}^{n}g_{i,h}(y)+o_{p}(n_{1}^{-1/2})$

uniformly in $y$, where the last step is by (S.7) and the definition of
$g_{i,h}(y)$. Then we further get

$\displaystyle\left|\frac{\partial b^{*\star}({\bf
x},\widehat{a}^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})-\frac{\partial b^{*\star}({\bf
x},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})\right|$ $\displaystyle=$
$\displaystyle\left|\frac{\partial}{\partial{\bm{\zeta}}^{\rm
T}}\left[\frac{\mathrm{E}_{p}^{\star}\\{(\widehat{a}^{*\star}-a^{*\star})(Y)\rho^{*}(Y)\mid{\bf
x},{\bm{\zeta}}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x},{\bm{\zeta}}\\}}\right](\widehat{\bm{\zeta}}-{\bm{\zeta}})\right|$
$\displaystyle=$
$\displaystyle\left|\frac{\mathrm{E}_{p}^{\star}\\{(\widehat{a}^{*\star}-a^{*\star})(Y)\rho^{*}(Y){\bf
S}_{\bm{\zeta}}^{\star\rm T}(Y,{\bf x},{\bm{\zeta}})\mid{\bf
x}\\}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})\right.$
$\displaystyle\left.-\mathrm{E}_{p}^{\star}\\{(\widehat{a}^{*\star}-a^{*\star})(Y)\rho^{*}(Y)\mid{\bf
x}\\}\frac{\mathrm{E}_{p}[\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\\}{\bf
S}_{\bm{\zeta}}^{\star\rm T}(Y,{\bf x},{\bm{\zeta}})\mid{\bf
x}]}{[\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}]^{2}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})\right|$ $\displaystyle\leq$
$\displaystyle\left\|\widehat{a}^{*\star}-a^{*\star}\right\|_{\infty}\frac{\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\|{\bf
S}_{\bm{\zeta}}^{\star}(Y,{\bf x},{\bm{\zeta}})\|_{2}\mid{\bf
x}\\}\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}}{\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}$
$\displaystyle+\left\|\widehat{a}^{*\star}-a^{*\star}\right\|_{\infty}\mathrm{E}_{p}^{\star}\\{\rho^{*}(Y)\mid{\bf
x}\\}\frac{\mathrm{E}_{p}^{\star}[\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\\}\|{\bf
S}_{\bm{\zeta}}^{\star}(Y,{\bf x},{\bm{\zeta}})\|_{2}\mid{\bf
x}]\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}}{[\mathrm{E}_{p}^{\star}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}]^{2}}$ $\displaystyle=$ $\displaystyle
o_{p}(n_{1}^{-1/4})O_{p}(n_{1}^{-1/2})=o_{p}(n_{1}^{-1/2}),$

because $\|\widehat{\bm{\zeta}}-{\bm{\zeta}}\|_{2}=O_{p}(n_{1}^{-1/2})$,
$\mathrm{E}_{p}^{\star}\\{\|{\bf S}_{\bm{\zeta}}^{\star}(Y,{\bf
x},{\bm{\zeta}})\|_{2}\mid{\bf x}\\}$ is bounded, and
$\|\widehat{a}^{*\star}-a^{*\star}\|_{\infty}=O_{p}\left\\{(n_{1}h)^{-1/2}\log
n_{1}+h^{2}\right\\}=o_{p}(n^{-1/4})$ by (S.8) and Condition (A6). Hence,
noting that $b^{*\star}({\bf x},a,{\bm{\zeta}})$ is linear with respect to
$a$, we get

$\displaystyle b^{*\star}({\bf x},\widehat{a}^{*\star},\widehat{\bm{\zeta}})$
$\displaystyle=$ $\displaystyle b^{*\star}({\bf
x},\widehat{a}^{*\star},{\bm{\zeta}})+\frac{\partial b^{*\star}({\bf
x},\widehat{a}^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(n_{1}^{-1/2})$ (S.9)
$\displaystyle=$ $\displaystyle b^{*\star}({\bf
x},a^{*\star},{\bm{\zeta}})+b^{*\star}({\bf
x},\widehat{a}^{*\star}-a^{*\star},{\bm{\zeta}})+\frac{\partial
b^{*\star}({\bf x},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(n_{1}^{-1/2})$ $\displaystyle=$
$\displaystyle b^{*\star}({\bf
x},a^{*\star},{\bm{\zeta}})+n_{1}^{-1}\sum_{i=1}^{n}b^{*\star}({\bf
x},g_{i,h},{\bm{\zeta}})+\frac{\partial b^{*\star}({\bf
x},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(n_{1}^{-1/2})$

uniformly in ${\bf x}$ by Condition (A4). Now we analyze $\widehat{\theta}$.
By the definition of $\widehat{\theta}$,

$\displaystyle\sqrt{n_{1}}(\widehat{\theta}-\theta)$ $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{y_{i}-b^{*\star}({\bf
x}_{i},\widehat{a}^{*\star},\widehat{\bm{\zeta}})\\}+\frac{1-r_{i}}{1-\pi}\\{b^{*\star}({\bf
x}_{i},\widehat{a}^{*\star},\widehat{\bm{\zeta}})-\theta\\}\right]$ (S.10)
$\displaystyle=$ $\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\phi_{\rm
eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})+\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-b^{*\star}({\bf
x}_{i},\widehat{a}^{*\star},\widehat{\bm{\zeta}})\\}$ $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\phi_{\rm eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})-T_{1}-T_{2}+o_{p}(1),$

where

$\displaystyle T_{1}$ $\displaystyle\equiv$ $\displaystyle
n_{1}^{-1/2}n^{-1}\sum_{i=1}^{n}\sum_{j=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}b^{*\star}({\bf
x}_{i},g_{j,h},{\bm{\zeta}}),$ $\displaystyle T_{2}$ $\displaystyle\equiv$
$\displaystyle
n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}\frac{\partial
b^{*\star}({\bf x}_{i},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}\sqrt{n_{1}}(\widehat{\bm{\zeta}}-{\bm{\zeta}}).$

Note that since $\|E(G_{j,h})\|_{\infty}=O(h^{2})$ by (S.8) and
$b^{*\star}({\bf x},a,{\bm{\zeta}})$ is linear with respect to $a$,

$\displaystyle\|E\\{b^{*\star}({\bf
x},G_{j,h},{\bm{\zeta}})\\}\|_{\infty}=\|b^{*\star}\\{{\bf
x},E(G_{j,h}),{\bm{\zeta}}\\}\|_{\infty}=O(h^{2}).$ (S.11)

Hence using the property of the U-statistic and Condition (A6), we can rewrite
$T_{1}$ as

$\displaystyle T_{1}$ $\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}E\\{b^{*\star}({\bf
x}_{i},G_{j,h},{\bm{\zeta}})\mid{\bf x}_{i},r_{i},r_{i}y_{i}\\}$
$\displaystyle+n_{1}^{-1/2}\sum_{j=1}^{n}\mathrm{E}\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*\star}({\bf
X}_{i},g_{j,h},{\bm{\zeta}})\mid{\bf x}_{j},r_{j},r_{j}y_{j}\right]$
$\displaystyle-n^{1/2}n_{1}^{-1/2}\mathrm{E}\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*\star}({\bf
X}_{i},G_{j,h},{\bm{\zeta}})\right]+O_{p}(n_{1}^{-1/2})$ $\displaystyle=$
$\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}\mathrm{E}\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*\star}({\bf
X}_{i},g_{j,h},{\bm{\zeta}})\mid{\bf
x}_{j},r_{j},r_{j}y_{j}\right]+O_{p}(n_{1}^{-1/2}+nn_{1}^{-1/2}h^{2})$
$\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}\int\left[\rho^{*}(y)\mathrm{E}\\{b^{*\star}({\bf
X},g_{j,h},{\bm{\zeta}})\mid y\\}p_{Y}(y)-\mathrm{E}\\{b^{*\star}({\bf
X},g_{j,h},{\bm{\zeta}})\mid y\\}q_{Y}(y)\right]dy+o_{p}(1)$ $\displaystyle=$
$\displaystyle n_{1}^{-1/2}\sum_{j=1}^{n}\int{\cal
L}^{*\star}(g_{j,h})(y)\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1)$
$\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}\int\\{v_{j,h}(y)-{\cal
L}_{j,h}^{*\star}(a^{*\star})(y)\\}\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1)$
$\displaystyle=$ $\displaystyle n_{1}^{-1/2}\sum_{j=1}^{n}r_{j}\int
K_{h}(y-y_{j})\\{y-b^{*\star}({\bf
x}_{j},a^{*\star},{\bm{\zeta}})\\}\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1).$

On the other hand, using $\mathrm{E}\\{b^{*\star}({\bf
X},a^{*\star},{\bm{\zeta}})\mid y\\}=y$, $T_{2}$ can be written as

$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\mathrm{E}\left[\left\\{\frac{R}{\pi}\rho^{*}(Y)-\frac{1-R}{1-\pi}\right\\}\frac{\partial
b^{*\star}({\bf X},a^{*\star},{\bm{\zeta}})}{\partial{\bm{\zeta}}^{\rm
T}}\right]\sqrt{n_{1}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(1)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial{\bm{\zeta}}^{\rm
T}}\mathrm{E}\left[\left\\{\frac{R}{\pi}\rho^{*}(Y)-\frac{1-R}{1-\pi}\right\\}Y\right]\sqrt{n_{1}}(\widehat{\bm{\zeta}}-{\bm{\zeta}})+o_{p}(1)$
$\displaystyle=$ $\displaystyle o_{p}(1).$

Therefore, (S.10) leads to

$\displaystyle\sqrt{n_{1}}(\widehat{\theta}-\theta)$ $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\phi_{\rm eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})+n_{1}^{-1/2}\sum_{i=1}^{n}r_{i}h({\bf
x}_{i},y_{i})+o_{p}(1)$ $\displaystyle=$ $\displaystyle
n^{-1/2}\sum_{i=1}^{n}\left\\{\sqrt{\pi}\phi_{\rm eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})+\frac{r_{i}}{\sqrt{\pi}}h({\bf
x}_{i},y_{i})\right\\}+o_{p}(1),$

where

$\displaystyle h({\bf x}_{i},y_{i})$ $\displaystyle\equiv$ $\displaystyle\int
K_{h}(y-y_{i})\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-y\\}\\{\rho^{*}(y)-\rho(y)\\}dy$
$\displaystyle=$ $\displaystyle\int K(t)\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-y_{i}-ht\\}\\{\rho^{*}(y_{i}+ht)-\rho(y_{i}+ht)\\}dt$
$\displaystyle=$ $\displaystyle\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-y_{i}\\}\\{\rho^{*}(y_{i})-\rho(y_{i})\\}+[\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-y_{i}\\}\\{{\rho^{*}}^{\prime\prime}(y_{i})-\rho^{\prime\prime}(y_{i})\\}-2\\{{\rho^{*}}^{\prime}(y_{i})-\rho^{\prime}(y_{i})\\}]$
$\displaystyle\times\frac{h^{2}}{2}\int t^{2}K(t)dt+O(h^{4})$

by Conditions (A2) and (A5). Hence, by Condition (A6), we have

$\displaystyle\sqrt{n_{1}}(\widehat{\theta}-\theta)$ $\displaystyle=$
$\displaystyle n^{-1/2}\sum_{i=1}^{n}\left[\sqrt{\pi}\phi_{\rm
eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})+\frac{r_{i}}{\sqrt{\pi}}\\{b^{*\star}({\bf
x}_{i},a^{*\star},{\bm{\zeta}})-y_{i}\\}\\{\rho^{*}(y_{i})-\rho(y_{i})\\}\right]+o_{p}(1).$

∎

###  S.8 Proof of Theorem 2

First let

$\displaystyle b^{*}({\bf x},a,\mathrm{E}_{p})$ $\displaystyle\equiv$
$\displaystyle\frac{\mathrm{E}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf x}\\}},$

and for any function $g(\cdot,\mathrm{E}_{p})$, its $k$th Gateaux derivative
with respect to $\mathrm{E}_{p}$ at $\mu_{1}$ in the direction $\mu_{2}$ as

$\displaystyle\frac{\partial^{k}g(\cdot,\mu_{1})}{\partial{\mathrm{E}_{p}}^{k}}(\mu_{2})$
$\displaystyle\equiv$
$\displaystyle\frac{\partial^{k}g(\cdot,\mu_{1}+h\mu_{2})}{\partial
h^{k}}\bigg{|}_{h=0}.$

We have

$\displaystyle\frac{\partial b^{*}({\bf
x},a,\mathrm{E}_{p})}{\partial\mathrm{E}_{p}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$
$\displaystyle\frac{(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}{[\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}]^{2}}$
$\displaystyle-\frac{(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}\mathrm{E}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{[\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}]^{2}}$ $\displaystyle=$ $\displaystyle b^{*}({\bf
x},a,\mathrm{E}_{p})\left[\frac{(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}}-\frac{(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}\right]$ $\displaystyle=$ $\displaystyle b^{*}({\bf
x},a,\mathrm{E}_{p})o_{p}(n_{1}^{-1/4}),$
$\displaystyle\frac{\partial^{2}b^{*}({\bf
x},a,\mu)}{\partial{\mathrm{E}_{p}}^{2}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$
$\displaystyle\frac{-2(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}{[\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}]^{3}}$
$\displaystyle\times\left[(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{a(Y)\rho^{*}(Y)\mid{\bf
x}\\}\mu\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf x}\\}\right.$
$\displaystyle\left.-(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}\mu\\{a(Y)\rho^{*}(Y)\mid{\bf x}\\}\right]$ $\displaystyle=$
$\displaystyle o_{p}(n_{1}^{-1/2})$

for any bounded $a(y)$ and $\mu(\cdot\mid{\bf x})$ since
$|(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})(\cdot\mid{\bf
x})|=o_{p}(n_{1}^{-1/4})$ by the assumption, and these hold uniformly with
respect to ${\bf x}$ by Condition (A4). Then by the Taylor expansion and mean
value theorem, for any bounded $a(y)$ and some $\alpha\in(0,1)$,

$\displaystyle b^{*}({\bf x},a,\widehat{\mathrm{E}}_{p})$ $\displaystyle=$
$\displaystyle b^{*}({\bf x},a,\mathrm{E}_{p})+\frac{\partial b^{*}({\bf
x},a,\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})+\frac{1}{2}\frac{\partial^{2}b^{*}\\{{\bf
x},a,\mathrm{E}_{p}+\alpha(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\\}}{\partial{\mathrm{E}_{p}}^{2}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
(S.12) $\displaystyle=$ $\displaystyle b^{*}({\bf
x},a,\mathrm{E}_{p})+\frac{\partial b^{*}({\bf
x},a,\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})+o_{p}(n_{1}^{-1/2}).$

Noting that $a^{*}(y)$ is bounded under Condition (A7), we further get

$\displaystyle\widehat{\cal L}^{*}(a^{*})(y)$ $\displaystyle=$ $\displaystyle
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})b^{*}({\bf
x}_{i},a^{*},\widehat{\mathrm{E}}_{p})$ $\displaystyle=$ $\displaystyle
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\left\\{b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})+\frac{\partial b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})+o_{p}(n_{1}^{-1/2})\right\\}$
$\displaystyle=$ $\displaystyle n_{1}^{-1}\sum_{i=1}^{n}{\cal
L}_{i,h}^{*}(a^{*})(y)+o_{p}(n_{1}^{-1/2})$

uniformly in $y$ by Condition (A4), where we define

$\displaystyle{\cal L}_{i,h}^{*}(a)(y)\equiv
r_{i}K_{h}(y-y_{i})\frac{\mathrm{E}_{p}\\{a(Y)\rho^{*}(Y)\mid{\bf
x}_{i}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)\mid{\bf
x}_{i}\\}+\pi/(1-\pi)\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf x}_{i}\\}}.$

The last equality above is because $E\\{b^{*}({\bf
X},a^{*},\mathrm{E}_{p})\mid y\\}=y$ from the definition of $a^{*}$, hence

$\displaystyle n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\frac{\partial
b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$ $\displaystyle
n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})\frac{\partial b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})-\frac{\partial[p_{Y}(y)E\\{b^{*}({\bf
X},a^{*},\mathrm{E}_{p})\mid
y\\}]}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$ $\displaystyle
n_{1}^{-1/4}\left[n_{1}^{-1}\sum_{i=1}^{n}r_{i}K_{h}(y-y_{i})n_{1}^{1/4}\frac{\partial
b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})-p_{Y}(y)E\left\\{n_{1}^{1/4}\frac{\partial
b^{*}({\bf
X},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\mid
y\right\\}\right]$ $\displaystyle=$ $\displaystyle
n_{1}^{-1/4}O_{p}\left\\{(n_{1}h)^{-1/2}\log{n_{1}}+h^{2}\right\\}$
$\displaystyle=$ $\displaystyle o_{p}(n_{1}^{-1/2})$

under Conditions (A4)-(A7). In addition, for any bounded function $a(y)$,

$\displaystyle\|(\widehat{\cal L}^{*}-{\cal
L}^{*})(a)\|_{\infty}=O_{p}\left\\{(n_{1}h)^{-1/2}\log{n_{1}}+h^{2}\right\\}+o_{p}(n_{1}^{-1/4})=o_{p}(n_{1}^{-1/4})$

under Conditions (A4)-(A7), and the assumption
$|\widehat{\mathrm{E}}_{p}\\{a(Y)\mid{\bf x}\\}-\mathrm{E}_{p}\\{a(Y)\mid{\bf
x}\\}|=o_{p}(n_{1}^{-1/4})$. Similarly,
$\|\widehat{v}-v\|_{\infty}=o_{p}(n_{1}^{-1/4})$. Then using that ${\cal
L}^{*-1}$ is a bounded linear operator by Lemma 2, $\widehat{a}^{*}(y)$ can be
expressed as

$\displaystyle\widehat{a}^{*}(y)$ $\displaystyle=$ $\displaystyle\\{{\cal
L}^{*}+(\widehat{\cal L}^{*}-{\cal L}^{*})\\}^{-1}(\widehat{v})(y)$
$\displaystyle=$ $\displaystyle\\{{\cal L}^{*-1}-{\cal L}^{*-1}(\widehat{\cal
L}^{*}-{\cal L}^{*}){\cal
L}^{*-1}\\}\\{v+(\widehat{v}-v)\\}(y)+o_{p}(n_{1}^{-1/2})$ $\displaystyle=$
$\displaystyle a^{*}(y)+{\cal L}^{*-1}(\widehat{v}-v)(y)-\\{{\cal
L}^{*-1}(\widehat{\cal L}^{*}-{\cal L}^{*})\\}(a^{*})(y)+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle a^{*}(y)+{\cal
L}^{*-1}\\{\widehat{v}-\widehat{\cal L}^{*}(a^{*})\\}(y)+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle
a^{*}(y)+n_{1}^{-1}\sum_{i=1}^{n}g_{i,h}(y)+o_{p}(n_{1}^{-1/2})$

uniformly in $y$ where $g_{i,h}(y)\equiv{\cal L}^{*-1}\\{v_{i,h}-{\cal
L}_{i,h}^{*}(a^{*})\\}(y)$. Hence, noting that $b^{*}({\bf
x},a,\mathrm{E}_{p})$ is linear with respect to $a$ and using (S.12), we get

$\displaystyle b^{*}({\bf x},\widehat{a}^{*},\widehat{\mathrm{E}}_{p})$
$\displaystyle=$ $\displaystyle b^{*}({\bf x},a^{*},\mathrm{E}_{p})+b^{*}({\bf
x},\widehat{a}^{*}-a^{*},\mathrm{E}_{p})+\frac{\partial b^{*}({\bf
x},\widehat{a}^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})+o_{p}(n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle b^{*}({\bf
x},a^{*},\mathrm{E}_{p})+n_{1}^{-1}\sum_{i=1}^{n}b^{*}({\bf
x},g_{i,h},\mathrm{E}_{p})+\frac{\partial b^{*}({\bf
x},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})+o_{p}(n_{1}^{-1/2})$

uniformly in ${\bf x}$ by Condition (A4). The last equality holds because for
any bounded $a(y)$,

$\displaystyle|b^{*}({\bf
x},a,\mathrm{E}_{p})|\leq\|a\|_{\infty}\frac{\mathrm{E}_{p}\\{\rho^{*}(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{*2}(Y)+\pi/(1-\pi)\rho^{*}(Y)\mid{\bf
x}\\}}=O(\|a\|_{\infty}),$

and

$\displaystyle\frac{\partial b^{*}({\bf
x},\widehat{a}^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})-\frac{\partial
b^{*}({\bf
x},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$ $\displaystyle\frac{\partial b^{*}({\bf
x},\widehat{a}^{*}-a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$ $\displaystyle b^{*}({\bf
x},\widehat{a}^{*}-a^{*},\mathrm{E}_{p})o_{p}(n_{1}^{-1/4})$ $\displaystyle=$
$\displaystyle\left[O_{p}(h^{2})+O_{p}\\{(n_{1}h)^{-1/2}\log(n_{1})\\}\right]o_{p}(n_{1}^{-1/4})$
$\displaystyle=$ $\displaystyle o_{p}(n_{1}^{-1/2}),$

where the third equality is due to $\|E(G_{i,h})\|_{\infty}=O(h^{2})$ by
(S.8), and the last equality holds because $h^{2}=o_{p}(n_{1}^{-1/2})$ and
$n_{1}\\{\log(n_{1})\\}^{-4}h^{2}\to\infty$ by Condition (A6). Now, by the
definition of $\widetilde{\theta}$,

$\displaystyle\sqrt{n_{1}}(\widetilde{\theta}-\theta)$ (S.13) $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{y_{i}-b^{*}({\bf
x}_{i},\widehat{a}^{*},\widehat{\mathrm{E}}_{p})\\}+\frac{1-r_{i}}{1-\pi}\\{b^{*}({\bf
x}_{i},\widehat{a}^{*},\widehat{\mathrm{E}}_{p})-\theta\\}\right]$
$\displaystyle=$ $\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\phi_{\rm
eff}^{*}({\bf
x}_{i},r_{i},r_{i}y_{i})+\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}\left\\{b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})-b^{*}({\bf
x}_{i},\widehat{a}^{*},\widehat{\mathrm{E}}_{p})\right\\}$ $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\phi_{\rm eff}^{*}({\bf
x}_{i},r_{i},r_{i}y_{i})-T_{1}-T_{2}+o_{p}(1),$

where

$\displaystyle T_{1}$ $\displaystyle\equiv$ $\displaystyle
n_{1}^{-1/2}n^{-1}\sum_{i=1}^{n}\sum_{j=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}b^{*}({\bf
x}_{i},g_{j,h},\mathrm{E}_{p}),$ $\displaystyle T_{2}$ $\displaystyle\equiv$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}\frac{\partial
b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p}).$

Note that $E\\{b^{*}({\bf x}_{i},G_{j,h},\mathrm{E}_{p})\mid{\bf
x}_{i}\\}=b^{*}\\{{\bf x}_{i},E(G_{j,h}),\mathrm{E}_{p}\\}=O(h^{2})$ by (S.8).
Then using the property of the U-statistic and Condition (A6), we can express
$T_{1}$ as

$\displaystyle T_{1}$ $\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}E\\{b^{*}({\bf
x}_{i},G_{j,h},\mathrm{E}_{p})\mid{\bf x}_{i},r_{i},r_{i}y_{i}\\}$
$\displaystyle+n_{1}^{-1/2}\sum_{j=1}^{n}E\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*}({\bf
X}_{i},g_{j,h},\mathrm{E}_{p})\mid{\bf x}_{j},r_{j},r_{j}y_{j}\right]$
$\displaystyle-n^{1/2}n_{1}^{-1/2}E\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*}({\bf
X}_{i},G_{j,h},\mathrm{E}_{p})\right]+O_{p}(n_{1}^{-1/2})$ $\displaystyle=$
$\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}E\left[\left\\{\frac{R_{i}}{\pi}\rho^{*}(Y_{i})-\frac{1-R_{i}}{1-\pi}\right\\}b^{*}({\bf
X}_{i},g_{j,h},\mathrm{E}_{p})\mid{\bf
x}_{j},r_{j},r_{j}y_{j}\right]+O_{p}(nn_{1}^{-1/2}h^{2}+n_{1}^{-1/2})$
$\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}\int\left[\rho^{*}(y)E\\{b^{*}({\bf
X},g_{j,h},\mathrm{E}_{p})\mid y\\}p_{Y}(y)-E\\{b^{*}({\bf
X},g_{j,h},\mathrm{E}_{p})\mid y\\}q_{Y}(y)\right]dy+o_{p}(1)$
$\displaystyle=$ $\displaystyle n_{1}^{-1/2}\sum_{j=1}^{n}\int{\cal
L}^{*}(g_{j,h})(y)\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1)$
$\displaystyle=$ $\displaystyle
n_{1}^{-1/2}\sum_{j=1}^{n}\int\\{v_{j,h}(y)-{\cal
L}_{j,h}^{*}(a^{*})(y)\\}\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1)$
$\displaystyle=$ $\displaystyle n_{1}^{-1/2}\sum_{j=1}^{n}r_{j}\int
K_{h}(y-y_{j})\\{y-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\left\\{\rho^{*}(y)-\rho(y)\right\\}dy+o_{p}(1),$

where

$\displaystyle\int K_{h}(y-y_{j})\\{y-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\left\\{\rho^{*}(y)-\rho(y)\right\\}dy$
$\displaystyle=$ $\displaystyle\int K(t)\\{y_{j}+ht-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\\{\rho^{*}(y_{j}+ht)-\rho(y_{j}+ht)\\}dt$
$\displaystyle=$ $\displaystyle\\{y_{j}-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\\{\rho^{*}(y_{j})-\rho(y_{j})\\}$
$\displaystyle+\left[\\{y_{j}-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\\{{\rho^{*}}^{\prime\prime}(y_{j})-\rho^{\prime\prime}(y_{j})\\}+2\\{{\rho^{*}}^{\prime}(y_{j})-\rho^{\prime}(y_{j})\\}\right]\frac{h^{2}}{2}\int
t^{2}K(t)dt+O(h^{4})$ $\displaystyle=$ $\displaystyle\\{y_{j}-b^{*}({\bf
x}_{j},a^{*},\mathrm{E}_{p})\\}\\{\rho^{*}(y_{j})-\rho(y_{j})\\}+o(n_{1}^{-1/2})$

under Conditions (A2), (A5), and (A6). On the other hand, using
$E\\{b^{*}({\bf X},a^{*},\mathrm{E}_{p})\mid y\\}=y$, $T_{2}$ can be written
as

$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\sqrt{n_{1}}n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}\frac{\partial
b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle-\sqrt{n_{1}}\frac{\partial}{\partial\mathrm{E}_{p}}E\left[\left\\{\frac{R}{\pi}\rho^{*}(Y)-\frac{1-R}{1-\pi}\right\\}b^{*}({\bf
X},a^{*},\mathrm{E}_{p})\right](\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})$
$\displaystyle=$ $\displaystyle
n_{1}^{1/4}\left(n^{-1}\sum_{i=1}^{n}\left\\{\frac{r_{i}}{\pi}\rho^{*}(y_{i})-\frac{1-r_{i}}{1-\pi}\right\\}n_{1}^{1/4}\frac{\partial
b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})}{\partial{\mathrm{E}_{p}}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\right.$
$\displaystyle\left.-E\left[\left\\{\frac{R}{\pi}\rho^{*}(Y)-\frac{1-R}{1-\pi}\right\\}n_{1}^{1/4}\frac{\partial
b^{*}({\bf
X},a^{*},\mathrm{E}_{p})}{\partial\mathrm{E}_{p}}(\widehat{\mathrm{E}}_{p}-\mathrm{E}_{p})\right]\right)$
$\displaystyle=$ $\displaystyle n_{1}^{1/4}O_{p}(n^{-1/2})$ $\displaystyle=$
$\displaystyle o_{p}(1).$

Therefore, (S.13) leads to

$\displaystyle\sqrt{n_{1}}(\widetilde{\theta}-\theta)$ $\displaystyle=$
$\displaystyle n^{-1/2}\sum_{i=1}^{n}\left[\sqrt{\pi}\phi_{\rm eff}^{*}({\bf
x}_{i},r_{i},r_{i}y_{i})+\frac{r_{i}}{\sqrt{\pi}}\\{b^{*}({\bf
x}_{i},a^{*},\mathrm{E}_{p})-y_{i}\\}\\{\rho^{*}(y_{i})-\rho(y_{i})\\}\right]+o_{p}(1).$

∎

###  S.9 Proposed Doubly Flexible Estimation for ${\bm{\theta}}$ such that
$\mathrm{E}_{q}\\{{\bf U}({\bf X},Y,{\bm{\theta}})\\}={\bf 0}$

####  S.9.1 General Approach to Estimating ${\bm{\theta}}$

Let $\dim({\bf U})=\dim({\bm{\theta}})$, $\mathrm{E}_{q}\\{\partial{\bf
U}({\bf X},Y,{\bm{\theta}})/\partial{\bm{\theta}}^{\rm T}\\}$ be invertible,
and ${\bf A}\equiv[\mathrm{E}_{q}\\{\partial{\bf U}({\bf
X},Y,{\bm{\theta}})/\partial{\bm{\theta}}^{\rm T}\\}]^{-1}$. We first
establish ${\cal F}$, the family of all influence functions for estimating
${\bm{\theta}}$. In Section S.10.1, we show that

$\displaystyle{\cal F}\equiv\left[\frac{r}{\pi}[\rho(y){\bf A}\\{{\bf U}({\bf
x},y,{\bm{\theta}})-{\bf b}({\bf x})\\}+{\bf c}]+\frac{1-r}{1-\pi}\\{{\bf
A}{\bf b}({\bf x})-{\bf c}\\}:\mathrm{E}\\{{\bf b}({\bf X})\mid
y\\}=\mathrm{E}\\{{\bf U}({\bf X},y,{\bm{\theta}})\mid y\\},\forall{\bf
c}\right].$

We now derive the efficient influence function ${\bm{\phi}}_{\rm eff}({\bf
x},r,ry)$ which corresponds to the semiparametric efficiency bound and also
provides guidance on constructing flexible estimators for ${\bm{\theta}}$. The
derivation is provided in Section S.10.2.

######  Proposition S.3.

The efficient influence function ${\bm{\phi}}_{\rm eff}({\bf x},r,ry)$ for
${\bm{\theta}}$ is

$\displaystyle{\bm{\phi}}_{\rm eff}({\bf x},r,ry)$ $\displaystyle=$
$\displaystyle\frac{r}{\pi}\rho(y){\bf A}\left[{\bf U}({\bf
x},y,{\bm{\theta}})-\frac{\mathrm{E}_{p}\\{{\bf U}({\bf
x},Y,{\bm{\theta}})\rho^{2}(Y)+{\bf a}(Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)+\pi/(1-\pi)\rho(Y)\mid{\bf x}\\}}\right]$
$\displaystyle+\frac{1-r}{1-\pi}\frac{{\bf A}\mathrm{E}_{p}\\{{\bf U}({\bf
x},Y,{\bm{\theta}})\rho^{2}(Y)+{\bf a}(Y)\rho(Y)\mid{\bf
x}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)+\pi/(1-\pi)\rho(Y)\mid{\bf x}\\}},$

where ${\bf a}(y)$ satisfies

$\displaystyle\mathrm{E}\left[\frac{\mathrm{E}_{p}\\{{\bf U}({\bf
X},Y,{\bm{\theta}})\rho^{2}(Y)+{\bf a}(Y)\rho(Y)\mid{\bf
X}\\}}{\mathrm{E}_{p}\\{\rho^{2}(Y)+\pi/(1-\pi)\rho(Y)\mid{\bf X}\\}}\mid
y\right]=\mathrm{E}\\{{\bf U}({\bf X},y,{\bm{\theta}})\mid y\\}.$

The efficient estimator can be obtained by solving the estimating equation
$\sum_{i=1}^{n}{\bm{\phi}}_{\rm eff}({\bf x}_{i},r_{i},r_{i}y_{i})={\bf 0}$,
where the efficient influence function ${\bm{\phi}}_{\rm eff}({\bf x},r,ry)$
is given in Proposition S.3. However, constructing ${\bm{\phi}}_{\rm eff}({\bf
x},r,ry)$ requires two possibly unknown functions $p_{Y\mid{\bf X}}(y,{\bf
x})$ and $\rho(y)$. We will consider replacing the two functions by working
models $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ and $\rho^{*}(y)$.

####  S.9.2 Proposed Estimator $\widehat{\bm{\theta}}$: Doubly Flexible in
$\rho^{*}(y)$ and $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$

Under the adopted working models $p_{Y\mid{\bf X}}^{\star}(y,{\bf x})$ and
$\rho^{*}(y)$, we obtain an estimator $\widehat{\bm{\theta}}_{t}$ by solving
the estimating equation $\sum_{i=1}^{n}{\bm{\phi}}_{\rm eff}^{*\star}({\bf
x}_{i},r_{i},r_{i}y_{i})={\bf 0}$. Given that ${\bf A}$ is invertible, this is
equivalent to solving

$\displaystyle\sum_{i=1}^{n}\left[\frac{r_{i}}{\pi}\rho^{*}(y_{i})\\{{\bf
U}({\bf x}_{i},y_{i},{\bm{\theta}})-{\bf b}^{*\star}({\bf