# Review of Quasi-Randomization Approaches for Estimation from Non-probability
Samples

Vladislav Beresovsky<EMAIL_ADDRESS>[ Julie
<EMAIL_ADDRESS>[ Terrance D.
<EMAIL_ADDRESS>[ Vladislav Beresovsky is Research
Mathematical Statistician, U.S. Bureau of Labor Statisticspresep= ]e1.  Julie
Gershunskaya is Mathematical Statistician, U.S. Bureau of Labor
Statisticspresep= ]e2.  Terrance D. Savitsky is Senior Research Mathematical
Statistician, U.S. Bureau of Labor Statisticspresep= ]e3.

###### Abstract

The recent proliferation of computers and the internet have opened new
opportunities for collecting and processing data. However, such data are often
obtained without a well-planned probability survey design. Such non-
probability based samples cannot be automatically regarded as representative
of the population of interest. Several classes of methods for estimation and
inferences from non-probability samples have been developed in recent years.
The quasi-randomization methods assume that non-probability sample selection
is governed by an underlying latent random mechanism. The basic idea is to use
information collected from a probability (“reference”) sample to uncover
latent non-probability survey participation probabilities (also known as
“propensity scores”) and use them in estimation of target finite population
parameters. In this paper, we review and compare theoretical properties of
recently developed methods of estimation survey participation probabilities
and study their relative performances in simulations.

Data combining,

non-probability sample,

reference sample,

participation probabilities,

sample likelihood,

variance estimation,

###### keywords:

, and

##  1 Introduction

Probability based sampling is a widely accepted approach to estimation of
target population parameters. However, traditional well-planned probability
surveys involve substantial costs and the labor intensive and time consuming
data collection. In contrast, recent proliferation of computers and the
internet created new ways for acquiring data. A lot of information now can be
collected seamlessly and is often available as a simple by-product of other
business activities. Such data sometimes are called “convenience” samples. The
demand for exploiting these resources is steadily growing. However, the
downside is that, no matter how big the samples are, data obtained without a
well-planned probability survey design cannot be automatically regarded as
being representative of the population of interest. Bethlehem (2010), Meng
(2018), VanderWeele and Shpitser (2011), among others, warn that naïve
estimation from non-probability based sources may lead to substantial biases.

A number of methods for proper estimation and inferences from non-probability
samples has been developed. The methodology can be broadly grouped into model-
based, quasi-randomization, and doubly-robust approaches. Elliott and Valliant
(2017), Valliant (2020), Beaumont and Rao (2021), Rao (2021), Wu (2022)
provide a comprehensive account of this rapidly developing field.

In the current paper, we focus on quasi-randomization approaches. These
methods assume that non-probability sample selection is governed by an
underlying latent random mechanism. The basic idea is to use information
collected from a probability sample to uncover latent non-probability survey
participation probabilities and use them in the estimation of target finite
population parameters. In this context, the probability sample is often
regarded as “reference” sample.

It is assumed that both the non-probability and reference samples contain a
common set of covariates and the non-probability sample participation
probabilities are governed by these covariates and could be estimated from a
model. Valliant (2020) also notes the importance of the _common support_
requirement, which means that the target population units are assumed to have
positive probabilities to be included into both samples conditional on the
shared set of covariates. For a broad range of typically used covariates, such
as large demographic groups, there is usually enough overlap in the
characteristics of units, so that the common support requirement is satisfied
for such covariates. It becomes less tenable, however, for more detailed
classifications and for covariates identifying smaller population groups. In
this case, we are faced with a potentially insufficient representation of some
population groups in either of the two samples. Application of an efficient
and robust methodology would lead to better utilization of available detailed
covariates for a more accurate estimation of participation probabilities.

Since estimation of participation probabilities depends on a model, it is
important to use a methodology that allows efficient model development,
checking, and estimation. Beaumont et al. (2023) proposed several methods of
variables selection based on the pseudo-likelihood approach of Chen, Li and Wu
(2020), employing the Akaike Information Criteria (AIC) and classification and
regression trees (CART). Machine Learning (ML) techniques are especially
promising when there is a large number of covariates. However, ML methods may
lead to terminal nodes where it is difficult to satisfy the common support
requirement, thus using an efficient objective function becomes critical for
model selection. In this paper, we review several methods for estimation of
survey participation probabilities and discuss our likelihood-based approach
that is more efficient than pseudo-likelihood based approaches for practically
relevant scenarios. This paper can be considered as a sequel to Savitsky et
al. (2022) who used a Bayesian framework for the modeling. In the current
paper, we derive asymptotic results and compare theoretical properties of the
discussed methods, and explore their relative performances using a numerical
study and simulations.

Our method is based on the so-called “stacked” samples setup, where the non-
probability and probability samples are put together and the overlapping
units, if any, appear twice in the stacked set. There is a common perception
in the literature that such a construction is only viable if inclusion
probabilities are small, in which case the probability of the overlap could be
ignored. See, for example, Elliott (2009) and Elliott and Valliant (2017).
Liu, Scholtus and De Waal (2023) and Kim and Tam (2021) advise to remove the
overlapping units from the stacked set. In practice, we often do not know
which units belong to both samples. Liu, Scholtus and De Waal (2023) proposed
to identify the common units via probabilistic record linkage under the
assumption that it can be done with high probability, and Kim and Tam (2021)
used a nonparametric classification method. In this paper, we provide the
proof that, under usual assumptions, inclusion indicators in the staked
samples are independent conditional on covariates even when the probability of
the overlap is large. This allows us to use the likelihood for independent
Bernoulli variables formulated over all units in the stacked sample and avoid
the need to identify and remove the overlapping units.

In Section 2.1, we introduce the setup and notation for the modeling of
convenience sample participation probabilities. In Section 2.2, we review the
method of Chen, Li and Wu (2020) that directly models the unknown convenience
sample inclusion probabilities, $\pi_{c}$, in a likelihood defined on the
population. Section 2.3 introduces an alternative class of methods that
_indirectly_ models $\pi_{c}$ in their likelihood statements through composing
a functional relationship between probabilities of inclusion from the
population to probabilities of being in a given part of the observed “stacked”
convenience and reference samples. We state the relationship and provide its
proof that utilizes conditional independence of inclusions in the convenience
and reference samples. In Section 2.4, we describe methods for modeling
convenience sample participation probabilities using this relationship.
Theoretical properties of the estimators of model parameters and respective
inverse propensity weighted (IPW) estimators of the population mean are given
in Section 3 for the $3$ methods on which we focus. Section 4 presents a
simulation study that compares methods under various scenarios that include
varying reference and convenience sample sizes and sampling fractions, as well
as different degrees of the common support. We provide concluding remarks in
Section 5.

##  2 Estimation of convenience sample participation probabilities

###  2.1 Setup and notation

In this paper we focus on estimation of the finite population mean
${\mu}={N^{-1}}\mathop{\sum}_{i\in U}y_{i}$, for units $i=1,\ldots,N$ in
target population $U$, in the situation where variables of interest $y_{i}$
are observed only on the non-probability (convenience) sample $S_{c}$ of size
$n_{c}$. In the sequel we often use the term “convenience” sample to generally
refer to a non-probability sample. We assume that the observed convenience
sample is a realization of a latent (unknown) random sampling design and
define _unknown_ probability
$\pi_{ci}=\pi_{c}(\mathbf{x}_{i})=P\\{I_{ci}=1|i\in U,\mathbf{x}_{i}\\}$ for
units in population $U$ to be included into the convenience sample. Here,
$I_{ci}$ is the inclusion indicator of population unit $i$, taking on the
value of 1 if unit $i$ is included into $S_{c}$, and $0$ if it is in
$U\backslash S_{c}$; $\mathbf{x}_{i}$ is a set of known covariates. We refer
to $\pi_{ci}$ as a “participation” probability (as differentiated from an
inclusion probability) to recognize a typical realization of a convenience
sample by self-selection, rather than solicitation.

A probability (reference) sample $S_{r}$ of size $n_{r}$ is also assumed to be
available. It is selected from the same finite population $U$ under a known
probability survey design with inclusion probabilities
$\pi_{ri}=P\\{I_{ri}=1\mid i\in U,\mathbf{x}^{{}^{\prime}}_{i}\\}$, where
$I_{ri}$ is the reference sample inclusion indicator of unit $i$. Suppose
samples $S_{c}$ and $S_{r}$ share a common set of covariates
$\mathbf{x}_{i}=\mathbf{x}^{{}^{\prime}}_{i}$.

We assume the following about the selection process:

  1. A1:

Each unit in the population has a positive participation probability:
$\pi_{ci}>0$ for all $i\in U.$

  2. A2:

Sample $S_{c}$ selection mechanism is ignorable given $\mathbf{x}_{i}$:
$P\\{I_{ci}=1|i\in U,\mathbf{x}_{i},y_{i}\\}=P\\{I_{ci}=1|i\in
U,\mathbf{x}_{i}\\}$ for all $i\in U$.

  3. A3:

Indicators $I_{ci}$ and $I_{cj}$ are independent given $\mathbf{x}_{i}$ and
$\mathbf{x}_{j}$ for $i\neq j.$

  4. A4:

Indicators $I_{ci}$ and $I_{ri}$ are independent given $\mathbf{x}_{i}$ for
any $i\in U.$

When the randomized sampling design and associated inclusion probabilties are
known, a condition like Assumption A1: ensures that the resulting design
distribution (that governs the taking of repeated samples from the target
population) will produce a collection of samples that together represent all
population subgroups. In our setting where the $\pi_{c}$ are estimated,
Assumption A1: may be interpreted as a “common support” assumption that
requires a complete overlap of the support of predictor values (or the
distributions of predictor values) between the convenience and reference
samples. See also discussion on a common support requirement in Valliant
(2020).

Quasi-randomization approaches propose methods for estimation of non-
probability sample participation probabilities $\pi_{ci}$ by combining non-
probability and probability samples. These estimated probabilities are then
used in the IPW estimator of the finite population mean.

###  2.2 Direct estimation of participation probabilities

We review the method of Chen, Li and Wu (2020) for combining information from
both samples to estimate the non-probability / convenience sample
participation probabilities. Their method is in the class of pseudo-likelihood
approaches where a log-likelihood is formulated over the population units, as
if the whole finite population was observed. Next, since only a probability
sample is observed, the log-likelihood is replaced by a pseudo log-likelihood
formulation, using the probability sample weights from the observed reference
sample. Model parameters are obtained as a solution to respective pseudo-
likelihood based estimating equations. Finally, the estimated non-probability
sample participation probabilities are used to form the IPW estimator of the
finite population mean.

Chen, Li and Wu (2020) (hereafter CLW) write a log-likelihood over units in
$U$, with respect to Bernoulli variable $I_{ci}$:

$\displaystyle\ell^{CLW}(\boldsymbol{\beta})$
$\displaystyle=\mathop{\sum}_{i\in
U}\left\\{I_{ci}\log\left[{\pi_{ci}(\boldsymbol{\beta})}\right]+(1-I_{ci})\log\left[1-\pi_{ci}\left(\boldsymbol{\beta}\right)\right]\right\\}$
(2) $\displaystyle=\mathop{\sum}_{i\in
S_{c}}\log\left[\frac{\pi_{ci}(\boldsymbol{\beta})}{1-\pi_{ci}\left(\boldsymbol{\beta}\right)}\right]+\mathop{\sum}_{i\in
U}\log\left[1-\pi_{ci}\left(\boldsymbol{\beta}\right)\right],$

where $\boldsymbol{\beta}$ is the parameter vector in a logistic regression
model
$\text{logit}\left[\pi_{ci}(\boldsymbol{\beta})\right]=\boldsymbol{\beta^{T}}\mathbf{x}_{i}$.

We refer to this method as “direct” estimation because $\pi_{ci}$ is directly
used in the Bernoulli likelihood (2) formulated for the _unobserved_ non-
probabilty sample indicators $I_{ci}$.

Since population $U$ is not available, CLW employ a pseudo-likelihood approach
by, first, regrouping units, as in (2), and then replacing the sum over the
finite population with its probability sample based estimate:

$\displaystyle\hat{\ell}^{CLW}(\boldsymbol{\beta})$
$\displaystyle=\mathop{\sum}_{i\in
S_{c}}{\log\left[\frac{\pi_{ci}(\boldsymbol{\beta})}{1-\pi_{ci}\left(\boldsymbol{\beta}\right)}\right]}+\mathop{\sum}_{i\in
S_{r}}{w_{ri}\log\left[1-\pi_{ci}\left(\boldsymbol{\beta}\right)\right]},$ (3)

where weights $w_{ri}=\pi_{ri}^{-1}$ are inverse values of the reference
sample inclusion probabilities $\pi_{ri}$. Estimates are obtained by solving
respective pseudo-likelihood based estimating equations.

Indicator $I_{ci}$ is typically not observed because the finite population is
not generally available; in particular, one does not know which units from the
finite population are selected into the convenience sample. In the sequel we
define an indicator for whether a unit from a pooled reference and convenience
sample is included in the convenience sample, which is observed.

