# Data integration of non-probability and probability samples with predictive
mean matching

Piotr Chlebicki111Adam Mickiewicz University in Poznań, Poland; Statistical
Office in Poznań, ul. Wojska Polskiego 27/29 60-624 Poznań, Poland., Łukasz
Chrostowski222Adam Mickiewicz University in Poznań, Poland, Maciej
Beręsewicz333Poznań University of Economics and Business, Institute of
Informatics and Quantitative Economics, Department of Statistics, Al.
Niepodległości 10, 61-875 Poznań, Poland, E-mail:
<EMAIL_ADDRESS>Statistical Office in Poznań, ul. Wojska
Polskiego 27/29 60-624 Poznań, Poland.

###### Abstract

In this paper we study predictive mean matching mass imputation estimators to
integrate data from probability and non-probability samples. We consider two
approaches: matching predicted to observed ($\hat{y}-y$ matching) or predicted
to predicted ($\hat{y}-\hat{y}$ matching) values. We prove the consistency of
two semi-parametric mass imputation estimators based on these approaches and
derive their variance and estimators of variance. Our approach can be employed
with non-parametric regression techniques, such as kernel regression, and the
analytical expression for variance can also be applied in nearest neighbour
matching for non-probability samples. We conduct extensive simulation studies
in order to compare the properties of this estimator with existing approaches,
discuss the selection of $k$-nearest neighbours, and study the effects of
model mis-specification. The paper finishes with empirical study in
integration of job vacancy survey and vacancies submitted to public employment
offices (admin and online data). Open source software is available for the
proposed approaches.

Keywords: mass imputation, non-probability surveys, variance estimation,
bootstrap, job vacancy survey

###  1 Introduction

With the availability of large sets of administrative data, voluntary internet
panels, social media and big data, inference with non-probability samples is
being heavily studied in the statistical literature [21, 29, 22, 28]. Because
of their non-statistical character and unknown sampling mechanism, these
sources cannot be used directly for estimating population characteristics.

Several inference approaches have been proposed in the literature with respect
to data from non-probability samples, which either involve data integration
with population level data or probability samples from the same population
[37, see for recent review]. Main techniques for inference based on such data
can broadly be classified into inverse probability weighting (IPW) estimators,
prediction estimators (PE) and doubly robust (DR) estimators, which include
both IPW and PE. In this paper we focus on mass imputation (MI) estimators,
which are a variant of prediction estimators used for integrating probability
and non-probability samples [29, cf.]. The general idea is to impute values of
the target variable $Y$ for all units in a probability sample based on values
observed in non-probability samples.

In recent papers two distinct MI estimators have been proposed: (1) the
nearest neighbour (NN; cf. [40]) imputation estimator and (2) the parametric
(PAR) and non-parametric (NPAR) imputation estimator [31, 25, cf.]. The latter
two consist in using a model for $\mathbb{E}[Y|\boldsymbol{X}]$ with
parameters estimated from a non-probability sample and then using predictions
from this model as imputed values for a probability sample. In contrast, NN
imputation involves finding $k$ nearest neighbours for each unit in a
probability sample from a set of units from a non-probability sample (imputed
values are observed values of $Y$ from the non-probability sample).

In this paper we focus on the predictive mean matching (PMM) MI estimator in
the context of data integration, which can be seen as a technique combining NN
and parametric (or non-parametric) imputation estimators. PMM combined with
multiple imputation was introduced by [32, 34] for imputing missing data in
surveys. [35] and [30] advocated PMM imputation because of its robustness to
model mis-specification. The asymptotic properties of the PMM estimator were
recently discussed in the context of survey non-response by [38]. Initial work
on multiply robust estimation was done by [24]. However, the PMM estimator has
not yet been discussed or studied in the context of data integration involving
non-probability and probability samples.

Our contribution can be summarized as follows:

  1. 1.

We study properties of two variants of semi-parametric PMM estimators for the
data integration problem in two settings: matching either by predicted-
observed ($\hat{y}-y$) or predicted-predicted ($\hat{y}-\hat{y}$) values. The
latter option is the standard PMM estimator, which is used in most of the
existing studies, while the first one is less commonly discussed.

  2. 2.

We prove the consistency of these estimators under suitable assumptions where
population size is either known or estimated. The theorem can easily be
extended to non-parametric models, such as kernel or local linear regression.
We also discuss consistency when the model is mis-specified.

  3. 3.

We derive analytical and bootstrap variance estimators and study their
performance for finite populations. Furthermore, the proposed closed form
expression for analytical variance can also be used for the NN estimator. We
also prove a result analogous to [40] concerning asymptotic and large non-
probability sample behaviour of the variance estimator.

The structure of the paper is as follows. In section 2 we provide basic
notation, describe two PMM algorithms and discuss the assumptions underlying
the PMM estimators. Section 3 provides proof of consistency of the PMM
estimators and the exact variance and its estimators. Section 4 describes one
of the simulation studies, whose results verify the proposed approach under
simple random sampling where two linear and non-linear models are considered.
Section 5 presents our empirical study. The article ends with conclusions and
the Appendix including details concerning the proofs as well as additional
simulation studies. All calculations in the article can be reproduced using R
[33] and the nonprobsvy package [27] developed by the authors. All codes are
available at https://github.com/ncn-foreigners/paper-nonprob-pmm.

###  2 Basic setup

####  2.1 Notation and motivation

Let $U=\\{1,...,N\\}$ denote the target population consisting of $N$ labelled
units. Each unit $i$ has an associated vector of auxiliary variables
$\boldsymbol{x}_{i}$ (a realisation of the random vector $\boldsymbol{X}$ in
the super-population) and the study variable $y_{i}$ (a realisation of the
random variable $Y$ in the super-population). Let
$\\{(y_{i},\boldsymbol{x}_{i}),i\in S_{A}\\}$ be a dataset of a non-
probability sample of size $n_{A}$ and let
$\\{\left(\boldsymbol{x}_{i},\pi_{i}\right),i\in S_{B}\\}$ be a dataset of a
probability sample of size $n_{B}$, where only information about variables
$\boldsymbol{X}$ and inclusion probabilities $\pi$ are available. Let $\delta$
be an indicator of inclusion into non-probability sample. Each unit in the
sample $S_{B}$ has been assigned a design-based weight given by
$d_{i}=1/\pi_{i}$. There is no overlap between the two samples. The setting is
summarised in Table 2.1.

Table 2.1: Two sample setting. Sample | ID |  | Sample weight
---
$d=\pi^{-1}$
Covariates $\boldsymbol{x}$ | Study variable $y$
Non-probability sample ($S_{A})$ | 1 | ? | $\checkmark$ | $\checkmark$
| $\vdots$ | ? | $\vdots$ | $\vdots$
| $n_{A}$ | ? | $\checkmark$ | $\checkmark$
Probability sample ($S_{B}$) | 1 | $\checkmark$ | $\checkmark$ | ?
| $\vdots$ | $\vdots$ | $\vdots$ | ?
| $n_{B}$ | $\checkmark$ | $\checkmark$ | ?

The goal is to estimate a finite population mean
$\displaystyle\mu_{y}=\frac{1}{N}\sum_{i=1}^{N}y_{i}$ of the target variable
$y$. As values of $y_{i}$ are not observed in the probability sample, it
cannot be used to estimate the target quantity. Instead, one could try
combining the non-probability and probability samples to estimate $\mu_{y}$.

In our paper we focus on prediction, which is an appealing solution as it
allows us to make available micro-level data with fully imputed target
variables and design specification. Moreover, we impute observed values of
$Y$, which makes it easier to explain them to users of the micro-data (since
observations are not generated from an unknown model).

####  2.2 Predictive mean matching algorithms

In this section we describe two general algorithms for semi-parametric PMM
estimators. Let $\hat{y}_{i}$ be a predicted value for unit $i$ from the
assumed regression model
$m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)$. Algorithm 1
(without step 2B) is a $\hat{y}-y$ matching imputation estimator, in which
units from the probability sample $S_{B}$ are matched based on $\hat{y}$ to
those from the non-probability sample $S_{A}$ based on $y$. This estimator is
given by (3) (throughout the paper this estimator is denoted as $\hat{y}-y$ or
the PMM A estimator). Algorithm 1 (without step 2A) represents the standard
$\hat{y}-\hat{y}$ matching (throughout the paper this estimator is denoted as
$\hat{y}-\hat{y}$ or as the PMM B estimator). The resulting estimator is given
by the same expression (3) but the matching is done differently.

1 Estimate regression model
$\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}]=m(\boldsymbol{x},\boldsymbol{\beta})$.

2A ($\hat{y}-y$ imputation; PMM A) Impute
$\hat{y}_{i}=m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)$ for
$i\in S_{B}$ and map each $i\in S_{B}$ to $\hat{\nu}(i)$ where:
$\displaystyle\hat{\nu}(i)=\operatorname*{arg\,min}_{j\in
S_{A}}\lVert\hat{y}_{i}-y_{j}\rVert$ or
$\displaystyle\hat{\nu}(i)=\operatorname*{arg\,min}_{j\in
S_{A}}d\left(\hat{y}_{i},y_{j}\right)$ if $d$ is not norm-induced.

If $k>1$ instead let:  $\hat{\nu}(i,z)=\operatorname*{arg\,min}_{\displaystyle
j\in
S_{A}\setminus\bigcup_{t=1}^{z-1}\\{\hat{\nu}(i,t)\\}}d\left(\hat{y}_{i},y_{j}\right).$
(1)

2B ($\hat{y}-\hat{y}$ imputation; PMM B) Impute
$\hat{y}_{i}=m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right),\hat{y}_{j}=m\left(\boldsymbol{x}_{j},\hat{\boldsymbol{\beta}}\right)$
for $i\in S_{B},j\in S_{A}$ and map each $i\in S_{B}$ to $\hat{\nu}(i)$ where:
$\displaystyle\hat{\nu}(i)=\operatorname*{arg\,min}_{j\in
S_{A}}\lVert\hat{y}_{i}-\hat{y}_{j}\rVert$ or
$\displaystyle\hat{\nu}(i)=\operatorname*{arg\,min}_{j\in
S_{A}}d\left(\hat{y}_{i},\hat{y}_{j}\right)$ if $d$ is not norm-induced. If
$k>1$ let instead let:
$\hat{\nu}(i,z)=\operatorname*{arg\,min}_{\displaystyle j\in
S_{A}\setminus\bigcup_{t=1}^{z-1}\\{\hat{\nu}(i,t)\\}}d\left(\hat{y}_{i},\hat{y}_{j}\right),$
(2) i.e. the $\hat{\nu}(i,z)$ is the $z$’th nearest neighbour from the sample.

3 Use imputed value in estimation:
$\displaystyle\hat{\mu}=\frac{1}{N}\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}y_{\hat{\nu}(i)}$, or if $k>1$ then instead:
$\hat{\mu}=\frac{1}{N}\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}.$ (3)

Algorithm 1 Two predictive mean matching imputation algorithms

######  Remark 1.

(Regarding $\hat{y}-y$ matching) If the study variable $Y$ is a binary
variable, then $\hat{y}-y$ matching estimator (3) is just a modification of
the GLM estimator, in which, instead of predicted values, we impute a
classification for each $i\in S_{B}$. The resulting estimate may still be
consistent according to Theorem 1 but may have comparatively poor finite
sample properties in practice.

######  Remark 2.

(Regarding nearest neighbour) If we let
$m(\boldsymbol{x}_{i},\boldsymbol{\beta})=\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}$
and use only one covariate, the resulting estimate (3) for $\hat{y}-\hat{y}$
matching is exactly the nearest neighbour estimate, since linear scaling
preserves partial linear order $\leq$ on $\mathbb{R}_{+}$.

When $m$ is a linear regression and $d$ is a Euclidean distance
$\hat{y}-\hat{y}$ matching is a "pseudo-nearest neighbour" matching with a
weighted (by estimated regression coefficients) distance function
$d^{\prime}:(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\mapsto\left|(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})^{T}\hat{\boldsymbol{\beta}}\right|$.
If the number of covariates is greater than one, the partial order generated
by the distance function is not deterministic in the general case, since
$\hat{\boldsymbol{\beta}}$ is a random variable, and therefore one cannot
simply assert that the $\hat{y}-\hat{y}$ matching estimator inherits
properties from the NN estimator.

####  2.3 Assumptions

We state the main assumptions used for studying properties of estimators
defined in (3) with a distance function being the Euclidean metric. See [31,
38] for a comparison with assumptions for other mass imputation estimators.

Throughout the paper we assume that in the super-population model every $Y$
variable has a finite variance, i.e. $\forall
i:\mathbb{E}\left[Y_{i}^{2}\right]<\infty$.

  * (A1)

$\mathbb{E}[Y|\boldsymbol{X}]=m\left(\boldsymbol{X},\boldsymbol{\beta}_{0}\right)$
almost surely for some $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$ and
continuous $m$ (as a function from $\mathbb{R}^{p}\times\mathbb{R}^{p}$ to
$\mathbb{R}$), with $\hat{\boldsymbol{\beta}}$ being a consistent estimator
for $\boldsymbol{\beta}_{0}$.

  * (A2)

