#  Biased-sample empirical likelihood weighting: an alternative to inverse
probability weighting

Yukun Liu  KLATASDS-MOE, School of Statistics, East China Normal University,
Shanghai 200062, China  Yan Fan Corresponding author<EMAIL_ADDRESS>School of Statistics and Information, Shanghai University of International
Business and Economics, Shanghai 201620, China

###### Abstract

Inverse probability weighting (IPW) is widely used in many areas when data are
subject to unrepresentativeness, missingness, or selection bias. An inevitable
challenge with the use of IPW is that the IPW estimator can be remarkably
unstable if some probabilities are very close to zero. To overcome this
problem, at least three remedies have been developed in the literature:
stabilizing, thresholding, and trimming. However the final estimators are
still IPW type estimators, and inevitably inherit certain weaknesses of the
naive IPW estimator: they may still be unstable or biased. We propose a
biased-sample empirical likelihood weighting (ELW) method to serve the same
general purpose as IPW, while completely overcoming the instability of IPW-
type estimators by circumventing the use of inverse probabilities. The ELW
weights are always well defined and easy to implement. We show theoretically
that the ELW estimator is asymptotically normal and more efficient than the
IPW estimator and its stabilized version for missing data problems and unequal
probability sampling without replacement. Its asymptotic normality is also
established under unequal probability sampling with replacement. Our
simulation results and a real data analysis indicate that the ELW estimator is
shift-equivariant, nearly unbiased, and usually outperforms the IPW-type
estimators in terms of mean square error.

Keywords: causal inference, empirical likelihood, inverse probability
weighting, missing data, unequal probability sampling

##  1 Introduction

Inverse probability weighting (IPW) has long been accepted as the standard
estimation procedure under unequal probability samplings with and without
replacement ever since the work of Hansen and Hurwitz (1943) and Horvitz and
Thompson (1952). IPW always produces an unbiased or asymptotically unbiased
estimator with an elegant expression, regardless of the complexity of the
underlying sampling plan, and this method therefore enjoys great popularity.
As well as survey sampling, it has been widely used in many other areas,
including missing data problems (Robins et al., 1994; Wooldridge, 2007; Tan,
2010; Kim and Shao, 2014), treatment effect estimation or program evaluation
(Rosenbaum and Rubin, 1983; Rosenbaum, 2002; Imbens and Wooldridge, 2009;
Hirano et al., 2003; Cattaneo, 2010; Young et al., 2019; Zhao, 2019; Tan,
2020), personalized medicine (Zhang et al., 2012; Jiang et al., 2017), and
survival data analysis (Robins and Rotnitzky, 1992; Robins, 1993; Bang and
Tsiatis, 2000; Ma and Yin, 2011; Dong et al., 2020), where IPW is renamed
inverse probability of censoring weighting. In recent years, accompanied by
optimal subsampling, the IPW method has also proved to be an effective
approach to validate statistical inferences for big data (Wang et al., 2018,
2019; Yu et al., 2020).

Roughly speaking, these applications of the IPW method can be categorized into
the following three representative types.

######  Problem A (Missing data problem).

Let $Z=(Y,X),$ with $Y$ being a response variable that is subject to
missingness and $X$ an always-observed covariate. Denote by $D$ a non-
missingness indicator, with $D=1$ if $Y$ is observed and 0 otherwise. Suppose
$(D_{i},D_{i}Y_{i},X_{i}),$ $i=1,2,\ldots,N,$ are $N$ independent and
identical copies of $(Y,X,D)$. The goal is to estimate
$\theta={\mathbb{E}}(Y)$.

######  Problem B (Unequal probability sampling without replacement, UPS-WOR).

Let $\\{Z_{1},\ldots,Z_{N}\\}$ be a population of interest with known size $N$
and let $\\{(D_{i},D_{i}Z_{i}):i=1,2,\ldots,N\\}$ be a sample drawn by an UPS-
WOR. The goal is to estimate the population mean
$\theta=(1/N)\sum_{i=1}^{N}Z_{i}$.

######  Problem C (Unequal probability sampling with replacement, UPS-WR).

The population and the parameter of interest are the same as in Problem B,
except that the sample is drawn by a UPS-WR.

Special cases of Problem A include treatment effect estimation in the
framework of Rubin’s potential outcomes (Rubin, 1974), as well as program
evaluation in economics and other social sciences. In a broad sense,
estimation of optimal treatment regime and survival data analysis both belong
to Problem A. Poisson sampling, a popular unequal probability sampling method,
can also be regarded as a special case of Problem A. See Section 3. Optimal
subsamplings in the context of big data are special cases of Problems B and C.

Through weighting the observations by the reciprocal of a certain probability
of inclusion in the sample, the IPW estimator is able to account for
unrepresentativeness, missingness or selection bias caused by non-random lack
of information or non-random selection of observations. However, the IPW
estimator can be highly unstable if there are extremely small probabilities,
which can result in biased estimation or poor finite-sample performance of the
accompanying asymptotic-normality-based inference (Busso et al., 2014; Kang
and Schafer, 2007; Robins et al., 2007; Imbens and Wooldridge, 2009; Cao et
al., 2009; Han et al., 2019). As pointed out by Robins et al. (2007) with
regard to double-robust estimators (which are IPW-type estimators) in missing
data problems, ‘Whenever the “inverse probability” weights are highly
variable, …, a small subset of the sample will have extremely large weights
relative to the remainder of the sample. In this setting, no estimator of the
marginal mean $\mu={\mathbb{E}}(Y)$ can be guaranteed to perform well.’ In
casual inference with observational studies, this is the well-known limited-
or non-overlap problem in covariate distributions in different treatment
groups (Crump et al., 2009; Khan and Tamer, 2010; Yang and Ding, 2018). The
IPW estimator becomes inflated disproportionately or even breaks down in
survival analysis when the number of patients at risk in the tails of the
survival curves of censoring times is too small (Robins and Finkelstein, 2000;
Dong et al., 2020). To guarantee that the IPW estimator possesses consistency,
asymptotic normality, and satisfactory finite-sample performance, it is usual
to impose an unnatural lower boundedness assumption on the probabilities
(Rosenbaum and Rubin, 1983; Mccaffrey et al., 2013; Sun and Tchetgen Tchetgen,
2018), although tiny probabilities are frequently encountered in practice,
especially when the propensity scores are estimated from data (Yang and Ding,
2018; Ma and Wang, 2020).

To overcome this notorious problem, at least three remedies have been proposed
in the literature: stabilizing, thresholding, and trimming. The stabilizing
method (Hájek, 1971) rescales the IPW estimator so that the weights sum to 1
(Kang and Schafer, 2007). Although straightforward, it can often sharply
reduce the instability of the IPW estimator. The thresholding method, proposed
by Zong et al. (2019) for solving Problem B, replaces those probabilities that
are less than a given threshold by that threshold while keeping others
unchanged. The parameter of interest is then estimated by IPW with the
modified probabilities. Zong et al. (2019) proposed an easy-to-use threshold
determining procedure and showed that, in general, the resulting IPW estimator
works better than the naive IPW estimator. This method can reduce the negative
effect of highly heterogeneous inclusion probabilities, and hence leads to
improved estimation efficiency, although at the cost of an estimation bias.

The trimming method excludes those observations with probabilities less than a
given threshold or, equivalently, sets their weights to zero (Crump et al.,
2009; Sasaki and Ura, 2018). However, the exact amount of trimming is usually
ad hoc and can affect the performance of the IPW estimator and the
corresponding confidence interval in nontrivial ways. Ma and Wang (2020)
systematically investigated the large-sample behaviour of the IPW estimator
after trimming and found it to be sensitive to the choice of trimming
threshold and subject to a non-negligible bias. They proposed a bias-corrected
and trimmed IPW estimator with an adaptively trimming threshold, and a
Studentized bootstrap procedure for interval estimation. Their estimator was
shown to be insensitive to small probability weights as well as being able to
correct the bias caused by trimming. The bias correction technique is built on
local polynomial regressions (Fan and Gijbels, 1996), which further require a
bandwidth. Inappropriate choices of the trimmed threshold and the bandwidth
may affect the performance of their estimator. More importantly, the bias
correction technique depends on the target quantity (e.g. $Y$ in Problem A) to
be weighted, which makes their method inapplicable to weighted optimization
problems, such as optimal treatment regime estimation (Zhang et al., 2012).

The final point estimators of the stabilizing, trimming, and thresholding
methods are all based on IPW, although they adopt different strategies to
reduce the detrimental effect of extremely small probabilities. These IPW-type
estimators inevitably inherit certain weaknesses of the naive IPW estimator:
they are either still unstable or biased. Also, the accompanying intervals,
regardless of whether they are asymptotic-normality-based or resampling-based,
often exhibit much undercoverage. See our simulation results in Section 4.

In this paper, we propose a biased-sample empirical likelihood weighting (ELW)
estimation method to serve the same general purpose as IPW in handling
incomplete or biased data while overcoming its instability. We systematically
investigate its finite- and large-sample properties in the context of missing
data problems and unequal probability samplings with and without replacement.
The proposed ELW estimation method has several advantages over the IPW-type
methods and the usual EL method.

  1. 1.

The ELW method circumvents the use of inverse probabilities and therefore
never suffers from extremely small or even zero selection probabilities. It
takes the maximum empirical likelihood estimates of the probability masses of
a multinomial distribution as weights, which always range from 0 to 1. This is
the most significant advantage of the ELW method over IPW and its variants.

  2. 2.

The ELW weights are always well defined. By contrast, the usual empirical
likelihood weights suffer from the well-known convex hull constraint or the
empty-set problem: they are undefined if the origin lies outside the convex
hull of certain transformed data points (Tsao, 2004; Chen et al., 2008; Liu
and Chen, 2010).

  3. 3.

Like the stabilized IPW estimator, the ELW weights always sum to 1, which
gives the ELW estimator the nice property of shift-equivariance.
Unfortunately, the naive IPW estimator, the trimmed IPW estimator of Zong et
al. (2019), and the IPW estimator of Ma and Wang (2020) are all sensitive to a
location shift in the response or the parameter of interest.

  4. 4.

The ELW weights are very convenient to calculate. Their calculation involves
only solving a univariate rational equation, which can be done efficiently by
the commonly used bisection algorithm. In contrast to the IPW estimator of Ma
and Wang (2020), the ELW estimator is free of any tuning parameter and is
hence more computationally efficient. The ELW weights depend only on the
propensity scores and the total sample size, and therefore the ELW method is
directly applicable to weighted optimization problems.

  5. 5.

As we shall show, the ELW estimator is theoretically more efficient than the
IPW estimator in at least two scenarios: missing data problems and unequal
probability samplings without replacement. Treatment effect estimation in
observational studies under the potential outcome framework of Rubin (1974)
can be regarded as a special case of missing data problems. This is a bonus of
ELW, since the construction of the ELW weights makes use of side information.
Under unequal probability sampling with replacement, we cannot tell which of
the ELW and IPW estimators wins theoretically. Nevertheless, our simulation
results indicate that the ELW estimator often has smaller mean square errors
and the accompanying interval has better coverage accuracy in most cases.

A crucial requirement of ELW is knowledge of the size of the finite population
of interest or a larger independent and identically distributed sample that
includes the observed data as a subsample. This is also required by the
original IPW method and some of its variants, and is available in most
situations. For example, in missing data problems, the size of the overall
dataset is clearly known, and in unequal probability sampling problems, the
size of the finite population from which the sample was drawn is usually known
a priori, since we need to construct a sampling frame before sampling. This
mild requirement implies that the ELW method has many potential applications
beyond missing data problems, sample surveys and casual inference.

The remainder of this article is organized as follows. In Section 2, we
introduce the ELW method by studying Problem A, namely estimating the
population mean of a response when data are subject to missingness. In Section
3, we extend the ELW method to unequal probability sampling by solving
Problems B and C. A simulation study and a real-life data analysis are
conducted in Sections 4 and 5 to demonstrate the usefulness and advantage of
the ELW method. Section 6 concludes with some discussion. All technical proofs
can be found in the supplementary material.

##  2 Empirical likelihood weighting

For ease of exposition, for the time being, we assume that the inclusion
probability or propensity score $\pi(Z)={\rm pr}(D=1|Z)$ is completely known
and always positive, although our method allows $\pi(Z)$ to take zero values.
The case with unknown propensity score is considered in Section 2.4. We take
the parameter of interest to be $\theta={\mathbb{E}}\\{g(Z)\\}$ for a user-
specific function $g(\cdot)$, for example $g(Z)=Y$ in Problem A.

Denote the data by $\\{(D_{i},D_{i}Z_{i}),i=1,2,\ldots,N\\}$, with
$Z_{i}=(Y_{i},X_{i})$ or simply $\\{z_{i},i=1,2,\ldots,n\\}$, where
$z_{i}=(y_{i},x_{i})$ and $n=\sum_{j=1}^{N}D_{j}$; the covariates $X_{i}$ with
$D_{i}=0$ do not come into play in most of this paper. The data
$\\{z_{i},i=1,2,\ldots,n\\}$ is in fact a biased sample of the underlying
population if all $\pi(z_{i})$ are not equal. The standard IPW estimator of
$\theta$ is

$\displaystyle\hat{\theta}_{\rm
IPW}=\frac{1}{N}\sum_{i=1}^{N}\frac{D_{i}}{\pi(Z_{i})}g(Z_{i})=\frac{1}{N}\sum_{i=1}^{n}\frac{g(z_{i})}{\pi(z_{i})}.$
(1)

This expression for the IPW estimator indicates that it becomes extremely
unstable when some of the $\pi(Z_{i})$ with $D_{i}=1$ are close to zero, and
that the terms with $D_{i}=0$ actually contribute nothing to it. The Hájek
estimator, or stabilized IPW (SIPW) estimator,