We must tie the associated probability of membership in the pooled sample to
the convenience sample participation probability that we seek to estimate
using a model likelihood. In the next Section, we present such an expression.

###  2.3 Indirect estimation of participation probabilities

We next derive a formula relating the probability of membership in the pooled
sample to the convenience sample participation probability that we seek to
estimate using a model likelihood. This formula is used by the two methods
introduced in Section 2.4 to estimate the latent convenience sample
participation probabilities. The formula, named Core Relationship for
Independent Sampling Probabilities (CRISP), is based on the assumption that,
after conditioning on covariates, inclusions into the non-probability and
reference samples are independent (see Assumption A4: in Section 2.1).

Let us consider stacking together samples $S_{c}$ and $S_{r}$ and denote the
resulting set by $S$. We define indicator variable $I_{zi}$ on such a combined
sample $S$ as follows: $I_{zi}=1$ if unit $i$ belongs to the non-probability
sample (part 1), and $I_{zi}=0$ if unit $i$ belongs to the probability sample
(part 2). Under this construction, if there is an overlap between the two
samples, $S_{c}$ and $S_{r}$, then the overlapping units are included into the
stacked set, $S$, twice: as part of the non-probability sample (with
$I_{zi}=1$) and as part of the reference probability sample (with $I_{zi}=0$).
Savitsky et al. (2022) use first principles to prove the relationship between
probabilities $\pi_{zi}=P\\{I_{zi}=1|i\in S,\mathbf{x}_{i}\\}$ of being in
part 1 of the stacked set, on the one hand, and inclusion probabilities,
$\pi_{ci}$ and $\pi_{ri}$, on the other hand. For completeness, we reproduce
this result in Theorem 2.1 of this Section. Furthermore, in Theorem 2.2 below,
we show that, regardless of the overlap, under certain usual assumptions,
indicators $I_{zi}$ and $I_{zj}$ are conditionally independent for any two
units $i$ and $j$ in stacked set $S$.

We note that the stacked samples setup should be contrasted with the combined
sample approach of Elliott (2009) and Elliott and Valliant (2017), where they
use the assumption that the overlap is negligible between the samples. The
stacked samples setup allows for an overlap of any size – from no overlap to a
situation where one sample is a subset of another; moreover, there is no need
to know the size of the overlap or the identities of the units that belong to
both samples.

We generalize the set-up in Section 2.1 to allow each of the convenience and
reference samples to participate or be selected from respective frames that
are _subsets_ of target population $U$. Namely, let the non-probability /
convenience sample be realized from frame $U_{c}\subseteq U$ and probability /
reference sample be selected from frame $U_{r}\subseteq U$ with respective
inclusion probabilities $\pi_{ci}=P\\{i\in S_{c}|i\in U_{c}\\}$ and
$\pi_{ri}=P\\{i\in S_{r}|i\in U_{r}\\}$. Define
$p_{ci}=p_{c}(\mathbf{x}_{i})=P\\{i\in U_{c}|i\in U,\mathbf{x}_{i}\\}$ and
$p_{ri}=p_{r}(\mathbf{x}_{i})=P\\{i\in U_{r}|i\in U,\mathbf{x}_{i}\\}$ to be
coverage probabilities. A coverage probability denotes the probability that a
unit from the target population is included in the frame from which sample
units participate. In order to ensure that realized frames are representative
of target population, $U$, we assume that both $(p_{ci}>0,p_{ri}>0)$,
strictly. This assumption of positive coverage probabilities may be violated
in the case of hard-to-reach subgroups in the population. Practitioners should
carefully examine the correlation between population characteristics of
interest, on the one hand, and covariates, $\mathbf{x}_{i}$, used to define
coverage probabilities, on the other hand. If some subgroups are known, by
design, to have $0$ probability to participate in the convenience sample, then
Valliant (2020) recommend to exclude these subgroups from the studied
population.

For example, suppose household units in $U_{c}\subseteq U$ have broadband
internet access and can potentially participate in a web survey. Assume, each
household in the target population has positive probability $p_{ci}=P\\{i\in
U_{c}|i\in U,\mathbf{x}_{i}\\}$ to have broadband internet access and let this
probability be _known_ for some broadly defined demographic groups
$\mathbf{x}_{i}.$ At the same time, survey participation probabilities
$\pi_{ci}=P\\{i\in S_{c}|i\in U_{c},\mathbf{x}_{i}\\}$ are _not known_.
Considering the reference sample case, suppose each person in the target
population has positive probability to have a mailable address and let this
probability $p_{ri}=P\\{i\in U_{r}|i\in U,\mathbf{x}_{i}\\}$ be known for a
given demography. Frame $U_{r}$ is a subset of $U$ consisting of persons who
actually have mailable addresses and thus could be selected into the reference
sample. In this example, $p_{ci}$ and $p_{ri}$ play the role of known coverage
probabilities.

Figure 1: Gridded area represents observed convenience $S_{c}$ and reference
$S_{r}$ samples stacked together to form combined sample $S$: $S=S_{c}+S_{r}$;
under this scheme, if samples $S_{c}$ and $S_{r}$ overlap, the overlapping
units are included in $S$ twice. Convenience sample $S_{c}$ is selected from
frame $U_{c}$, and reference sample $S_{r}$ is selected from frame $U_{r}$,
where $U_{c}$ and $U_{r}$ are subsets of target population $U$:
$U_{c}\subseteq U$ and $U_{r}\subseteq U$. In this setup, two identical copies
of target population $U$ are stacked together, so that $2U=U+U$.

######  Theorem 2.1 (CRISP).

Assuming samples $S_{c}$ and $S_{r}$ are independently selected from,
respectively, $U_{c}\subseteq U$ and $U_{r}\subseteq U$ with coverage
probabilities, $(p_{ci}>0,p_{ri}>0)$, the following relationship holds:

$\displaystyle{\pi_{zi}}=\frac{{\pi_{ci}}p_{ci}}{{\pi_{ci}}p_{ci}+{\pi_{ri}}p_{ri}}.$
(4)

###### Proof.

Since the sampling selection into the samples is performed independently, the
combined set $S$ can be viewed as emerging from the following imaginary
stacked populations scheme displayed in Figure 1. We consider two copies of
target population $U$, where one copy of the population includes frame
$U_{c}$, the other copy includes $U_{r}$. We _stack_ the two copies of $U$
together and denote the result by $2U$: $2U=U+U$. We next use the scheme of
Figure 1 in our proof.

For such a setup, by the Law of Total Probability (LTP), we have:

$\displaystyle P\left\\{i\in S_{c}|i\in 2U\right\\}=P\left\\{i\in S_{c}|i\in
U_{c},i\in U\right\\}\times$ (5) $\displaystyle P\left\\{i\in U_{c}|i\in
U\right\\}P\\{i\in U\mid i\in 2U\\}=\frac{1}{2}\pi_{ci}p_{ci}.$

In (5), we used the fact that $i\in S_{c}$ implies that $i\in U_{c}$ since we
draw the convenience sample from its associated frame, $U_{c}$, and thus
$P\left\\{i\in S_{c}|i\in U_{c},i\in U\right\\}=P\left\\{i\in S_{c}|i\in
U_{c}\right\\}=\pi_{ci}.$

Similarly,

$\displaystyle P\left\\{i\in S_{r}|i\in 2U\right\\}=P\left\\{i\in S_{r}|i\in
U_{r},i\in U\right\\}\times$ $\displaystyle P\left\\{i\in U_{r}|i\in
U\right\\}P\\{i\in U\mid i\in 2U\\}=\frac{1}{2}\pi_{ri}p_{ri}.$

The total probability of being included into the pooled sample is simply the
sum of probabilities for the non-overlapping sets:

$\displaystyle P\left\\{i\in S|i\in 2U\right\\}$ $\displaystyle=P\left\\{i\in
S_{c}|i\in 2U\right\\}+P\left\\{i\in S_{r}|i\in 2U\right\\}$
$\displaystyle=\frac{1}{2}\pi_{ci}p_{ci}+\frac{1}{2}\pi_{ri}p_{ri}.$

Finally, by the definition of conditional probability,

$\displaystyle P\left\\{i\in S_{c}|i\in S,i\in 2U\right\\}$
$\displaystyle=\frac{P\left\\{i\in S_{c}|i\in 2U\right\\}}{P\left\\{i\in
S|i\in 2U\right\\}}.$ (8)

Formula (4) directly follows from putting (5) and (2.3) in formula (8).
∎

In the special case where $U_{c}=U_{r}$ and $p_{ci}=p_{ri}$, the CRISP formula
(4) reduces to:

$\pi_{zi}=\frac{\pi_{ci}}{\pi_{ci}+\pi_{ri}}.$ (9)

######  Remark 1.

In the approach of Elliott (2009) and Elliott and Valliant (2017) expression
(9) is an _approximation_ under the assumption of a negligible overlap between
samples $S_{c}$ and $S_{r}.$ Note, that unless sampling inclusion
probabilities are small, the probability of a unit being included in both
samples is not trivial and should be accounted for. Under our setup, this
expression _exactly_ holds without assuming a small overlap.

######  Remark 2.

The stacking devise given in Figure 1 is useful not only as a “trick” for the
derivation of formula (4). In usual practical situations, samples come from
different sources and we _do not observe_ which units belong to both the
reference and non-probability samples and thus we would simply stack the two
samples together without matching. Accounting for the probability of the
overlap would contradict with such a setup. However, if information on the
identity of the overlapping units is available to researches, then it would be
advisable to use this additional information and consider the union,
$S_{c}\cup S_{r},$ of the two samples instead of the stacked sample setup. The
corresponding formula for $\pi_{zi}$ would require a simple adjustment to
account for the probability of the overlap, namely, by subtracting
$\pi_{ci}\pi_{ri}$ in the denominator of (9):

$\pi_{zi}=\frac{\pi_{ci}}{\pi_{ci}+\pi_{ri}-\pi_{ci}\pi_{ri}}.$ (10)

We do not consider this setup in the current paper.

######  Remark 3.

We focus inference on estimation of participation probabilities, $\pi_{ci}$,
under known coverage probabilities, $(p_{ri},p_{ci})$, using our CRISP result
of Equation (4). We make brief mention that our CRISP formula may also be used
in other sampling situations. For example, we may define $(p_{ri},p_{ci})$ as
first stage participation probabilities in a multi-stage sampling design
(rather than viewing them as coverage probabilities) to construct a sampling
frame. In this set-up, $(\pi_{ri},\pi_{ci})$, are second stage inclusion
probabilities for units selected from the realized frame. The first stage
probabilities may be _unknown_ , but once the sampling frame is defined, the
second stage inclusion probabilities, by contrast, can be set and _known_. In
this set-up, the focus for inference would be the first stage participation
probabilities in the sampling frame, $(p_{ri},p_{ci})$, rather than
$(\pi_{ri},\pi_{ci})$.

One class of examples of such a multi-stage participation and selection design
belongs to the opt-in panels recruited over the internet. The probability of
signing up for a web-based panel (“the first stage”) is not known. On “the
second stage”, once the panel has been recruited, a sample – with _known_
inclusion probabilities – would be selected for a given survey.

See also related discussion in Elliott and Valliant (2017) and Valliant
(2020).

Let us now consider joint inclusion probabilities of units $i,j\in U$, $i\neq
j$, conditional on covariates. Denote

$\displaystyle p_{c,ij}=p_{c}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle=P\\{i,j\in U_{c}|i,j\in U,\mathbf{x}_{i},\mathbf{x}_{j}\\},$
$\displaystyle p_{r,ij}=p_{r}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle=P\\{i,j\in U_{r}|i,j\in U,\mathbf{x}_{i},\mathbf{x}_{j}\\},$
$\displaystyle\pi_{c,ij}=\pi_{c}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle=P\\{i,j\in S_{c}|i,j\in
U_{c},\mathbf{x}_{i},\mathbf{x}_{j}\\},$
$\displaystyle\pi_{r,ij}=\pi_{r}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle=P\\{i,j\in S_{r}|i,j\in
U_{r},\mathbf{x}_{i},\mathbf{x}_{j}\\}.$

For $i,j\in S$, $i\neq j,$ denote
$\pi_{z,ij}=\pi_{z}(\mathbf{x}_{i},\mathbf{x}_{j})=P\\{i,j\in S_{c}|i,j\in
S,\mathbf{x}_{i},\mathbf{x}_{j}\\}.$

######  Theorem 2.2 (Independence of $I_{zi}$ under the stacked sample setup).

Assume, inclusion indicator variables in each sampling stage are independent
given $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ for $i\neq j$:
$p_{c,ij}=p_{ci}p_{cj},$ $p_{r,ij}=p_{ri}p_{rj},$
$\pi_{c,ij}=\pi_{ci}\pi_{cj},$ $\pi_{r,ij}=\pi_{ri}\pi_{rj}.$

Then indicator variables $I_{zi}$, $I_{zj}$ are independent given
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$: $\pi_{z,ij}=\pi_{zi}\pi_{zj}$ for
$i,j\in S$, $i\neq j.$