$\mathbb{P}[Y\in dy|\boldsymbol{X},\delta]=\mathbb{P}[Y\in dy|\boldsymbol{X}]$
almost surely, i.e. $Y$ is independent of the sampling process for the non-
probability sample conditional on $\boldsymbol{X}$, and for every event
$M_{j,i_{1},\dots,i_{k}}$ event depending only on units $j,i_{1},\dots,i_{k}$
in the population, if $\mathbb{P}(M_{j,i_{1},\dots,i_{k}})>0$ for infinitely
many tuples $(i_{1},\dots,i_{k})$ and we set $j$ from the super-population
then:

$\lim_{n_{A}\rightarrow\infty}\mathbb{P}[1\leq\\#\\{j,i_{1},\dots,i_{k}\in
S_{A}:M_{i_{1},\dots,i_{k}}\\}|j\in S_{B}]=1.$

In this case we consider event $M$ to depend only on units
$j,i_{1},\dots,i_{k}$ if and only if it belongs to $\sigma$-algebra generated
by $\\{\mathmybb{1}_{\\{i_{1}\in S_{A}\\}},\mathmybb{1}_{\\{i_{1}\in
S_{B}\\}},\boldsymbol{x}_{i_{1}},y_{i_{1}},\ldots,\mathmybb{1}_{\\{j\in
S_{A}\\}},\mathmybb{1}_{\\{j\in S_{B}\\}},\boldsymbol{x}_{j},y_{j}\\}$.

  * (A3)

The variable of interest is square-integradable (in the super-population
model). We study the asymptotic properties of (3) in the setting of a sequence
of populations such that:

$\exists C_{1},C_{2}\forall i:C_{1}\leq\frac{N\pi_{i}}{n_{B}}\leq C_{2},$

and each population in the sequence is a random sample from some the super-
population and where $\displaystyle\lim_{n_{A}\rightarrow\infty}n_{B}=\infty$
a.s.

  * (A4)

We require that:

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{j=1}^{N}\left(m(\boldsymbol{x}_{i},\boldsymbol{\beta}_{0})-y_{j}\right)=0\text{
in probability.}$

Note that by virtue of (A1) the same holds after replacing
$\boldsymbol{\beta}_{0}$ with $\hat{\boldsymbol{\beta}}$.

  * (A5)

Matching between $S_{A}$ and $S_{B}$ is asymptotically weakly uniquely
identifiable i.e. for each $i\in S_{B}$ asymptotically there exists at most 1
set of $y_{\hat{\nu}(i,1)},\dots,y_{\hat{\nu}(i,k)}$, meaning there is only
one choice for $\\{\hat{\nu}(i,1),\dots,\hat{\nu}(i,k)\\}$ or for every choice
of $\\{\hat{\nu}(i,1),\dots,\hat{\nu}(i,k)\\}$ we get the same
$y_{\hat{\nu}(i,1)},\dots,y_{\hat{\nu}(i,k)}$. In other words, with
probability 1 the following property holds:

$\forall i\in S_{B},j\in
S_{A}:\left(d\left(y_{i},y_{\hat{\nu}(i,k)}\right)=d\left(y_{i},y_{j}\right)\implies
y_{j}=y_{\hat{\nu}(i,k)}\right).$

This assumption is met, for instance, if at least one covariate is continuous
and with support on $\mathbb{R}$ and the regression coefficient $\beta$
associated with this variable is non-zero.

  * (B1)

(For $\hat{y}-\hat{y}$ matching only) For any $\varepsilon>0$ we assume that:

$\lim_{n_{A}\rightarrow\infty}\mathbb{P}\left[\left.\exists
j_{1},j_{2},\dots,j_{k}\in
S_{A}:\max_{t}\lVert\boldsymbol{x}_{i}-\boldsymbol{x}_{j_{t}}\rVert<\varepsilon\right|i\in
S_{B}\right]=1,$

i.e. asymptotically we will cover at least $k$ units that have similar
$\boldsymbol{X}$ values to $i\in S_{B}$.

  * (C1)

(For $\hat{y}-y$ matching only) The support of $\boldsymbol{X}$ in $S_{B}$ (or
equivalently in $U$) is such that:

$\forall\varepsilon>0,y:\left(\mathbb{P}\left[Y\in\left(y-\varepsilon,y+\varepsilon\right)\right]>0\implies\mathbb{P}\left[m(\boldsymbol{X},\boldsymbol{\beta})\in\left(y-\varepsilon,y+\varepsilon\right)\right]>0\right),$

for $\boldsymbol{\beta}$ in some arbitrarily small open neighbourhood of
$\boldsymbol{\beta}_{0}$.

#####  2.3.1 Discussion regarding the assumptions

Assumption (A3) describes a setting in which we explore consistency and is not
important for finite sample performance. Assumptions (B1) and (C1) are not
strictly necessary for reasonable performance in finite samples, but the PMM
estimator will perform better if the non-probability sample is reasonably
diverse, which is basically what (B1) and (C1) are getting at. Assumption (C1)
is satisfied if, for example, at least one covariate has support on
$\mathbb{R}$ and $m(\cdot,\boldsymbol{\beta})$ depends on this covariate and
$Y$ is continuous under reasonable link functions, or if $Y$ is discrete and
$m(\boldsymbol{X},\boldsymbol{\beta})=\exp\left(\boldsymbol{X}^{T}\boldsymbol{\beta}\right)$.

A reasonable question is whether assumption (A4) can be proved from other
assumptions. The answer is negative since we can take $Y$ to be independent
from $\boldsymbol{X}$ and $\forall t=1,\ldots,N:y_{t}\sim\mathcal{N}(0,t)$
independently, then $m\equiv 0$ and
$\displaystyle\frac{1}{N}\sum_{j=1}^{N}y_{j}\sim\mathcal{N}\left(0,\frac{N+1}{2}\right)$
which does not converge to $0$ in distribution (and therefore does not
converge in probability). Nevertheless, it is trivially met by almost all
regression models used in practice.

In proving consistency we do not assume positivity, contrary to [31] or [40].
Instead we have some very weak assumption as part of (A2) and assumption (B1)/
(C1) depending on the PMM estimator. (B1)/ (C1) is less restrictive and allows
discrete $\boldsymbol{x}$, as in [31]. In the case of discrete
$\boldsymbol{X}$ (e.g. age, gender, etc.) we do not require that for all units
$\displaystyle\mathbb{P}(\delta=1|\boldsymbol{X}=\boldsymbol{x}_{i})>0$, as
may happen in quota sampling or Internet surveys, where certain units may be
missing due to the lack of access to the Internet; instead, we only require
that each level of $\boldsymbol{X}$ will eventually be represented (with at
least $k$ units) as sample size increases to infinity. The case of continuous
variables is analogous, i.e. if we set $i\in U$ and $\varepsilon>0$ we require
that eventually sample $S_{A}$ will contain at least $k$ units such that
distance between a $\boldsymbol{x}_{i}$ and $\boldsymbol{x}$ corresponding to
these $k$ units will be no greater than $\varepsilon$. This convergence may in
principal be arbitrarily slow (although this will have its impact on speed of
convergence for $\hat{\mu}$).

###  3 Main results

####  3.1 Consistency

The proof strategy, which is somewhat similar to the proof of consistency
presented in [31], can essentially be reduced to first proving that (3) and
the hypothetical Horvitz-Thompson (HT) estimator, i.e. what would be the HT
estimator if we knew values of $Y$ for units in $S_{B}$, have the same limit,
and then inferring the consistency of (3) from the consistency of the
hypothetical HT in the setting described in section 2.

######  Theorem 1.

Under assumptions (A1)-(A5) and (C1) for $\hat{y}-y$ matching and (B1) for
$\hat{y}-\hat{y}$ matching, estimators in (3) are consistent. Additionally,
the modified estimators with unknown population size:

$\hat{\mu}=\frac{1}{\displaystyle\sum_{r\in S_{B}}\pi_{r}^{-1}}\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}.$

are also consistent.

###### Proof.

To establish consistency first notice that for:

$\displaystyle|\hat{\mu}-\mu|$ $\displaystyle=\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j=1}^{N}y_{j}\right|$
(4) $\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right|,$ (5)

both (4) and (5) have the same limit in probability almost surely, since by
inverse triangle inequality:

$\displaystyle\left|\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j=1}^{N}y_{j}\right|-\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right|\right|\leq$
$\displaystyle\left|\frac{1}{N}\left(\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j=1}^{N}y_{j}\right)-\frac{1}{N}\left(\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right)\right|=$
$\displaystyle\left|-\frac{1}{N}\sum_{j=1}^{N}y_{j}+\frac{1}{N}\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right|,$

and the HT estimator of the population mean is consistent in this setting, so
it is sufficient to check whether (3) is consistent for the hypothetical HT
estimator. By the same argument for

$\displaystyle|\hat{\mu}-\mu|$
$\displaystyle=\left|\frac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\frac{1}{N}\sum_{j=1}^{N}y_{j}\right|$
$\displaystyle\frac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right|,$

an analogous property holds if $N$ is unknown.

By inverse triangle inequality and (A4) for the expression in (5) the
following occurs:

$\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}y_{j}\right|\leq$ (6)
$\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|+\underbrace{\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\left(m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)-y_{i}\right)\right|}_{\displaystyle\overset{n_{A}\rightarrow\infty}{\longrightarrow}0\text{
in probability}},$

again we can replace $\dfrac{1}{N}$ by $\dfrac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}$ and obtain an analogous property. Thus we have reduced
the problem of proving the consistency of (3) to proving that:

$\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\left(\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|\longrightarrow
0\text{ in probability.}$ (7)

(Consistency of the PMM A estimator with $\hat{y}-y$) Notice that from the
continuity of $m$ we have the following implications:

$\displaystyle\forall\varepsilon^{\prime}>0:\mathbb{P}\left(\left.|\boldsymbol{x}_{i}-\boldsymbol{X}|<\varepsilon^{\prime}\right|\\{i\in
S_{B}\\}\right)>0,$ (8) $\displaystyle\implies$
$\displaystyle\forall\varepsilon>0,\boldsymbol{\beta}:\mathbb{P}\left(\left.|m(\boldsymbol{x}_{i},\boldsymbol{\beta})-m(\boldsymbol{X},\boldsymbol{\beta})|<\varepsilon\right|\\{i\in
S_{B}\\}\right)>0.$

and:

$\forall\varepsilon>0,\boldsymbol{\beta}:\\{\exists j\in
S_{B}:|y_{j}-m(\boldsymbol{x}_{i},\boldsymbol{\beta})|<\varepsilon\\}\cap\\{i\in
S_{B}\\},$

happens with positive probability. Notice for $\hat{y}-y$ matching
$\displaystyle\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)$
is exactly the minimized quantity from Algorithm 1. By virtue of the above
reasoning and assumptions (A2), (C1) we have:

$\displaystyle\forall\varepsilon>0,i:$
$\displaystyle\lim_{n_{A}\rightarrow\infty}\mathbb{P}\left[\left.k\leq\\#\left\\{j\in
S_{A}:\left|y_{j}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|<\varepsilon\right\\}\right|\\{i\in
S_{B}\\}\right]=1,$ (9)

from which one infers that:

$\displaystyle\lim_{n_{A}\rightarrow\infty}\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\left(\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|\leq$
$\displaystyle\lim_{n_{A}\rightarrow\infty}\max_{i\in
S_{B}}\left\\{\frac{n_{B}}{N\pi_{i}}\right\\}\max_{i\in
S_{B}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\leq\frac{\varepsilon}{C_{1}},$

so the $\hat{y}-y$ matching estimator from Algorithm 1 is consistent.

(Consistency of the PMM B estimator with $\hat{y}-\hat{y}$ matching) In the
case of $\hat{y}-\hat{y}$ matching by (A1), (A2), (B1) and (8) we have:

$\forall\varepsilon>0,i:\lim_{n_{A}\rightarrow\infty}\mathbb{P}\left[\left.k\leq\\#\left\\{j\in
S_{A}:\left|m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{j},\hat{\boldsymbol{\beta}}\right)\right|<\varepsilon\right\\}\right|\\{i\in
S_{B}\\}\right]=1.$ (10)

In other words, there will be at least $k$ units in $S_{A}$ that are "similar"
to $i\in S_{B}$ and we can asymptotically replace
$m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)$ with
$m\left(\boldsymbol{x},\hat{\boldsymbol{\beta}}\right)$ value of any of those
$k$-units.

Notice that the term in (7) has the same limit as:

$\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}\left(m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|\leq$
(11) $\displaystyle\max_{i\in
S_{B}}\left\\{\frac{n_{B}}{N\pi_{i}}\right\\}\max_{i\in
S_{B}}\left|\frac{1}{k}\sum_{t=1}^{k}\left(m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|,$

since by inverse triangle inequality

$\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|-\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\left|\frac{1}{k}\sum_{t=1}^{k}m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\right|\leq$
$\displaystyle\frac{1}{N}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}\left(y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)\right)\right|,$

which tends to $0$ in probability with $n_{A}\rightarrow\infty$ by the
consistency of the Horvitz-Thompson estimator for the random variable
$\displaystyle\frac{1}{k}\sum_{t=1}^{k}\left(y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)\right)$,
whose mean is $0$ since:

$\displaystyle\mathbb{E}\left[\left.\frac{1}{k}\sum_{t=1}^{k}\left(y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)\right)\right|\sigma\left(\\{\hat{\nu}(i,t)\\}_{i,t=1,\dots
0,k}\right)\right]$ $\displaystyle=0\text{ a.s.}$
$\displaystyle\implies\mathbb{E}\left[\frac{1}{k}\sum_{t=1}^{k}\left(y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)\right)\right]$
$\displaystyle=0$