$\displaystyle\hat{\theta}_{\rm
SIPW}=\frac{\sum_{i=1}^{N}D_{i}g(Z_{i})/\pi(Z_{i})}{\sum_{j=1}^{N}D_{j}/\pi(Z_{j})}=\frac{\sum_{i=1}^{n}g(z_{i})/\pi(z_{i})}{\sum_{j=1}^{n}1/\pi(z_{j})}$

replaces $N$ by $\sum_{j=1}^{N}D_{j}/\pi(Z_{j})$, so that the weights of
$g(Z_{i})$ sum to 1 and the estimator becomes more stable. Since the size $N$
is known, the zero-value $D_{i}$ together with the other single-value $D_{i}$
contain information about
${\mathbb{E}}(D)={\mathbb{E}}\\{{\mathbb{E}}(D|Z)\\}={\mathbb{E}}\\{\pi(Z)\\}$.
The IPW estimator and its variants ignore such side information, and are not
able to utilize it as well, and they consequently have potential losses of
efficiency. As a popular and flexible non-parametric technique, empirical
likelihood (EL) (Owen, 1988, 1990, 2001) can conveniently and efficiently make
use of side information to achieve improvements in efficiency. This motivates
us to develop a biased-sample empirical likelihood weighting (ELW) estimation
method to serve the same purpose as the IPW estimator, while overcoming its
instability and improving its estimation efficiency.

###  2.1 ELW estimator

Let the distribution function of $Z$ be $F(z)={\rm pr}(Z\leq z)$, where the
inequality holds element-wise for vector-valued $Z$. To estimate
$\theta={\mathbb{E}}\\{g(Z)\\}=\int g(z)\,dF(z)$, it suffices to estimate
$F(z)$. We consider the problem of estimating $F$ by discarding those $Z_{i}$
with $D_{i}=0$, although these quantities may be partially accessible. The
likelihood based on the remaining data is

$\displaystyle\tilde{L}=\prod_{i=1}^{N}\\{\pi(Z_{i}){\rm
pr}(Z_{i})\\}^{D_{i}}(1-\alpha)^{1-D_{i}},$ (2)

where $\alpha={\rm pr}(D=1)={\mathbb{E}}\\{\pi(Z)\\}$. In the principle of EL,
we model $F(z)$ by a discrete distribution function
$F(z)=\sum_{i=1}^{N}p_{i}I(Z_{i}\leq z)$. This function is not well defined,
because those $Z_{i}$ with $D_{i}=0$ are not fully available. Nevertheless,
our final estimator, say $\hat{p}_{i}$, of $p_{i}$ satisfies
$\hat{p}_{i}=\hat{p}_{i}D_{i}$, making the resulting maximum EL estimators of
$F$ and $\theta$ well-defined statistics.

With $p_{i}$ in place of ${\rm pr}(Z_{i})$ in (2) and taking logarithms, we
have the empirical log-likelihood

$\displaystyle\tilde{\ell}=\sum_{i=1}^{N}[D_{i}\log(p_{i})+D_{i}\log\\{\pi(Z_{i})\\}+(1-D_{i})\log(1-\alpha)].$
(3)

Those $p_{i}$ that are feasible satisfy

$\displaystyle p_{i}\geq
0,\quad\sum_{i=1}^{N}p_{i}=1,\quad\sum_{i=1}^{N}p_{i}\\{\pi(Z_{i})-\alpha\\}=0.$
(4)

We emphasize that although those $Z_{i}$ with $D_{i}=0$ appear in
$\tilde{\ell}$ and $\sum_{i=1}^{N}p_{i}\\{\pi(Z_{i})-\alpha\\}=0$, they have
no likelihood contribution or any influence on the resulting EL method.

The proposed EL estimator of $F(z)$, or equivalently of the $p_{i}$, is
obtained by maximizing the empirical log-likelihood (3) subject to (4). For
fixed $\alpha$, the maximum of the log-EL in (3) subject to (4) is attained at

$\displaystyle
p_{i}=\frac{1}{n}\frac{D_{i}}{1+\lambda(\alpha)\\{\pi(Z_{i})-\alpha\\}},$ (5)

where $\lambda(\alpha)$ satisfies

$\displaystyle\frac{1}{n}\sum_{i=1}^{N}\frac{D_{i}}{1+\lambda(\alpha)\\{\pi(Z_{i})-\alpha\\}}\\{\pi(Z_{i})-\alpha\\}=0.$
(6)

Putting (5) into (3) gives the profile log-EL of $\alpha$ (up to a constant
that is independent of $\alpha$)

$\displaystyle\ell(\alpha)=\sum_{i=1}^{N}\\{-D_{i}\log[1+\lambda(\alpha)\\{\pi(Z_{i})-\alpha\\}]+(1-D_{i})\log(1-\alpha)\\}.$

This immediately gives $\hat{\alpha}=\arg\max\ell(\alpha)$, the EL estimator
of $\alpha$. Accordingly, the EL estimators of $p_{i}$ and $F(z)$ are

$\displaystyle\hat{p}_{i}=\frac{1}{n}\frac{D_{i}}{1+\lambda(\hat{\alpha})\\{\pi(Z_{i})-\hat{\alpha}\\}}$
(7)

and $\hat{F}(z)=\sum_{i=1}^{N}\hat{p}_{i}I(Z_{i}\leq z)$. Finally, the EL
estimator or the ELW estimator of $\theta$ is

$\displaystyle\hat{\theta}_{{\rm ELW}}=\int
g(z)d\hat{F}(z)=\sum_{i=1}^{N}\hat{p}_{i}g(Z_{i}).$ (8)

Obviously, both $\hat{F}(z)$ and $\hat{\theta}$ are well-defined statistics
because $\hat{p}_{i}=D_{i}\hat{p}_{i}$.

######  Remark 1.

The focus of this paper is the derivation of a better weighting method than
IPW. We achieve this by taking the probability masses of the maximum empirical
likelihood estimator of $F(z)$ as weights. Our ELW method requires the total
sample size $N$ and the (estimated) selection probabilities. After
multiplication by ${N\choose n},$ the likelihood in (2) is a degenerate case
(with only one capture occasion) of the full likelihood for capture–recapture
data of Liu et al. (2017), in which $N$ is unknown and the focus of
estimation. The case with known $N$ has been studied by Li and Qin (1998) for
biased and truncated data, although their interest was not in better
weighting, but rather in estimating the distributions of the target and
truncation variables, together with the unknown parameters involved.

######  Remark 2.

When the propensity scores are available for those with $D_{i}=1,$ the ELW
method works even if the responses $Y_{i}$ and the covariates $X_{i}$ are
missing simultaneously for all data $Z_{i}=(Y_{i},X_{i})$ with $D_{i}=0$.
However, it may be subject to loss of efficiency when all covariates are
observed. To improve efficiency in this case, we may incorporate estimating
equations such as $\sum_{i=1}^{N}p_{i}\\{h(X_{i})-\bar{h}\\}=0$ in the
definition of the empirical log-likelihood (3), where $h$ is a user-specific
function and $\bar{h}=(1/N)\sum_{i=1}^{N}h(X_{i})$. This strategy was also
adopted by Qin et al. (2008) to improve efficiency and by Han (2013) and Han
and Wang (2014) to construct multiple robust estimators. The price paid for
this modification, however, is that the resulting ELW weights are undefined if
$\bar{h}$ lies outside the convex hull of
$\\{h(X_{i}):D_{i}=1,i=1,2,\ldots,N\\}$ (Tsao, 2004; Chen et al., 2008; Liu
and Chen, 2010). The probability that the ELW weights are undefined can be
large when the sample size is small and/or the dimension of $h$ is high (Tsao,
2004); the calculational burden of the ELW weights will also be heavier. As
the focus of this paper is on the development of an estimation method that is
robust to small probabilities, we do not incorporate side information in the
proposed ELW method for the time being.

###  2.2 Practical implementation

The key to calculating the proposed EL estimators, including the EL estimator
$\hat{F}$ of $F$ and the ELW estimator $\hat{\theta}_{{\rm ELW}}$, is to
calculate $\hat{\alpha}$ by maximizing $\ell(\alpha)$. This necessitates a
double iterative algorithm because $\ell(\alpha)$ involves an implicit
function $\lambda(\alpha)$, and thus it seems to be rather a difficult task.
We find a more convenient solution, in which we need only solve a univariate
equation.

Mathematically, $\hat{\alpha}=\arg\max\ell(\alpha)$ is a solution to

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\left[\frac{D_{i}\lambda}{1+\lambda\\{\pi(Z_{i})-\alpha\\}}-\frac{1-D_{i}}{1-\alpha}\right].$
(9)

Combining (6) and (9) gives

$\displaystyle\lambda=\frac{N-n}{n(1-\alpha)}.$ (10)

Putting this expression into (6) leads to an equivalent equation for $\alpha$:

$\displaystyle\sum_{i=1}^{N}\frac{D_{i}(\pi(Z_{i})-\alpha)}{n/N+(1-n/N)\pi(Z_{i})-\alpha}$
$\displaystyle=$ $\displaystyle 0.$ (11)

As (11) has multiple roots, it is necessary to identify the interval
containing the desired root. Denote the observed $Z_{i}$ by
$z_{1},\ldots,z_{n}$ and define $\xi_{i}=n/N+(1-n/N)\pi(z_{i})$ for
$i=1,2,\ldots,n$. Equation (11) is further equivalent to $K(\alpha)=0$, where
$K(\alpha)=\sum_{i=1}^{n}\\{\pi(z_{i})-\alpha\\}/(\xi_{i}-\alpha)$. Because
$\alpha\in(0,1)$, $\xi_{i}\geq n/N$, and $n/N$ is a consistent estimator of
$\alpha$, the desired root of $K(\alpha)=0$ should lie between $0$ and
$\min\xi_{i}$. Actually, there must exist one and only one solution to
$K(\alpha)=0$ between $0$ and $\min\xi_{i}$. Because $\xi_{i}\geq\pi(z_{i})$,
it follows that $K\\{\min\pi(z_{i})\\}\geq 0$,
$\lim_{\alpha\uparrow\min\xi_{i}}K(\alpha)=-\infty$, and that $K(\alpha)$ is
strictly decreasing between $0$ and $\min\xi_{i}$. By the intermediate value
theorem, there must exist one and only one solution, denoted by
$\hat{\alpha}$, in $[\min\pi(z_{i}),\min\xi_{i})$ such that
$K(\hat{\alpha})=0$. It is worth noting that if all the $\pi(z_{i})$ are equal
and equal to $\alpha_{0}$, then $\hat{\alpha}=\alpha_{0}$ and the resulting
$\hat{p}_{i}$ are all equal to $1/n$, and the ELW estimator reduces to the
sample mean, i.e., $\hat{\theta}_{{\rm ELW}}=(1/n)\sum_{i=1}^{n}g(z_{i})$.
Otherwise, all $\pi(z_{i})$ ($i=1,2,\ldots,n$) are not equal to each other,
and $\hat{\alpha}$, $\hat{p}_{i}$, and $\hat{\theta}_{{\rm ELW}}$ are all non-
degenerate.

The proposed ELW estimation procedure can be implemented by Algorithm 1. The
second step can be efficiently achieved by a bi-section search algorithm, and
the remaining steps all involve only closed-form calculations, making the ELW
procedure extremely easy to implement.

Algorithm 1 ELW estimation procedure

0: The missing dataset
$\\{(D_{i},D_{i}g(Z_{i}),\pi(Z_{i})):i=1,2,\ldots,N\\}$.

0: The ELW estimate, $\hat{\theta}_{{\rm ELW}}$, of
$\theta={\mathbb{E}}\\{g(Z)\\}$.

Step 1. Calculate $n=\sum_{i=1}^{N}D_{i}$,
$\zeta_{l}=\min\\{\pi(Z_{i}):D_{i}=1,i=1,2,\ldots,N\\}$ and
$\zeta_{u}=n/N+(1-n/N)\zeta_{l}$.

Step 2. Calculate $\hat{\alpha}$ by solving (11) in the interval
$[\zeta_{l},\zeta_{u})$, and calculate
$\lambda(\hat{\alpha})=(N-n)/\\{n(1-\hat{\alpha})\\}$.

Step 3. Calculate
$\hat{p}_{i}=D_{i}n^{-1}[1+\lambda(\hat{\alpha})\\{\pi(Z_{i})-\hat{\alpha}\\}]^{-1}$
for $i=1,2,\ldots,N$.

Step 4. Calculate $\hat{\theta}_{{\rm
ELW}}=\sum_{i=1}^{N}\hat{p}_{i}D_{i}g(Z_{i})$.

###  2.3 Finite- and large-sample properties

The non-zero EL weights are $(1-\hat{\alpha})/\\{N(\xi_{i}-\hat{\alpha})\\}$
for $1\leq i\leq n$. We use the maximum weight ratio $\kappa=(\max_{1\leq
i\leq n}\xi_{i}-\hat{\alpha})/(\min_{1\leq i\leq n}\xi_{i}-\hat{\alpha})$
among the non-zero EL weights to quantify the dispersion between the EL
weights. The following proposition establishes an upper bound on $\kappa$.

######  Proposition 1.

Suppose $\pi(z_{i})$ $(1\leq i\leq n)$ take $m$ $\geq 2$ distinct values
$\pi_{(1)}<\ldots<\pi_{(m)}$ $(m\geq 2)$. If there exists
$\varepsilon\in(0,1)$ such that $\pi_{(m)}-\pi_{(1)}>\varepsilon$ and
$n_{1}/n<1-\varepsilon,$ then $\kappa\leq n/\varepsilon^{3}$.

Proposition 1 indicates that the ELW method works even if the smallest
$\pi(z_{i})$ is as small as zero. However, the maximum weight ratio of the IPW
estimator, namely $\max_{1\leq i\leq n}\pi(z_{i})/\min_{1\leq i\leq
n}\pi(z_{i})$, has no such a guarantee, and the IPW estimator becomes
extremely unstable when some of the $\pi(z_{i})$ are close to zero. In
particular, it fails to work when $\min_{1\leq i\leq n}\pi(z_{i})$ is exactly
zero. Our ELW estimator successfully and completely overcomes this issue,
which is its most significant advantage over the traditional IPW in finite-
sample performance.