The proof is analogous to Theorem 2.1, it is given in Appendix A.

######  Remark 4.

The joint conditional independence of inclusion indicators $I_{ci}$ and
$I_{cj}$ stated in Assumption A3: is a serious supposition that, nevertheless,
is reasonable given conditioning on a rich set of covariates. The same
reasoning can be applied to the assumption of conditional independence of the
reference sample inclusion indicators $I_{ri}$ and $I_{rj}$ given
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}.$

###  2.4 Modeling using CRISP based on valid likelihood for the observed
sample

Beresovsky (2019) proposed using CRISP expression (4) to obtain estimates of
$\pi_{ci}$ directly from a likelihood formulated on the combined sample. This
approach is labeled “Implicit Logistic Regression (ILR)” because parameters
are entering the logistic link function “implicitly” through a composite
function.

For simplicity and coherence with competitor approaches, in what follows we
consider the special case scenario corresponding to the formula of (9).

Participation probabilities $\pi_{ci}=\pi_{ci}(\boldsymbol{\beta})$ are
parameterized as in the CLW approach by
$\text{logit}(\pi_{ci}(\boldsymbol{\beta}))=\boldsymbol{\beta^{T}}\mathbf{x}_{i},$
and identity (9) is used to present $\pi_{zi}$ as a composite function of
$\boldsymbol{\beta}$; that is,
$\pi_{zi}=\pi_{zi}(\pi_{ci}(\boldsymbol{\beta}))=\pi_{ci}(\boldsymbol{\beta})/(\pi_{ri}+\pi_{ci}(\boldsymbol{\beta})).$
The log-likelihood for observed Bernoulli variable $I_{zi}$ is

$\displaystyle\ell^{ILR}(\boldsymbol{\beta})=\mathop{\sum}_{i\in
S_{c}}\log\left[{\pi_{zi}(\pi_{ci}(\boldsymbol{\beta}))}\right]+\mathop{\sum}_{i\in
S_{r}}\log\left[1-\pi_{zi}\left(\pi_{ci}(\boldsymbol{\beta})\right)\right].$
(11)

This approach may be labeled “indirect” estimation because $\pi_{c}$ is
indirectly used through the CRISP in Bernoulli likelihood formulated for the
_observed_ indicators $I_{zi}$ of convenience sample on the “stacked” set
$S=S_{c}+S_{r}$.

######  Remark 5.

Note, by Theorem 2.2, indicators $I_{zi}$ are _independent_ for units on $S$,
even if samples $S_{c}$ and $S_{r}$ have common units. We use exact likelihood
for independent Bernoulli variables. Thus, standard statistical properties and
techniques can be applied. For example, we can apply usual likelihood based
model selection methods such as AIC or BIC.

The score equations are obtained from (11) by taking the derivatives, with
respect to $\boldsymbol{\beta},$ of composite function
$\pi_{zi}=\pi_{zi}(\pi_{ci}(\boldsymbol{\beta})).$ This way, the estimates of
$\pi_{ci}$ are obtained directly from (11) in a single step. The resulting
score equations are

$\displaystyle\boldsymbol{0}=S^{ILR}(\boldsymbol{\beta})=\frac{\partial\ell^{ILR}(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}}=$
(12) $\displaystyle\mathop{\sum}_{i\in
S_{c}}(1-\pi_{zi})(1-\pi_{ci})\mathbf{x}_{i}-\mathop{\sum}_{i\in
S_{r}}\pi_{zi}(1-\pi_{ci})\mathbf{x}_{i}.$

###  2.5 Pseudo-ILR method

Wang, Valliant and Li (2021) also developed an indirect method for estimation
of $\pi_{ci}$ as a function of $\pi_{zi}$, which we reveal is a special case
of (9). For their Adjusted Logistic Propensity (ALP) weighting method, Wang,
Valliant and Li (2021) use a stacking approach as do we. They introduce an
imaginary construct consisting of two parts: they _stack_ together non-
probability sample $S_{c}$ (part 1) and finite population $U$ (part 2). Since
non-probability sample units belong to the finite population, they appear in
the stacked set twice. They formulate a Bernoulli likelihood for variable
$I_{{\delta}i}$, where $I_{{\delta}i}=1$ if unit $i$ belongs to part 1 of the
stacked set; and $I_{{\delta}i}=0$, otherwise:

$\displaystyle\ell^{ALP}(\boldsymbol{\gamma})=\mathop{\sum}_{i\in
S_{c}}\log\left[{\pi_{{\delta}i}(\boldsymbol{\gamma})}\right]+\mathop{\sum}_{i\in
U}\log\left[1-\pi_{{\delta}i}\left(\boldsymbol{\gamma}\right)\right],$ (13)

where $\boldsymbol{\gamma}$ is the parameter vector in a logistic regression
model for
$\pi_{{\delta}i}(\mathbf{x}_{i})=P\\{I_{{\delta}i}=1|\mathbf{x}_{i}\\}.$ At
the second step of the ALP method, estimates of $\pi_{ci}$ are derived from
identity

$\pi_{{\delta}i}=\frac{\pi_{ci}}{\pi_{ci}+1}.$ (14)

Their identity may also be derived from our result in (9) by replacing
$\pi_{ri}$ with $1$ because it is assumed that the entire finite population is
observed.

Since the finite population is not available, they apply a pseudo-likelihood
approach:

$\displaystyle\hat{\ell}^{ALP}(\boldsymbol{\gamma})=\mathop{\sum}_{i\in
S_{c}}\log\left[{\pi_{{\delta}i}(\boldsymbol{\gamma})}\right]+\mathop{\sum}_{i\in
S_{r}}w_{ri}\log\left[1-\pi_{{\delta}i}\left(\boldsymbol{\gamma}\right)\right],$
(15)

leading to an estimate of $\pi_{{\delta}i}.$ However, the actual goal is to
find probabilities $\pi_{ci}$ rather than $\pi_{{\delta}i}.$ On the second
step of the ALP approach, Wang, Valliant and Li (2021) obtain the estimate of
$\pi_{ci}$ using (14).

Wang, Valliant and Li (2021) noted that in their simulations the ALP estimator
was more efficient than the CLW, especially in the case of large non-
probability and small probability sample sizes.

The estimation method of Wang, Valliant and Li (2021) can be modified to a
one-step estimation procedure similar to ILR: $\pi_{ci}$ can be parameterized
using the logistic link function as
$\text{logit}(\pi_{ci}(\boldsymbol{\beta}))=\boldsymbol{\beta^{T}}\mathbf{x}_{i}$,
while probabilities $\pi_{{\delta}i}$ in (14) could be viewed as a composite
function,
$\pi_{{\delta}i}=\pi_{{\delta}i}(\pi_{ci}(\boldsymbol{\beta}))=\pi_{ci}(\boldsymbol{\beta})/(1+\pi_{ci}(\boldsymbol{\beta}))$.
We call this method pseudo-ILR (PILR), and the pseudo-likelihood is

$\displaystyle\hat{\ell}^{PILR}(\boldsymbol{\beta})=\mathop{\sum}_{i\in
S_{c}}\log\left[{\pi_{{\delta}i}(\pi_{ci}(\boldsymbol{\beta}))}\right]+\mathop{\sum}_{i\in
S_{r}}w_{ri}\log\left[1-\pi_{{\delta}i}(\pi_{ci}(\boldsymbol{\beta}))\right].$
(16)

The score equations are:

$\displaystyle\boldsymbol{0}=S^{PILR}(\boldsymbol{\beta})=\frac{\partial\hat{\ell}^{PILR}(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}}=$
$\displaystyle\mathop{\sum}_{i\in
S_{c}}(1-\pi_{{\delta}i})(1-\pi_{ci})\mathbf{x}_{i}-\mathop{\sum}_{i\in
S_{r}}w_{ri}\pi_{{\delta}i}(1-\pi_{ci})\mathbf{x}_{i}.$

This change in estimation of model parameters makes the approach more
efficient and it avoids cases where estimates of $\pi_{ci}$ become greater
than $1$, as may occur under the ALP approach where the estimation is
performed in two steps.

######  Remark 6.

Similar to Remark 5, by Theorem 2.2, indicators $I_{\delta i}$ and $I_{\delta
j}$ are independent given $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ for $i\neq j$
in the stacked $S_{c}$ and $U$ set, even though each $S_{c}$ unit is included
into the stacked set twice, once with $I_{\delta i}=1$ and once with
$I_{\delta i}=0.$ Thus, the likelihood for independent Bernoulli variables is
well justified for the stacked $S_{c}$ and $U$ set.

######  Remark 7.

In any of these approaches, there is no need to restrict the parameterization
to a linear form and one could parameterize the likelihood as
$\text{logit}(\pi_{ci}(\boldsymbol{\beta}))=\psi(\mathbf{x}_{i};\boldsymbol{\beta})$
for some general function $\psi.$ For example, Savitsky et al. (2022) used
splines to estimate the model parameters.

###  2.6 Inverse probability weighted estimator

In the sequel we consider the following Hájek form of the inverse probability
weighted estimator of finite population mean $\mu$:

$\hat{\mu}=\frac{1}{\hat{N}}\mathop{\sum}_{i\in
S_{c}}\frac{y_{i}}{\hat{\pi}_{ci}},$ (18)

where $\hat{N}=\sum_{i\in S_{c}}(\hat{\pi}_{ci})^{-1}$ and $\hat{\pi}_{ci}$ is
an estimate of ${\pi}_{ci}$.

Note that the Hájek estimator of parameter $\mu$ can be presented as a
solution of the following estimating equation:

$U(\mu)=\mathop{\sum}_{i\in U}I_{ci}\hat{\pi}_{ci}^{-1}(y_{i}-\mu)=0.$ (19)

##  3 Theoretical properties

Theoretical properties of the CLW estimators are given in Chen, Li and Wu
(2020). In Section 3.1, we present asymptotic variances of the estimated
parameters and rates of convergence to true values for both of the ILR and
PILR. We also include the CLW variance to facilitate our comparison of the
relative efficiencies of these methods. In Section 3.2, we discuss the
expected behavior of the asymptotic variances for ILR, PLR and CLW. Section
3.3 presents a numerical study of relative efficiency of the estimators for
the methods over a sequence of sampling fractions.

###  3.1 Asymptotic variances of CLW, ILR, and PILR estimators

We study asymptotic properties of the estimators using the following
framework. Consider a sequence of finite populations $U_{\nu}$, indexed by
$\nu$, of increasing sizes $N_{\nu}$, so that $N_{\nu}\to\infty$ as
$\nu\to\infty,$ and let non-probability sample $S_{c\nu}$ of size $n_{c\nu}$
and reference probability sample $S_{r\nu}$ of size $n_{r\nu}$ be
independently drawn from population $U_{\nu}.$ The probability sample design
is assumed to be the same for each $U_{\nu}$ and assume that both sample sizes
increase $n_{c\nu}\to\infty$ and $n_{r\nu}\to\infty$ as $\nu\to\infty$ in such
a way that sampling fractions $f_{c}=\lim_{\nu\to\infty}n_{c\nu}/N_{\nu}$ and
$f_{r}=\lim_{\nu\to\infty}n_{r\nu}/N_{\nu}.$ The same limiting process in
considered by Chen, Li and Wu (2020). As is customary, to simplify notation,
we drop index $\nu$ in the rest of the paper.

We first present the asymptotic variances of the estimates of parameters $\mu$
and $\boldsymbol{\beta}$ in a general form applicable for all three methods.
Such a representation aids in comparing methods:

$\displaystyle Var(\hat{\mu})$ $\displaystyle\doteq$
$\displaystyle\frac{1}{N^{2}}\\{Var[U(\mu)]-2\boldsymbol{b}^{T}Cov[U(\mu),{S}(\boldsymbol{\beta})]+\boldsymbol{b}^{T}Var[{S}(\boldsymbol{\beta})]\boldsymbol{b}\\},$
(20) $\displaystyle Var(\boldsymbol{\hat{\beta}})$ $\displaystyle\doteq$
$\displaystyle\boldsymbol{H}^{-1}Var[{S}(\boldsymbol{\beta})]\boldsymbol{H}^{-1},$
(21)

where $U(\mu)$ is the score function for $\mu$, ${S}(\boldsymbol{\beta})$ are
the score functions for parameters $\boldsymbol{\beta}$ in a given model for
$\pi_{ci}$, and
$\boldsymbol{b}=S_{\boldsymbol{\beta}}^{-1}U_{\boldsymbol{\beta}}^{T}$,
$U_{\boldsymbol{\beta}}=E[{\partial
U(\mu)}/{\partial\boldsymbol{\beta}^{T}}]$,
$S_{\boldsymbol{\beta}}=E[{\partial
S(\boldsymbol{\beta})}/{\partial\boldsymbol{\beta}}]$,
$\boldsymbol{H}=-S_{\boldsymbol{\beta}}.$ See Appendix B for the derivation.