by (A1). By virtue of (10) we have:

$\forall\varepsilon>0:\lim_{n_{A}\rightarrow\infty}\max_{i\in
S_{B}}\left|\frac{1}{k}\sum_{t=1}^{k}m\left(\boldsymbol{x}_{\hat{\nu}(i,t)},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|<\varepsilon,$

so by virtue of (11), $\hat{\mu}$ estimators from (3) are consistent.
Additionally, for unknown $N$:

$\displaystyle\frac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\left|\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right)\right|\leq\underbrace{\max_{i\in
S_{B}}\left\\{\frac{n_{B}\pi_{i}^{-1}}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\right\\}}_{\overset{{n_{A}\rightarrow\infty}}{\longrightarrow}0\text{
in probability}},$

since

$\displaystyle\max_{i\in
S_{B}}\left\\{\left|\frac{n_{B}\pi_{i}^{-1}}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|-\frac{n_{B}}{\pi_{i}N}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\right|\right\\}$
$\displaystyle\leq n_{A}\max_{i\in
S_{B}}\left\\{\frac{1}{\pi_{i}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\left|\frac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}-\frac{1}{N}\right|\right\\}$ $\displaystyle\leq\max_{i\in
S_{B}}\left\\{\frac{n_{B}}{N\pi_{i}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)\right|\right\\}\frac{1}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\left|\sum_{r\in
S_{B}}\pi_{r}^{-1}-N\right|\overset{n_{A}\rightarrow\infty}{\longrightarrow}0\text{
in probability},$

so if $N$ is unknown and estimated by $\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}$, corrected estimates for $\hat{\mu}$ from (3) are also
consistent. ∎

Since $\hat{\mu}$ and the hypothetical HT estimator have the same limit in
order for (3) to be asymptotically normal with
$\left(\mu^{(HT)}_{n_{B}},V^{(HT)}_{n_{B}}\right)$, it is both sufficient and
necessary for the HT estimator to be asymptotically normal with
$\left(\mu^{(HT)}_{n_{B}},V^{(HT)}_{n_{B}}\right)$.

######  Remark 3.

(Regarding the proof of Theorem 1) If $Y$ is multivariate we can, under
analogous assumptions (though these assumptions may be somewhat less realistic
in a multivariate case), replace $|\cdot|$ norm by the Euclidean norm
$\lVert\cdot\rVert$ and get the consistency of the PMM estimators in a
multivariate case.

Additionally, Proof 1 can be rewritten to accommodate distance functions other
than the Euclidean metric (under reasonable restrictions) almost in an exactly
analogous manner.

######  Remark 4.

(Regarding robustness of $\hat{y}-\hat{y}$ matching) If the assumption that
model $\mathbb{E}[Y|\boldsymbol{X}]=m(\boldsymbol{X},\boldsymbol{\beta}_{0})$
is correctly specified is violated and instead we have chosen the
$m^{\prime}(\boldsymbol{X},\cdot)$ model with $m$ and $m^{\prime}$ being
different, then in order for (3) to be consistent $m^{\prime}$ only needs to
satisfy (A1) and:

$\forall i\in
S_{B}:\bigcup_{j=1}^{k}\left\\{\hat{\nu}(i,j)\right\\}=\bigcup_{j=1}^{k}\left\\{\hat{\nu}^{\prime}(i,j)\right\\},$

where:

$\displaystyle\hat{\nu}(i,z)$
$\displaystyle=\operatorname*{arg\,min}_{\displaystyle j\in
S_{A}\setminus\bigcup_{t=1}^{z-1}\\{\hat{\nu}(i,t)\\}}\left|m\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}\right)-m\left(\boldsymbol{x}_{j},\hat{\boldsymbol{\beta}}\right)\right|,$
$\displaystyle\hat{\nu}^{\prime}(i,z)$
$\displaystyle=\operatorname*{arg\,min}_{\displaystyle j\in
S_{A}\setminus\bigcup_{t=1}^{z-1}\\{\hat{\nu}^{\prime}(i,t)\\}}\left|m^{\prime}\left(\boldsymbol{x}_{i},\hat{\boldsymbol{\beta}}^{\prime}\right)-m^{\prime}\left(\boldsymbol{x}_{j},\hat{\boldsymbol{\beta}}^{\prime}\right)\right|.$

In other words, the linear order created by $\hat{m}$ and $\hat{m}^{\prime}$
asymptotically agrees (at least approximately), which is realistically
possible.

For example, [38] discussed robustness of the $\hat{y}-\hat{y}$ PMM estimator
in cases when the regression model $m$ has a Lipschitz property i.e. when:

$\exists L\in\mathbb{R}\forall
i,j:\left|m\left(\boldsymbol{x}_{i},\boldsymbol{\beta}\right)-m\left(\boldsymbol{x}_{j},\boldsymbol{\beta}\right)\right|\leq
L\lVert\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\rVert.$

######  Remark 5.

(Regarding assumption (A1)) We can relax assumption (A1) by assuming that for
continuous $m$: $\mathbb{E}[Y|\boldsymbol{X}]=m(\boldsymbol{X})$ almost
surely, i.e. a continuous function is a parameter, not just a vector in the
Euclidean space. Then, under modified assumption (A1): $\hat{m}$ is consistent
for $m$ in the sense that:

$\forall\varepsilon>0:\lim_{n_{A}\rightarrow\infty}\mathbb{P}\left[\lVert
m-\hat{m}\rVert_{\infty}<\varepsilon\right]=\lim_{n_{A}\rightarrow\infty}\mathbb{P}\left[\sup_{\boldsymbol{x}\in\mathbb{R}^{p}}\lVert
m(\boldsymbol{x})-\hat{m}(\boldsymbol{x})\rVert<\varepsilon\right]=1,$

and almost surely continuous, we can establish the consistency of (3), where
matching is done using $\hat{m}$ with almost exactly the same argument.

This means that certain non-parametric regression techniques, such as kernel
regression via the Nadaraya–Watson estimator or even local linear regression
techniques, could potentially be used in (3) under suitable assumptions (to
ensure their consistency).

In the Appendix we present simulation study that covers these examples. We
have not determined whether convergence in $\lVert\cdot\rVert_{\infty}$ can be
replaced by convergence in weaker topology.

####  3.2 Exact variance

Let us now study the variance of (3) choosing the case when $Y$ is univariate
for the sake of simplicity. The results can easily be generalized to
multivariate $Y$ by replacing products of scalars with outer products of
vectors.

######  Theorem 2.

If the population size $N$ is known, the variance of estimators (3) is given
by

$\mathbb{V}\left[\hat{\mu}\right]=V_{1}+V_{2},$ (12)

where

$\displaystyle V_{1}$
$\displaystyle=\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\frac{1}{\pi_{i}}\left(\frac{1}{\pi_{i}}-1\right)\mathbb{E}\left[\left(\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)^{2}\right]$
$\displaystyle+\frac{1}{N^{2}k^{2}}\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{N}\pi_{i}^{-1}\pi_{j}^{-1}\pi_{ij}^{-1}\mathbb{E}\left[\text{cov}(\mathmybb{1}_{\\{i\in
S_{B}\\}},\mathmybb{1}_{\\{j\in
S_{B}\\}})\sum_{t^{\prime},t=1}^{k}y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right]$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{i}\pi_{j}\pi_{ij}}\mathbb{E}\left[\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)\left(\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right)\right],$
(13)

and

$V_{2}=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\text{cov}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right)=\mathbb{V}\left[\frac{1}{N}\sum_{i=1}^{N}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right].$
(14)

The proof of Theorem 2 is given in Appendix A. The $V_{1}$ term is just a
variance of the HT estimator for the mean of imputed values via Algorithm 1
and the $V_{2}$ term can be seen as a compensation for inherent randomness
(induced into $S_{A}$ by unknown sampling, regression estimation etc.)
resulting from PMM imputation.

######  Remark 6.

(Regarding Theorem 2) When deriving the expression for
$\mathbb{V}\left[\hat{\mu}\right]$ we conditioned variance on knowing values
of $y_{i}$ and $\hat{\nu}(i,t)$ for $i=1,\dotso,N$, $t=1,\dotso,k$. We can
choose to do the same for nearest neighbour imputation (if we treat $k$ as a
random variable and not a hyper-parameter we additionally have to condition on
$k$) in the spirit of [40]. This gives us an exactly analogous expression for
$V_{1},V_{2}$ with matching variables $\hat{\nu}(i,t)$ being constructed by
the kNN algorithm from which exact variance for kNN (conditional on $k$ if it
is a random variable) imputation can be obtained (a finite population
correction $V_{2}$ was not considered in [40]).

As a corollary of Theorem 2 we have the following approximation for (3) with
unknown population size estimated via the usual
$\displaystyle\hat{N}=\sum_{r\in S_{B}}\pi_{r}^{-1}$:

$\displaystyle\mathbb{V}\left[\frac{N}{\displaystyle\sum_{r\in
S_{B}}\pi_{r}^{-1}}\hat{\mu}\right]\approx
V_{1}+V_{2}+\mathbb{E}\left[\hat{\mu}\right]^{2}\frac{\mathbb{V}\left[\hat{N}\right]}{N^{2}}-2\frac{\text{cov}\left(\hat{N},\hat{\mu}\right)}{N}$
$\displaystyle=V_{1}+V_{2}+\frac{\mathbb{E}\left[\hat{\mu}\right]^{2}}{N^{2}}\sum_{i,j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{i}\pi_{j}}-2\frac{1}{N^{2}}\sum_{i,j=1}^{N}\pi_{i}^{-1}\pi_{j}^{-1}\text{cov}\left(\mathmybb{1}_{i\in
S_{B}},\mathmybb{1}_{j\in
S_{B}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(j,t)}\right)$ (15)
$\displaystyle=V_{1}+V_{2}+\frac{\mathbb{E}\left[\hat{\mu}\right]^{2}}{N^{2}}\sum_{i,j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{i}\pi_{j}}-2\frac{1}{N^{2}}\sum_{i,j=1}^{N}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{i}\pi_{j}}\mathbb{E}\left[\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(j,t)}\right]$

Since $\displaystyle\hat{N}$ is unbiased for $N$ we obtain (3.2) by Taylor
series expansion.

####  3.3 Variance estimator

#####  3.3.1 "Analytic" variance estimator

Theorem 3 provides a consistent estimator of the variance estimator provided
in Theorem 2.

######  Theorem 3.

If second order inclusion probabilities are known, assumptions of Theorem 1
are satisfied, the standard variance estimator of the "hypothetical" HT
estimator is consistent and
$\displaystyle\left|\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\right|$ is
bounded for all $i,j$, then:

$\widehat{\mathbb{V}\left[\hat{\mu}\right]}=\hat{V}_{1}=\frac{1}{N^{2}}\sum_{i\in
S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\frac{\displaystyle\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}}{\pi_{i}}\frac{\displaystyle\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}}{\pi_{j}},$
(16)

is a consistent estimator for $\mathbb{V}\left[\hat{\mu}\right]$.

The full proof of Theorem 3 is presented in Appendix B.

The result $\displaystyle\lim_{n_{A}\rightarrow\infty}V_{2}=0$, which is a
corollary to Theorem 3, is intuitively expected since
$\displaystyle\mathbb{V}\left[\frac{1}{N}\sum_{i=1}^{N}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right]$
is an alternative expression for $V_{2}$ and under (B1) or (C1) asymptotically
we will almost surely find $k$ units that will not differ “too much” from each
$i\in S_{B}$ and therefore in large samples $y_{\hat{\nu}(i,1)}\approx
y_{\hat{\nu}(i,2)}\approx\dotso\approx y_{\hat{\nu}(i,k)}$, which implies that
$\displaystyle\forall
i:\mathbb{V}\left[\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(i,t^{\prime})}\right]\approx
0$, therefore variance of the sum of
$\displaystyle\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}$
will not grow faster than $N^{2}$.

The first component $\hat{V}_{1}$ is the Horvitz-Thomspon variance estimator
for the mean computed via imputed values. Computation of the term
$\hat{V}_{2}$, which in the light of Theorem 3 can be seen as a finite sample
correction, given by:

$\hat{V}_{2}=\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\pi_{ij}^{-1}\widehat{\text{cov}}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right),$

is more challenging as it requires an estimator for covariance of imputed
values. Fortunately, as we demonstrate in the simulation study we can omit the
$\hat{V}_{2}$ term for large $S_{A}$ samples.

If the population size $N$ is unknown and estimated from the probability
sample we can use the plug-in estimator for (3.2) given by

$\displaystyle\frac{1}{\hat{N}^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\frac{\displaystyle\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}}{\pi_{i}}\frac{\displaystyle\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}}{\pi_{j}}$
$\displaystyle+\frac{1}{\hat{N}^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\pi_{ij}^{-1}\widehat{\text{cov}}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right),$
$\displaystyle+\frac{\hat{\mu}^{2}}{\hat{N}^{2}}\sum_{i,j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{i}\pi_{j}}-2\frac{\hat{\mu}}{\hat{N}^{2}}\sum_{i,j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}$

where the first two lines refer to estimators for $V_{1}$ and $V_{2}$ under
estimated population size and the last line is a correction for randomness of
$\hat{N}$.