Next, we show that asymptotically our ELW estimator is unbiased and more
efficient than the IPW estimator. This is a bonus of using ELW, and also a
significant advantage that it has over the conventional IPW in large-sample
performance. For ease of presentation, let
$B_{11}={\mathbb{E}}\\{1/\pi(Z)\\}$, and, for two generic functions $f$ and
$h$, define
$B_{fh}={\mathbb{E}}\\{f(Z)h^{\mathrm{\scriptscriptstyle\top}}(Z)/\pi(Z)\\}$,
$B_{f1}={\mathbb{E}}\\{f(Z)/\pi(Z)\\}$, and
$G=(B_{g1}-\theta_{0})/(B_{11}-1)$. We denote $A^{\otimes
2}=AA^{\mathrm{\scriptscriptstyle\top}}$ for a vector or matrix $A$.

######  Theorem 1.

Let $\alpha_{0}\in(0,1)$ be the truth of $\alpha$. Suppose that
${\mathbb{V}\rm ar}\\{\pi(Z)|D=1\\}>0$ and that $B_{11}$ and $B_{gg}$ are both
finite. As $N$ goes to infinity,

  1. (a)

$\sqrt{N}(\hat{\theta}_{{\rm
ELW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm ELW}}),$
where

$\displaystyle\Sigma_{{\rm ELW}}$ $\displaystyle=$
$\displaystyle{\mathbb{V}\rm
ar}\left\\{\left(\frac{Dg(Z)}{\pi(Z)}-\theta_{0}\right)-G\frac{D-\pi(Z)}{\pi(Z)}\right\\}=(B_{gg}-\theta_{0}^{\otimes
2})-G^{\otimes 2}(B_{11}-1);$

  2. (b)

$\sqrt{N}(\hat{\theta}_{{\rm
IPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm IPW}}),$
where

$\displaystyle\Sigma_{{\rm IPW}}={\mathbb{V}\rm
ar}\left\\{\frac{Dg(Z)}{\pi(Z,\beta_{0})}-\theta_{0}\right\\}=B_{gg}-\theta_{0}^{\otimes
2};$

  3. (c)

$\sqrt{N}(\hat{\theta}_{{\rm
SIPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm SIPW}}),$
where

$\displaystyle\Sigma_{{\rm SIPW}}$ $\displaystyle=$
$\displaystyle{\mathbb{V}\rm
ar}\left\\{\left(\frac{Dg(Z)}{\pi(Z,\beta_{0})}-\theta_{0}\right)-\theta_{0}\left(\frac{D}{\pi(Z,\theta_{0})}-1\right)\right\\}$
$\displaystyle=$ $\displaystyle(B_{gg}-\theta_{0}^{\otimes
2})-\theta_{0}(B_{g1}-\theta_{0})^{\mathrm{\scriptscriptstyle\top}}-(B_{g1}-\theta_{0})\theta_{0}^{\mathrm{\scriptscriptstyle\top}}+\theta_{0}^{\otimes
2}(B_{11}-1);$

  4. (d)

the ELW estimator $\hat{\theta}_{{\rm ELW}}$ is more efficient than the IPW
estimator $\hat{\theta}_{{\rm IPW}}$ and the SIPW estimator
$\hat{\theta}_{{\rm SIPW}},$ i.e. $\Sigma_{{\rm ELW}}\leq\Sigma_{{\rm IPW}}$
and $\Sigma_{{\rm ELW}}\leq\Sigma_{{\rm SIPW}},$ where the equalities hold
only if $\pi(Z)$ is degenerate.

The assumption ${\mathbb{V}\rm ar}\\{\pi(Z)|D=1\\}>0$ guarantees that with
probability tending to 1, the propensity scores are not all equal to each
other, and so the ELW estimator is non-degenerate. A reasonable estimator of
$\Sigma_{{\rm ELW}}$ is required in the construction of Wald-type confidence
intervals for $\theta$. Inspired by the fact that $\hat{p}_{i}\approx
D_{i}/\\{N\pi(Z_{i})\\}$, we propose to estimate $\Sigma_{{\rm ELW}}$ with the
ELW method by

$\displaystyle\widehat{\Sigma}_{{\rm ELW}}=(\hat{B}_{gg}-\hat{\theta}_{{\rm
ELW}}^{\otimes 2})-(\hat{B}_{g1}-\hat{\theta}_{{\rm ELW}})^{\otimes
2}/(\hat{B}_{11}-1),$ (12)

where $\hat{B}_{11}=N\sum_{i=1}^{N}(\hat{p}_{i})^{2}$,
$\hat{B}_{g1}=N\sum_{i=1}^{N}g(Z_{i})(\hat{p}_{i})^{2}$, and
$\hat{B}_{gg}=N\sum_{i=1}^{N}\\{g(Z_{i})\\}^{\otimes 2}(\hat{p}_{i})^{2}.$ It
is worth stressing that the ELW-based variance estimator is again insensitive
to small probabilities, since it circumvents the use of inverse probabilities.

###  2.4 Estimated propensity score

In many situations, such as missing data problems and causal inference, the
propensity score is unknown and needs be estimated from the observed data. The
ELW and IPW estimators have different large-sample behaviours if we take the
variability of the estimated propensity score into account. Suppose that
$\pi(\cdot)$ is parametrically modelled by $\pi(Z,\beta)$ with true parameter
value $\beta_{0}$.

######  Condition 1.

The following conditions are satisfied:

  1. (a)

The function $\pi(Z,\beta)$ is continuous in $\beta$ and
$\pi(Z,\beta_{0})=\pi(Z)$ for all $Z$.

  2. (b)

There exists a constant $\varepsilon>0$ such that (1)
${\mathbb{E}}\\{\sup_{\beta:|\beta-\beta_{0}|<\varepsilon}\|\nabla_{\beta}\pi(Z,\beta)\|\\}<\infty$
and (2)
${\mathbb{E}}\\{\sup_{\beta:|\beta-\beta_{0}|<\varepsilon}1/\pi(Z,\beta)\\}<\infty$.

  3. (c)

The truth $\alpha_{0}$ of $\alpha$ satisfies $0<\alpha_{0}<1$ and
${\mathbb{V}\rm ar}\\{\pi(Z,\beta_{0})|D=1\\}>0$.

Condition 1(a) holds when the non-missingness indicator $D$ follows logistic
and probit models, which are commonly used in the literature. Condition 1(b)
together with the other conditions guarantee the consistency of $\hat{\alpha}$
and hence the consistency of the ELW estimator. Under Condition 1(c), with
probability tending to 1, the observed propensity scores are not all equal to
each other, and so the ELW estimator is non-degenerate.

Under Condition 1(a), $\pi(Z,\beta_{0})=\pi(Z)$ and therefore
$B_{11}={\mathbb{E}}\\{1/\pi(Z,\beta_{0})\\}$,
$B_{g1}={\mathbb{E}}\\{g(Z)/\pi(Z,\beta_{0})\\}$, and
$B_{gg}={\mathbb{E}}[\\{g(Z)\\}^{\otimes 2}/\pi(Z,\beta_{0})].$ We still use
$G=(B_{g1}-\theta_{0})/(B_{11}-1)$. In addition, we denote
$B_{\dot{\pi}1}={\mathbb{E}}\\{\nabla_{\beta}\pi(Z,\beta_{0})/\pi(Z,\beta_{0})\\}$
and
$B_{g\dot{\pi}}={\mathbb{E}}[g(Z)\\{\nabla_{\beta}\pi(Z,\beta_{0})\\}^{\mathrm{\scriptscriptstyle\top}}/\pi(Z,\beta_{0})].$

######  Theorem 2.

Assume Condition 1 and that $\hat{\beta}$ satisfies
$\hat{\beta}-\beta_{0}=N^{-1}\sum_{i=1}^{N}h(D_{i},Z_{i})+o_{p}(N^{-1/2}),$
where the influence function $h(D,Z)$ has zero mean. Suppose that
$B_{11},B_{gg},B_{hh},$ and $B_{\dot{\pi}\dot{\pi}}$ are all finite. As $N$
goes to infinity,

  1. (a)

$\sqrt{N}(\hat{\theta}_{{\rm
ELW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm ELW},e}),$
where

$\displaystyle\Sigma_{{\rm ELW},e}={\mathbb{V}\rm
ar}\left\\{\left(\frac{Dg(Z)}{\pi(Z,\beta_{0})}-\theta_{0}\right)-G\frac{D-\pi(Z,\beta_{0})}{\pi(Z,\beta_{0})}+(GB_{1\dot{\pi}}-B_{g\dot{\pi}})h(D,Z)\right\\};$

  2. (b)

$\sqrt{N}(\hat{\theta}_{{\rm
IPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm IPW},e}),$
where

$\displaystyle\Sigma_{{\rm IPW},e}={\mathbb{V}\rm
ar}\left\\{\left(\frac{Dg(Z)}{\pi(Z,\beta_{0})}-\theta_{0}\right)-B_{g\dot{\pi}}h(D,Z)\right\\};$

  3. (c)

$\sqrt{N}(\hat{\theta}_{{\rm
SIPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm SIPW},e}),$
where

$\displaystyle\Sigma_{{\rm SIPW},e}={\mathbb{V}\rm
ar}\left\\{\frac{D_{i}g(Z_{i})-D_{i}\theta_{0}}{\pi(Z_{i},\beta_{0})}+(\theta_{0}B_{1\dot{\pi}}-B_{g\dot{\pi}})h(D_{i},Z_{i})\right\\}.$

According to Theorem 2, the ELW, IPW and SIPW estimators still follow
asymptotic normal distributions when the propensity score involves a finite-
dimensional unknown parameter. However, in general, the efficiency gain of the
ELW estimator over the IPW and SIPW estimators is no longer guaranteed.

The most popular choice for $\hat{\beta}$ is the maximum likelihood estimator.
For example, if data are missing at random (Rubin, 1974) and the covariates
$X_{i}$ are all observed, the maximum likelihood estimator $\hat{\beta}$ is
the maximizer of
$\sum_{i=1}^{N}[D_{i}\log\pi(Z_{i},\beta)+(1-D_{i})\pi\\{1-\pi(Z_{i},\beta)\\}]$.
In this case,

$\displaystyle h(D,Z)$ $\displaystyle=$
$\displaystyle\frac{D-\pi(Z,\beta_{0})}{\pi(Z,\beta_{0})\\{1-\pi(Z,\beta_{0})\\}}(\tilde{B}_{\dot{\pi}\dot{\pi}})^{-1}\nabla_{\beta}\pi(Z,\beta_{0}),$

where
$\tilde{B}_{\dot{\pi}\dot{\pi}}={\mathbb{E}}[\\{\nabla_{\beta}\pi(Z,\beta_{0})\\}^{\otimes
2}/\\{\pi(Z,\beta_{0})(1-\pi(Z,\beta_{0}))\\}].$ The asymptotic variance of
the ELW estimator is

$\displaystyle\Sigma_{{\rm ELW},e}$ $\displaystyle=$ $\displaystyle
B_{gg}-\theta_{0}^{\otimes
2}-G(B_{11}-1)G^{\mathrm{\scriptscriptstyle\top}}-(GB_{1\dot{\pi}}-B_{g\dot{\pi}})(\tilde{B}_{\dot{\pi}\dot{\pi}})^{-1}(GB_{1\dot{\pi}}-B_{g\dot{\pi}})^{\mathrm{\scriptscriptstyle\top}}.$

Again, an ELW estimator can be constructed for $\Sigma_{{\rm ELW},e}$.

###  2.5 Resampling-based interval estimation

Based on Theorems 1 and 2, we can construct Wald-type confidence intervals for
$\theta$ once a consistent estimator for the asymptotic variance is available.
The asymptotic normality of the ELW, IPW, and SIPW estimators requires
$B_{11}={\mathbb{E}}[\\{\pi(Z_{i})\\}^{-1}]<\infty$ and $B_{gg}<\infty$. If
this is violated, the Wald-type confidence intervals may not have the promised
coverage probability. This dilemma can be overcome by resampling. We propose
to construct confidence intervals for $\theta$ by the resampling method in
Algorithm 2.

Algorithm 2 Wald confidence region based on resampling and ELW

0: The missing dataset
$\\{(D_{i},D_{i}g(Z_{i}),\pi(Z_{i})):i=1,2,\ldots,N\\}$. Calculate the ELW
estimator $\hat{\theta}_{{\rm ELW}}$ and the proposed variance estimator
$\widehat{\Sigma}_{{\rm ELW}}$, and define
$T_{N}=\sqrt{N}(\widehat{\Sigma}_{{\rm ELW}})^{-1/2}(\hat{\theta}_{{\rm
ELW}}-\theta_{0})$.

0: Wald confidence region for $\theta={\mathbb{E}}\\{g(Z)\\}$ based on
resampling and ELW

Step 1. Draw $M\ll N$ (e.g. $M=\sqrt{N}$) observations, say
$(D_{i}^{*},D_{i}^{*}g(Z_{i}^{*}),\pi(Z_{i}^{*}))$ ($1\leq i\leq M$), from the
original sample by simple random sampling without replacement.

Step 2. Calculate the counterparts of $\hat{\theta}_{{\rm IPW}}$ and
$\widehat{\Sigma}_{{\rm ELW}}$ based on the subsample, denoted by
$\hat{\theta}_{{\rm ELW}}^{*}$ and $\widehat{\Sigma}_{{\rm ELW}}^{*}$ .
Construct $T_{M}^{*}=\sqrt{M}(\widehat{\Sigma}_{{\rm
ELW}}^{*})^{-1/2}(\hat{\theta}_{{\rm ELW}}^{*}-\hat{\theta}_{{\rm ELW}})$.