The first term in (20), $Var[U(\mu)],$ is the design-based variance under the
probability sampling design for the convenience sample, and the other two
terms are related to the estimation of unknown participation probabilities
$\pi_{ci}$ involved in the probability weighted estimator for $\mu.$
Importantly, (20) contains variance $Var[{S}(\boldsymbol{\beta})]$ due to
estimation of parameters $\boldsymbol{\beta},$ which shows that efficiency of
estimation of $\mu$ depends on the quality of the model for $\pi_{ci}.$ The
score function for $\boldsymbol{\beta}$ can be further decomposed into two
independent parts, corresponding to the non-probability and reference samples:

$\displaystyle{S}(\boldsymbol{\beta})={S_{c}}(\boldsymbol{\beta})+{S_{r}}(\boldsymbol{\beta}),$
(22)

and the corresponding variance is

$\displaystyle
Var[S(\boldsymbol{\beta})]=Var[S_{c}(\boldsymbol{\beta})]+Var[S_{r}(\boldsymbol{\beta})]=:\mathbf{A}+\mathbf{D}.$
(23)

Unless the reference sample is large, the design variance
$\mathbf{D}=Var[S_{r}(\boldsymbol{\beta})]$ is a significant contributor into
the overall variance. In Section 3.2 below, we discuss and compare the effect
of the design-based variance for each of the three methods.

We now follow with a formal statement that spells out detailed formula for
each of the methods, CLW, ILR, and PILR. Assume the following regularity
conditions:

  1. C1

Matrix $\sum_{i\in U}\mathbf{x}_{i}\mathbf{x}_{i}^{T}$ is positive-definite.

  2. C2

There exist $c_{2}\geq c_{1}>0$ such that $c_{1}\leq\pi_{ci}N/n_{c}\leq c_{2}$
and $c_{1}\leq\pi_{ri}N/n_{r}\leq c_{2}$ for all units $i.$

######  Theorem 3.1 (Consistency of the estimators).

Assume logistic regression parameterization for the non-probability sample
participation probabilities,
$\text{logit}\left[\pi_{ci}(\boldsymbol{\beta})\right]=\boldsymbol{\beta^{T}}\mathbf{x}_{i}.$
Then, under regularity conditions C1-C2, estimators
$\hat{\boldsymbol{\eta}}_{m}=(\hat{\mu}_{m},\hat{\boldsymbol{\beta}}_{m})$ for
$\boldsymbol{\eta}=({\mu},{\boldsymbol{\beta}})$ based on CLW, ILR, or PILR
methods are design consistent,
$\hat{\boldsymbol{\eta}}_{m}-\boldsymbol{\eta}=O_{p}(\text{min}(n_{c},n_{r})^{-1/2})$,
and variances are

$\displaystyle Var(\hat{\mu}_{m})$ $\displaystyle\doteq$ $\displaystyle
N^{-2}\mathop{\sum}_{i\in U}\frac{1-\pi_{ci}}{\pi_{ci}}(y_{i}-\mu)^{2}$
$\displaystyle-2N^{-2}\boldsymbol{b}_{m}^{T}\boldsymbol{C}_{m}+N^{-2}\boldsymbol{b}_{m}^{T}(\boldsymbol{A}_{m}+\boldsymbol{D}_{m})\boldsymbol{b}_{m},$
$\displaystyle Var(\boldsymbol{\hat{\beta}}_{m})$ $\displaystyle\doteq$
$\displaystyle\boldsymbol{H}_{m}^{-1}(\boldsymbol{A}_{m}+\boldsymbol{D}_{m})\boldsymbol{H}_{m}^{-1},$
(25)

where subscript $m$ stands for methods (CLW, ILR, or PILR),
$\boldsymbol{b}_{m}=\boldsymbol{H}_{m}^{-1}\sum_{i\in
U}(1-\pi_{ci})(y_{i}-\mu)\mathbf{x}_{i},$ and the other terms are spelled out
as follows:

For the CLW method,

$\displaystyle\boldsymbol{C}_{CLW}$ $\displaystyle=\sum_{i\in
U}(1-\pi_{ci})(y_{i}-\mu)\mathbf{x}_{i},$ $\displaystyle\boldsymbol{H}_{CLW}$
$\displaystyle=\boldsymbol{A}_{CLW}=\sum_{i\in
U}\pi_{ci}(1-\pi_{ci})\mathbf{x}_{i}\mathbf{x}_{i}^{T},$

and

$\displaystyle\boldsymbol{D}_{CLW}=Var_{d}\left[\sum_{i\in
S_{r}}w_{ri}\pi_{ci}\mathbf{x}_{i}\right]$

is the design-based variance-covariance matrix under the probability sampling
design for $S_{r}$.

For the ILR method,

$\displaystyle\boldsymbol{C}_{ILR}=\mathop{\sum}_{i\in
U}(y_{i}-\mu)(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i},$
$\displaystyle\boldsymbol{H}_{ILR}=\sum_{i\in
U}(\pi_{ci}+\pi_{ri})\pi_{zi}(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T},$
$\displaystyle\boldsymbol{A}_{ILR}=\boldsymbol{H}_{ILR}-\sum_{i\in
U}(\pi_{ri}+1)\pi_{ci}\pi_{zi}(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T},$

and

$\displaystyle\boldsymbol{D}_{ILR}=Var_{d}\left[\mathop{\sum}_{i\in
S_{r}}\pi_{zi}(1-\pi_{ci})\mathbf{x}_{i}\right]$

is the design-based variance-covariance matrix under the probability sampling
design for $S_{r}$.

For the PILR method,

$\displaystyle\boldsymbol{C}_{PILR}=\mathop{\sum}_{i\in
U}(y_{i}-\mu)(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i},$
$\displaystyle\boldsymbol{H}_{PILR}=\sum_{i\in
U}(\pi_{ci}+1)\pi_{{\delta}i}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T},$
$\displaystyle\boldsymbol{A}_{PILR}=\boldsymbol{H}_{PILR}-2\mathop{\sum}_{i\in
U}\pi_{ci}\pi_{{\delta}i}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T},$

and

$\displaystyle\boldsymbol{D}_{PILR}=Var_{d}\left[\sum_{i\in
S_{r}}w_{ri}\pi_{ci}(1-\pi_{ci})(1-\pi_{{\delta}i})\mathbf{x}_{i}\right]$

is the design-based variance-covariance matrix under the probability sampling
design for $S_{r}$.

The details of the proof for ILR and PILR are given in Appendix B.

###  3.2 Analysis of the relative efficiency of the methods

In this Section, we take a closer look at the components of variances
presented in Section 3.1.

Consider a situation where the whole population $U$ is available in place of a
reference sample (that is, sampling fraction $f_{r}=n_{r}/N=1$). Since
Bernoulli likelihood (2) used in the CLW method describes the true data-
generating process at the population level, the CLW approach is the preferable
one for estimation in the case where the whole $U$ is observed. Under this
scenario, design-based variance-covariance matrices $\boldsymbol{D}_{CLW},$
$\boldsymbol{D}_{ILR}$, and $\boldsymbol{D}_{PILR}$ are reduced to zero
matrices, reflecting the fact that variability due to the probability sampling
does not contribute to respective total variances of parameters ${\mu}$ and
${\boldsymbol{\beta}}.$ In particular, the variance of parameter
$\boldsymbol{\beta}$ given by (25) would be reduced to
$\boldsymbol{H}_{CLW}^{-1}$ and equal to the variance of the maximum
likelihood estimator from Bernoulli likelihood (2). In this special case, the
ILR and PILR methods are identical and yield estimators that are less
efficient than CLW estimators.

Let us now turn to a more common situation where probability sample $S_{r}$
with sampling fraction $f_{r}=n_{r}/N<1$ is available rather than the whole
population $U.$ To a large degree, respective contributions of design
variances $\boldsymbol{D}_{CLW},$ $\boldsymbol{D}_{ILR}$, and
$\boldsymbol{D}_{PILR}$ determine the size of the total variances. Consider

$\displaystyle\boldsymbol{D}_{CLW}$
$\displaystyle={Var_{d}}\left[\sum\nolimits_{S_{r}}{{g_{i}}{{\mathbf{x}}_{i}}}\right],$
(26) $\displaystyle\boldsymbol{D}_{ILR}$
$\displaystyle={Var_{d}}\left[\sum\nolimits_{S_{r}}{{\frac{g_{i}}{1+g_{i}}}\left(1-{{\pi}_{ci}}\right){{\mathbf{x}}_{i}}}\right],$
(27) $\displaystyle\boldsymbol{D}_{PILR}$
$\displaystyle={Var_{d}}\left[\sum\nolimits_{S_{r}}{{g_{i}}\frac{1-{{\pi}_{ci}}}{1+{{\pi}_{ci}}}{{\mathbf{x}}_{i}}}\right],$
(28)

where
$g_{i}=g\left(\mathbf{x}_{i}\right)={\pi_{c}\left(\mathbf{x}_{i}\right)}/{\pi_{r}\left(\mathbf{x}_{i}\right)}\in\left(0,\infty\right)$.
Variation in ratio $g_{i}$ depends on the degree of common support, or
overlap, between two samples in the covariate-defined domains. In the case of
“high overlap”, $\pi_{ci}$ and $\pi_{ri}$ would have “similar dependence” on
$\mathbf{x}_{i}$, and there will be relatively minor fluctuations of ratios
$g_{i}$. However, in the case of “low overlap”, $\pi_{ci}$ and $\pi_{ri}$
would have an “opposite dependence” on $\mathbf{x}_{i}$ and there would be a
significant variation of $g_{i}.$

In context of the observational studies, Firth (1993) and Heinze and Schemper
(2002) considered estimation of parameters of generalized linear models for
the case of small overlap, or “separation”, between treated and control units
in variable-defined domains. They showed that separation in variable-defined
domains may result in biased and unstable estimates of model parameters.

Expressions (26) - (28) indicate that susceptibility to the effect of a low
overlap in variable-defined domains varies between the three methods.
Expression (26) for $\boldsymbol{D}_{CLW}$ includes $g_{i}$ as a factor that
potentially can vary in the whole interval $\left(0,\infty\right)$. Expression
(28) defining design variance $\boldsymbol{D}_{PILR}$ differs from the CLW
expression by factor $(1-\pi_{ci})/(1+\pi_{ci})<1$, that could soften the
effect from the variability of $g_{i}$. In expression (27) for
$\boldsymbol{D}_{ILR},$ factor ${g_{i}}/(1+g_{i})$ varies between 0 and 1,
thus indicating lower susceptibility of the ILR method to the low overlap
problem compared to CLW or PILR.

This quick analysis demonstrates that low overlap in covariate-defined domains
may affect stability of the estimates, and the impact varies between the three
methods.

The differences in estimated variances between the three methods are expected
to be more pronounced for the situation of a large non-probability sampling
fraction coupled with a small probability sampling fraction. This situation is
encountered in the practical setting of a large-sized administrative data for
the non-probability / convenience sample and a small reference sample obtained
from a probability survey. In the sequel, this expectation is confirmed in
numerical and simulation studies.

###  3.3 Numerical study using theoretical quantities

We present a numerical study of this influence of sampling fractions for
reference and convenience samples, on the one hand, on the quality of
estimates for the $3$ methods, on the other hand, by computing variances for
CLW, ILR, and PILR on a grid of values of sampling fractions
$\left(f_{r},f_{c}\right)=\left(n_{r}/N,n_{c}/N\right)$ under the low and high
overlap scenarios.

For each unit $i=1\ldots,N$ in finite population $U$ of size $N=100,000$, we
generate covariate $x_{i}$ as independent standard normal variable, $x_{i}\sim
N\left(0,1\right)$. Outcome $y_{i}$ is generated as a normal variable with
mean depending on $x_{i}$ and constant variance, $y_{i}\sim
N\left(1+x_{i},1.5^{2}\right)$.

We use a Poisson sampling design with participation probabilities $\pi_{ci}$
to select convenience sample $S_{c}$ from population $U$. Probabilities
$\pi_{ci}$ are generated depending on $x_{i}$ as

$\text{logit}(\pi_{ci})=\beta_{c0}+\beta_{c1}x_{i}.$ (29)

We vary values of intercept $\beta_{c0}$ to scan across a range of sampling
fractions $f_{c}\in(0,1)$. The value of slope $\beta_{c1}$ is set to $1.$

We select reference sample $S_{r}$ using the probability proportional to size
(PPS) sampling with the measure of size $m_{ri}$ generated as

$\text{logit}(m_{ri})=1+\beta_{r}x_{i}.$ (30)

For any size measures $m_{r}\left(x_{i}\right)$, overall sampling fraction
$f_{r}$ can be tuned to a desired value using functionality of sampling
R-package by Tille and Matei (2021).