Unfortunately, the $V_{2}$ term cannot be estimated directly from samples
$S_{A}$ and $S_{B}$; therefore we propose a simple "mini-bootstrap" to
estimate
$\displaystyle\text{cov}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right)$.
The procedure is explained in Algorithm 2.

1 Sample $n_{A}$ units from $S_{A}$ with replacement to create
$S_{A}^{\prime}$ (if pseudo-weights are present inclusion probabilities should
be proportional to their inverses).

2 Estimate regression model
$\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}^{\prime}]=m(\boldsymbol{x}^{\prime},\boldsymbol{\beta})$
based on $j\in S_{A}^{\prime}$ from step 1

3 Compute $\hat{\nu}^{\prime}(i,t)$ for $t=1,\dots,k,i\in S_{B}$ using
estimated $m(\boldsymbol{x}^{\prime},\cdot)$ and
$\\{(y_{j},\boldsymbol{x}_{j})|j\in S_{A}^{\prime}\\}$

4 Compute $\displaystyle\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}^{\prime}(i)}$
using $Y$ values from $S_{A}^{\prime}$.

5 Repeat steps 1-4 $M$ times (we set $M=50$ in our simulations)

6 Estimate covariance between imputed values for each pair $i,j\in S_{B}$
using constructed pseudo-sample with values from 4.

Algorithm 2 Non-parametric mini-bootstrap estimator for covariance terms

The estimates from Algorithm 2 are satisfactory for $\hat{y}-y$ matching and a
little pessimistic for $\hat{y}-\hat{y}$ matching judging by numerical
simulations from section 4. But the performance of Algorithm 2 is not
explained by any solid theoretical argument.

The intuitive reason as to why Algorithm 2 may be an appropriate bootstrap for
$\hat{V}_{2}$ is that as per equation (14), $V_{2}$ is just the variance of
the mean of imputed values for the whole population. Since imputation to the
whole population can be done with just the $S_{A}$ sample (if we know
$\boldsymbol{X}$ values for the entire population), the sampling for sample
$A$ is the only source of randomness that is relevant to $V_{2}$. Therefore,
we only need to account for sampling in $S_{A}$ and use the probability sample
$S_{B}$ to estimate the
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}$
by the usual weighted estimator and compute the variance of estimates obtained
via Algorithm 2. This also explains why results from Algorithm 2 are a little
pessimistic since we are essentially estimating

$\mathbb{V}\left[\frac{1}{N}\sum_{i=1}^{N}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right],$

by an estimator for:

$\mathbb{V}\left[\frac{1}{N}\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right],$

which induces more randomness.

####  3.4 Bootstrap variance estimator

In this section we provide a bootstrap variance estimator for the mean. Since
we estimate $V_{2}$ using the “mini-bootstrap” described by Algorithm 2, we
need to estimate $V_{1}$ using the appropriate bootstrap approach for
probability samples.

1 Sample with replacement $n_{A}$ units from $S_{A}$ to create
$S_{A}^{\prime}$ (if pseudo-weights are present, inclusion probabilities
should be proportional to their inverses).

2 Sample with replacement $n_{B}$ units from $S_{B}$ according to the sampling
design to create $S_{B}^{\prime}$.

3 Estimate regression model parameters
$\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x}^{\prime}]=m(\boldsymbol{x}^{\prime},\boldsymbol{\beta})$
based on $j\in S_{A}^{\prime}$ from step 1

4 Follow steps 2A/2B-3 from Algorithm 1 depending on the imputation method
($\hat{y}-\hat{y}$ or $\hat{y}-y$) to compute $\hat{\mu}$ in current
iteration.

5 Repeat steps 1-4 $L$ times (we set $L=500$ in our simulations).

6 Estimate $V_{2}$ term using Algorithm 2.

7 Estimate variance $V_{1}$ term as
$\displaystyle\hat{V}_{1}=\frac{1}{L-1}\sum_{i=1}^{L}\left(\hat{\mu}_{i}-\bar{\mu}\right)^{2}$
where $\displaystyle\bar{\mu}=\frac{1}{L}\sum_{i=1}^{L}\hat{\mu}_{i}$ and
final bootstrap variance as $\hat{V}=\hat{V}_{1}+\hat{V}_{2}$

Algorithm 3 Bootstrap variance estimator

In the simulation study we compare the analytic and the bootstrap approach and
show that both estimators have the assumed nominal coverage rate.

###  4 Simulation study

Our simulation study follows the procedure described in [31]. We generate a
finite population of size $N=10^{5},$ random variables $X\sim\mathcal{N}(2,1)$
and $\varepsilon\sim\mathcal{N}(0,1)$ and assume the following models: 1)
$y_{1i}=1+2x_{i}+\varepsilon_{i}$; 2) $y_{2i}=3+x_{i}+2\varepsilon_{i}$; and
3) $y_{3i}=2.5+0.5x_{i}^{2}+\varepsilon_{i}$. We consider two sample sizes for
the non-probability sample $n_{A}=500$ and $n_{A}=1,000$. Sample $S_{A}$ is
drawn from the population using simple random sampling within two strata
created by $x_{i}\leq 2$ and $x_{i}>2$ with the following sample sizes:
$n_{strata1}=0.7n_{A}$ an $n_{strata2}=0.3n_{A}$, respectively. For the
probability sample we assume a fixed sample size $n_{B}=500$ and simple random
sampling without replacement from the population. In each of $500$ simulation
runs and for every $y$ we calculated

  * •

The naive estimator (sample mean) from sample $A$ given by
$\displaystyle\hat{\mu}=n_{A}^{-1}\sum_{i\in S_{A}}y_{i}$,

  * •

The mass imputation estimator using a linear regression model
$\hat{y}_{i}=\hat{\beta}_{0}+\hat{\beta}_{1}x_{i}$, where
$(\hat{\beta}_{0},\hat{\beta}_{1})$ are regression coefficients estimated from
the non-probability sample. We denote this estimator as MI-GLM.

  * •

Two PMM mass imputation estimators given by (3) using $k=1$ denoted as PMM1A
for $\hat{y}-y$ and PMM1B for $\hat{y}-\hat{y}$ matching,

  * •

Two PMM mass imputation estimators given by (3) using $k=5$ denoted as PMM5A
for $\hat{y}-y$ and PMM5B for $\hat{y}-\hat{y}$ matching.

Three results are reported in this section. Firstly, we report Monte Carlo
bias (Bias), standard error (SE), relative mean square error (RMSE) and
empirical coverage rate (CR) of 95% confidence intervals using the analytical
variance estimator based on $R=500$ simulations for each $Y$ variables:
$\text{Bias}=\bar{\hat{\mu}}-\mu$,
$\text{SE}=\sqrt{\frac{\sum_{r=1}^{R}\left(\hat{\mu}^{(r)}-\bar{\hat{\mu}}\right)^{2}}{R-1}}$
and $\text{RMSE}=\sqrt{\text{Bias}^{2}+\text{SE}^{2}}$, where
$\bar{\hat{\mu}}=\sum_{r=1}^{R}\hat{\mu}^{(r)}/R$ and
$\hat{\mu}^{\left(r\right)}$ is an estimate of the mean in the $r$-th
replication. Secondly, we compare how the CI’s change if we drop $\hat{V}_{2}$
from the PMM variance estimator. Finally, we compare analytic and bootstrap
estimators for the PMM estimators discussed in this section.

Table 4.1 presents the main results of the simulation study. For the first
variable ($Y_{1}$), where the model is correctly specified, values of both the
bias and the RMSE are similar to the MI-GLE estimator proposed by [31]. The
main difference between the proposed PMM estimators is the standard error,
which is larger for $\hat{y}-\hat{y}$ imputation compared to the standard
$\hat{y}-y$ imputation. This relationship is particularly visible for the
larger sample ($n=1000$). In the case of $Y_{1}$, there is no significant
difference between the proposed PMM $\hat{y}-y$ and the MI-GLM estimator.

Table 4.1: Simulated Bias, SE, RMSE (multiplied by 100) and CI of 5 estimators Estimator | $Y_{1}$ | $Y_{2}$ | $Y_{3}$
---|---|---|---
| Bias | SE | RMSE | CR | Bias | SE | RMSE | CR | Bias | SE | RMSE | CR
$n_{A}=500$
Naive | -64.36 | 6.83 | 64.73 |  | -32.29 | 8.81 | 33.47 |  | -64.33 | 6.54 | 64.66 |
MI-GLM | -0.52 | 10.33 | 10.34 | 94.20 | -0.56 | 10.21 | 10.22 | 95.60 | -6.34 | 10.44 | 12.22 | 89.60
PMM1A | -0.59 | 10.32 | 10.34 | 95.00 | -0.56 | 10.20 | 10.22 | 96.20 | -6.07 | 10.43 | 12.06 | 89.00
PMM1B | -0.60 | 12.21 | 12.22 | 93.00 | -0.46 | 16.66 | 16.67 | 92.80 | -0.97 | 12.51 | 12.54 | 93.00
PMM5A | -0.82 | 10.25 | 10.28 | 95.00 | -0.57 | 10.20 | 10.22 | 95.80 | -5.90 | 10.40 | 11.95 | 89.60
PMM5B | -1.11 | 10.81 | 10.87 | 94.60 | -1.02 | 12.13 | 12.18 | 95.00 | -2.37 | 11.09 | 11.34 | 94.60
$n_{A}=1,000$
Naive | -64.01 | 5.21 | 64.22 |  | -31.89 | 6.88 | 32.62 |  | -64.01 | 4.97 | 64.20 |
MI-GLM | -0.32 | 9.88 | 9.89 | 94.80 | -0.14 | 8.27 | 8.28 | 93.80 | -6.01 | 9.50 | 11.24 | 88.20
PMM1A | -0.34 | 9.87 | 9.87 | 95.00 | -0.14 | 8.28 | 8.28 | 94.40 | -5.83 | 9.48 | 11.13 | 88.80
PMM1B | -0.66 | 11.25 | 11.27 | 95.60 | -0.74 | 13.34 | 13.36 | 96.00 | -0.79 | 11.70 | 11.73 | 95.60
PMM5A | -0.43 | 9.84 | 9.85 | 94.80 | -0.14 | 8.27 | 8.27 | 93.40 | -5.66 | 9.45 | 11.01 | 88.80
PMM5B | -0.59 | 10.11 | 10.13 | 95.00 | -0.37 | 9.49 | 9.50 | 96.20 | -1.19 | 10.54 | 10.61 | 94.80

As regards the second variable $Y_{2}$, a similar pattern can be observed as
for $Y_{2}$ with one exception. The larger the number of nearest neighbours,
the smaller the bias and the standard error for the $\hat{y}-\hat{y}$ matching
PMM estimator. In the case of $\hat{y}-y$ matching PMM, increasing the number
of nearest neighbours has no significant effect on the bias or the standard
error. The empirical coverage rates are close to the nominal 95%.

The main differences can be observed for the last variable $Y_{3}$, which
includes the non-linear part (i.e. $x_{i}^{2}$). The PMM B with
$\hat{y}-\hat{y}$ fitting is now characterised by the lowest bias of all
estimators under consideration and a CR close to the nominal 95%. This is in
line with Remark 4 regarding robustness of $\hat{y}-\hat{y}$ matching. MI-GLM
and the PMM with $\hat{y}-y$ are characterised with similar bias, SE, RMSE and
CR.

Table 4.2: Empirical coverage rate intervals for two cases: the effect of $V_{2}$ on the variance estimator and analytic and bootstrap variance estimators. | $Y_{1}$ | $Y_{2}$ | $Y_{3}$
---|---|---|---
Effect of $V_{2}$ on the variance estimator and its empirical coverage rate
| $V_{1}$ | $V_{1}+V_{2}$ | $V_{1}$ | $V_{1}+V_{2}$ | $V_{1}$ | $V_{1}+V_{2}$
$n_{A}=500$
PMM1A | 90.80 | 95.00 | 58.60 | 96.20 | 81.00 | 89.00
PMM1B | 89.40 | 93.00 | 78.20 | 92.80 | 89.60 | 93.00
PMM5A | 90.60 | 95.00 | 58.80 | 95.80 | 80.80 | 89.60
PMM5B | 91.20 | 94.60 | 71.00 | 95.00 | 90.20 | 94.60
$n_{A}=1,000$
PMM1A | 94.40 | 95.00 | 76.20 | 94.40 | 85.00 | 88.80
PMM1B | 93.80 | 95.60 | 87.20 | 96.00 | 93.80 | 95.60
PMM5A | 94.40 | 94.80 | 76.00 | 93.40 | 85.20 | 88.80
PMM5B | 94.00 | 95.00 | 83.20 | 96.20 | 92.60 | 94.80
Comparison of analytic and bootstrap variance estimators
| Analytic | Bootstrap | Analytic | Bootstrap | Analytic | Bootstrap
$n_{A}=500$
PMM1A | 95.00 | 94.40 | 96.20 | 95.80 | 89.00 | 89.40
PMM1B | 93.00 | 93.60 | 92.80 | 96.00 | 93.00 | 93.40
PMM5A | 95.00 | 94.80 | 95.80 | 96.00 | 89.60 | 89.80
PMM5B | 94.60 | 94.60 | 95.00 | 96.20 | 94.60 | 94.40
$n_{A}=1,000$
PMM1A | 95.00 | 94.60 | 94.40 | 94.00 | 88.80 | 89.20
PMM1B | 95.60 | 95.20 | 96.00 | 97.80 | 95.60 | 95.80
PMM5A | 94.80 | 95.40 | 93.40 | 93.60 | 88.80 | 88.80
PMM5B | 95.00 | 95.80 | 96.20 | 97.80 | 94.80 | 95.00