Step 3. Repeat Steps 1 and 2 $B=1000$ times and denote the resulting test
statistics by $\\{T_{M,i}^{*}:i=1,2,\ldots,B\\}$. Let
$t_{i}^{*}=\|T_{M,i}^{*}-\bar{T}^{*}\|$, where
$\bar{T}^{*}=(1/B)\sum_{i=1}^{B}T_{M,i}^{*}$. Denote the $(1-\alpha)$
empirical quantile of the $t_{i}^{*}$ by $q_{1-\alpha}^{*}$. Then a
$(1-\alpha)$-level confidence region for $\theta$ can be constructed as
$\\{\theta:\|\sqrt{N}(\widehat{\Sigma}_{{\rm ELW}})^{-1/2}(\hat{\theta}_{{\rm
ELW}}-\theta)-\bar{T}^{*}\|\leq q_{1-\alpha}^{*}\\}.$

In the case of the estimated propensity score $\widehat{\pi}(Z_{i})$, we
replace $\pi(Z_{i})$ and $\widehat{\Sigma}_{{\rm ELW}}$ by
$\widehat{\pi}(Z_{i})$ and $\widehat{\Sigma}_{{\rm ELW},e}$, respectively. The
ELW variance estimator $\widehat{\Sigma}_{{\rm ELW}}$ converges in probability
to ${\Sigma}_{{\rm ELW}}$, which is assumed to be positive definite. This,
together with Theorems 1 and 2, implies that $T_{N}$ converges in distribution
to the standard normal, an obviously continuous distribution. By Corollary 2.1
of Politis and Romano (1994), the empirical distribution of $T_{M}^{*}$ is a
uniformly consistent estimator of the distribution of $T_{N}$, which is
formally summarized in Theorem 3. This validates the interval estimator
produced by Algorithm 2.

######  Theorem 3.

Assume the conditions in Theorem 1 (for a known propensity score) or those in
Theorem 2 (for an estimated propensity score) are satisfied. As
$N\rightarrow\infty,$ if $M\rightarrow\infty$ and $M/N\rightarrow 0,$ then
$\sup_{t\geq 0}|P(T_{N}\leq t)-P^{*}(T_{M}^{*}\leq t)|=o_{p}(1),$ where
$P^{*}$ is the conditional probability given the original sample.

##  3 Extension to unequal probability samplings

We have shown that the ELW method works for missing data problems, where data
are independent and identically distributed. The results in Theorem 1 also
hold for Poisson sampling, a special case of UPS-WOR. Under Poisson sampling,
the indicators $D_{i}$ are all independent of each other and
$n=\sum_{i=1}^{N}D_{i}$ is random. Suppose there exists a function
$\pi(\cdot)$ such that $P(D_{i}=1)=\pi(Z_{i})$; then
$\\{(D_{i},D_{i}Z_{i}),i=1,2,\ldots,N\\}$ can be regarded as missing data and
the ELW method applies. However, Theorem 1 no longer holds under more general
UPS-WORs, where the sample size is fixed and the data are correlated, or UPS-
WRs, where a piece of data may be sampled multiple times. In this section, we
extend the ELW method to UPS-WOR and UPS-WR by solving Problems B and C.

###  3.1 ELW for UPS-WOR

With the notation of Problem B, denote the inclusion probability of $Z_{i}$ by
$\pi_{i}$. The population size $N$ and all the $\pi_{i}$ are known a priori.
The IPW estimator of $\theta$ is the famous Horvitz–Thompson estimator
(Horvitz and Thompson, 1952)

$\displaystyle\hat{\theta}_{\rm
IPW}=\frac{1}{N}\sum_{i=1}^{N}\frac{D_{i}}{\pi_{i}}g(Z_{i}).$ (13)

With $\pi_{i}$ in place of $\pi(Z_{i})$, we propose to directly apply the ELW
method from Section 2.1 and estimate $\theta_{0}=(1/N)\sum_{i=1}^{N}g(Z_{i})$
by $\hat{\theta}_{{\rm ELW}}$ in (8). Although the ELW estimators for UPS-WOR
and for missing data problems have the same form, their random behaviours are
totally differentm because the joint distributions of the observed data are
different.

The asymptotic normality of the IPW estimator, although difficult, has been
established for many commonly used samplings, including simple random sampling
with or without replacement (Erdös and Rényi, 1959; Hájek, 1960), rejective
sampling with unequal probabilities (H$\rm\acute{a}$jek, 1964), stratified
unequal probability sampling with or without replacement (Krewski and Rao,
1981; Bickel and Freedman, 1984), and two-phase sampling (Chen and Rao, 2007).
See Wu and Thompson (2020) for a comprehensive review.

We establish the asymptotic normality of the ELW estimator following the
results of Patterson et al. (2001) and Brändén and Jonasson (2012). A crucial
requirement in these papers for the series of random elements is linear
negative dependence or negative association, which are defined formally as
follows.

######  Definition 1 (Patterson et al. (2001)).

A sequence of random variables, $\\{X_{i}\\},$ is said to be linearly
negatively dependent (LIND) if for any disjoint subsets of indices $A,B$ and
positive constants $\lambda_{j},$

$\displaystyle P\left(\sum_{k\in A}\lambda_{k}X_{k}\leq s,\sum_{j\in
B}\lambda_{j}X_{j}\leq t\right)\leq P\left(\sum_{k\in A}\lambda_{k}X_{k}\leq
s\right)\cdot P\left(\sum_{j\in B}\lambda_{j}X_{j}\leq t\right)$

for any real numbers $s$ and $t$.

######  Definition 2 (Brändén and Jonasson (2012)).

A family $X=(X_{1},\ldots,X_{N})$ of random variables is said to be negatively
associated (NA) if
${\mathbb{E}}\\{f(X)g(X)\\}\leq{\mathbb{E}}\\{f(X)\\}{\mathbb{E}}\\{g(X)\\}$
for all increasing functions $f$ and $g$ such that there exists
$I\subset\\{1,2,\ldots,N\\}$ such that $f$ depends only on $\\{X_{i}:i\in
I\\}$ and $g$ depends only on $\\{X_{j}:i\neq I\\}$.

Patterson et al. (2001) proved that if a series of random variables is NA,
then it must be LIND, and that for LIND random series, asymptotic normality of
the IPW estimator holds under mild additional conditions. Brändén and Jonasson
(2012) showed that many commonly used $\pi ps$ samplings, including
conditional Poisson sampling, Sampford sampling, Pareto sampling, and pivotal
sampling satisfy the NA property. In summary, these results imply that, under
mild additional conditions, the IPW estimator for these commonly used $\pi ps$
samplings without replacement follows an asymptotic normal distribution.

Suppose that there is a sequence of large datasets, $\\{Z_{\nu
1},\ldots,Z_{\nu N_{\nu}}\\}$, of size $N_{\nu}$, $\nu=1,2,\ldots,$ all of
which are regarded as finite sets of non-random numbers. The large dataset
under study is one of them. For a particular large dataset $\\{Z_{\nu
1},\ldots,Z_{\nu N_{\nu}}\\}$, a sample of size $n_{\nu}$ is taken from it
with a prespecified UPS-WOR. Define $D_{\nu k}=1$ if $Z_{\nu k}$ is selected,
and 0 otherwise for $1\leq k\leq N_{\nu}$. Then,
$n_{\nu}=\sum_{k=1}^{N_{\nu}}D_{\nu k}$. Suppose that
$N_{\nu}\rightarrow\infty$ and $\alpha_{\nu 0}\equiv
n_{\nu}/N_{\nu}\rightarrow\alpha_{0}\in(0,1)$, as $\nu\rightarrow\infty$.

######  Theorem 4.

Let $\\{Z_{\nu k},1\leq k\leq N_{\nu},\nu\geq 1\\}$ be an array of
$q$-dimensional constant vectors and $\\{D_{\nu k},1\leq k\leq N_{\nu},\nu\geq
1\\}$ be an array of row-wise LIND random variables. Let $\pi_{\nu
k}={\mathbb{E}}(D_{\nu k})$ and $W_{\nu
k}=(g^{\mathrm{\scriptscriptstyle\top}}(Z_{\nu
k}),1)^{\mathrm{\scriptscriptstyle\top}}/\pi_{\nu k},$ and define
$S_{\nu}^{2}=N_{\nu}^{-1}{\mathbb{V}\rm ar}\\{\sum_{k=1}^{N_{\nu}}W_{\nu
k}(D_{\nu k}-\pi_{\nu k})\\}.$ Suppose that

$\displaystyle(1)$
$\displaystyle\lim_{\nu\rightarrow\infty}S_{\nu}^{2}=V_{0}\quad\mbox{ is a
positive definite matrix},$ $\displaystyle(2)$
$\displaystyle\lim_{\nu\rightarrow\infty}N_{\nu}^{-1}\sum_{1\leq k\neq l\leq
N_{\nu}}\|W_{\nu k}\|\|W_{\nu l}\|{\mathbb{C}\rm ov}(D_{\nu k},D_{\nu
l})=0,\quad\mbox{and}$ $\displaystyle(3)$
$\displaystyle\lim_{\nu\rightarrow\infty}N_{\nu}^{-1}\sum_{k=1}^{N_{\nu}}\|W_{\nu
k}\|^{2}P(\|W_{\nu k}\||D_{\nu k}-\pi_{\nu
k}|\geq\varepsilon\sqrt{N_{\nu}})=0\ \mbox{for any}\ \varepsilon>0.$

Then, as $\nu\rightarrow\infty,$
$S_{\nu}^{-1}N_{\nu}^{-1/2}\sum_{k=1}^{N_{\nu}}W_{\nu k}(D_{\nu k}-\pi_{\nu
k}){\overset{d}{\longrightarrow\;}}N(0,I_{q+1}).$

For notational simplicity, we shall suppress the subscript $\nu$ and use
$N\rightarrow\infty$ instead of $\nu\rightarrow\infty$. This theorem implies
that
$N^{-1/2}\sum_{k=1}^{N}W_{k}(D_{k}-\pi_{k}){\overset{d}{\longrightarrow\;}}N(0,V_{0})$.
This result is sufficient to establish the asymptotic normalities of the IPW
and ELW estimators. Define
$B_{2}=\lim_{N\rightarrow\infty}(1/N)\sum_{k=1}^{N}\\{g(Z_{k})\\}^{\otimes
2}$. We redefine the notation from Section 2. For two generic functions $f$
and $h$, define
$B_{fh}=\lim_{N\rightarrow\infty}N^{-1}\sum_{k=1}^{N}f(Z_{k})h^{\mathrm{\scriptscriptstyle\top}}(Z_{k})/\pi_{k}$
and write $B_{11}=\lim_{N\rightarrow\infty}N^{-1}\sum_{k=1}^{N}1/\pi_{k}$ and
$B_{f1}=\lim_{N\rightarrow\infty}N^{-1}\sum_{k=1}^{N}f(Z_{k})/\pi_{k}$.

######  Theorem 5.

Under the conditions in Theorem 4, as $\nu$ goes to infinity,

  1. (a)

$N^{1/2}(\hat{\theta}_{{\rm
IPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm IPW}}),$
where $\Sigma_{{\rm IPW}}=B_{gg}-B_{2};$

  2. (b)

$\sqrt{N_{\nu}}(\hat{\theta}_{{\rm
SIPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm SIPW}}),$
where

$\Sigma_{{\rm SIPW}}=\Sigma_{{\rm IPW}}-\frac{(B_{g1}-\theta_{0})^{\otimes
2}}{B_{11}-1}+\frac{(B_{g1}-\theta_{0}B_{11})^{\otimes 2}}{B_{11}-1};$

  3. (c)

$\sqrt{N_{\nu}}(\hat{\theta}_{{\rm
ELW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm ELW}}),$
where $\Sigma_{{\rm ELW}}=\Sigma_{{\rm IPW}}-(B_{g1}-\theta_{0})^{\otimes
2}/(B_{11}-1);$

  4. (d)

the ELW estimator $\hat{\theta}_{{\rm ELW}}$ is more efficient than
$\hat{\theta}_{{\rm IPW}}$ and $\hat{\theta}_{{\rm SIPW}},$ i.e. $\Sigma_{{\rm
ELW}}\leq\Sigma_{{\rm IPW}}$ and $\Sigma_{{\rm ELW}}\leq\Sigma_{{\rm SIPW}},$
where the equalities hold only if the $\pi_{k}$ are all equal.

By Theorem 5, again the ELW, IPW, and SIPW estimators are asymptotically
normal, and the ELW estimator is asymptotically the most efficient. A variance
estimator similar to (12) can be constructed as $\widehat{\Sigma}_{{\rm
ELW}}=\hat{B}_{gg}-\hat{B}_{2}-(\hat{B}_{g1}-\hat{\theta}_{{\rm
ELW}})^{\otimes 2}/(\hat{B}_{11}-1),$ where
$\hat{B}_{fh}=N\sum_{k=1}^{N}(\hat{p}_{k})^{2}f(Z_{k})h^{\mathrm{\scriptscriptstyle\top}}(Z_{k})$
and $\hat{B}_{2}=\sum_{k=1}^{N}\hat{p}_{k}\\{g(Z_{k})\\}^{\otimes 2}.$

######  Remark 3.

Under UPS-WOR, the $Z_{\nu k}$ are non-random and only the $D_{\nu k}$ are
random variables. Unlike the cases in missing data problems, the ELW weights
$\hat{p}_{\nu k}$ do not lead to an estimate of a certain distribution
function. They are simply taken as weights for parameter estimation.

###  3.2 ELW for UPS-WR

We still assume that the finite population of interest
$\\{Z_{1},\ldots,Z_{N}\\}$, where the subscript $\nu$ is suppressed, lies in a
hypothetical infinite sequence (indexed by $\nu$) of finite populations of
increasing size. The parameter to be estimated is
$\theta_{0}=(1/N)\sum_{i=1}^{N}g(Z_{i})$. Suppose a sample of size $n$ is
drawn by an UPS-WR from $\\{Z_{1},\ldots,Z_{N}\\}$. Let the sampling
probability of the $i$th individual at each sampling occasion be $q_{i}^{*}$
($i=1,2,\ldots,N$). Clearly, $q_{i}^{*}>0$ and $\sum_{i=1}^{N}q_{i}^{*}=1$.
Denote the observed data by $z_{i}$ and the corresponding sampling probability
by $q_{i}$ for $i=1,2,\ldots,n$.