For the low overlap scenario, we simulate the “opposite dependence” of
$\pi_{ci}$ and $\pi_{ri}$ on covariate $x_{i}$ by setting $\beta_{r}=-1.$ For
the high overlap scenario, we generate “similar dependence” by setting
$\beta_{r}=1,$ so that it is the same as slope $\beta_{c1}$ in the generating
model for $\pi_{ci}$. Figure 2 shows the resulting histogram of sample
densities in domains of covariate $x_{i}.$ Figure 2 helps to visualize and
clarify what is meant by low and high overlap, in a simple case of one
covariate. In the case of multivariate models with a large number of auxiliary
variables and interactions, it is likely to have very little overlap in some
of the variable-defined domains, even when samples are large. A similar
phenomenon in machine learning is called “the curse of dimensionality” and it
is recognized as a reason for unstable and not reproducible predictions.

Figure 2:  High overlap: $\beta_{r}=\beta_{c1}$; Low overlap:
$\beta_{r}=-\beta_{c1}$ Figure 3:  Relative standard error of ILR and CLW
estimates of propensity model parameter $\beta_{c1}$ Figure 4:  Relative
standard error of ILR and CLW estmates of population mean $\mu$ Figure 5:
Relative standard error of ILR and PILR estimates of propensity model
parameter $\beta_{c1}$ Figure 6:  Relative standard error of ILR and PILR
estmates of population mean $\mu$

We use formulas of Section 3.1 to compute variances of estimates ${\mu}$ and
${\boldsymbol{\beta}},$ for each considered method on the grid of sampling
fractions $f_{r}\in\left[0.02,1\right]$ and $f_{c}\in\left[0.01,0.85\right]$,
for high and low overlap in variable-defined domains. Figures 3 and 4 show
ratios of standard errors of ILR and CLW estimates of $\mu$ and $\beta_{c1}$.
Similarly, Figures 5 and 6 show ratios of standard errors of ILR and PILR
estimates of $\mu$ and $\beta_{c1}$. Notice that design variance
$\boldsymbol{D}_{m}$ enters variance formula (3.1) and (25) for both estimates
of $\mu$ and $\boldsymbol{\beta}$. This explains similar dependencies of
standard error ratios on sampling fractions $f_{c}$ and $f_{r}$ in Figures for
estimates of $\mu$ and $\beta_{c1}$.

Figure 3 demonstrates that CLW is the most efficient (lowest standard error)
method for estimating response propensity when sampling fractions $f_{r}$
approach $1$. In this limit CLW pseudo-likelihood (2) practically matches
exact Bernoulli likelihood (2) for non-probability sample indicator $I_{ci}$,
which explains CLW advantage over the other methods. This advantage is not
observed for the estimates of population means in Figure 4.

Relative efficiency of the CLW method disappears for smaller sampling fraction
$f_{r}$, which is the case for majority of probability surveys used in
statistical offices. This effect is particularly pronounced for the low
overlap scenario. The CLW method becomes increasingly inefficient relative to
ILR when non-probability fraction $f_{c}$ is large while probability sampling
fraction $f_{r}$ is small. This scenario is practically relevant when large
administrative data, usually viewed as convenience sample, collects
information which can complement and improve estimates from a relatively small
probability survey.

##  4 Simulation study

We conduct a Monte Carlo simulation study to compare the relative performances
of the CLW, ILR, and PILR estimates under various scenarios for sample sizes,
sampling fractions, and degrees of overlap of shared covariates of the
convenience and reference samples, that allow us to validate our theoretical
results in (26), (27) and (28). A secondary purpose of our simulation study is
to demonstrate the validity of our theory-derived, plug-in variance estimators
of (3.1) and (25) for variances of the estimated parameters $\mu$ and
$\boldsymbol{\beta}$. Our simulation study also demonstrates that inference
for ILR remains robust no matter the degree of overlap in units between the
reference and convenience samples in the case of independent sampling into
each.

###  4.1 Simulation setup

Table 1 outlines the various scenarios utilized in a Monte Carlo (MC)
simulation study. In each scenario, convenience and reference samples are
drawn $A=1,000$ times from the _same_ population. The convenience sample
selection probabilities were generated using expression (29) with true
coefficient value $\beta_{c1}=1.$ The outcome variable $y_{i}$ is generated
from the normal distribution as $y_{i}\sim N\left(1+x_{i},1.5^{2}\right)$,
where $x_{i}$ is generated as the independent standard normal variate, so that
mean $\mu=1$.

Each of the scenarios is implemented under low and high overlap settings, as
described in Section 3.2 and illustrated in Figure 2. Sampling fractions
$f_{c}$ are made equal to approximately $0.10$ or $0.01$ by setting intercept
$\beta_{c0}$ in the non-probability survey participation probabilities model
(29) to, respectively, $-2.5$ or $-5.0$, while keeping the model coefficient
$\beta_{c1}=1.$ Reference samples with fractions $f_{r}$ equal to $0.01$ and
$0.10$ are selected using the probability proportional to size (PPS) without
replacement sampling using sampling R-package by Tille and Matei (2021). The
large convenience sample coupled with selecting the whole finite population as
the reference sample in Scenario S7 in Table 1 allows us to study the
“ultimate overlap” case in which all units of the large non-probability sample
appear in the stacked set twice. This scenario is implemented by setting
intercept $\beta_{c0}=0$.

Table 1: Summary of considered simulation scenarios. | Sampling | Sample
---|---|---
| fraction | size
Scenario | $f_{c}$ | $f_{r}$ | $n_{c}$ | $n_{r}$
S1 | 0.01 | 0.01 | 600 | 600
S2 | 0.01 | 0.01 | 100 | 100
S3 | 0.10 | 0.10 | 600 | 600
S4 | 0.10 | 0.10 | 100 | 100
S5 | 0.01 | 0.10 | 100 | 1,000
S6 | 0.10 | 0.01 | 1,000 | 100
S7 | 0.50 | 1.00 | 500 | 1,000

###  4.2 Simulation results

Table 2 reports root mean squared errors (RMSE) for estimators of parameters
${\beta}_{c1}$ and ${\mu}$ for scenarios S1-S6 listed in Table 1. The top
three lines of the table report results for $\hat{\beta}_{c1}$ and the bottom
three lines show results for $\hat{\mu}$, for each of the methods. The left
half of the table provides results under the high overlap in covariate-defined
domains, and the right half of the table shows results under the low overlap.
The RMSEs are computed as
$\text{RMSE}(\hat{\mu})=\sqrt{{A^{-1}}\sum\nolimits_{\alpha=1}^{A}(\hat{\mu}^{(\alpha)}-\mu)^{2}}$
and
$\text{RMSE}(\hat{\beta}_{c1})=\sqrt{{A^{-1}}\sum\nolimits_{\alpha=1}^{A}(\hat{\beta}_{c1}^{(\alpha)}-{\beta_{c1}})^{2}}$,
and $\hat{\mu}^{(\alpha)},\hat{\beta}_{c1}^{(\alpha)}$ are estimates of
parameters at simulation run $\alpha$ based on respective methods.

The ILR approach provides better results compared to the other methods, while
the PILR approach gives better results compared to the CLW for both
$\hat{\mu}$ and $\hat{\beta}_{c1}$ in almost every scenario considered. The
results are closer for the high overlap cases, and the differences become more
pronounced under the low overlap. As expected, for all estimators, RMSEs are
generally lower when the sample size is larger. We also note relatively large
values of RMSE in the CLW approach under the low overlap cases of scenarios
S2, S4, and S6 where the reference sample is small.

Table 2:  Root mean square error (RMSE) of estimated parameters ${\beta}_{c1}$ and $\mu$ under high and low overlap for scenarios S1-S6 of Table 1 | High overlap | Low overlap
---|---|---
| S1 | S2 | S3 | S4 | S5 | S6 | S1 | S2 | S3 | S4 | S5 | S6
$\hat{\beta}_{ILR}$ | 0.07 | 0.17 | 0.08 | 0.22 | 0.12 | 0.14 | 0.08 | 0.22 | 0.10 | 0.25 | 0.14 | 0.15
$\hat{\beta}_{PILR}$ | 0.08 | 0.22 | 0.08 | 0.22 | 0.12 | 0.16 | 0.18 | 0.46 | 0.13 | 0.31 | 0.17 | 0.26
$\hat{\beta}_{CLW}$ | 0.09 | 0.26 | 0.09 | 0.29 | 0.12 | 0.23 | 0.22 | 1.50 | 0.21 | 1.22 | 0.19 | 0.68
$\hat{\mu}_{ILR}$ | 0.13 | 0.33 | 0.12 | 0.27 | 0.30 | 0.14 | 0.13 | 0.32 | 0.12 | 0.29 | 0.32 | 0.13
$\hat{\mu}_{PILR}$ | 0.14 | 0.38 | 0.12 | 0.29 | 0.31 | 0.17 | 0.22 | 0.54 | 0.15 | 0.37 | 0.35 | 0.25
$\hat{\mu}_{CLW}$ | 0.15 | 0.41 | 0.14 | 0.33 | 0.31 | 0.24 | 0.26 | 0.68 | 0.23 | 0.64 | 0.37 | 0.60

We next examine a more detailed summary for scenarios S5 and S6 of uneven
sample sizes between convenience and reference samples in Table 3 (for
$\beta_{c1}$) and Table 4 (for $\mu$). In each MC iteration, we estimated
variances by using a plug-in variance estimator based on the asymptotic
variance formulas presented in Section 3.1. Let $\widehat{V}^{(\alpha)}$
denote variance estimate for a parameter ($\mu$ or $\beta_{c1}$) obtained at
MC iteration $\alpha$. It can be calculated from expressions (3.1) or (25) by:
1) plugging-in estimated values of participation probabilities
$\hat{\pi}_{ci}$; and 2) replacing summation over population with weighted
sums over non-probability and probability samples. As usual, sampling weights
are equal to inverse of estimated participation $\hat{\pi}_{ci}^{-1}$ and
selection $\pi_{ri}^{-1}$ probabilities .

Table 3:  Estimate of propensity parameter ${\beta}_{c1}$ for sample fractions $\left(f_{c},f_{r}\right)=\left(0.01,0.1\right)$ and $\left(0.1,0.01\right)$ (scenarios S5 and S6) | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S5: $n_{c}=100,n_{r}=1000$
ILR | 1.00 | 0.12 | 0.11 | 0.95 | 1.01 | 0.14 | 0.13 | 0.95
PILR | 1.00 | 0.12 | 0.12 | 0.95 | 1.03 | 0.16 | 0.15 | 0.94
CLW | 1.01 | 0.12 | 0.12 | 0.95 | 1.05 | 0.19 | 0.16 | 0.93
| Scenario S6: $n_{c}=1000,n_{r}=100$
ILR | 1.00 | 0.14 | 0.14 | 0.94 | 1.01 | 0.15 | 0.15 | 0.95
PILR | 1.01 | 0.16 | 0.16 | 0.95 | 1.05 | 0.26 | 0.25 | 0.94
CLW | 1.05 | 0.23 | 0.22 | 0.95 | 1.29 | 0.62 | 0.51 | 0.92
Table 4:  Estimate of population mean ${\mu}$ for sample fractions $\left(f_{c},f_{r}\right)=\left(0.01,0.1\right)$ and $\left(0.1,0.01\right)$ (scenarios S5 and S6) | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S5: $n_{c}=100,n_{r}=1000$
ILR | 1.07 | 0.30 | 0.29 | 0.87 | 1.03 | 0.32 | 0.28 | 0.88
PILR | 1.07 | 0.30 | 0.30 | 0.88 | 1.00 | 0.35 | 0.34 | 0.90
CLW | 1.07 | 0.30 | 0.31 | 0.88 | 0.98 | 0.37 | 0.37 | 0.91
| Scenario S6: $n_{c}=1000,n_{r}=100$
ILR | 1.02 | 0.14 | 0.14 | 0.94 | 1.01 | 0.13 | 0.13 | 0.95
PILR | 1.00 | 0.17 | 0.18 | 0.96 | 0.96 | 0.24 | 0.27 | 0.96
CLW | 0.97 | 0.23 | 0.29 | 0.95 | 0.75 | 0.54 | Inf | 0.95

The columns in both tables are defined, as follows:

  * •

“$Mean$” is the mean estimate of a parameter ($\mu$ or $\beta_{c1}$) over the
MC simulations;

  * •

“$SE$” is the standard error of the estimates computed over MC simulations, so
it is a repeated sampling-derived variance estimator;

  * •

“$\overline{\widehat{SE}}$” is the average over MC simulations computed based
on the plug-in variance estimates $\widehat{V}^{(\alpha)}$ for parameters
($\mu$ or $\beta_{c1}$). We first compute $\widehat{SE}$ as the square root of
the plug-in variance estimate for each MC simulation iteration
$\sqrt{\widehat{V}^{(\alpha)}}$. $\widehat{SE}$ is then used to construct
$95\%$ confidence intervals and to assess coverage. Finally, we construct
$\overline{\widehat{SE}}$ as the average
$\overline{\widehat{SE}}=\sqrt{A^{-1}\sum\nolimits_{\alpha=1}^{A}\widehat{V}^{(\alpha)}}$
over the MC simulation iterations to facilitate the comparison with $SE$.

  * •

“$95\%CI$” shows the percentage of times, over the MC simulation runs, the
true parameter is covered by the $95$ percent normal confidence intervals
constructed based on $\sqrt{\widehat{V}^{(\alpha)}}$.