Table 4.2 shows two results. The first part of the table shows the effect of
the second component $V_{2}$ on the variance estimator of the mean. The second
compares the analytic and bootstrap variance estimators as discussed in the
previous section. As suggested by theory, the proportion of variance assigned
to the $V_{2}$ component decreases with sample size. Moreover, only for the
$Y_{2}$ variable the impact of $V_{2}$ is not negligible, the CR based only on
the probability part ($V_{1}$) is significantly smaller than the expected 95%,
with values around 60% for $\hat{y}-y$ and 80-70% for $\hat{y}-\hat{y}$
matching. The second part suggests that the proposed analytical and bootstrap
variance estimation approaches yield the same results for all study variables
and sample sizes. As expected, the number of nearest neighbours does not
affect the CR as only small differences are observed.

In the Appendix we provide additional results: 1) choosing $k$ number of
nearest neighbours does not improve coverage but lowers MSE (as was also shown
by [40] for the NN estimator); 2) variable selection for the MI-GLM and PMM
estimators using the SCAD penalty (as in [39]) improves estimates and the
coverage rate and should be considered in applications; 3) non-parametric
regression using a local linear model can be employed, although in practice
better non-parametric methods ought to be considered; 4) the proposed
estimators perform similarly to the MI-GLM estimator when the positivity
assumption is violated with some exceptions when a non-linear model is
considered; and finally 5) the initial result on the multiply robust procedure
proposed by [24] seems to work well but provides comparable results to
standard PMM B matching (i.e. without the property of being multiply robust).
This, however, requires further studies as no theory has been developed for
the multiply robust PMM estimator.

###  5 Empirical study

In this section we present an attempt to integrate administrative and survey
data about job vacancies for the end of 2022Q2 in Poland. The aim was to
estimate the share of vacancies aimed at Ukrainian workers. We defined our
target variable as follows: $Y$ indicates whether the vacancy has been
translated into Ukrainian. As the survey only provides information about the
number of vacancies and not specific vacancies, our models were based on the
independent variable $y$ defined as the share of vacancies translated into
Ukrainian calculated separately for each unit. We present a brief description
of the datasets, but we encourage readers to read the full description of the
dataset and the relationships between the target and auxiliary variables as
described in [23].

The first source we used is the Job Vacancy Survey (JVS, also known as the
Labour Demand Survey), which is a stratified sample of 100,000 units, with a
response rate of about 60% ($S_{B}$). The survey population consists of
companies and their local units with 1 or more employees. The sampling frame
includes information about NACE (19 levels), region (16 levels), sector (2
levels), size (3 levels) and the number of employees according to
administrative data integrated by Statistics Poland (RE). The JVS sample
contains 304 strata created separately for enterprises with up to 9 employees
and those with more than 10 employees [36, cf.]. Of the 60,000 responding
units, around 7,000 reported at least one vacancy. Our target population
included units with at least one vacancy; according to the survey, there were
between 38,000 and 43,000 of such units at the end of 2022Q2.

The second source is the Central Job Offers Database (CBOP), which is a
register of all vacancies submitted to Public Employment Offices (PEOs –
$S_{A}$). CBOP is available online and can be accessed via API. CBOP includes
all types of contracts and jobs outside Poland, so data cleaning was carried
out to align the definition of a vacancy with that used in the JVS. CBOP data
collected via API include information about whether a vacancy is outdated
(e.g. 17% of vacancies were outdated when downloaded at the end of 2022Q2).
CBOP also contains information about unit identifiers (REGON and NIP), so we
were able to link units to the sampling frame to obtain auxiliary variables
with the same definitions as those used in the survey (24% of records
contained no identifier because the employer has the right to withhold this
information). The final CBOP dataset contained about 8,500 units included in
the sampling frame.

The following variables were considered: Region (16 levels), NACE (19 levels),
Sector (2 levels), Size (3 levels), $\log$(RE), $\log$(# vacancies) (the
number of vacancies), $I$(# vacancies = 1) (whether employer seeks only one
person). We considered 5 estimators in the empirical study (for all models
linear regression is assumed)

  * •

Mass imputation (MI) estimators: MI-GLM [31], PMM A, PMM B ,

  * •

Inverse probability weighting (IPW) estimator with calibration constraints
[26, cf.],

  * •

Doubly robust (DR) estimator with IPW defined above [26, cf.].

Variance was estimated using the following bootstrap approach: 1) JVS sample
was resampled using a stratified bootstrap approach, 2) CBOP was resampled
using simple random sampling with replacement. This procedure was repeated
$B=500$ times. Table 5.1 shows point estimates (denoted as Point), bootstrap
standard errors (denoted as SE), the coefficients of variation (CV) and 95%
confidence intervals.

Table 5.1: Estimates of the share of job vacancies aimed at Ukrainians at the end of 2022Q2 in Poland Estimator | Point | SE | CV | 2.5% | 97.5%
---|---|---|---|---|---
Naïve | 20.51 | – | – | – | –
IPW | 21.79 | 0.60 | 2.75 | 20.57 | 22.96
MI-GLM | 22.72 | 0.57 | 2.51 | 21.49 | 23.79
PMM A | 22.79 | 0.57 | 2.51 | 21.58 | 23.84
PMM B | 22.79 | 0.57 | 2.51 | 21.58 | 23.84
DR | 21.79 | 0.60 | 2.75 | 20.57 | 22.96

As can be seen in Table 5.1, the point estimates produced by all the
estimators in the study are at a similar level and fluctuate around 22%, with
the exception of the Naïve estimator, for which the estimated share of job
vacancies aimed at Ukrainians is lower and amounts to 20.51%. CV values for
the discussed estimators are similar, with the lowest values for mass
imputation estimators. Results for MI-GLM, PMM A and PMM B in this empirical
study are not distinguishable.

###  6 Conclusions

We have discussed the asymptotic properties of two PMM estimators that can be
used to integrate data from non-probability and probability samples. We have
proved their consistency and derived analytical variance estimators. The
results obtained in the study hold not only for parametric but also for non-
parametric models used for matching. The variance estimators can be applied to
the NN estimator for non-probability samples.

In the simulation study, we have shown that the PMM estimator based on
$\hat{y}-\hat{y}$ matching is robust to model mis-specification and has a
lower bias than existing alternatives. This is in line with results on PMM
estimators for survey non-response or statistical matching. However, this gain
is achieved at the expense of an increase in variance, but only for linear
models. In additional simulation studies presented in the appendix, we
verified the use of PMM estimators for different settings.

A natural extension of the study is that we can maintain consistency if the
design weights $d_{i}$ are replaced with calibration weights, which is often
the case when working with survey sample datasets. The next interesting
question regards the choice of $k$ hyper-parameters, and some suggestions are
given in the appendix. Another problem worth investigating is whether we can
replace probability convergence in $\lVert\cdot\rVert_{\infty}$ with weaker
convergence for estimators of $\hat{m}$ from Remark 5\. In finite dimensional
spaces all norms and topologies induced by them are equivalent, so it only
matters for non-parametric regression models. Further studies may focus on the
best method for estimating $V_{2}$, for which we have proposed a "mini"
bootstrap.

### Acknowledgements

The authors’ work has been financed by the National Science Centre in Poland,
OPUS 22, grant no. 2020/39/B/HS4/00941. Codes to reproduce simulations from
the paper are available at https://github.com/ncn-foreigners/paper-nonprob-
pmm.

Piotr Chlebicki and Łukasz Chrostowski are master students at the Faculty of
Mathematics and Computer Science at Adam Mickiewicz University in Poznań,
Poland (as of March 2024).

Contribution: Piotr Chlebicki is the main author responsible for deriving
properties of the MI PMM estimators and proving Theorems 1-3. Łukasz
Chrostowski was responsible for implementing these and other methods for non-
probability samples in the nonprobsvy package in the R language. Maciej
Beręsewicz was responsible for the initial idea, selected simulation studies
and the empirical study and prepared the final version of the manuscript.

### References

  * [1] Jean-Francois Beaumont  “Are probability surveys bound to disappear for the production of official statistics”  In _Survey Methodology_ 46.1 Statistics Canada, 2020, pp. 1–28
  * [2] Maciej Beręsewicz  “A two-step procedure to measure representativeness of internet data sources”  In _International Statistical Review_ 85.3 Wiley Online Library, 2017, pp. 473–493
  * [3] Maciej Beręsewicz and Marcin Szymkowiak  “Inference for non-probability samples using the calibration approach for quantiles”, https://arxiv.org/abs/2403.09726, 2024  arXiv:2403.09726 [stat.ME]
  * [4] Sixia Chen, David Haziza and Alexander Stubblefield  “A note on multiply robust predictive mean matching imputation with complex survey data”  In _Survey Methodology_ 47.1, 2021, pp. 215–222
  * [5] Sixia Chen, Shu Yang and Jae Kwang Kim  “Nonparametric Mass Imputation for Data Integration”  In _Journal of Survey Statistics and Methodology_ 10.1, 2022, pp. 1–24  DOI: 10.1093/jssam/smaa036
  * [6] Yilin Chen, Pengfei Li and Changbao Wu  “Doubly robust inference with nonprobability survey samples”  In _Journal of the American Statistical Association_ 115.532 Taylor & Francis, 2020, pp. 2011–2021
  * [7] Łukasz Chrostowski, Maciej Beręsewicz and Piotr Chlebicki  “ncn-foreigners/nonprobsvy: Initial release”  Zenodo, 2023  DOI: 10.5281/zenodo.10280114
  * [8] Constance F Citro  “From multiple modes for surveys to multiple data sources for estimates”  In _Survey Methodology_ 40.2 Statistics Canada, 2014, pp. 137–162
  * [9] Michael R. Elliott and Richard Valliant  “Inference for Nonprobability Samples”  In _Statistical Science_ 32.2, 2017  DOI: 10.1214/16-STS598
  * [10] Nicholas J Horton and Stuart R Lipsitz  “Multiple imputation in practice: comparison of software packages for regression models with missing variables”  In _The American Statistician_ 55.3 Taylor & Francis, 2001, pp. 244–254
  * [11] Jae Kwang Kim, Seho Park, Yilin Chen and Changbao Wu  “Combining Non-Probability and Probability Survey Samples Through Mass Imputation”  In _Journal of the Royal Statistical Society Series A: Statistics in Society_ 184.3, 2021, pp. 941–963  DOI: 10.1111/rssa.12696
  * [12] Roderick JA Little  “Missing-data adjustments in large surveys”  In _Journal of Business & Economic Statistics_ 6.3 Taylor & Francis, 1988, pp. 287–296
  * [13] R Core Team  “R: A Language and Environment for Statistical Computing”, 2023  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [14] Donald B Rubin  “Statistical matching using file concatenation with adjusted weights and multiple imputations”  In _Journal of Business & Economic Statistics_ 4.1 Taylor & Francis, 1986, pp. 87–94
  * [15] Nathaniel Schenker and Jeremy MG Taylor  “Partially parametric techniques for multiple imputation”  In _Computational Statistics & Data analysis_ 22.4 Elsevier, 1996, pp. 425–446
  * [16] Statistics Poland  “Methodological report The demand for labour”, 2021  URL: https://stat.gov.pl/obszary-tematyczne/rynek-pracy/popyt-na-prace/zeszyt-metodologiczny-popyt-na-prace,3,1.html
  * [17] Changbao Wu  “Statistical inference with non-probability survey samples”  In _Survey Methodology_ 48, 2022, pp. 283–311
  * [18] Shu Yang and Jae Kwang Kim  “Asymptotic theory and inference of predictive mean matching imputation using a superpopulation model framework”  In _Scandinavian Journal of Statistics_ 47.3, 2020, pp. 839–861  DOI: 10.1111/sjos.12429
  * [19] Shu Yang, Jae Kwang Kim and Rui Song  “Doubly Robust Inference when Combining Probability and Non-Probability Samples with High Dimensional Data”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 82.2, 2020, pp. 445–465  DOI: 10.1111/rssb.12354
  * [20] Shu Yang, Jae-Kwang Kim and Youngdeok Hwang  “Integration of data from probability surveys and big found data for finite population inference using mass imputation”  In _Survey Methodology_ 47, 2021, pp. 29–58