To apply the IPW and ELW methods from Section 2, we define
$\pi_{i}^{*}=nq_{i}^{*}$ ($1\leq i\leq N$) and $\pi_{i}=nq_{i}$ ($1\leq i\leq
n$), which may be greater than 1. The IPW estimator of $\theta_{0}$ under UPW-
WR is the Hansen–Hurwitz estimator (Hansen and Hurwitz, 1943),

$\displaystyle\hat{\theta}_{{\rm
IPW}}=\frac{1}{nN}\sum_{i=1}^{n}\frac{g(z_{i})}{q_{i}}=\frac{1}{N}\sum_{i=1}^{n}\frac{g(z_{i})}{\pi_{i}}.$
(14)

The ELW method from Section 2 was introduced with data
$\\{(D_{i},D_{i}Z_{i},\pi(Z_{i})):1\leq i\leq N\\}$. However, the notation
$D_{i}$ is not appropriate under UPS-WR, because a unit from the population
may be drawn multiple times. As the ELW estimator depends only on $N$ and
$\\{(z_{i},\pi_{i}):1\leq i\leq n\\}$, it is applicable to UPS-WR through a
modified version of Algorithm 1.

Algorithm 3 ELW estimation procedure under UPS-WR

0: The sample drawn, $\\{(g(z_{i}),\pi_{i}):i=1,2,\ldots,n\\}$; The size of
the finite population, $N$.

0: The ELW estimate, $\hat{\theta}_{{\rm ELW}}$, of
$\theta={\mathbb{E}}\\{g(Z)\\}$.

Step 1. Calculate $\xi_{i}=n/N+(1-n/N)\pi_{i}$ ($1\leq i\leq n$),
$\zeta_{l}=\min_{1\leq i\leq n}\pi_{i}$, and $\zeta_{u}=\min_{1\leq i\leq
n}\xi_{i}$.

Step 2. Calculate $\hat{\alpha}$ by solving
$0=\sum_{i=1}^{n}(\pi_{i}-\alpha)/(\xi_{i}-\alpha)$ in the interval
$[\zeta_{l},\zeta_{u})$, and calculate
$\lambda(\hat{\alpha})=(N-n)/\\{n(1-\hat{\alpha})\\}$.

Step 3. Calculate
$\hat{p}_{i}=n^{-1}\\{1+\lambda(\hat{\alpha})(\pi_{i}-\hat{\alpha})\\}^{-1},i=1,2,\ldots,n$.

Step 4. Calculate the ELW estimate $\hat{\theta}_{{\rm
ELW}}=\sum_{i=1}^{n}\hat{p}_{i}g(z_{i})$.

Define $\alpha_{0}=n/N$ and
$B_{fh}=N^{-1}\sum_{k=1}^{N}f(Z_{k})h^{\mathrm{\scriptscriptstyle\top}}(Z_{k})/\pi_{k}^{*}$
for two generic functions $f$ and $h$. The asymptotic normalities of the IPW,
SIPW, and ELW estimators are established in the following theorem.

######  Theorem 6.

Suppose that there exist $\alpha_{0}^{*}>0,$ $\theta_{0}^{*},$
$B_{11}^{*}>1/\alpha_{0}^{*},$ $B_{g1}^{*},$ and a positive-definite matrix
$B_{gg}^{*}$ such that as $\nu\rightarrow\infty,$ we have
$\alpha_{0}\rightarrow\alpha_{0}^{*},$ $\theta_{0}\rightarrow\theta_{0}^{*},$
$B_{1}\rightarrow B_{11}^{*},$ $B_{g1}\rightarrow B_{g1}^{*},$ and
$B_{gg}\rightarrow B_{gg}^{*}$. As $\nu\rightarrow\infty,$

  1. (a)

$\sqrt{n}(\hat{\theta}_{{\rm
IPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm IPW}})$ with
$\Sigma_{{\rm IPW}}=\alpha_{0}^{*}B_{gg}^{*}-(\theta_{0}^{*})^{\otimes 2};$

  2. (b)

$\sqrt{n}(\hat{\theta}_{{\rm
SIPW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm SIPW}}),$
with $\Sigma_{{\rm
SIPW}}=\alpha_{0}^{*}\\{B_{gg}^{*}-\theta_{0}^{*}B_{g1}^{*}-B_{g1}^{*}\theta_{0}^{*}+(\theta_{0}^{*})^{\otimes
2}B_{11}^{*}\\};$

  3. (c)

$\sqrt{n}(\hat{\theta}_{{\rm
ELW}}-\theta_{0}){\overset{d}{\longrightarrow\;}}N(0,\Sigma_{{\rm ELW}}),$
where

$\displaystyle\Sigma_{{\rm ELW}}$ $\displaystyle=$
$\displaystyle\alpha_{0}^{*}B_{gg}^{*}-(\theta_{0}^{*})^{\otimes
2}+\frac{(1-\alpha_{0}^{*})^{2}(\theta_{0}^{*}B_{11}^{*}-B_{g1}^{*})^{\otimes
2}}{(\alpha_{0}^{*}B_{11}^{*}-1)(B_{11}^{*}-1)^{2}}-\frac{(\alpha_{0}^{*}B_{g1}^{*}-\theta_{0}^{*})^{\otimes
2}}{\alpha_{0}^{*}B_{11}^{*}-1}.$

As $\hat{p}_{i}\approx 1/(N\pi_{i})$ ($1\leq i\leq n$), similarly to (12), we
propose to estimate $\Sigma_{{\rm ELW}}$ by

$\widehat{\Sigma}_{{\rm ELW}}=\alpha_{0}\hat{B}_{gg}-(\hat{\theta}_{{\rm
ELW}})^{\otimes 2}+\frac{(1-\alpha_{0})^{2}(\hat{\theta}_{{\rm
ELW}}\hat{B}_{11}-\hat{B}_{g1})^{\otimes
2}}{(\alpha_{0}\hat{B}_{11}-1)(\hat{B}_{11}-1)^{2}}-\frac{(\alpha_{0}\hat{B}_{g1}-\hat{\theta}_{{\rm
ELW}})^{\otimes 2}}{\alpha_{0}\hat{B}_{11}-1},$

where
$\hat{B}_{fh}=N\sum_{k=1}^{n}(\hat{p}_{k})^{2}f(z_{k})h^{\mathrm{\scriptscriptstyle\top}}(z_{k}).$

##  4 Simulation study

We conduct simulations to investigate the finite-sample performance of the
proposed ELW estimator and the accompanying asymptotic-normality-based
interval estimator. For comparison, we also take into account the IPW
estimator, the SIPW estimator, and some popular variants of the IPW estimator:

  1. 1.

The modified IPW estimator of Zong et al. (2019) (ZZZ for short):
$\hat{\theta}_{\rm ZZZ}=N^{-1}\sum_{i=1}^{N}D_{i}Y_{i}/\tilde{\pi}_{i},$ where
$\tilde{\pi}_{i}=\max\\{\pi_{(K)},\pi(X_{i})\\}$, $K$ is the maximum $i$ such
that $\pi_{(i)}\leq 1/(i+1)$, and $\\{\pi_{(1)},\ldots,\pi_{(N)}\\}$ are the
propensity scores in increasing order.

  2. 2.

The trimmed IPW estimator of Crump et al. (2009) (CHIM for short):

$\displaystyle\hat{\theta}_{\rm
CHIM}=\left.\sum_{i=1}^{N}\frac{D_{i}Y_{i}}{\pi(X_{i})}\cdot
I\\{\alpha\leq\pi(X_{i})\leq
1-\alpha\\}\right/\sum_{i=1}^{N}I\\{\alpha\leq\pi(X_{i})\leq 1-\alpha\\},$

where $\alpha$ is obtained by minimizing a variance term and $I(\cdot)$ is the
indicator function.

  3. 3.

The IPW estimator of Ma and Wang (2020) with $s=1$ and $s=2$, denoted by MW1
and MW2, respectively. Following Ma and Wang (2020), we set the tuning
parameters $b_{N}$ and $h_{N}$ in MW1 and MW2 to the respective solutions of
$b_{N}^{s}N^{-1}\sum_{i=1}^{N}I\\{\pi(X_{i})\leq b_{N}\\}=1/(2N)$ and
$h_{N}^{5}\sum_{i=1}^{N}I\\{\pi(X_{i})\leq h_{N}\\}=1.$ For details, see the
discussion below Theorem 3 of Ma and Wang (2020) and Section III of their
supplementary material.

We simulate data from two examples, which represent missing data and unequal
probability sampling problems, respectively. All numbers reported in this
simulation study are calculated based on $M=5000$ simulated random samples.

######  Example 1 (Missing data).

Instead of generating $X,$ we generate the propensity score $\pi(X)$ from
$P(\pi(X)\leq u)=u^{\gamma-1}$ $(0\leq u\leq 1)$ with $\gamma=1.5$ or $2.5$.
Given $\pi(X),$ we generate $Y$ from
$Y=\mu\\{\pi(X)\\}+c\cdot(\eta-4)/\sqrt{8},$ where $c=1$ or $0.1,$ and
$\eta\sim\chi_{4}^{2},$ and the missingness status $D$ of $Y$ follows the
Bernoulli distribution with success probability $\pi(X)$. Four choices of
$\mu(t)$ are considered: $\mu(t)=\cos(2\pi t)$ (Model 1), $\mu(t)=1-t$ (Model
2), $\mu(t)=\cos(2\pi t)+5$ (Model 3), and $\mu(t)=6-t$ (Model 4). The sample
size is $N=2000$ and the parameter of interest is $\theta={\mathbb{E}}(Y)$.

This example is a modified version of Example 1 in Section III of the
supplementary material of Ma and Wang (2020), who considered the cases with
$\gamma=1.5$, $c=1$, and $N=2000$ for Models 1 and 2. The parameter $\gamma$
($\gamma>1$) controls the tail behaviour of $1/\pi(X)$. When $\gamma>2$, the
tail is light and ${\mathbb{E}}\\{1/\pi(X)\\}=(\gamma-1)/(\gamma-2)$ is
finite. In this case, if $g$ is bounded, then the conditions in Theorem 1 are
generally fulfilled, and the asymptotic normalities of the ELW, IPW, and SIPW
estimators are guaranteed. However, in the case of $1<\gamma\leq 2$, the tail
is heavy and ${\mathbb{E}}\\{1/\pi(X)\\}=\infty$, which violates the
conditions of Theorem 1: the ELW, IPW, and SIPW estimators no longer follow
asymptotically normal distributions. The constant $c$ controls the influence
of the random error on the response variable; a smaller $c$ leads to a smaller
noise. Models 3 and 4 are simply Models 1 and 2 with a mean shift.

Point estimation  As a measure of the finite-sample performance of a generic
estimator $\tilde{\theta}$, we define its scaled root mean square error (RMSE)
as ${\rm
RMSE}(\tilde{\theta})=\sqrt{N}\times\\{(1/M)\sum_{j=1}^{M}(\tilde{\theta}_{j}-\theta)^{2}\\}^{1/2}$,
where $\tilde{\theta}_{j}$ is the estimate $\tilde{\theta}$ based on the $j$th
simulated random sample. Table 1 presents a comparison of the RMSEs of the
seven estimators, and Figures 1 and 2 display the corresponding boxplots for
$\gamma=1.5$ and 2.5, respectively. For clearer presentation, we ignore the
boxplots of the IPW estimator, because it fluctuates too dramatically.

In terms of RMSE, ELW outperforms IPW, SIPW, ZZZ, and CHIM in almost all
scenarios, which is clearly confirmed by their boxplots. The only exception is
the scenario with $\gamma=2.5,c=1$ for Model 1, where the RMSE (1.17) of ELW
is slightly greater than the minimum RMSE (1.14) of IPW, SIPW, ZZZ, and CHIM.
The boxplots also indicate that ELW is always nearly unbiased in all
scenarios. ELW also outperforms MW1 and MW2 in most cases. The only exceptions
are the scenarios with $\gamma=2.5$ for Model 1 and those with
$\gamma=1.5,c=1$ for Models 1 and 2. In the least favourable scenario
($\gamma=1.5,c=1$, Model 2), the RMSE of ELW is greater than those of MW1 and
MW2 by at most $(5.13-3.78)/3.78\approx 35.7\%$. By contrast, the RMSEs of MW1
and MW2 can be more than 12 times that of ELW; see the scenario with
$\gamma=1.5$ and $c=0.1$ for Model 4. Although MW1 and MW2 often have smaller
RMSEs than ELW for Model 1, the boxplots in Figures 1 and 2 indicate that they
tend to have either non-ignorable biases or larger variances.

Models 3 and 4 are simply Models 1 and 2 with a mean shift. When we change
Models 1 and 2 to Models 3 and 4, respectively, and keep the remaining
settings unchanged, the boxplots demonstrate that ELW clearly performs the
best: it not only is nearly unbiased, but also has the smallest variance.
Meanwhile, ELW, CHIM, and SIPW have nearly unchanged RMSEs. This makes sense,
because their weights all sum to 1. Unfortunately, IPW, ZZZ, MW1, and MW2 are
all very sensitive to a mean shift in the data generating process, since their
weights do not sum to 1.

When $c$ decreases from 1 to 0.1, the influence of random error become
negligible and we expect all methods to exhibit better performance. Indeed,
all methods have decreasing RMSEs, except for IPW. ELW has the largest rates
of decline in RMSE: these rates are at least 69% and 42% when $\gamma=1.5$ and
$2.5$, respectively. However, the RMSEs of ZZZ, MW1, and MW2 have nearly no
reduction for Models 3 and 4. ELW performs in the most stable manner, whereas
the other methods have either extremely large fluctuations or remarkable
biases.