We generally observe relatively high values of $\overline{\widehat{SE}}$ for
CLW and to a lesser degree for PILR. So, though the methods may achieve
nominal coverage, the resulting average of interval widths over the MC
iterations (which are related to $\overline{\widehat{SE}}$) are wider,
indicating that CLW and PILR estimators are less efficient than that of ILR.
This relative efficiency difference is particularly pronounced in the low
overlap setting.

We note that in the high overlap case in scenario S6, where the reference
sample size is relatively small, there is about 7% bias in the estimates of
$\mu$ based on all methods. Respective normal confidence intervals do not
provide nominal coverage of the true mean. At the same time, estimates of
$\beta_{c1}$ are nearly unbiased and equally efficient for all approaches. In
smaller-sized convenience samples nonlinear dependence of an estimator of
$\mu$ on estimated participation probabilities may result in biased estimates
of the mean. Since $\hat{\mu}$ is a consistent estimator, the bias would
decrease in larger samples. Bias of estimates for high overlap and small
convenience samples also can be seen on boxplots in Figure 12 of Appendix C.

Lastly, the very similar values of $SE$ and $\overline{\widehat{SE}}$ across
the $3$ approaches in Tables 3 and Table 4 demonstrate the validity of the
theoretically-derived plug-in variance estimator. In some cases in Table 4,
estimates of standard error deviate from $SE$ over the MC simulations. We
attribute such discrepancies to biased estimates of participation
probabilities, i.e. CLW estimate for Scenario S6 under high overlap, or to
nonlinear dependence of population mean $\mu$ on $\pi_{c}$ and small non-
probability sample size $n_{c}$ (ILR estimate for Scenario 5 under low
overlap). These cases need to be further studied for better explanation.

Table 5:  Estimates of propensity parameter ${\beta}_{c1}$ and population mean
${\mu}$ for $S_{r}=U$, non-probability sampling fraction $f_{c}=0.5$,
population size $N=1,000$ (scenario S7)

| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
---|---|---|---|---
$\hat{\beta}_{ILR/PILR}$ | 1.01 | 0.11 | 0.11 | 0.95
$\hat{\beta}_{CLW}$ | 1.01 | 0.09 | 0.09 | 0.94
$\hat{\mu}_{ILR/PILR}$ | 1.00 | 0.07 | 0.07 | 0.94
$\hat{\mu}_{CLW}$ | 1.00 | 0.08 | 0.08 | 0.94

Table 5 shows estimates for the “ultimate overlap” scenario S7 where the whole
finite population is available (and used as the reference sample), so that
units of the large convenience sample appear in the stacked set to two samples
twice. In this case, sampling fractions and sizes are large and there is also
a high overlap in the covariate-defined domains. When the whole population is
used as the reference sample, CLW pseudo-likelihood (2) coincides with the
population-based likelihood (2), thus making CLW the optimal method. ILR and
PILR are idetical in this case since pseudo-likelihood (16) of the PILR method
coincides with the ILR likelihood (11). As expected, CLW based estimates of
${\beta}_{c1}$ are more efficient than ILR/PILR based estimates.
Interestingly, variances of the estimates of $\mu$ are still smaller for
ILR/PILR methods. These observations are also confirmed by the numerical study
presented in Figures 3 and 4 of Section 3.2.

Simulation results verify that variance expressions (3.1) and (25) work
accurately for all considered estimation methods.

As mentioned in Section 1, there is a perception that the existence of the
overlapping units in the stacked sample set would prohibit formulation of a
likelihood under the independence assumption for the ILR and PILR methods. The
closeness of the estimated $\overline{\widehat{SE}}$ and MC simulation based
$SE$ numbers verifies formulas (3.1) and (25) and, indirectly, confirms the
validity of the conditional independence assumption. At the same time, if
overlapping units could be identified, this information could be used to
modify the CRISP formula as in (10), potentially leading to more efficient
estimation.

Additional results for scenarios S1-S4 are reported in Tables 6 \- 9 of
Appendix C. Boxplots in Figures 7 \- 12 in Appendix C show the distribution of
relative biases of respective estimated parameters under scenarios S1-S6.

##  5 Concluding remarks

In this paper we reviewed $3$ recent approaches for estimation of a non-
probability survey participation probabilities. We showed that in many
practical settings the likelihood based ILR approach provides for more
efficient estimates of the model parameters and the IPW estimate of the finite
population mean.

The ILR and PILR methods are based on the CRISP formula relating inclusion
probabilities of unobserved sample indicators to the probability that can be
modeled over an observed stacked samples set. We proved conditional
independence of the indicator variables on the stacked set under the
assumption of conditional independence of the samples inclusion indicators.
This property allows the use of a likelihood for independent Bernoulli
variables formulated over the observed data. An advantage for the ILR approach
that results from this property is that the construction of a valid Bernoulli
likelihood on the observed pooled sample allows for application of standard
statistical estimation methods.

We derived asymptotic variances of the estimators and compared theoretical
properties of the CLW, ILR, and PILR methods. The formulas revealed that the
efficiency of the estimates is impacted by the design-based variance over the
repeated probability samples. Small reference samples induce increases in
variances of the target parameters in all methods; however, the degree of the
effect varies among the methods with the ILR approach being the least
susceptible to the effect of a small reference sample size.

The variances are also affected by the degree of the overlap in variable-
defined domains between the two samples. Low overlap may cause instability in
the estimates even with relatively large reference samples. Any samples
combining method depends on the existence of a set of good covariates shared
by both samples. Machine learning methods can be useful for the effective
modeling. However, ML methods such as classification and regression trees may
result in the low-overlap terminal nodes having a small number of units from
either sample. The ILR based likelihood could be used as an objective function
in these methods, making them more effective.

The methodology is general and can be extended in various ways. At present,
for the ILR approach, reference sample inclusion probabilities were assumed to
be known for units in both samples. Savitsky et al. (2022) co-modeled them
along with the non-probability sample participation probabilities and used
splines in place of the linear regression model. This approach could be
extended to include the outcome variables in the modeling as well.

## References

  * Beaumont and Rao (2021) [author] Beaumont, Jean-FrancoisJ.-F. and Rao, J. N. K.J. N. K. (2021). Pitfalls of making inferences from non-probability samples: Can data integration through probability samples provide remedies. The Survey Statistician 83 11 – 22.
  * Beaumont et al. (2023) [author] Beaumont, Jean-FrançoisJ.-F., Bosa, KevenK., Brennan, AndrewA., Charlebois, JoanneJ. and Chu, KennethK. (2023). Handling non-probability samples through inverse probability weighting with an application to Statistics Canada’s crowdsourcing data. Survey Methodology (accepted in 2023 and expected to appear in 2024).
  * Beresovsky (2019) [author] Beresovsky, VladislavV. (2019). On application of a response propensity model to estimation from web samples. In ResearchGate.
  * Bethlehem (2010) [author] Bethlehem, J.J. (2010). Selection Bias in Web Surveys. International Statistical Review 78 161 –- 188.
  * Chen, Li and Wu (2020) [author] Chen, YilinY., Li, PengfeiP. and Wu, ChangbaoC. (2020). Doubly Robust Inference With Nonprobability Survey Samples. Journal of the American Statistical Association 115 2011-2021. 10.1080/01621459.2019.1677241
  * Elliott (2009) [author] Elliott, Michael R.M. R. (2009). Combining Data from Probability and Non-Probability Samples Using Pseudo-Weights. Survey Practice 2 813–845.
  * Elliott and Valliant (2017) [author] Elliott, Michael R.M. R. and Valliant, RichardR. (2017). Inference for Nonprobability Samples. Statistical Science 32 249 – 264. 10.1214/16-STS598
  * Firth (1993) [author] Firth, DavidD. (1993). Bias reduction of maximum likelihood estimates. Biometrika 80 27–38.
  * Heinze and Schemper (2002) [author] Heinze, G.G. and Schemper, M.M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine 21 2409–2419.
  * Kim and Tam (2021) [author] Kim, J. K.J. K. and Tam, S. M.S. M. (2021). Data Integration by Combining Big Data and Survey Sample Data for Finite Population Inference. International Statistical Review 89 382–-401.
  * Liu, Scholtus and De Waal (2023) [author] Liu, An-ChiaoA.-C., Scholtus, SanderS. and De Waal, TonT. (2023). Correcting Selection bias in Big Data by pseudo-weighting. Journal of Survey Statistics and Methodology 11 1181–-1203.
  * Meng (2018) [author] Meng, X. L.X. L. (2018). Statistical Paradises and Paradoxes in Big Data (I): Law of Large Populations, Big Data Paradox, and the 2016 US Presidential Election. Annals of Applied Statistics 12 685 – 726.
  * Rao (2021) [author] Rao, J. N. K.J. N. K. (2021). On Making Valid Inferences by Integrating Data from Surveys and Other Sources. Sankhya 83 242 – 272.
  * Savitsky et al. (2022) [author] Savitsky, Terrance D.T. D., Williams, Matthew R.M. R., Gershunskaya, JulieJ., Beresovsky, VladislavV. and Johnson, Nels G.N. G. (2022). Methods for Combining Probability and Nonprobability Samples Under Unknown Overlaps.  https://doi.org/10.48550/arXiv.2208.14541.
  * Tille and Matei (2021) [author] Tille, YvesY. and Matei, AlinaA. (2021). sampling: Survey Sampling.  https://CRAN.R-project.org/package=sampling.
  * Valliant (2020) [author] Valliant, RichardR. (2020). Comparing Alternatives for Estimation from Nonprobability Samples. Journal of Survey Statistics and Methodology 8 231 – 263.
  * VanderWeele and Shpitser (2011) [author] VanderWeele, T. J.T. J. and Shpitser, I.I. (2011). A new criterion for confounder selection. Biometrics 67 1406 —- 1413.
  * Wang, Valliant and Li (2021) [author] Wang, L.L., Valliant, R.R. and Li, Y.Y. (2021). Adjusted logistic propensity weighting methods for population inference using nonprobability volunteer-based epidemiologic cohorts. Stat Med. 40 5237–5250. 10.1002/sim.9122
  * Wu (2022) [author] Wu, ChangbaoC. (2022). A new criterion for confounder selection. Survey Methodology 48 283 —- 311.

##  Appendix A Proof of Theorem 2.2

By LTP, we have:

$\displaystyle\begin{split}&P\left\\{i,j\in S_{c}|i,j\in 2U\right\\}\\\
&=P\left\\{i,j\in S_{c}|i,j\in U_{c},i,j\in U\right\\}P\left\\{i,j\in
U_{c}|i,j\in U\right\\}\\\ &\times P\\{i,j\in U\mid i,j\in
2U\\}=\frac{1}{4}\pi_{c,ij}p_{c,ij}=\frac{1}{4}\pi_{ci}\pi_{cj}p_{ci}p_{cj}.\end{split}$
(31)

Similarly,

$\displaystyle\begin{split}&P\left\\{i,j\in S_{r}|i,j\in 2U\right\\}\\\
&=P\left\\{i,j\in S_{r}|i,j\in U_{r},i,j\in U\right\\}P\left\\{i,j\in
U_{r}|i,j\in U\right\\}\\\ &\times P\\{i,j\in U\mid i,j\in
2U\\}=\frac{1}{4}\pi_{r,ij}p_{r,ij}=\frac{1}{4}\pi_{ri}\pi_{rj}p_{ri}p_{rj};\end{split}$

$\displaystyle\begin{split}&P\left\\{i\in S_{c},j\in S_{r}|i,j\in
2U\right\\}\\\ &=P\left\\{i\in S_{c},j\in S_{r}|i\in U_{c},j\in U_{r},i,j\in
U\right\\}\\\ &\times P\left\\{i\in U_{c},j\in U_{r}|i,j\in
U\right\\}P\\{i,j\in U\mid i,j\in 2U\\}\\\
&=\frac{1}{4}\pi_{ci}\pi_{rj}p_{ci}p_{rj};\end{split}$

$\displaystyle\begin{split}&P\left\\{i\in S_{r},j\in S_{c}|i,j\in
2U\right\\}\\\ &=P\left\\{i\in S_{r},j\in S_{c}|i\in U_{r},j\in U_{c},i,j\in
U\right\\}\\\ &\times P\left\\{i\in U_{r},j\in U_{c}|i,j\in
U\right\\}P\\{i,j\in U\mid i,j\in 2U\\}\\\
&=\frac{1}{4}\pi_{ri}\pi_{cj}p_{ri}p_{cj}.\end{split}$

The total probability of being included into the pooled sample is the sum of
probabilities for the non-overlapping sets:

$\displaystyle\begin{split}&P\\{i,j\in S|i,j\in 2U\\}=\\\ &P\\{i,j\in
S_{c}|i,j\in 2U\\}+P\\{i,j\in S_{r}|i,j\in 2U\\}\\\ &+P\\{i\in S_{c},j\in
S_{r}|i,j\in 2U\\}\\\ &+P\\{i\in S_{r},j\in S_{c}|i,j\in 2U\\}\\\
&=\frac{1}{4}(\pi_{ci}\pi_{cj}p_{ci}p_{cj}+\pi_{ri}\pi_{rj}p_{ri}p_{rj}+\pi_{ci}\pi_{rj}p_{ci}p_{rj}\\\
&+\pi_{ri}\pi_{cj}p_{ri}p_{cj})\\\
&=\frac{1}{4}(\pi_{ci}p_{ci}+\pi_{ri}p_{ri})(\pi_{cj}p_{cj}+\pi_{rj}p_{rj})\end{split}$
(32)

Finally, by the definition of conditional probability,

$\displaystyle P\left\\{i,j\in S_{c}|i,j\in S,i,j\in 2U\right\\}$
$\displaystyle=\frac{P\left\\{i,j\in S_{c}|i,j\in 2U\right\\}}{P\left\\{i,j\in
S|i,j\in 2U\right\\}}.$ (33)

$\pi_{z,ij}=\pi_{zi}\pi_{zj}$ follows from putting (31) and (32) in formula
(33).

##  Appendix B Proof of Theorem 3.1

Consider vector of parameters
$\boldsymbol{\eta^{T}}=(\mu,\boldsymbol{\beta^{T}})$. The full set of
estimating equations is

$\boldsymbol{\Phi}(\boldsymbol{\eta})=\begin{pmatrix}U(\mu)\\\
S(\boldsymbol{\beta})\end{pmatrix}=\boldsymbol{0},$ (34)

where $U(\mu)=\mathop{\sum}_{i\in U}I_{ci}\pi_{ci}^{-1}(y_{i}-\mu)$ is a score
function for $\mu$ and $S(\boldsymbol{\beta})$ are score function for
respective models for $\beta.$ Let
$\boldsymbol{\eta}_{0}^{T}=(\mu_{0},\boldsymbol{\beta}_{0}^{T})$ be a vector
of true parameters satisfying
$E\boldsymbol{\Phi}(\boldsymbol{\eta}_{0})=\boldsymbol{0}$ and
$\hat{\boldsymbol{\eta}}$ a solution of (34). From first-order Taylor
expansion of $\boldsymbol{\Phi}(\hat{\boldsymbol{\eta}})$ around
$\boldsymbol{\eta}_{0}$, we have

$\displaystyle\hat{\boldsymbol{\eta}}-\boldsymbol{\eta}_{0}\doteq\left[E\\{\boldsymbol{\phi}(\boldsymbol{\eta}_{0})\\}\right]^{-1}\boldsymbol{\Phi}(\boldsymbol{\eta}_{0}),$
(35)

where

$E\\{\boldsymbol{\phi}(\boldsymbol{\eta})\\}=E\frac{\partial\boldsymbol{\Phi}(\boldsymbol{\eta})}{\partial\boldsymbol{\eta}}=\begin{pmatrix}U_{\mu}&U_{\boldsymbol{\beta}}\\\
\boldsymbol{0}&S_{\boldsymbol{\beta}}\end{pmatrix},$ (36) $\displaystyle
U_{\mu}$ $\displaystyle=E\frac{\partial U}{\partial\mu}=-\mathop{\sum}_{i\in
U}EI_{ci}\pi_{ci}^{-1}=-N,$ $\displaystyle U_{\boldsymbol{\beta}}$
$\displaystyle=E\frac{\partial
U}{\partial\boldsymbol{\beta^{T}}}=\mathop{\sum}_{i\in
U}(y_{i}-\mu)EI_{ci}\frac{\partial}{\partial\boldsymbol{\beta^{T}}}\pi_{ci}^{-1}$
$\displaystyle=-\mathop{\sum}_{i\in
U}(y_{i}-\mu)(1-\pi_{ci})\mathbf{x}_{i}^{T}.$

From (35), the approximate variance is

$\displaystyle
Var(\hat{\boldsymbol{\eta}})\doteq\left[E\\{\boldsymbol{\phi}(\boldsymbol{\eta}_{0})\\}\right]^{-1}Var\\{\boldsymbol{\Phi}(\boldsymbol{\eta}_{0})\\}\left[E\\{\boldsymbol{\phi}(\boldsymbol{\eta}_{0})\\}\right]^{-1},$
(37)

where the inverse of matrix (36) comes out to

$\left[E\\{\boldsymbol{\phi}(\boldsymbol{\eta})\\}\right]^{-1}=N^{-1}\begin{pmatrix}-1&\boldsymbol{b}^{T}\\\
\boldsymbol{0}&NS_{\boldsymbol{\beta}}^{-1}\end{pmatrix},$ (38)

where $\boldsymbol{b}=S_{\boldsymbol{\beta}}^{-1}U_{\boldsymbol{\beta}}^{T}$
and

$\displaystyle
Var[\boldsymbol{\Phi}(\boldsymbol{\eta}_{0})]=\begin{pmatrix}Var[U(\mu)]&\boldsymbol{C}^{T}\\\
\boldsymbol{C}&Var[S(\boldsymbol{\beta})]\end{pmatrix},$ (39)

where $\boldsymbol{C}=Cov[U(\mu),S(\boldsymbol{\beta})].$

Finally, putting (38) and (39) into (37), we find that

$\displaystyle Var(\hat{\boldsymbol{\eta}})\doteq$ (40)
$\displaystyle\begin{pmatrix}N^{-2}\\{Var[U(\mu)]-2\boldsymbol{b}^{T}\boldsymbol{C}+\boldsymbol{b}^{T}Var[S(\boldsymbol{\beta})]\boldsymbol{b}\\}&\;0\\\
-N^{-1}S_{\boldsymbol{\beta}}^{-1}[\boldsymbol{C}-Var[S(\boldsymbol{\beta})]\boldsymbol{b}]&\;0\end{pmatrix}$
$\displaystyle+\begin{pmatrix}0\;&-N^{-1}[\boldsymbol{C^{T}}-\boldsymbol{b}^{T}Var[S(\boldsymbol{\beta})]S_{\boldsymbol{\beta}}^{-1}\\\
0\;&S_{\boldsymbol{\beta}}^{-1}Var[S(\boldsymbol{\beta})]S_{\boldsymbol{\beta}}^{-1}\end{pmatrix}.$

where the design variance of $U(\mu)$ under the Poisson sampling is

$\displaystyle Var\\{U(\mu)\\}$ $\displaystyle=\mathop{\sum}_{i\in
U}Var(I_{ci})\left(\frac{y_{i}-\mu}{\pi_{ci}}\right)^{2}$
$\displaystyle=\mathop{\sum}_{i\in
U}(1-\pi_{ci})\pi_{ci}^{-1}(y_{i}-\mu)^{2}.$

If we present $S(\boldsymbol{\beta})$ as a sum of two independent parts
$S(\boldsymbol{\beta})=S_{c}(\boldsymbol{\beta})+S_{r}(\boldsymbol{\beta}),$
the variance can be written as

$\displaystyle Var[S(\boldsymbol{\beta})]$
$\displaystyle=Var[S_{c}(\boldsymbol{\beta})]+Var[S_{r}(\boldsymbol{\beta})].$

Let us now spell out (40) for the ILR and PILR methods.

The ILR case:

The score function is

$\displaystyle S(\boldsymbol{\beta})$ $\displaystyle=\mathop{\sum}_{i\in
S_{c}}(1-\pi_{zi})(1-\pi_{ci})\mathbf{x}_{i}-\mathop{\sum}_{i\in
S_{r}}\pi_{zi}(1-\pi_{ci})\mathbf{x}_{i}$ $\displaystyle=\mathop{\sum}_{i\in
U}[I_{ci}(1-\pi_{zi})-I_{ri}\pi_{zi}](1-\pi_{ci})\mathbf{x}_{i},$

that can be decomposed into:

$\displaystyle S_{c}(\boldsymbol{\beta})$ $\displaystyle=\mathop{\sum}_{i\in
U}I_{ci}(1-\pi_{zi})(1-\pi_{ci})\mathbf{x}_{i},$ $\displaystyle
S_{r}(\boldsymbol{\beta})$ $\displaystyle=-\mathop{\sum}_{i\in
U}I_{ri}\pi_{zi}(1-\pi_{ci})\mathbf{x}_{i}.$

Next,

$\displaystyle\boldsymbol{H}_{ILR}$
$\displaystyle=-S_{\boldsymbol{\beta}}=-E\frac{\partial
S(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}}$ $\displaystyle=\sum_{i\in
U}(\pi_{ci}+\pi_{ri})\pi_{zi}(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}.$

The variance of $S_{c}(\boldsymbol{\beta})$ is

$\displaystyle\boldsymbol{A}_{ILR}=Var[S_{c}(\boldsymbol{\beta})]$
$\displaystyle=\mathop{\sum}_{i\in
U}\pi_{ci}(1-\pi_{ci})(1-\pi_{zi})^{2}(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}$
$\displaystyle=\boldsymbol{H}_{ILR}-\sum_{i\in
U}(\pi_{ri}+1)\pi_{ci}\pi_{zi}(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}$

The variance of $S_{r}(\boldsymbol{\beta})$ is
$\boldsymbol{D}_{ILR}=Var_{d}[\mathop{\sum}_{i\in
U}I_{ri}\pi_{zi}(1-\pi_{ci})\mathbf{x}_{i}]$ is the variance-covariance matrix
under the probability sample design.

For the covariance, using the independence of indicators $I_{ci}$ and
$I_{ri}$, we have

$\displaystyle\boldsymbol{C}_{ILR}=Cov[U(\mu),S(\boldsymbol{\beta})]$
$\displaystyle=\mathop{\sum}_{i\in
U}Cov(I_{ci},I_{ci})\pi_{ci}^{-1}(y_{i}-\mu)(1-\pi_{zi})(1-\pi_{ci})\mathbf{x}_{i}$
$\displaystyle=\mathop{\sum}_{i\in
U}(y_{i}-\mu)(1-\pi_{zi})(1-\pi_{ci})^{2}\mathbf{x}_{i}.$

Variance formulas for the parameters in the ILR method follow by putting
together the above formulas.

It is easy to see, using regularity condition C2, that

$\displaystyle Var({U\mu})$
$\displaystyle=N^{2}O\left(\frac{1}{n_{c}}\right),$
$\displaystyle\boldsymbol{b}_{ILR}^{T}\boldsymbol{C}_{ILR}$
$\displaystyle=N^{2}O\left(\frac{n_{c}+n_{r}}{n_{c}n_{r}}\frac{n_{r}}{n_{c}+n_{r}}\right),$
$\displaystyle\boldsymbol{b}_{ILR}^{T}(\boldsymbol{A}_{ILR}+\boldsymbol{D}_{ILR})\boldsymbol{b}_{ILR}$
$\displaystyle=N^{2}O\left(\frac{(n_{c}+n_{r})^{2}}{n_{c}^{2}n_{r}^{2}}\left\\{\frac{n_{c}n_{r}^{2}}{(n_{c}+n_{r})^{2}}+\frac{n_{c}^{2}n_{r}}{(n_{c}+n_{r})^{2}}\right\\}\right),$
$\displaystyle\boldsymbol{H}_{ILR}^{-1}\left(\boldsymbol{A}_{ILR}+\boldsymbol{D}_{ILR}\right)\boldsymbol{H}_{ILR}$
$\displaystyle=O\left(\frac{(n_{c}+n_{r})^{2}}{n_{c}^{2}n_{r}^{2}}\left\\{\frac{n_{c}n_{r}^{2}}{(n_{c}+n_{r})^{2}}+\frac{n_{c}^{2}n_{r}}{(n_{c}+n_{r})^{2}}\right\\}\right).$

After putting it into the variance formulas, we have

$\displaystyle Var(\hat{\mu}_{ILR})$
$\displaystyle=O\Big{(}\frac{1}{n_{c}}+\frac{1}{n_{r}}\Big{)},$ $\displaystyle
Var(\boldsymbol{\hat{\beta}}_{ILR})$
$\displaystyle=O\Big{(}\frac{1}{n_{c}}+\frac{1}{n_{r}}\Big{)}.$

Hence, $\hat{\mu}_{ILR}$ and $\boldsymbol{\hat{\beta}}_{ILR}$ are consistent
estimates of, respectively, $\mu$ and $\boldsymbol{\beta}$ and
$\hat{\mu}_{ILR}-\mu=O_{p}(\text{min}(n_{c},n_{r})^{-1/2})$ and
$\boldsymbol{\hat{\beta}}_{ILR}-\boldsymbol{{\beta}}=O_{p}(\text{min}(n_{c},n_{r})^{-1/2})$.

The PILR case:

The score function is

$\displaystyle\begin{split}S(\boldsymbol{\beta})&=\mathop{\sum}_{i\in
S_{c}}(1-\pi_{{\delta}i})(1-\pi_{ci})\mathbf{x}_{i}\\\ &-\mathop{\sum}_{i\in
S_{r}}w_{ri}\pi_{{\delta}i}(1-\pi_{ci})\mathbf{x}_{i}\\\ \end{split}$
$\displaystyle=\mathop{\sum}_{i\in
U}\left\\{I_{ci}(1-\pi_{{\delta}i})-I_{ri}w_{ri}\pi_{{\delta}i}\right\\}(1-\pi_{ci})\mathbf{x}_{i},$

that can be decomposed into:

$\displaystyle S_{c}({\boldsymbol{\beta}})$ $\displaystyle=\mathop{\sum}_{i\in
U}(I_{ci}(1-\pi_{{\delta}i})-\pi_{{\delta}i})(1-\pi_{ci})\mathbf{x}_{i},$
$\displaystyle S_{r}({\boldsymbol{\beta}})$ $\displaystyle=\mathop{\sum}_{i\in
U}(1-I_{ri}w_{ri})\pi_{{\delta}i}(1-\pi_{ci})\mathbf{x}_{i}.$

Next,

$\displaystyle\boldsymbol{H}_{PILR}$
$\displaystyle=-S_{\boldsymbol{\beta}}=-E\frac{\partial
S(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}}$ $\displaystyle=\sum_{i\in
U}(\pi_{ci}+1)\pi_{{\delta}i}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}.$