### References

  * [21] Jean-Francois Beaumont  “Are probability surveys bound to disappear for the production of official statistics”  In _Survey Methodology_ 46.1 Statistics Canada, 2020, pp. 1–28
  * [22] Maciej Beręsewicz  “A two-step procedure to measure representativeness of internet data sources”  In _International Statistical Review_ 85.3 Wiley Online Library, 2017, pp. 473–493
  * [23] Maciej Beręsewicz and Marcin Szymkowiak  “Inference for non-probability samples using the calibration approach for quantiles”, https://arxiv.org/abs/2403.09726, 2024  arXiv:2403.09726 [stat.ME]
  * [24] Sixia Chen, David Haziza and Alexander Stubblefield  “A note on multiply robust predictive mean matching imputation with complex survey data”  In _Survey Methodology_ 47.1, 2021, pp. 215–222
  * [25] Sixia Chen, Shu Yang and Jae Kwang Kim  “Nonparametric Mass Imputation for Data Integration”  In _Journal of Survey Statistics and Methodology_ 10.1, 2022, pp. 1–24  DOI: 10.1093/jssam/smaa036
  * [26] Yilin Chen, Pengfei Li and Changbao Wu  “Doubly robust inference with nonprobability survey samples”  In _Journal of the American Statistical Association_ 115.532 Taylor & Francis, 2020, pp. 2011–2021
  * [27] Łukasz Chrostowski, Maciej Beręsewicz and Piotr Chlebicki  “ncn-foreigners/nonprobsvy: Initial release”  Zenodo, 2023  DOI: 10.5281/zenodo.10280114
  * [28] Constance F Citro  “From multiple modes for surveys to multiple data sources for estimates”  In _Survey Methodology_ 40.2 Statistics Canada, 2014, pp. 137–162
  * [29] Michael R. Elliott and Richard Valliant  “Inference for Nonprobability Samples”  In _Statistical Science_ 32.2, 2017  DOI: 10.1214/16-STS598
  * [30] Nicholas J Horton and Stuart R Lipsitz  “Multiple imputation in practice: comparison of software packages for regression models with missing variables”  In _The American Statistician_ 55.3 Taylor & Francis, 2001, pp. 244–254
  * [31] Jae Kwang Kim, Seho Park, Yilin Chen and Changbao Wu  “Combining Non-Probability and Probability Survey Samples Through Mass Imputation”  In _Journal of the Royal Statistical Society Series A: Statistics in Society_ 184.3, 2021, pp. 941–963  DOI: 10.1111/rssa.12696
  * [32] Roderick JA Little  “Missing-data adjustments in large surveys”  In _Journal of Business & Economic Statistics_ 6.3 Taylor & Francis, 1988, pp. 287–296
  * [33] R Core Team  “R: A Language and Environment for Statistical Computing”, 2023  R Foundation for Statistical Computing  URL: https://www.R-project.org/
  * [34] Donald B Rubin  “Statistical matching using file concatenation with adjusted weights and multiple imputations”  In _Journal of Business & Economic Statistics_ 4.1 Taylor & Francis, 1986, pp. 87–94
  * [35] Nathaniel Schenker and Jeremy MG Taylor  “Partially parametric techniques for multiple imputation”  In _Computational Statistics & Data analysis_ 22.4 Elsevier, 1996, pp. 425–446
  * [36] Statistics Poland  “Methodological report The demand for labour”, 2021  URL: https://stat.gov.pl/obszary-tematyczne/rynek-pracy/popyt-na-prace/zeszyt-metodologiczny-popyt-na-prace,3,1.html
  * [37] Changbao Wu  “Statistical inference with non-probability survey samples”  In _Survey Methodology_ 48, 2022, pp. 283–311
  * [38] Shu Yang and Jae Kwang Kim  “Asymptotic theory and inference of predictive mean matching imputation using a superpopulation model framework”  In _Scandinavian Journal of Statistics_ 47.3, 2020, pp. 839–861  DOI: 10.1111/sjos.12429
  * [39] Shu Yang, Jae Kwang Kim and Rui Song  “Doubly Robust Inference when Combining Probability and Non-Probability Samples with High Dimensional Data”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 82.2, 2020, pp. 445–465  DOI: 10.1111/rssb.12354
  * [40] Shu Yang, Jae-Kwang Kim and Youngdeok Hwang  “Integration of data from probability surveys and big found data for finite population inference using mass imputation”  In _Survey Methodology_ 47, 2021, pp. 29–58

##  Part I Appendix

###  Appendix A Proof of Theorem 2

###### Proof.

Using law of total variance we have:

$\displaystyle\mathbb{V}\left[\hat{\mu}\right]$
$\displaystyle=\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\mathbb{V}\left[\left.\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]$ (17)
$\displaystyle+\frac{1}{N^{2}k^{2}}\mathbb{V}\left[\mathbb{E}\left[\left.\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right].$ (18)

by Cauchy product (18) can be expressed as:

$\displaystyle\frac{1}{N^{2}k^{2}}\mathbb{V}\left[\mathbb{E}\left[\left.\sum_{i=1}^{N}\frac{1}{\pi_{i}}\mathmybb{1}_{\\{i\in
S_{B}\\}}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]=\frac{1}{N^{2}k^{2}}\mathbb{V}\left[\sum_{i=1}^{N}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right]$
$\displaystyle=\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\left(\sum_{i=1}^{N}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)^{2}\right]-\frac{1}{N^{2}k^{2}}\left(\mathbb{E}\left[\sum_{i=1}^{N}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right]\right)^{2}$
$\displaystyle=\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\left(\sum_{i=1}^{N}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)\left(\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right)\right]$
$\displaystyle-\frac{1}{N^{2}k^{2}}\left(\sum_{i=1}^{N}\mathbb{E}\left[\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right]\right)\left(\sum_{j=1}^{N}\mathbb{E}\left[\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right]\right)$
$\displaystyle=\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right]$
$\displaystyle-\frac{1}{N^{2}k^{2}}\left(\sum_{i=N}^{N}\sum_{j=1}^{N}\mathbb{E}\left[\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right]\mathbb{E}\left[\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right]\right)$
$\displaystyle=\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{k}\sum_{t=1}^{k}\mathbb{E}\left[y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right]$
$\displaystyle-\frac{1}{N^{2}k^{2}}\left(\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{k}\sum_{t^{\prime}=1}^{k}\mathbb{E}\left[y_{\hat{\nu}(i,t)}\right]\mathbb{E}\left[y_{\hat{\nu}(j,t^{\prime})}\right]\right)$
$\displaystyle=\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t^{\prime}=1}^{k}\sum_{t=1}^{k}\left(\mathbb{E}\left[y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right]-\mathbb{E}\left[y_{\hat{\nu}(i,t)}\right]\mathbb{E}\left[y_{\hat{\nu}(j,t^{\prime})}\right]\right)$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\text{cov}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right)$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\text{cov}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right),$

which can be estimated by:

$\frac{1}{N^{2}}\sum_{i=1}^{n_{B}}\sum_{j=1}^{n_{B}}\frac{1}{\pi_{ij}}\widehat{\text{cov}}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right),$
(19)

assuming we have an estimate for:

$\text{cov}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)},\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}\right).$

The remaining term in (17) on the other hand is:

$\displaystyle\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\mathbb{V}\left[\left.\sum_{i\in
S_{B}}\frac{1}{\pi_{i}}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]=$
$\displaystyle\frac{1}{N^{2}k^{2}}\mathbb{E}\left[\mathbb{V}\left[\left.\sum_{i=1}^{N}\mathmybb{1}_{\\{i\in
S_{B}\\}}\frac{1}{\pi_{i}}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]=$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\frac{1}{\pi_{i}^{2}}\mathbb{E}\left[\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\mathbb{V}\left[\left.\mathmybb{1}_{\\{i\in
S_{A}\\}}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]+$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{N}\frac{1}{\pi_{i}\pi_{j}}\mathbb{E}\left[\sum_{t^{\prime},t=1}^{k}y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\text{cov}\left[\left.\mathmybb{1}_{\\{i\in
S_{A}\\}},\mathmybb{1}_{\\{j\in
S_{A}\\}}\right|\sigma\left(\\{y_{i},\hat{\nu}(i,t)\\}_{i\in
U,t\in\\{1,\dots,k\\}}\right)\right]\right]=$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\frac{1}{\pi_{i}^{2}}\mathbb{E}\left[\left(\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)^{2}\pi_{i}\left(1-\pi_{i}\right)\right]+$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{N}\frac{1}{\pi_{i}\pi_{j}}\text{cov}(\mathmybb{1}_{\\{i\in
S_{A}\\}},\mathmybb{1}_{\\{j\in
S_{A}\\}})\mathbb{E}\left[\sum_{t^{\prime},t=1}^{k}y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right]=$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{i=1}^{N}\frac{1-\pi_{i}}{\pi_{i}}\mathbb{E}\left[\left(\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)^{2}\right]+$
$\displaystyle\frac{1}{N^{2}k^{2}}\sum_{\begin{subarray}{c}i,j=1\\\ i\neq
j\end{subarray}}^{N}\frac{1}{\pi_{i}\pi_{j}}\text{cov}(\mathmybb{1}_{\\{i\in
S_{B}\\}},\mathmybb{1}_{\\{j\in
S_{B}\\}})\mathbb{E}\left[\sum_{t^{\prime},t=1}^{k}y_{\hat{\nu}(i,t)}y_{\hat{\nu}(j,t^{\prime})}\right],$

which (if second order inclusion probabilities are known) can be estimated by:

$\displaystyle\frac{1}{N^{2}}\sum_{i\in
S_{B}}\left(1-\pi_{i}\right)\frac{\displaystyle\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\right)^{2}}{\pi_{i}^{2}}+\frac{1}{N^{2}}\sum_{i\in
S_{B}}\sum_{\begin{subarray}{c}j\in S_{B}\\\ j\neq
i\end{subarray}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\frac{\displaystyle\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}}{\pi_{i}}\frac{\displaystyle\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}}{\pi_{j}}.$

∎

###  Appendix B Proof of Theorem 3

###### Proof.

Let $\hat{\mu}_{HT}$ denote the hypothetical Horvitz-Thompson estimate of
$\mu$. Since by virtue of Theorem 1 $|\hat{\mu}-\hat{\mu}_{HT}|\rightarrow 0$
in probability, it also follows that $\hat{\mu}\rightarrow\hat{\mu}_{HT}$ in
distribution and therefore:

$\lim\left|\mathbb{V}\left[\hat{\mu}\right]-\mathbb{V}\left[\hat{\mu}_{HT}\right]\right|=\lim\left|V_{1}+V_{2}-\mathbb{V}\left[\hat{\mu}_{HT}\right]\right|=0$

in the sequence of samples described in the assumptions underlying Theorem 1.
Furthermore:

$\displaystyle\left|\hat{V}_{1}-\widehat{\mathbb{V}\left[\hat{\mu}_{HT}\right]}\right|$
$\displaystyle\leq\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|$
$\displaystyle\leq\frac{n_{B}^{2}}{N^{2}}\max_{i,j\in
S_{B}}\left\\{\frac{1}{\pi_{i}\pi_{j}}\left|1-\frac{\pi_{i}\pi_{j}}{\pi_{ij}}\right|\right\\}\max_{i,j\in
S_{B}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|$
$\displaystyle\leq\frac{n_{B}^{2}}{N^{2}}\max_{i,j\in
S_{B}}\left\\{\pi_{ij}^{-1}\right\\}\max_{i,j\in
S_{B}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|$

By elementary inequalities:

$\displaystyle\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|=$
$\displaystyle\Bigg{|}\left(\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right)\left(\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right)-2y_{i}y_{j}+y_{i}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}+y_{j}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\Bigg{|}\leq$
$\displaystyle\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|+\left|y_{i}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}+y_{j}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-2y_{i}y_{j}\right|$

continuing this sequence:

$\displaystyle\left|y_{i}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}+y_{j}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-2y_{i}y_{j}\right|\leq$
$\displaystyle\left|y_{i}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|+\left|y_{j}\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}y_{j}\right|\leq\left|y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|+\left|y_{j}\right|\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|$

so if $\displaystyle\max_{i,j\in
S_{B}}\left\\{\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\right\\}$ is bounded:

$\displaystyle\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{i}y_{j}\right|\leq$
$\displaystyle\underbrace{\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|}_{\overset{n_{A}\rightarrow\infty}{\longrightarrow}0}+$
$\displaystyle\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|+\frac{1}{N^{2}}\sum_{i\in
S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|y_{j}\right|\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|$

the second line follows from:

$\displaystyle\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|\leq$
$\displaystyle\max_{i,j\in
S_{B}}\left\\{\left|\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\right|\right\\}\frac{1}{N}\sum_{i\in
S_{B}}\frac{1}{N}\frac{1}{\pi_{i}}\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|$

which is a product of a essentially bounded random variable and two random
variables converging to 0 in probability. By consistency of (3) and
consistency of HT estimator for
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}|y_{i}|$:

$\displaystyle\frac{1}{N^{2}}\sum_{i\in S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|y_{i}\right|\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|+\frac{1}{N^{2}}\sum_{i\in
S_{B}}\sum_{j\in
S_{B}}\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}\pi_{i}\pi_{j}}\left|y_{j}\right|\left|\frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i,t)}-y_{i}\right|\leq$
$\displaystyle 2\max_{i,j\in
S_{B}}\left\\{\frac{\pi_{ij}-\pi_{i}\pi_{j}}{\pi_{ij}}\right\\}\underbrace{\frac{1}{N}\sum_{i\in
S_{B}}\frac{\left|y_{i}\right|}{\pi_{i}}}_{\displaystyle\overset{n_{A}\rightarrow\infty}{\longrightarrow}\frac{1}{N}\sum_{i=1}^{N}|y_{i}|}\underbrace{\frac{1}{N}\sum_{j\in
S_{B}}\frac{1}{\pi_{j}}\left|\frac{1}{k}\sum_{t^{\prime}=1}^{k}y_{\hat{\nu}(j,t^{\prime})}-y_{j}\right|}_{\displaystyle\overset{n_{A}\rightarrow\infty}{\longrightarrow}0}$