When $\gamma$ increases from $1.5$ to 2.5, ELW clearly outperforms the
competitors in all scenarios except those for Model 1. All methods exhibit
similar performance for Models 1 and 2. However for Models 3 and 4, IPW, ZZZ,
MW1, and MW2 have much larger fluctuations, compared with their performance
for both Models 1 and 2. This indicates that they are sensitive to a mean
shift, which is undesirable.

Roughly speaking, among the seven estimators under comparison, the ELW
estimator is the most reliable in almost all scenarios. Both the RMSE results
and the boxplots indicate that MW1 and MW2 can exhibit very different
performances. In other words, the performance of the method of Ma and Wang
(2020) can be affected by the choice of underlying tuning parameters. We have
also conducted simulations for $N=50$ and 500, and $\gamma=1.3$ and 1.9. See
Section 8 of the supplementary material for the RMSE results and the
corresponding boxplots. The general findings are similar.

Interval estimation  Two confidence intervals for $\theta$ can be constructed
based on the ELW estimator $\hat{\theta}_{{\rm ELW}}$. One is the Wald
confidence interval (ELW-an for short) based on the asymptotic normality of
$\hat{\theta}_{{\rm ELW}}$, where the asymptotic variance is estimated using
the ELW method. The other is the resampling-based interval estimator (ELW-re)
given in Section 2.5. Similar intervals (IPW-an, IPW-re, SIPW-an, and SIPW-re)
can be constructed when the IPW and SIPW estimators take the place of the ELW
estimator in the estimations of both $\theta$ and the asymptotic variances. We
compare these six confidence intervals with those of Ma and Wang (2020) based
on their resampling method and the MW1 and MW2 point estimators, which are
denoted by MW1-re and MW2-re, respectively.

We generate random samples of size $N=2000$ from Example 1, and calculate the
coverage probabilities and average lengths of the eight confidence intervals
at the 95% confidence level. The results are tabulated in Table 2. MW2-re has
the most accurate coverage accuracy, followed by ELW-re when $\gamma=1.5$ and
$c=1.0$ for Models 1 and 2. When Models 1 and 2 are replaced by Models 3 and
4, the coverage probabilities of ELW-re remain nearly unchanged; however,
those for MW1-re and MW2-re decrease sharply by more than 5% and 10%,
respectively. When $c$ decreases from 1.0 to 0.1, the coverage accuracy of
ELW-re becomes better or is still acceptable, although both MW1-re and MW2-re
perform much more poorly. With different tuning parameters, MW1-re and MW2-re
often have quite different coverage probabilities and average lengths, which
again shows that the performance of the method of Ma and Wang (2020) can be
greatly affected by different choices of tuning parameters. SIPW-re has very
close coverage probabilities to ELW-re in most cases, whereas its average
lengths are generally much greater than those of the latter.

As expected, all asymptotic-normality-based Wald intervals exhibit remarkable
undercoverage when $\gamma=1.5$, because the asymptotic normalities are
generally violated. In the meantime, all resampling-based intervals have
improved performance. When $\gamma$ increases to 2.5, all intervals except
MW1-re and MW2-re have very desirable coverage accuracy, and the asymptotic-
normality-based intervals have close or even better coverage probabilities
compared with the resampling-based intervals.

######  Example 2 (Unequal probability sampling).

The finite population consists of $N=3000$ observations
$\\{(x_{i},y_{i}):1\leq i\leq N\\}$. We generate $x_{i}$ from the uniform
distribution on $[0,2],$ and $y_{i}=\mu(x_{i})+\sqrt{3(1-\rho^{2})}e_{i},$
where $e_{i}\sim N(0,1)$ are independent of each other and of $x_{i}$. We
consider four models for $\mu(x)$, namely $\sqrt{3}\rho x$ (Model 1),
$\sqrt{3}\rho(x+x^{2})$ (Model 2), $\sqrt{3}\rho x+5$ (Model 3), and
$\sqrt{3}\rho(x+x^{2})+5$ (Model 4), and three unequal probability samplings,
namely Poisson sampling, pivotal sampling, and PPS (probability proportional
to pize) sampling. We fix $n=500,$ take $\pi_{i}=nx_{i}/\sum_{j=1}^{N}x_{j}$
for Poisson sampling and pivotal sampling, and take $x_{i}$ as the size for
PPS sampling. The goal is to estimate $\theta=(1/N)\sum_{i=1}^{N}y_{i}$.

Example 2 is a modified version of Example 2 from Zong et al. (2019) and is
designed to represent unequal probability sampling. Poisson sampling and
pivotal sampling are UPS-WORs, with the sample size being random in the former
and fixed in the latter. PPS sampling is a popular UPS-WR. In particular,
pivotal sampling satisfies the NA property (Brändén and Jonasson, 2012) and
therefore must satisfy the LIND condition (Patterson et al., 2001) required by
Theorem 4.

We exclude CHIM, MW1, and MW2 from the subsequent comparison because they are
all designed not for unequal probability sampling but for missing data
problems. Table 4 presents the simulated RMSEs for the IPW, SIPW, ZZZ, and ELW
estimators, and Figure 3 displays the accompanying boxplots under pivotal
sampling. The boxplots under Poisson and PPS samplings are similar, and can be
found in the supplementary material.

ELW always has the smallest RMSEs under Poisson sampling. It still perform the
best under both pivotal and PPS samplings, except for the scenarios with
$\rho=0.8$ for Models 1 and 2. When the RMSE of ELW is not the smallest, the
boxplots of ZZZ and ELW are close to each other. In the remaining scenarios,
the competitors of ELW have either much larger variances or much larger biases
than ELW. Similar to Example 1, Models 3 and 4 in Example 2 are also Models 1
and 2 with a mean shift. Again, the performance of ELW and SIPW is equivariant
under such a mean shift; however, the IPW and ZZZ estimators are very
sensitive, and their performances become much worse when the mean of $y_{i}$
moves away from the origin. The performances of the four estimators seem to be
insensitive to the choice of $\rho$.

For interval estimation under unequal probability sampling, no valid
resampling-based confidence intervals are available as far as we know.
Therefore, we compare only IPW-an, SIPW-a, and ELW-an, which are the
confidence intervals based on the asymptotic normalities of the IPW, SIPW, and
ELW estimators. Table 3 presents their simulated coverage probabilities and
average lengths. Under both Poisson and pivotal samplings, overall ELW-an has
not only more accurate coverage probabilities but also shorter average lengths
than SIPW-an. Under PPS sampling, when $\rho=0.8$, ELW-an still wins, although
both ELW-an and SIPW-an have undercoverages of 3% or more. However, when
$\rho=0.2$, SIPW-an has better and more desirable coverage probabilities. This
implies that ELW-an still has room for improvement.

In summary, the ELW point estimator has the most reliable overall performance,
is shift-equivariant, and is nearly unbiased in all cases. The proposed
resampling-based ELW interval estimator often has desirable coverage accuracy
and short lengths in missing data problems, whether the proportion of
extremely small propensity scores is small or large. The asymptotic-normality-
based ELW interval estimator has desirable performance under UPW-WORs such as
Poisson and pivotal samplings, and acceptable performance under UPS-WRs such
as PPS sampling.

##  5 Real data analysis

LaLonde (1986) estimated the impact of the National Supported Work
Demonstration, a labour training programme, on post-intervention income
levels, using data from a randomized evaluation of the programme. To further
demonstrate the superiority of the proposed ELW method, we analyse the
LLvsPSID data from the R package cem, which is the Lalonde set of treated
units versus PSID (Panel Study of Income Dynamics) control individuals. The
data consist of 2787 observations (297 from treated units and 2490 from
control units) on 10 variables: treated (treatment indicator), age (age),
education (years of education), black (race, indicator variable), married
(marital status, indicator variable), nodegree (indicator variable of not
possessing a degree), re74 (real earnings in 1974), re75 (real earnings in
1975), re78 (real earnings in 1978), hispanic (ethnic, indicator variable),
u74 (unemployment in 1974, indicator variable), and u75 (unemployment in 1975,
indicator variable). The variable re78 is the post-treatment outcome.

Let $Y=$re78/10 000 be the response, let $D=$ treated, and let $Y(d)$ denote
the response of an individual whose treatment status is $D=d$. We shall not
address the original treatment effect estimation problem. Instead, we take the
data as missing data and wish to estimate the average earnings of the treated
in 1978. In other words, the parameter of interest is
$\theta={\mathbb{E}}\\{Y(1)\\}$. We first estimate the propensity scores by
fitting a linear logistic regression model of the treatment indicator $D$ on
the remaining eight variables (excluding $D$, $Y$, and re78). With the fitted
propensity scores, the IPW, SIPW, MW1, MW2, and ELW point estimates are 0.65,
0.92, 0.72, 0.70, and 1.11, respectively, and the corresponding resampling-
based interval estimates at the 95% level are $[-12.68,8.00]$, $[-4.16,3.15]$,
$[-3.04,1.26]$, $[-0.27,1.04]$, and $[-8.86,6.27]$, respectively. If we
replace all $Y$ by $Y+5$, the point estimates become 4.16, 5.92, 5.81, 5.56,
and 6.11 with interval estimates being $[-27.31,22.76]$, $[0.90,8.16]$,
$[0.92,7.16]$, $[5.05,6.22]$, and $[-3.87,11.44]$, respectively. As expected,
the SIPW and ELW point estimates are shift-equivariant, but the IPW estimator
and the MW estimators are not.

Figure 4 displays the fitted propensity scores of both the treated and control
groups. A clump of near-zero propensity scores in the treated group implies
that the standard IPW estimator is dramatically unstable. The excessive number
of near-zero propensity scores in both groups indicates that the distribution
of the inverse propensity score has a very heavy right tail similar to that in
the simulation scenario with $\gamma=1.5$ in Example 1. According to our
simulation experience in the case of $\gamma=1.5$, the ELW point estimator is
always unbiased or nearly unbiased, and its performance is the most stable in
most cases. By contrast, the other estimators SIPW, MW1, and MW2 may have
either much larger RMSEs or large biases. The ELW-re interval has the most
desirable and much better coverage accuracy than the other intervals. These
observations makes it reasonable to believe that the ELW point and interval
estimates, 1.11 and $[-8.86,6.27]$, are the most preferable for the estimation
of $\theta={\mathbb{E}}\\{Y(1)\\}$.

In addition, we regard the observations in the LLvsPSID data with non-zero
re75 as a finite population, and conduct Poisson, pivotal, and PPS samplings
with inclusion probabilities proportional to re75. We take the parameter of
interest to be the mean of $Y={\it re78}/10\,000+a$, with $a=0$ or 2. Table 5
presents the simulated RMSE results based on 5000 simulation repetitions with
a sample size (in pivotal and PPS samplings) or ideal sample size (in Poisson
sampling) of $200$. The ELW estimator has the smallest RMSEs under Poisson
sampling, regardless of whether $a=0$ or 2, and under pivotal and PPS
samplings when $a=2$. It also uniformly outperforms SIPW under all three
samplings. When $a=0$, its performance can be inferior to those of IPW and
ZZZ, which, however, are highly sensitive to a location shift in $Y$. The ELW
estimator again has the best overall performance under unequal probability
sampling.

##  6 Discussion

When introducing the ELW method, we assume that the parameter of interest has
a simple form $\theta={\mathbb{E}}\\{g(Z)\\}$, for convenience of
presentation. Our method can be directly extended to the parameters defined by
a group of estimation equations (Qin and Lawless, 1994). Much more effort
needs to be applied to the investigation of the large-sample properties of the
corresponding ELW estimator.

In the development of the ELW method, we regard $\alpha$ as unknown to
guarantee that the EL weights are always well defined. In some situations,
such as Poisson and PPS samplings, the true parameter value $\alpha_{0}$ of
$\alpha$ may be known a priori. Under PPS sampling, the coverage accuracy of
the ELW-based Wald interval is not quite satisfactory. Hopefully, replacing
$\alpha$ with $\alpha_{0}$ or $n/N$ may help, although the EL weights may be
undefined. This poses the interesting problem of how to overcome the non-
definition problem of the EL method while simultaneously improving its
estimation efficiency.

When constructing confidence regions for the parameter $\theta$, we consider
only the Wald method in this paper. In missing data problems or under Poisson
sampling, an alternative is to construct confidence intervals based on the
empirical likelihood ratio function, which necessitates study of the limiting
distribution of the empirical likelihood ratio statistic. We believe that this
statistic asymptotically follows a central chi-square distribution, making it
convenient for the construction of empirical likelihood ratio regions with
asymptotically correct confidence level.

## Supplementary material

The supplementary material contains proofs of Proposition 1, Theorems 1, 2, 4,
5, and 6, and more simulation results.

##  Acknowledgements

The research was supported by the National Natural Science Foundation of China
(12171157 and 11971300), the Natural Science Foundation of Shanghai
(19ZR1420900), the State Key Program of the National Natural Science
Foundation of China (71931004 and 32030063), the Development Fund for Shanghai
Talents, and the 111 Project (B14019).

## References