The variance of $S_{c}(\boldsymbol{\beta})$ is

$\displaystyle\boldsymbol{A}_{PILR}$
$\displaystyle=Var[S_{c}(\boldsymbol{\beta})]$
$\displaystyle=\mathop{\sum}_{i\in
U}\pi_{ci}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}$
$\displaystyle-2\mathop{\sum}_{i\in
U}\pi_{ci}\pi_{{\delta}i}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}$
$\displaystyle=\boldsymbol{H}_{PILR}-2\mathop{\sum}_{i\in
U}\pi_{ci}\pi_{{\delta}i}(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}\mathbf{x}_{i}^{T}.$

The variance of $S_{r}(\boldsymbol{\beta})$ is
$\boldsymbol{D}_{PILR}=Var_{d}[\mathop{\sum}_{i\in
U}I_{ri}w_{ri}\pi_{{\delta}i}(1-\pi_{ci})\mathbf{x}_{i}]$, the variance-
covariance matrix under the probability sample design.

For the covariance, using the independence of indicators $I_{ci}$ and
$I_{ri}$, we have

$\displaystyle\boldsymbol{C}_{PILR}=Cov[U(\mu),S(\boldsymbol{\beta})]$
$\displaystyle=\mathop{\sum}_{i\in
U}Cov(I_{ci},I_{ci})\pi_{ci}^{-1}(y_{i}-\mu)(1-\pi_{{\delta}i})(1-\pi_{ci})\mathbf{x}_{i}$
$\displaystyle=\mathop{\sum}_{i\in
U}(y_{i}-\mu)(1-\pi_{{\delta}i})(1-\pi_{ci})^{2}\mathbf{x}_{i}.$

Variance formulas for the parameters in the PILR method follow by putting
together the above formulas.

We can see, using regularity condition C2, that

$\displaystyle Var({U\mu})$
$\displaystyle=N^{2}O\left(\frac{1}{n_{c}}\right),$
$\displaystyle\boldsymbol{b}_{PILR}^{T}\boldsymbol{C}_{PILR}$
$\displaystyle=N^{2}O\left({\frac{1}{n_{c}}}\right),$
$\displaystyle\boldsymbol{b}_{PILR}^{T}(\boldsymbol{A}_{PILR}+\boldsymbol{D}_{PILR})\boldsymbol{b}_{PILR}$
$\displaystyle=N^{2}O\Big{[}({\frac{1}{n_{c}^{2}}})(n_{c}+\frac{n_{c}^{2}}{n_{r}})\Big{]},$
$\displaystyle\boldsymbol{H}_{PILR}^{-1}(\boldsymbol{A}_{PILR}+\boldsymbol{D}_{PILR})\boldsymbol{H}_{PILR}$
$\displaystyle=O\Big{[}\frac{1}{n_{c}^{2}}(n_{c}+\frac{n_{c}^{2}}{n_{r}})\Big{]}.$

After putting it into the variance formulas, we have

$\displaystyle Var(\hat{\mu}_{PILR})$
$\displaystyle=O\Big{(}\frac{1}{n_{c}}+\frac{1}{n_{r}}\Big{)},$ $\displaystyle
Var(\boldsymbol{\hat{\beta}}_{PILR})$
$\displaystyle=O\Big{(}\frac{1}{n_{c}}+\frac{1}{n_{r}}\Big{)}.$

Hence, $\hat{\mu}_{ILR}$ and $\boldsymbol{\hat{\beta}}_{ILR}$ are consistent
estimates of, respectively, $\mu$ and $\boldsymbol{\beta}$ and
$\hat{\mu}_{ILR}-\mu=O_{p}(\text{min}(n_{c},n_{r})^{-1/2})$ and
$\boldsymbol{\hat{\beta}}_{ILR}-\boldsymbol{{\beta}}=O_{p}(\text{min}(n_{c},n_{r})^{-1/2})$.

##  Appendix C Simulation Results: Tables and plots

The Tables and Boxplots of this Section show the bulk of simulations results
discussed in the main part of the paper. Additional results for scenarios
S1-S4 are reported in Tables 6 \- 9. Boxplots in Figures 7 \- 12 show the
distribution of relative biases of respective estimated parameters under
scenarios S1-S6.

Under the high overlap and relatively large size of the reference sample
(upper-left quadrant in each Table), estimates perform similarly well for all
methods: estimates of both $\beta_{c1}$ and $\mu$ are nearly unbiased and
variance estimates perform well providing for the nominal coverage by the
normal confidence intervals. The only exception to this satisfactory
performance are the estimates of ${\mu}$ in Table 4 that are positively biased
by about $7\%$ for all the methods. Corresponding estimates for ${\beta}_{c1}$
(Table 3) are essentially unbiased.

Results in the lower-left quadrants of the Tables and boxplots correspond to
reference samples of smaller size $n_{r}=100$. Table 6 presents results for
small sampling fractions $f_{c}=f_{r}=0.01$. It shows that ILR estimates of
${\beta}_{c1}$ are almost unbiased and about $30\%$ more efficient than the
CLW based estimates. Also, the CLW based estimates are biased by about $7\%$.
PILR estimates are of intermediate quality. Estimates corresponding to a
larger sampling fractions $f_{c}=f_{r}=0.10$ in the lower-left quadrant of
Table 8 show significantly smaller variation between different estimation
methods. Numbers in the upper-right quadrants of these Tables correspond to
relatively large reference samples and low overlap in covariate-defined
domains. For small sampling fractions of $0.01$, ILR estimates of
${\beta}_{c1}$ have a significantly smaller bias and variance compared to the
estimates from the CLW and PILR methods. This bias in the estimate of
$\beta_{c1}$ also affects the variance estimates for the PILR and CLW methods,
while variance estimates for the ILR perform satisfactory. For the larger
sampling fraction of $0.10$, estimates based on the PILR and ILR methods are
comparable and better than the CLW estimates. A similar pattern is observed in
the estimates of population mean ${\mu}$.

Results in the lower-right quadrants of the Tables correspond to small
reference samples sizes and low overlap in covariates domains between the non-
probability and reference samples. For smaller sampling fractions of the non-
probability sample, estimates by the ILR method have significantly smaller
bias and variance comparing to estimates by other methods. For a larger
sampling fractions of $0.10$, PILR and ILR estimates in Tables 8 and 9 are
comparable and significantly better than CLW estimates. Variances of ILR
estimates are accurately estimated and produce close to nominal coverage of
population parameters by the confidence intervals; PILR variances are
reasonably accurate as well. However, for the CLW method, large biases of
point estimates result in very inaccurate variance estimates.

Table 6:  Estimate of propensity parameter ${\beta}_{c1}$ for sample fractions $f_{c}=f_{r}=0.01$ | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S1: $n_{c}=600,n_{r}=600$
ILR | 1.01 | 0.07 | 0.07 | 0.95 | 1.01 | 0.08 | 0.08 | 0.96
PILR | 1.01 | 0.08 | 0.08 | 0.95 | 1.05 | 0.18 | 0.14 | 0.88
CLW | 1.01 | 0.09 | 0.09 | 0.94 | 1.07 | 0.21 | 0.16 | 0.86
| Scenario S2: $n_{c}=100,n_{r}=100$
ILR | 1.01 | 0.17 | 0.17 | 0.95 | 1.04 | 0.22 | 0.21 | 0.95
PILR | 1.04 | 0.22 | 0.20 | 0.94 | 1.22 | 0.40 | 0.32 | 0.87
CLW | 1.07 | 0.25 | 0.22 | 0.93 | 1.44 | 1.44 | 5.88 | 0.86
Table 7:  Estimate of population mean ${\mu}$ for sample fractions $f_{c}=f_{r}=0.01$ | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S1: $n_{c}=600,n_{r}=600$
ILR | 1.00 | 0.13 | 0.13 | 0.93 | 1.00 | 0.13 | 0.13 | 0.93
PILR | 1.00 | 0.14 | 0.15 | 0.93 | 0.95 | 0.22 | 0.21 | 0.94
CLW | 0.99 | 0.15 | 0.15 | 0.94 | 0.93 | 0.25 | 0.24 | 0.95
| Scenario S2: $n_{c}=100,n_{r}=100$
ILR | 1.05 | 0.32 | 0.30 | 0.89 | 1.03 | 0.32 | 0.31 | 0.91
PILR | 1.01 | 0.38 | 0.42 | 0.91 | 0.87 | 0.52 | 1.03 | 0.93
CLW | 0.99 | 0.41 | 0.53 | 0.91 | 0.77 | 0.64 | Inf | 0.95
Table 8:  Estimate of propensity parameter ${\beta}_{c1}$ for sample fractions $f_{c}=f_{r}=0.1$ | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S3: $n_{c}=600,n_{r}=600$
ILR | 1.00 | 0.08 | 0.08 | 0.94 | 1.00 | 0.10 | 0.09 | 0.94
PILR | 1.00 | 0.08 | 0.08 | 0.95 | 1.01 | 0.13 | 0.12 | 0.93
CLW | 1.01 | 0.09 | 0.10 | 0.95 | 1.05 | 0.20 | 0.17 | 0.91
| Scenario S4: $n_{c}=100,n_{r}=100$
ILR | 1.04 | 0.22 | 0.20 | 0.95 | 1.04 | 0.25 | 0.24 | 0.94
PILR | 1.04 | 0.22 | 0.21 | 0.95 | 1.05 | 0.30 | 0.29 | 0.94
CLW | 1.07 | 0.28 | 0.27 | 0.96 | 1.32 | 1.18 | 7.83 | 0.93
Table 9:  Estimate of population mean ${\mu}$ for sample fractions $f_{c}=f_{r}=0.1$ | High overlap | Low overlap
---|---|---
| $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$ | $Mean$ | $SE$ | $\overline{\widehat{SE}}$ | $95\%CI$
| Scenario S3: $n_{c}=600,n_{r}=600$
ILR | 1.02 | 0.12 | 0.11 | 0.91 | 0.97 | 0.12 | 0.11 | 0.93
PILR | 1.01 | 0.12 | 0.12 | 0.91 | 0.96 | 0.15 | 0.14 | 0.95
CLW | 1.01 | 0.14 | 0.14 | 0.91 | 0.92 | 0.21 | 0.21 | 0.95
| Scenario S4: $n_{c}=100,n_{r}=100$
ILR | 1.03 | 0.27 | 0.26 | 0.90 | 1.00 | 0.29 | 0.25 | 0.89
PILR | 1.02 | 0.29 | 0.31 | 0.90 | 0.97 | 0.37 | 0.37 | 0.91
CLW | 1.00 | 0.33 | 0.41 | 0.91 | 0.81 | 0.61 | Inf | 0.93
Figure 7: Relative bias of the estimated propensity parameters ${\beta}_{c1}$
for sample fraction 0.01 over the MC simulations, scenarios S1 and S2 Figure
8: Relative bias of the estimated population mean ${\mu}$ for sample fraction
0.01 over the MC simulations, scenarios S1 and S2 Figure 9: Relative bias of
the estimated propensity parameters ${\beta}_{c1}$ for sample fraction 0.1
over the MC simulations, scenarios S3 and S4 Figure 10: Relative bias of the
estimated population mean ${\mu}$ for sample fraction 0.1 over the MC
simulations, scenarios S3 and S4 Figure 11: Relative bias of the estimated
propensity parameters ${\beta}_{c1}$ for sampling fractions
$\left(f_{c},f_{r}\right)$ equal to $\left(0.01,0.1\right)$ (upper row) and to
$\left(0.1,0.01\right)$ (lower row) over the MC simulations, scenarios S5 and
S6 Figure 12: Relative bias of the estimated population mean ${\mu}$ for
sampling fractions $\left(f_{c},f_{r}\right)$ equal to $\left(0.01,0.1\right)$
(upper row) and to $\left(0.1,0.01\right)$ (lower row) over the MC
simulations, scenarios S5 and S6