since by assumption
$\widehat{\mathbb{V}\left[\hat{\mu}_{HT}\right]}\rightarrow\mathbb{V}\left[\hat{\mu}_{HT}\right],\left|V_{1}+V_{2}-\mathbb{V}\left[\hat{\mu}_{HT}\right]\right|\rightarrow
0$ and by inequalities above
$\left|\hat{V}_{1}-\widehat{\mathbb{V}\left[\hat{\mu}_{HT}\right]}\right|\rightarrow
0$ it follows that $\left|V_{1}+V_{2}-\hat{V}_{1}\right|\rightarrow 0$. ∎

###  Appendix C Additional simulations

####  C.1 Simple procedure for choosing $k$ hyper-parameter

In this simulation study we provide evidence for choosing the $k$ hyper-
parameter. We consider the following study design:

  * •

Population size: $N=100,000$,

  * •

Probability sample size: $n_{B}=500$,

  * •

Expected non-probability sample size: $\mathbb{E}\left[n_{A}/N\right]\approx
0.20$,

  * •

$X_{1}$ was generated from the multivariate normal distribution transformed to
1d vector (the reason for that was to introduce correlation between
observations):

$\begin{pmatrix}\boldsymbol{Z}_{1}^{T}&\boldsymbol{Z}_{2}^{T}&\dots&\boldsymbol{Z}_{N/5}^{T}\end{pmatrix}$
where
$\boldsymbol{Z}_{k}\sim\mathcal{N}_{5}\left(\begin{pmatrix}1&1&1&1&1\end{pmatrix}^{T},\boldsymbol{\Sigma}\right)$
($\boldsymbol{Z}_{k}$ being i.i.d) with

$\boldsymbol{\Sigma}=\begin{pmatrix}1&U_{12}&U_{13}&U_{14}&U_{15}\\\
U_{12}&1&U_{23}&U_{24}&U_{25}\\\ U_{13}&U_{23}&1&U_{34}&U_{35}\\\
U_{14}&U_{24}&U_{34}&1&U_{45}\\\ U_{15}&U_{25}&U_{35}&U_{45}&1\\\
\end{pmatrix}$

in setting where
$U_{12},\dotso\sim\mathcal{U}\left(-\tfrac{1}{2},\tfrac{1}{2}\right)$ i.i.d.,

  * •

$X_{2}\sim\text{Exp}(1)$ i.i.d,

  * •

Vector of $\epsilon$’s analogously to $x_{1}$ being
$\begin{pmatrix}\boldsymbol{R}_{1}^{T}&\boldsymbol{R}_{2}^{T}&\dots&\boldsymbol{R}_{N/5}^{T}\end{pmatrix}^{T}$
where
$\boldsymbol{R}_{k}\sim\mathcal{N}_{5}\left(\boldsymbol{0},\boldsymbol{\Sigma}\right)$
and $U_{12},\dotso\sim\mathcal{U}\left(-0.7,1\right)$ i.i.d.,

  * •

We generated two vector of probabilities for selection into the non-
probability sample $p_{1}=\dfrac{\exp(x_{2})}{1+\exp(x_{2})}$ and
$p_{2}=\dfrac{\exp(x_{1})}{1+\exp(x_{1})}$,

  * •

The inclusion into non-probability sample $S_{A}$ was generated as

$\displaystyle\delta=\min\left\\{\text{Bernoulli}(p_{1}),\text{Bernoulli}(p_{2}),\mathmybb{1}_{\\{\epsilon>\epsilon_{.8}\\}}+\mathmybb{1}_{\\{\epsilon<\epsilon_{.2}\\}}\right\\}$

where $\epsilon_{t}$ is an empirical quantile of $\epsilon$ of order $t$.

  * •

Two target variables were generated according to the following formulas:

$\displaystyle Y_{1}$
$\displaystyle=1+\tfrac{1}{2}x_{1}+\tfrac{35}{100}x_{2}+\epsilon$
$\displaystyle Y_{2}$
$\displaystyle=1.2+(x_{1}-0.5)^{2}+\arctan(x_{2})^{3+\sin(x_{1}+x_{2})}+\sin(x_{1})\cos(x_{2})+\epsilon$

and $k$ being chosen by minimising
$\widehat{\mathbb{V}\left[\hat{\mu}\right]}$ in each simulation iteration.

We consider three estimators discussed in the paper: MI-GLM, PMM A
($\hat{y}-y$ matching) and PMM B ($\hat{y}-\hat{y}$ matching) with fixed $k=5$
and dynamic selection of $k$.

Table C.1 shows the results of the simulation study. As regards the first
variable ($Y_{1}$), slight differences are observed between the estimators
under consideration. The dynamic selection is characterised by a smaller bias
and a slightly better RMSE for the PMM A estimator. In the case of the highly
non-linear variable $Y_{2}$, both MI-GLM and PMM A estimators are
significantly biased, while PMM B with either dynamic or fixed $k$ is
characterised by negligible bias. In addition, PMM B with dynamic selection of
$k$ is characterised by a smaller RMSE.

Table C.1: Results for the simulation study on choosing of $k$ hyper-parameter (all measures multiplied by 100) Variable | Estimator |  $k$-selection | Bias | SE | RMSE | CR
---|---|---|---|---|---|---
Y1 | MI-GLM | – | 0.2132 | 3.1015 | 3.1088 | 95.40
| PMM A | Dynamic | 0.2117 | 3.1022 | 3.1094 | 90.80
| PMM A | Fixed | 0.2133 | 3.1018 | 3.1092 | 95.40
| PMM B | Dynamic | 0.4676 | 5.2867 | 5.3073 | 99.60
| PMM B | Fixed | 0.5688 | 5.2856 | 5.3161 | 99.40
Y2 | MI-GLM | – | -20.2760 | 8.0450 | 21.8137 | 26.00
| PMM A | Dynamic | -20.2724 | 8.0422 | 21.8094 | 26.20
| PMM A | Fixed | -20.2698 | 8.0420 | 21.8069 | 26.20
| PMM B | Dynamic | -0.7583 | 9.3356 | 9.3664 | 96.20
| PMM B | Fixed | -0.7613 | 9.4749 | 9.5054 | 97.20

####  C.2 Simulation with variable selection

In this case our study designs follows that proposed by [39]:

  * •

finite population size was set to $N=10000$,

  * •

for each unit we generate
$X_{i}=\begin{pmatrix}1&X_{1,i}&\ldots&X_{p-1,i}\end{pmatrix}^{T}$ where
$p=50$ and $X_{i,t}\sim\mathcal{N}(0,1)$ i.i.d for each $i$ and
$t=0,\dots,p-1$,

  * •

from the finite population, we select a non-probability sample $\mathcal{A}$
of size $n_{\mathrm{A}}\approx 2000$, according to the selection indicator
$I_{\mathrm{A},i}\sim$
$\operatorname{Bernoulli}\left(\pi_{\mathrm{A},i}\right)$,

  * •

we select a probability sample $\mathcal{A}$ of the average size
$n_{\mathrm{B}}=500$ under Poisson sampling with
$\pi_{\mathrm{B},i}\propto\left(0.25+\left|X_{1i}\right|+0.03\left|Y_{i}\right|\right)$
(we have 4 set of $\pi_{\mathrm{B}}$ as we have 4 $Y$ variables defined
below),

  * •

for the non-probability sample inclusion probability, we consider both linear
and non-linear sampling score models:

    * –

$\operatorname{logit}\left(\pi_{\mathrm{B},i}\right)=\alpha_{0}^{\mathrm{T}}X_{i}$,
where
$\alpha_{0}=\begin{pmatrix}-2&1&1&1&1&0&0&\ldots&0\end{pmatrix}^{\mathrm{T}}$
(model PSM I)

    * –

$\operatorname{logit}\left(\pi_{\mathrm{B},i}\right)=3.5+\alpha_{0}^{T}\log\left(X_{i}^{2}\right)-\sin\left(X_{3,i}+X_{4,i}\right)-X_{5,i}-X_{6,i}$,
where $\alpha_{0}=\begin{pmatrix}0&0&0&3&3&3&3&0&\ldots&0\end{pmatrix}^{T}$
(model PSM II)

  * •

to generate a continuous outcome variable $Y_{i}$, we consider both linear and
non-linear outcome models with
$\beta_{0}=\begin{pmatrix}1&0&0&1&1&1&1&0&\ldots&0\end{pmatrix}^{T}$

    * –

$Y_{i}=\beta_{0}^{\mathrm{T}}X_{i}+\epsilon_{i},\epsilon_{i}\sim\mathcal{N}(0,1)$
i.i.d (model OM I, denoted as $Y_{11}$),

    * –

$Y_{i}=1+\exp\left(3\sin\left(\beta_{0}^{T}X_{i}\right)\right)+X_{5,i}+X_{6,i}+\epsilon_{i},\epsilon_{i}\sim\mathcal{N}(0,1)$
i.i.d (model OM II, denoted as $Y_{12}$),

  * •

to generate a binary outcome variable $Y_{i}$, we consider both linear and
non-linear outcome models with
$\beta_{0}=\begin{pmatrix}1&0&0&3&3&3&3&0&\ldots&0\end{pmatrix}^{T}$:

    * –

$Y_{i}\sim\operatorname{Bernoulli}\left(\pi_{Y}(X)\right)$ with logit
$\left(\pi_{Y}(X)\right)=\beta_{0}^{T}X$ (model OM III, denoted as $Y_{21}$),

    * –

$Y_{i}\sim\operatorname{Bernoulli}\left(\pi_{Y}(X)\right)$ with logit
$\left(\pi_{Y}(X)\right)=2-\log\left(\left(\beta_{0}^{T}X\right)^{2}\right)+2X_{5,i}+2X_{6,i}$
(model OM IV, denoted as $Y_{22}$).

  * •

variable selection was done using SCAD penalty.

Table C.2 shows the results of the simulation study with variable selection.
As expected, when the model is linear (OM I and OM III), MI-GLM and PMM B
($\hat{y}-\hat{y}$ matching) are almost unbiased with coverage close to the
nominal rate. The only difference is the PMM A ($\hat{y}-y$ matching)
estimator, which is characterised by greater bias for both linear and non-
linear selection mechanisms. Variable selection reduces the bias and slightly
reduces the standard error for the estimators discussed. As a result, the RMSE
for OM I and OM III is slightly lower than without variable selection.

Table C.2: Results of the simulation with variable selection (all numbers multiplied by 100) |  | Linear PS (PSM I) | Non-linear PS (PSM II)
---|---|---|---
Estimator | Selection | Bias | SE | RMSE | CR | Bias | SE | RMSE | CR
$Y_{11}$ (OM I)
MI-GLM | No | -0.24 | 10.18 | 10.18 | 95.00 | -0.61 | 9.75 | 9.77 | 96.20
| Yes | -1.03 | 9.97 | 10.02 | 95.60 | -0.80 | 9.78 | 9.81 | 96.20
PMM A | No | 0.03 | 10.12 | 10.12 | 93.80 | -0.61 | 9.75 | 9.77 | 95.20
| Yes | -0.74 | 9.93 | 9.96 | 94.60 | -0.79 | 9.77 | 9.80 | 95.60
PMM B | No | 0.55 | 10.22 | 10.24 | 93.40 | -0.92 | 10.04 | 10.08 | 95.20
| Yes | -0.11 | 9.91 | 9.92 | 94.40 | -0.77 | 10.11 | 10.14 | 94.80
$Y_{12}$ (OM II)
MI-GLM | No | 82.55 | 23.62 | 85.87 | 8.40 | -122.37 | 11.40 | 122.90 | 0.00
| Yes | 87.80 | 20.49 | 90.16 | 2.20 | -119.18 | 10.25 | 119.62 | 0.00
PMM A | No | 82.55 | 23.63 | 85.86 | 0.00 | -122.38 | 11.40 | 122.91 | 0.00
| Yes | 87.69 | 20.56 | 90.07 | 0.00 | -119.21 | 10.27 | 119.65 | 0.00
PMM B | No | 84.41 | 29.09 | 89.28 | 4.40 | -121.89 | 20.25 | 123.56 | 0.00
| Yes | 88.87 | 27.74 | 93.10 | 2.40 | -111.56 | 20.14 | 113.36 | 0.00
$Y_{21}$ (OM III)
MI-GLM | No | 0.27 | 1.97 | 1.99 | 96.60 | 0.29 | 1.54 | 1.57 | 97.00
| Yes | 0.38 | 1.75 | 1.79 | 96.60 | 0.33 | 1.52 | 1.55 | 97.40
PMM A | No | 4.60 | 2.91 | 5.44 | 52.20 | 4.76 | 2.40 | 5.33 | 52.00
| Yes | 5.19 | 2.73 | 5.87 | 43.40 | 5.47 | 2.40 | 5.97 | 38.60
PMM B | No | 0.20 | 2.18 | 2.18 | 89.80 | 0.53 | 1.78 | 1.86 | 93.40
| Yes | 0.45 | 1.96 | 2.01 | 92.60 | 0.64 | 1.76 | 1.87 | 92.20
$Y_{22}$ (OM IV)
MI-GLM | No | 0.84 | 2.06 | 2.23 | 94.20 | -7.39 | 1.67 | 7.57 | 2.40
| Yes | 1.47 | 1.94 | 2.44 | 87.80 | -7.42 | 1.65 | 7.60 | 1.40
PMM A | No | 3.36 | 2.66 | 4.28 | 69.60 | -5.52 | 2.57 | 6.09 | 44.60
| Yes | 4.18 | 2.52 | 4.88 | 55.60 | -5.32 | 2.60 | 5.92 | 48.00
PMM B | No | 0.82 | 2.20 | 2.35 | 89.00 | -8.16 | 1.85 | 8.36 | 0.40
| Yes | 1.45 | 2.10 | 2.56 | 83.20 | -8.23 | 1.93 | 8.45 | 1.40