## References

  * Bang and Tsiatis (2000) Bang, H. and Tsiatis, A. A. (2000)  Estimating medical costs with censored data.  Biometrika, 87, 329–343.
  * Bickel and Freedman (1984) Bickel, P. J. and Freedman, D. A. (1984)  Asymptotic normality and the bootstrap in stratified sampling.  Annals of Statistics, 12, 470–482.
  * Brändén and Jonasson (2012) Brändén , P. and Jonasson, J. (2012)  Negative dependence in sampling.  Scandinavian Journal of Statistics, 39(4), 830–838.
  * Busso et al. (2014) Busso, M., DiNardo, J. E. and McCrary, J. (2014)  New evidence on the finite sample properties of propensity score reweighting and matching estimators.  Review of Economics and Statistics, 96, 885–897.
  * Cao et al. (2009) Cao, W., Tsiatis, A. A. and Davidian, M. (2009)  Improving efficiency and robustness of the doubly robust estimator for a population mean with incomplete data.  Biometrika, 96(3), 723–734.
  * Cattaneo (2010) Cattaneo, M. D. (2010)  Efficient semiparametric estimation of multi-valued treatment effects under ignorability.  Journal of Econometrics, 155, 138–154.
  * Chen and Rao (2007) Chen, J. and Rao, J. N. K. (2007) Asymptotic normality under two-phase sampling designs. Statistica Sinica, 17, 1047–1064.
  * Chen et al. (2008) Chen, J., Variyath, A. M. and Abraham, B. (2008)  Adjusted empirical likelihood and its properties.  Journal of Computational and Graphical Statistics, 17, 426–443.
  * Crump et al. (2009) Crump, R. K., Hotz, V. J., Imbens, G. W. and Mitnik, O. A. (2009)  Dealing with limited overlap in estimation of average treatment effects.  Biometrika, 96, 187–199.
  * Dong et al. (2020) Dong, G., Mao, L., Huang, B., Gamalo-Siebers, M., Wang, J., Yu, G. and Hoaglin, D. C. (2020)  The inverse-probability-of-censoring weighting (IPCW) adjusted win ratio statistic: an unbiased estimator in the presence of independent censoring.  Journal of Biopharmaceutical Statistics, 30(5), 882–899.
  * Erdös and Rényi (1959) Erdös, P. and Rényi, A. (1959)  On the central limit theorem for samples from a finite population.  Publication of the Mathematical Institute of the Hungarian Academy of Sciences, 4, 49–61.
  * Fan and Gijbels (1996) Fan, J., and Gijbels, I. (1996)  Local Polynomial Modelling and Its Applications. New York: Chapman and Hall.
  * Hájek (1960) Hájek, J. (1960)  Limiting distributions in simple random sampling from a finite population.  Publication of the Mathematical Institute of the Hungarian Academy of Sciences, 5, 361–374.
  * Hájek (1964) Hájek, J. (1964)  Asymptotic theory of rejective sampling with varying probabilities from a finite population.  Annals of Mathematical Statistics, 35, 1491–1523.
  * Hájek (1971) Hájek, J. (1971)  Discussion of ‘An essay on the logical foundations of survey sampling, Part One’ by D. Basu. In  Foundations of Statistical Inference (eds. V. P. Godambe and D. A. Sprott), vol. 236. Toronto: Holt, Rinehart and Winston.
  * Han (2013) Han, P. (2014)  Multiply robust estimation in regression analysis with missing data.  Journal of the American Statistical Association, 109, 1159–1173.
  * Han and Wang (2014) Han, P. and Wang, L. (2013)  Estimation with missing data: beyond double robustness.  Biometrika, 100, 417–430.
  * Han et al. (2019) Han, P., Kong, L., Zhao, J. and Zhou, X. (2019)  A general framework for quantile estimation with incomplete data.  Journal of the Royal Statistical Society, Series B, 82(2), 305–333.
  * Hansen and Hurwitz (1943) Hansen, M. and Hurwitz, W. (1943) On the theory of sampling from finite populations. Annals of Mathematical Statistics, 14, 333–362.
  * Hirano et al. (2003) Hirano, K., Imbens, G. W. and Ridder, G. (2003)  Efficient estimation of average treatment effects using the estimated propensity score.  Econometrica, 71, 1161–1189.
  * Horvitz and Thompson (1952) Horvitz, D. G. and Thompson, D. J. (1952)  A generalization of sampling without replacement from a finite universe.  Journal of the American Statistical Association, 47, 663–685.
  * Imbens and Wooldridge (2009) Imbens, G. W. and Wooldridge, J. (2009)  Recent developments in the econometrics of program evaluation.  Journal of Economic Literature, 47(1), 5–86.
  * Jiang et al. (2017) Jiang, R., Lu, W., Song, R. and Davidian, M. (2017)  On estimation of optimal treatment regimes for maximizing $t$-year survival probability.  Journal of the Royal Statistical Society, Series B, 79, 1165–1185.
  * Kang and Schafer (2007) Kang, J. D. and Schafer, J. L. (2007)  A comparison of alternative strategies for estimating a population mean from incomplete data.  Statistical Science, 22, 523–539.
  * Khan and Tamer (2010) Khan, S. and Tamer, E. (2010)  Irregular identification, support conditions, and inverse weight estimation.  Econometrica, 78(6), 2021–2042.
  * Kim and Shao (2014) Kim, J. K., and Shao, J. (2014)  Statistical Methods for Handling Incomplete Data.  Boca Raton, FL: CRC Press.
  * Krewski and Rao (1981) Krewski, D. and Rao, J. N. K. (1981)  Inference from stratified samples: properties of linearization, jackknife and balanced repeated replication methods.  Annals of Statistics, 9, 1010–1019.
  * LaLonde (1986) LaLonde, R. J. (1986)  Evaluating the econometric evaluations of training programs with experimental data.  American Economic Review, 76, 604–620.
  * Li and Qin (1998) Li, G. and Qin, J. (1998)  Semiparametric likelihood-based inference for biased and truncated data when the total sample size is known.  Journal of the Royal Statistical Society, Series, B, 60(1), 243–254.
  * Liu and Chen (2010) Liu, Y. and Chen, J. (2010)  Adjusted empirical likelihood with high-order precision.  Annals of Statistics, 38 (3), 1341–1362.
  * Liu et al. (2017) Liu, Y., Li, P. and Qin, J. (2017)  Maximum empirical likelihood estimation for abundance in a closed population from capture–recapture data.  Biometrika, 104, 527–543.
  * Ma and Yin (2011) Ma, Y. and Yin, G. (2011)  Censored quantile regression with covariate measurement errors.  Statistica Sinica, 21, 949–971.
  * Mccaffrey et al. (2013) Mccaffrey, D. F., Lockwood, J. R. and Setodji, C. M. (2013)  Inverse probability weighting with error-prone covariates.  Biometrika, 100(3), 671–680.
  * Ma and Wang (2020) Ma, X. and Wang, J. (2020)  Robust inference using inverse probability weighting.  Journal of the American Statistical Association, 115(532), 1851–1860.
  * Owen (1988) Owen, A. B. (1988)  Empirical likelihood ratio confidence intervals for a single functional.  Biometrika, 75, 237–249.
  * Owen (1990) Owen, A. B. (1990)  Empirical likelihood ratio confidence regions.  Annals of Statistics, 18, 90–120.
  * Owen (2001) Owen, A. B. (2001)  Empirical Likelihood.  New York: Chapman and Hall.
  * Patterson et al. (2001) Patterson, R. F., Smith, W. D., Taylor, R. L. and Bozorgnia, A. (2001)  Limit theorems for negatively dependent random variables.  Nonlinear Analysis, 47, 1283–1295.
  * Politis and Romano (1994) Politis, D. N. and Romano, J. P. (1994)  Large sample confidence regions based on subsamples under minimal assumptions.  Annals of Statistics, 22 (4), 2031–2050.
  * Qin and Lawless (1994) Qin, J. and Lawless, J. (1994)  Empirical likelihood and general equations.  Annals of Statistics, 22, 300–325.
  * Qin et al. (2008) Qin, J., Shao, J. and Zhang, B. (2008)  Efficient and doubly robust imputation for covariate-dependent missing responses.  Journal of the American Statistical Association, 103, 797–810.
  * Robins (1993) Robins, J. M. (1993)  Information recovery and bias adjustment in proportional hazards regression analysis of randomized trials using surrogate markers. In  Proceedings of the Biopharmaceutical Section, American Statistical Association, pp. 24–33. Alexandria, VA: American Statistical Association.
  * Robins and Finkelstein (2000) Robins, M. and Finkelstein, D. M. (2000)  Correcting for noncompliance and dependent censoring in an AIDS clinical trial with inverse probability of censoring weighted (IPCW) log-rank tests.  Biometrics, 56, 779–788.
  * Robins and Rotnitzky (1992) Robins, J. M. and Rotnitzky, A. (1992)  Recovery of information and adjustment for dependent censoring using surrogate markers.  In AIDS Epidemiology – Methodological Issues (eds. N. Jewell, K. Dietz and V. Farewell), pp. 297–331. Boston: Birkhäuser.
  * Robins et al. (1994) Robins, J. M., Rotnitzky, A. and Zhao, L. P. (1994)  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American Statistical Association, 89, 846–866.
  * Robins et al. (2007) Robins, J., Sued, M., Lei-Gomez, Q, and Rotnitzky, A. (2007)  Comment: Performance of double-robust estimators when ‘inverse probability’ weights are highly variable.  Statistical Science, 22(4), 544–559.
  * Rosenbaum (2002) Rosenbaum, P. R. (2002) Observational Studies. New York: Springer.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983)  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70, 41–55.
  * Rubin (1974) Rubin, D. B. (1974)  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of Educational Psychology, 66, 688–701.
  * Rubin (1974) Rubin, D. B. (1976)  Inference and missing data (with discussion).  Biometrika, 63, 581–592.
  * Sasaki and Ura (2018) Sasaki, Y. and Ura, T. (2018)  Estimation and inference for moments of ratios with robustness against large trimming bias.  arXiv:1709.00981.
  * Sun and Tchetgen Tchetgen (2018) Sun, B., and Tchetgen Tchetgen, E. J. (2018)  On inverse probability weighting for nonmonotone missing at random data.  Journal of the American Statistical Association, 113, 369–379.
  * Tan (2010) Tan, Z. (2010)  Bounded, efficient and doubly robust estimation with inverse weighting.  Biometrika, 97(3), 661–682.
  * Tan (2020) Tan, Z. (2020)  Model-assisted inference for treatment effects using regularized calibrated estimation with high-dimensional data.  Annals of Statistics, 48(2), 811–837.
  * Tsao (2004) Tsao, M. (2004)  Bounds on coverage probabilities of the empirical likelihood ratio confidence regions.  Annals of Statistics 32 1215–1221.
  * Wang et al. (2019) Wang, H., Yang, M. and Stufken, J. (2019)  Information-based optimal subdata selection for big data linear regression,  Journal of the American Statistical Association, 114, 393–405.
  * Wang et al. (2018) Wang, H., Zhu, R. and Ma, P. (2018)  Optimal subsampling for large sample logistic regression.  Journal of the American Statistical Association, 113, 829–844.
  * Wooldridge (2007) Wooldridge, J. M. (2007)  Inverse probability weighted estimation for general missing data problems.  Journal of Econometrics, 141, 1281–1301.
  * Wu and Thompson (2020) Wu, C. and Thompson, M. E. (2020)  Sampling Theory and Practice. Switzerland: Springer.
  * Yang and Ding (2018) Yang, S. and Ding, P. (2018)  Asymptotic inference of causal effects with observational studies trimmed by the estimated propensity scores.  Biometrika, 105(2), 487–493.
  * Young et al. (2019) Young, J. G., Logan, R. W., Robins, J. M. and Hernán, M. A. (2019)  Inverse probability weighted estimation of risk under representative interventions in observational studies.  Journal of the American Statistical Association, 114, 938–947.
  * Yu et al. (2020) Yu, J., Wang, H., Ai, M. and Zhang, H. (2020)  Optimal distributed subsampling for maximum quasi-likelihood estimators with massive data.  Journal of the American Statistical Association. DOI: 10.1080/01621459.2020.1773832.
  * Zhang et al. (2012) Zhang, B., Tsiatis, A. A., Laber, E. B. and Davidian, M. (2012)  A robust method for estimating optimal treatment regimes.  Biometrics, 68, 1010–1018.
  * Zhao (2019) Zhao, Q. (2019)  Covariate balancing propensity score by tailored loss functions.  Annals of Statistics, 47(2), 965–993.
  * Zong et al. (2019) Zong, X., Zhu, R. and Zou, G. (2019)  Improved Horvitz–Thompson estimator in survey sampling.  Survey Methodology, 45(1), 165–184.

Table 1:  Simulated RMSEs of the estimators under comparison when data are generated from Example 1 and $N=2000$. Smallest RMSEs are highlighted in bold.  $\gamma$ | $c$ | Model | IPW | SIPW | ZZZ | CHIM | MW1 | MW2 | ELW
---|---|---|---|---|---|---|---|---|---
1.5 | 1.0 | 1 | 24.72 | 8.05 | 6.05 | 8.00 | 4.96 | 4.35 | 5.51
1.5 | 1.0 | 2 | 17.89 | 6.17 | 5.95 | 6.17 | 4.84 | 3.78 | 5.13
1.5 | 1.0 | 3 | 69.08 | 7.49 | 27.27 | 7.49 | 18.29 | 9.59 | 5.21
1.5 | 1.0 | 4 | 110.80 | 6.49 | 27.11 | 6.49 | 18.31 | 9.83 | 5.21
1.5 | 0.1 | 1 | 14.76 | 4.89 | 4.48 | 4.87 | 3.00 | 2.52 | 1.60
1.5 | 0.1 | 2 | 26.23 | 2.16 | 4.44 | 2.15 | 2.88 | 1.36 | 0.71
1.5 | 0.1 | 3 | 68.12 | 4.74 | 27.04 | 4.73 | 17.81 | 8.94 | 1.61
1.5 | 0.1 | 4 | 140.05 | 2.21 | 26.86 | 2.19 | 18.02 | 9.03 | 0.74
2.5 | 1.0 | 1 | 2.11 | 2.11 | 1.97 | 2.11 | 1.93 | 1.87 | 2.02
2.5 | 1.0 | 2 | 2.06 | 1.81 | 1.90 | 1.81 | 1.89 | 1.82 | 1.72
2.5 | 1.0 | 3 | 7.64 | 2.15 | 6.77 | 2.15 | 6.33 | 5.32 | 2.05
2.5 | 1.0 | 4 | 8.14 | 1.85 | 7.31 | 1.85 | 6.87 | 6.01 | 1.70
2.5 | 0.1 | 1 | 1.49 | 1.33 | 1.14 | 1.33 | 1.07 | 0.98 | 1.17
2.5 | 0.1 | 2 | 1.22 | 0.69 | 1.01 | 0.69 | 0.93 | 0.77 | 0.42
2.5 | 0.1 | 3 | 7.63 | 1.31 | 6.60 | 1.31 | 6.21 | 5.14 | 1.18
2.5 | 0.1 | 4 | 8.26 | 0.68 | 7.13 | 0.68 | 6.80 | 5.85 | 0.42
Table 2:  Simulated coverage probabilities (%) of the interval estimators under comparison when data were generated from Example 1 with total sample size $N=2000$. The numbers in parentheses are average lengths.  $\gamma$ | $c$ | Model | IPW-an | IPW-re | SIPW-an | SIPW-re | MW1-re | MW2-re | ELW-an | ELW-re
---|---|---|---|---|---|---|---|---|---|---
1.5 | 1.0 | 1 | $\underset{(0.556)}{76.14}$ | $\underset{(3.029)}{84.52}$ | $\underset{(0.436)}{78.42}$ | $\underset{(1.455)}{88.86}$ | $\underset{(1.060)}{89.44}$ | $\underset{(0.780)}{92.86}$ | $\underset{(0.329)}{82.58}$ | $\underset{(1.044)}{91.48}$
1.5 | 1.0 | 2 | $\underset{(1.593)}{78.24}$ | $\underset{(6397.733)}{87.24}$ | $\underset{(0.366)}{85.08}$ | $\underset{(1.159)}{91.14}$ | $\underset{(1.165)}{89.68}$ | $\underset{(0.831)}{94.52}$ | $\underset{(0.312)}{81.62}$ | $\underset{(1.047)}{91.04}$
1.5 | 1.0 | 3 | $\underset{(2.948)}{76.76}$ | $\underset{(18.447)}{86.40}$ | $\underset{(0.446)}{77.84}$ | $\underset{(1.575)}{88.94}$ | $\underset{(3.325)}{82.74}$ | $\underset{(1.507)}{81.44}$ | $\underset{(0.333)}{82.02}$ | $\underset{(1.061)}{90.94}$
1.5 | 1.0 | 4 | $\underset{(3.282)}{78.74}$ | $\underset{(121.859)}{86.32}$ | $\underset{(0.370)}{85.36}$ | $\underset{(1.160)}{91.60}$ | $\underset{(3.276)}{82.32}$ | $\underset{(1.551)}{83.32}$ | $\underset{(0.312)}{82.54}$ | $\underset{(1.046)}{92.38}$
1.5 | 0.1 | 1 | $\underset{(0.491)}{76.94}$ | $\underset{(6.361)}{88.86}$ | $\underset{(0.304)}{79.58}$ | $\underset{(1.031)}{91.54}$ | $\underset{(0.543)}{79.76}$ | $\underset{(0.261)}{68.98}$ | $\underset{(0.126)}{92.86}$ | $\underset{(0.287)}{95.72}$
1.5 | 0.1 | 2 | $\underset{(0.477)}{76.22}$ | $\underset{(5.052)}{87.68}$ | $\underset{(0.130)}{76.54}$ | $\underset{(0.436)}{89.20}$ | $\underset{(0.552)}{81.92}$ | $\underset{(0.224)}{82.40}$ | $\underset{(0.049)}{87.62}$ | $\underset{(0.122)}{91.72}$
1.5 | 0.1 | 3 | $\underset{(2.952)}{77.40}$ | $\underset{(22.300)}{86.18}$ | $\underset{(0.304)}{78.00}$ | $\underset{(1.036)}{90.44}$ | $\underset{(3.177)}{80.50}$ | $\underset{(1.293)}{78.90}$ | $\underset{(0.125)}{91.84}$ | $\underset{(0.287)}{94.84}$
1.5 | 0.1 | 4 | $\underset{(2.958)}{77.88}$ | $\underset{(31.359)}{85.96}$ | $\underset{(0.131)}{76.22}$ | $\underset{(0.441)}{89.00}$ | $\underset{(3.125)}{81.76}$ | $\underset{(1.361)}{79.12}$ | $\underset{(0.049)}{85.38}$ | $\underset{(0.122)}{89.28}$
2.5 | 1.0 | 1 | $\underset{(0.178)}{94.02}$ | $\underset{(0.341)}{91.40}$ | $\underset{(0.177)}{93.60}$ | $\underset{(0.351)}{92.46}$ | $\underset{(0.323)}{91.92}$ | $\underset{(0.323)}{93.14}$ | $\underset{(0.169)}{93.32}$ | $\underset{(0.343)}{93.20}$
2.5 | 1.0 | 2 | $\underset{(0.169)}{93.52}$ | $\underset{(0.346)}{93.48}$ | $\underset{(0.150)}{93.72}$ | $\underset{(0.303)}{93.54}$ | $\underset{(0.324)}{92.38}$ | $\underset{(0.321)}{93.40}$ | $\underset{(0.141)}{93.72}$ | $\underset{(0.287)}{93.60}$
2.5 | 1.0 | 3 | $\underset{(0.630)}{94.14}$ | $\underset{(1.261)}{93.68}$ | $\underset{(0.177)}{93.86}$ | $\underset{(0.348)}{92.74}$ | $\underset{(1.029)}{89.40}$ | $\underset{(0.837)}{88.16}$ | $\underset{(0.169)}{93.56}$ | $\underset{(0.344)}{93.40}$
2.5 | 1.0 | 4 | $\underset{(0.682)}{94.26}$ | $\underset{(1.346)}{93.34}$ | $\underset{(0.150)}{93.80}$ | $\underset{(0.299)}{94.06}$ | $\underset{(1.118)}{90.24}$ | $\underset{(0.948)}{88.96}$ | $\underset{(0.141)}{93.46}$ | $\underset{(0.286)}{93.94}$
2.5 | 0.1 | 1 | $\underset{(0.107)}{95.18}$ | $\underset{(0.203)}{92.88}$ | $\underset{(0.107)}{94.94}$ | $\underset{(0.213)}{93.92}$ | $\underset{(0.176)}{91.18}$ | $\underset{(0.157)}{90.90}$ | $\underset{(0.098)}{94.04}$ | $\underset{(0.208)}{94.72}$
2.5 | 0.1 | 2 | $\underset{(0.095)}{93.76}$ | $\underset{(0.208)}{94.42}$ | $\underset{(0.054)}{94.02}$ | $\underset{(0.119)}{94.90}$ | $\underset{(0.159)}{89.40}$ | $\underset{(0.127)}{88.28}$ | $\underset{(0.035)}{94.50}$ | $\underset{(0.071)}{94.04}$
2.5 | 0.1 | 3 | $\underset{(0.616)}{94.42}$ | $\underset{(1.231)}{93.92}$ | $\underset{(0.108)}{94.64}$ | $\underset{(0.214)}{93.48}$ | $\underset{(0.986)}{89.36}$ | $\underset{(0.784)}{87.96}$ | $\underset{(0.099)}{94.02}$ | $\underset{(0.208)}{94.80}$
2.5 | 0.1 | 4 | $\underset{(0.664)}{94.08}$ | $\underset{(1.306)}{92.96}$ | $\underset{(0.054)}{94.34}$ | $\underset{(0.117)}{95.08}$ | $\underset{(1.081)}{89.90}$ | $\underset{(0.901)}{88.70}$ | $\underset{(0.035)}{94.38}$ | $\underset{(0.070)}{94.32}$
Table 3:  Coverage probabilities and average length of Wald intervals based on the IPW, SIPW, and ELW estimators when data were generated from Example 2 when $n=500$.  |  | Coverage probability | Average length
---|---|---|---
$\rho$ | Model | IPW-an | SIPW-an | ELW-an | IPW-an | SIPW-an | ELW-an
|  | Poisson sampling
0.2 | 1 | 93.28 | 95.34 | 93.86 | 0.476 | 0.280 | 0.251
0.2 | 2 | 94.26 | 93.18 | 93.52 | 0.507 | 0.333 | 0.275
0.2 | 3 | 93.06 | 95.16 | 93.96 | 1.771 | 0.278 | 0.252
0.2 | 4 | 93.94 | 93.78 | 93.58 | 1.811 | 0.334 | 0.274
0.8 | 1 | 94.72 | 90.90 | 93.52 | 0.431 | 0.343 | 0.232
0.8 | 2 | 96.04 | 91.84 | 95.34 | 0.684 | 0.776 | 0.479
0.8 | 3 | 93.80 | 90.26 | 93.54 | 1.752 | 0.340 | 0.233
0.8 | 4 | 94.66 | 91.78 | 94.92 | 1.928 | 0.778 | 0.479
|  | Pivotal sampling
0.2 | 1 | 96.82 | 94.34 | 92.84 | 0.494 | 0.270 | 0.247
0.2 | 2 | 98.66 | 92.52 | 93.68 | 0.528 | 0.313 | 0.262
0.2 | 3 | 97.86 | 94.84 | 93.40 | 1.789 | 0.267 | 0.246
0.2 | 4 | 98.52 | 92.48 | 93.56 | 1.828 | 0.311 | 0.263
0.8 | 1 | 99.98 | 89.40 | 94.90 | 0.438 | 0.317 | 0.214
0.8 | 2 | 100.00 | 91.46 | 97.48 | 0.698 | 0.729 | 0.436
0.8 | 3 | 98.98 | 88.80 | 94.46 | 1.734 | 0.310 | 0.214
0.8 | 4 | 99.88 | 89.78 | 97.10 | 1.936 | 0.724 | 0.435
|  | PPS sampling
0.2 | 1 | 92.34 | 94.48 | 91.34 | 0.433 | 0.289 | 0.253
0.2 | 2 | 93.54 | 92.96 | 91.80 | 0.451 | 0.333 | 0.262
0.2 | 3 | 91.60 | 94.52 | 91.28 | 1.560 | 0.288 | 0.253
0.2 | 4 | 91.70 | 93.06 | 92.28 | 1.558 | 0.334 | 0.263
0.8 | 1 | 97.52 | 89.46 | 90.64 | 0.311 | 0.333 | 0.188
0.8 | 2 | 100.00 | 90.86 | 91.94 | 0.419 | 0.762 | 0.330
0.8 | 3 | 92.18 | 89.58 | 91.32 | 1.471 | 0.335 | 0.187
0.8 | 4 | 94.34 | 90.56 | 91.04 | 1.468 | 0.775 | 0.330
Table 4:  Simulated RMSEs of the IPW, SIPW, ZZZ, and ELW estimators where data were generated from Example 2 with $n=500$.  |  | IPW | SIPW | ZZZ | ELW | IPW | SIPW | ZZZ | ELW | IPW | SIPW | ZZZ | ELW
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\rho$ | Model | Poisson sampling | Pivotal sampling | PPS sampling
0.2 | 1 | 9.04 | 4.34 | 6.28 | 3.93 | 40.67 | 4.38 | 5.03 | 3.78 | 7.41 | 4.39 | 5.11 | 3.87
0.2 | 2 | 11.38 | 5.36 | 6.95 | 4.16 | 40.64 | 5.36 | 4.89 | 3.92 | 7.28 | 5.31 | 4.94 | 4.03
0.2 | 3 | 36.21 | 4.59 | 23.56 | 3.89 | 179.21 | 4.38 | 17.48 | 3.78 | 27.46 | 4.39 | 18.01 | 3.87
0.2 | 4 | 36.55 | 5.36 | 24.06 | 4.16 | 179.18 | 5.36 | 17.27 | 3.92 | 27.29 | 5.31 | 17.76 | 4.03
0.8 | 1 | 8.07 | 5.66 | 5.85 | 3.39 | 24.90 | 5.88 | 3.07 | 2.91 | 4.53 | 5.74 | 3.13 | 3.00
0.8 | 2 | 10.79 | 12.84 | 9.24 | 6.56 | 24.84 | 13.29 | 2.79 | 5.03 | 4.27 | 13.14 | 2.79 | 5.28
0.8 | 3 | 33.58 | 5.66 | 23.09 | 3.39 | 163.52 | 5.88 | 15.86 | 2.91 | 24.97 | 5.74 | 16.36 | 3.00
0.8 | 4 | 35.35 | 12.84 | 25.53 | 6.56 | 163.41 | 13.29 | 15.03 | 5.03 | 24.32 | 13.14 | 15.40 | 5.28
Table 5:  Simulated RMSEs of the IPW, SIPW, ZZZ and ELW estimators when data were generated from the LLvsPSID dataset with $n=200$.  | IPW | ZZZ | IPW | ZZZ | SIPW | ELW
---|---|---|---|---|---|---
| $a=0$ | $a=2$ |  |
Poisson sampling | 9.35 | 8.44 | 19.27 | 16.33 | 8.41 | 6.14
Pivotal sampling | 5.07 | 3.91 | 12.17 | 7.63 | 7.15 | 4.66
PPS sampling | 5.46 | 4.13 | 13.86 | 8.19 | 8.70 | 5.51

Figure 1: Boxplots of the SIPW, ZZZ, CHIM, MW1, MW2, and ELW estimators when
data were generated from Example 1 with $N=2000$ and $\gamma=1.5$. The solid
horizontal line corresponds to the target parameter value.

Figure 2: Boxplots of the SIPW, ZZZ, CHIM, MW1, MW2, and ELW estimators when
data were generated from Example 1 with $N=2000$ and $\gamma=2.5$. The solid
horizontal line corresponds to the target parameter value.

Figure 3:  Boxplots of the SIPW, ZZZ, and ELW estimators when data were
generated by pivotal sampling from Example 2 with $n=500$. The solid
horizontal line corresponds to the target parameter value.

Figure 4:  Histograms of the variable $Y$ and the fitted propensity score in
the treated and control groups, based on the LLvsPSID data.