In the case of non-linear OM II and OM IV, all the proposed estimators are
characterised by large bias and variance. Furthermore, in this case variable
selection leads to inconclusive results. For example, for the linear PSM I and
OM II, variable selection increases the bias, but for some models it reduces
the standard error. For the non-linear PSM II, bias and standard error are
slightly reduced. In the case of binary (OM IV) and linear selection (PSM I),
the MI-GLM and PMM B estimators have a coverage close to nominal, but variable
selection increases bias while decreasing standard error. For this
combination, the RMSE is larger for the estimators with variable selection.
For non-linear selection (PSM II), all estimators are significantly biased and
the coverage is close to zero, except for the PMM A estimator.

####  C.3 Non-parametric regression methods

In this section we check the performance of the PMM estimators with non-
parametric regression through LOESS via stats::loess function in R with the
following specification to boost performance:

  * •

loess.control(surface="interpolate") – surface approximation using a KD-tree ,

  * •

loess.control(trace.hat="interpolate") – approximated trace of the smoother
matrix,

  * •

span set to 0.2.

We adapt the same setting as in section C.1 and define additional target
variable $Y_{3}=x_{1}x_{2}\epsilon$. Table C.3 shows the results of this
simulation study. As expected, the proposed PMM estimators based on non-
parametric regression are characterised by negligible bias, but the variance
is significantly larger than that of the MI-GLM estimator for $Y_{1}$. As
expected, the bias and variance for MI-GLM are large compared to the non-
parametric PMM estimators for non-linear variables ($Y_{2}$ and $Y_{3}$). Note
that the PMM B ($\hat{y}-\hat{y}$ matching) estimator is almost unbiased even
when the target variables are highly non-linear. Owing to large standard
errors for the last variable, the coverage is close to 100%.

Table C.3: Results for the simulation study with non-parametric predictive mean matching models (all numbers multiplied by 100) Y | Estimator | Non-parametric? | Bias | SE | RMSE | CR
---|---|---|---|---|---|---
Y1 | MI-GLM | No | 0.31 | 2.93 | 2.95 | 97.60
| PMM A | No | 0.31 | 2.93 | 2.95 | 97.23
|  | Yes | -0.36 | 3.08 | 3.10 | 97.23
| PMM B | No | 0.46 | 5.18 | 5.20 | 99.82
|  | Yes | 0.32 | 5.19 | 5.20 | 99.63
Y2 | MI-GLM | No | -19.93 | 7.61 | 21.33 | 25.65
| PMM A | No | -19.93 | 7.61 | 21.33 | 25.28
|  | Yes | 2.41 | 8.64 | 8.97 | 94.10
| PMM B | No | -1.05 | 9.02 | 9.08 | 97.97
|  | Yes | -1.57 | 9.15 | 9.28 | 97.05
Y3 | MI-GLM | No | 59.41 | 7.61 | 59.90 | 100.00
| PMM A | No | 0.83 | 1.28 | 1.52 | 100.00
|  | Yes | -2.24 | 5.08 | 5.55 | 99.45
| PMM B | No | 1.22 | 9.61 | 9.69 | 99.82
|  | Yes | 1.51 | 9.93 | 10.04 | 99.63

####  C.4 Violation of the positivity assumption

In this section we check the performance of the proposed estimators when the
positivity assumption is violated. We consider the same settings as in section
C.1 with additional $X_{3}$ variable and effects connected with this variable.

  * •

$X_{3}\sim\text{NegBin}(\mu=10,r=4)$ and then numbers over 40 were grouped
into one category "40",

  * •

we created two target variables:

    * –

$Y_{1}=-7+6X_{1}-5X_{2}+X_{3}^{T}\gamma+15\epsilon$,

    * –

$Y_{2}=-2+0.37\cdot(X_{1}-0.5)^{2}+X_{2}^{2}+X_{3}^{T}\gamma+5\epsilon$,

where $\gamma\sim\mathcal{U}\left(-6,10\right)$ are coefficients associated
with each level of $X_{3}$ variable (40 possible levels),

  * •

we consider two scenarios for under-coverage:

    * –

Stochastic: where the samplable sub-population is generated from
$\text{Binomial}(\pi_{u})$ where
$\operatorname{logit}(\pi_{u})\allowbreak=X_{3}^{T}\lambda$ and
$\lambda=\begin{pmatrix}0.75&0.72&\ldots&0.10\end{pmatrix}^{T}$, so the larger
the $X_{3}$ the less likely to be observed in the sub-population.

    * –

Deterministic: where the samplable sub-population is generated in the
following process:

      * *

target variable $\pi_{u}$ is defined as
$\operatorname{logit}(\pi_{u})=0.84+0.32X_{1}+0.68X_{2}+X_{3}^{T}\lambda$
where $\lambda$ is defined above,

      * *

inclusion indicator was generated as
$\mathmybb{1}_{\\{\pi_{u}>Q_{\pi_{u},0.25}\\}}$, where $Q_{\pi_{u},0.25}$ is
the 25% percentile of $\pi_{u}$.

  * •

probabilities of inclusion in the non-probability sample from the under-
represented populations were defined as follows:

    * –

Stochastic: $\operatorname{logit}(p_{i})=X_{2}-X_{1}-2$, so the expected size
is about 17% of the population.

    * –

Deterministic: $\operatorname{logit}(p_{i})=0.6X_{1}-X_{2}-2$, so the expected
size is about 20% of the population.

Table C.4 shows the results of the simulation study when the positivity
assumption is violated. For the linear $Y_{1}$ all the estimators yield the
same results with slightly lower bias for the PMM B ($\hat{y}-\hat{y}$ fit).
However, the PMM B suffers from larger standard errors and, as a result, the
standard error and RMSE are larger for deterministic and stochastic results
and the CR is close to 100%. For the non-linear $Y_{2}$, the main difference
exists for deterministic under-coverage. The bias of PMM B is lower than that
of the MI-GLM and PMM A estimators and the standard error is comparable. Under
this scenario, coverage is the highest (about 67%), but still does not reach
the nominal level of 95%. In the case of stochastic under-coverage, the
proposed estimators are comparable to the MI-GLM estimator, with a slightly
smaller bias for the PMM B estimator.

Table C.4: Results of the simulation study when the positivity assumption is violated (results multiplied by 100) Variable | Under-coverage | Estimator | Bias | SE | RMSE | CR
---|---|---|---|---|---|---
Y1 | Deterministic | MI-GLM | 0.36 | 43.17 | 43.17 | 96.20
|  | PMM A | 0.36 | 43.17 | 43.17 | 96.20
|  | PMM B | 2.44 | 62.25 | 62.30 | 99.20
| Stochastic | MI-GLM | 1.96 | 43.36 | 43.40 | 96.40
|  | PMM A | 1.96 | 43.36 | 43.40 | 96.40
|  | PMM B | 2.92 | 60.49 | 60.56 | 98.80
Y2 | Deterministic | MI-GLM | -84.66 | 32.60 | 90.72 | 29.40
|  | PMM A | -84.66 | 32.60 | 90.72 | 29.40
|  | PMM B | -60.20 | 33.37 | 68.83 | 71.20
| Stochastic | MI-GLM | 11.21 | 30.93 | 32.90 | 94.20
|  | PMM A | 11.21 | 30.93 | 32.90 | 94.00
|  | PMM | 5.89 | 31.54 | 32.09 | 98.20

The deterministic case is supposed to mirror online sampling (or sampling via
telephone connection) where a unit can only be sampled if it has the internet
connection (a telephone number). Such cases usually depend on age, mirrored by
$X_{3}$, and, to a smaller extent, on social status, which may be correlated
with covariates. Note that despite the violation of the positivity assumption
in the deterministic case
$\mathbb{P}\left[\left.\delta=1\right|\boldsymbol{X},Y\right]=\mathbb{P}\left[\left.\delta=1\right|\boldsymbol{X}\right]$
a.s is maintained. In the case of stochastic under-coverage, the selection of
units to the non-probability sample can be viewed as a two-stage sampling
process; in this case, positivity is not violated, which explains good results
of the MI-GLM estimator.

####  C.5 Multiply robust imputation

In this section we test whether the approach proposed by [24] is suitable for
non-probability samples. They proposed multiply robust imputation using
predictive mean matching, where predictions from multiple models are used as
independent variables. In this section we follow the same design as that
described in section C.1 but we create two target variables as follows:

  * •

$Y_{1}=1+\tfrac{1}{5}x_{1}+5x_{2}+\epsilon$,

  * •

$Y_{2}=-2+5\left(x_{1}-0.5\right)^{5}+x_{2}^{3}+\epsilon$.

For both models we consider two regression models:

  * •

Model 1: linear regression with $x_{1}$ and $x_{2}$,

  * •

Model 2: linear regression with $(x_{1}-0.5)^{5}$ and $x_{2}^{3}$,

so that for each model, one of the models considered is correctly specified.
For both PMM estimators without the multirobust property, we used a linear
model with $x_{1}$ and $x_{2}$.

Table C.5: Results for the simulation study with multiply robust predictive mean matching Variable | Estimator | MultiRobust | Bias | SE | RMSE | CR
---|---|---|---|---|---|---
$Y_{1}$ | MI-GLM | – | 0.00 | 0.23 | 0.23 | 93.60
| PMM A | No | 0.00 | 0.23 | 0.23 | 93.60
|  | Yes | 0.00 | 0.23 | 0.23 | 93.60
| PMM B | No | 0.00 | 0.23 | 0.23 | 94.00
|  | Yes | 0.00 | 0.23 | 0.23 | 94.40
$Y_{2}$ | MI-GLM | – | -17.40 | 6.37 | 18.53 | 22.60
| PMM A | No | -17.30 | 6.36 | 18.43 | 23.60
|  | Yes | 0.40 | 9.49 | 9.50 | 92.20
| PMM B | No | -1.33 | 8.77 | 8.87 | 90.80
|  | Yes | 0.40 | 9.49 | 9.50 | 92.20

Table C.5 contains the results for multiply robust predictive mean matching.
All linear estimators are characterised by the same bias, SE, RMSE and CR.
Differences can be observed for the non-linear $Y_{2}$. MI-GLM and the
standard PMM A are characterised with the same bias and RMSE, while the
addition of multiplicative robustness to the PMM A estimator significantly
reduces its bias. For the PMM B estimator, it does not matter whether the
fitting is done with a linear model with $x_{1}$ and $x_{2}$ or with
predictions from the two models considered. In this limited simulation study
none of the models for $Y_{2}$ reached the nominal coverage rate of 95%.

###  Appendix D Implementation in the R language

The proposed methods are implemented in the R language in the nonprobsvy
package. One should install the package from the github repository using the
following code:

    install.packages("remotes")
    remotes::install_github("ncn-foreigners/nonprobsvy")

Below we present functions for the PMM estimator under specific settings. In
general, the user should change the following arguments:

  * •

method_outcome = "pmm" – specification that the PMM estimator is used

  * •

controlOut – controls for MI estimator

    * –

k – number of nearest neighbours,

    * –

predictive_match – whether $\hat{y}-y$ (=1; PMM A) or $\hat{y}-\hat{y}$ (=2;
PMM B) matching should be used.

  * •

controlInf – controls for inference (variance estimation)

    * –

pmm_exact_se – whether $\hat{V}_{2}$ should be estimated using "mini"
bootstrap (=TRUE) or not (=FALSE).

Example codes are shown below.

  * •

Predictive mean matching with $\hat{y}-y$ matching (PMM A) with $\hat{V}_{2}$
estimated using "mini"-bootstrap presented in Section 3.

        PMM_A <- nonprob(
      outcome = y1 ~ x1 + x2,
      data = sample_nonprob,
      svydesign = sample_prob,
      method_outcome = "pmm",
      family_outcome = "gaussian",
      control_outcome = controlOut(k = 5,
                                   predictive_match = 1),
      control_inference = controlInf(pmm_exact_se = TRUE)
    )

  * •

Predictive mean matching with $\hat{y}-\hat{y}$ matching (PMM B) with
$\hat{V}_{2}$ estimated using "mini"-bootstrap presented in Section 3.

        PMM_B <- nonprob(
      outcome = y1 ~ x1 + x2,
      data = sample_nonprob,
      svydesign = sample_prob,
      method_outcome = "pmm",
      family_outcome = "gaussian",
      control_outcome = controlOut(k = 5,
                                   predictive_match = 2),
      control_inference = controlInf(pmm_exact_se = TRUE)
    )

  * •

Predictive mean matching with $\hat{y}-\hat{y}$ matching with local linear
regression (PMM B) with $\hat{V}_{2}$ estimated using "mini"-bootstrap
presented in Section 3.

        PMM_B <- nonprob(
      outcome = y1 ~ x1 + x2,
      data = sample_nonprob,
      svydesign = sample_prob,
      method_outcome = "pmm",
      family_outcome = "gaussian",
      control_outcome = controlOut(k = 5,
                                   predictive_match = 2,
                                   pmm_reg_engine = "loess"),
      control_inference = controlInf(pmm_exact_se = TRUE)
    )