--- title: "nonprobsvy: An R Package for Modern Methods for Non-Probability Surveys" author: - Łukasz Chrostowski (Analytic Partners) - Piotr Chlebicki (Stockholm University) - Maciej Beręsewicz (Poznań University of Economics and Business; Statistical Office in Poznań) date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 bibliography: nonprobsvy.bib link-citations: true vignette: > %\VignetteIndexEntry{nonprobsvy: An R Package for Modern Methods for Non-Probability Surveys} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} Sys.setenv(OMP_NUM_THREADS = "1", OPENBLAS_NUM_THREADS = "1", MKL_NUM_THREADS = "1", VECLIB_MAXIMUM_THREADS = "1", OMP_THREAD_LIMIT = "1") options(prompt = "R> ", continue = "+ ", width = 76) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, fig.align = "center" ) ``` > **Abstract.** This vignette is based on the paper Chrostowski, Ł., > Chlebicki, P., & Beręsewicz, M. (2026). nonprobsvy: An R package for modern > methods for non-probability surveys. Journal of Statistical Software, 117(2), > 1-37. https://doi.org/10.18637/jss.v117.i02. > The paper presents **nonprobsvy** -- an R package for inference > based on non-probability samples. The package implements various approaches > that can be categorized into three groups: prediction-based approach, inverse > probability weighting and doubly robust approach. In the package, we assume > the existence of either population-level data or probability-based population > information and leverage the **survey** package for inference. The package > implements both analytical and bootstrap variance estimation for the proposed > estimators. In the paper we present the theory behind the package, its > functionalities and a case study that showcases the usage of the package. The > package is aimed at scientists and researchers who would like to use > non-probability samples (e.g., big data, opt-in web panels, social media) to > accurately estimate population characteristics. > > *Keywords*: data integration, doubly robust estimation, propensity score > estimation, mass imputation, **survey**. # Introduction In official statistics, information about the target population and its characteristics is mainly collected through probability surveys, censuses or is obtained from administrative registers and covers all (or nearly all) units of the population. However, owing to increasing non-response rates, particularly unit non-response and non-contact, which result from the growing respondent burden as well as rising costs of surveys conducted by National Statistical Institutes, non-probability data sources are becoming more popular [@berkesewicz2017two; @beaumont2020probability; @biffignandi2021handbook]. Non-probability surveys, such as opt-in web panels, social media, scanner data, mobile phone data or voluntary register data, are currently being explored for use in the production of official statistics [@citro2014multiple; @daas2015big], public opinion studies [@Schonlau2017] or market research [cf. @Grow2022]. Since the selection mechanism underlying these sources is unknown, standard design-based inference methods cannot be directly applied and, in the case of large datasets, can lead to the big data paradox, i.e., the larger the sample, the larger the bias, as described by @meng2018statistical. Table 1 compares basic characteristics of probability and non-probability samples. In particular, it shows the advantages and disadvantages of each type with respect to the selection mechanism, the population coverage, bias, variance, costs and timeliness. In general, the quality of non-probability samples suffers from an unknown selection mechanism (i.e., unknown probabilities of inclusion) and under-coverage of certain groups from the population (e.g.,  older people). As a result, direct estimates based on non-probability samples are biased and, in most cases, are characterized by small variance owing to their size. Certainly, the costs and timeliness of these surveys are significantly smaller than those of probability samples. Table: **Table 1.** A comparison of probability and non-probability samples and their characteristics. | *Factor* | *Probability sample* | *Non-probability sample* | |:----------------------|:------------------------|:----------------------------------------| | Selection | Known probabilities | Unknown probabilities (self-selection) | | Coverage | Complete | May be incomplete | | Estimation bias | Unbiased under design | Potential systematic bias | | Variance of estimates | Typically high | Typically low | | Cost | High | Low | | Timeliness | Long delay | Very short delay | To address this problem, several approaches have been proposed, which rely on the estimation of propensity scores (PS; i.e., inclusion probabilities) for deriving inverse probability weights (IPW; also known as propensity score weighting/adjustment, cf. @lee2006propensity [@lee2009estimation]), on model-based prediction (in particular, mass imputation estimators; MI) and on the doubly robust (DR) approach involving IPW and MI estimators. Two main scenarios are usually considered: 1. Only population-level means or totals are available (e.g., from census, registers or sample surveys). 2. Unit-level data are available either in the form of registers covering the whole population or in the form of probability surveys [cf. @elliott_inference_2017]. @wu2022statistical classified these approaches into three groups that require a joint randomization framework involving a probability sampling design (denoted as $p$) and an outcome regression model (denoted as $\xi$) or a PS model (denoted as $q$). According to this classification, IPW estimators represent the $qp$ framework, MI estimators represent the $\xi p$ framework, and DR estimators can represent either the $qp$ or the $\xi p$ framework. Most approaches assume that population data is available in order to reduce the bias of non-probability sampling by reweighting to reproduce known population totals/means (i.e., IPW estimators); by modeling the target variable using various techniques (i.e., MI estimators); or by combining both approaches (e.g., DR estimators, cf. @chen2020doubly; see also multilevel regression and post-stratification, MRP; Mister-P, cf. @gelman1997poststratification). This topic has become very popular and a number of new methods have been proposed; for instance non-parametric approaches based on nearest neighbors [@yang2021integration], kernel density estimation [@chen_nonparametric_2022], empirical likelihood [@kim2023empirical], model-calibration with LASSO [@chen2018] or quantile balanced IPW [@beresewicz2025] to name a few. It should be highlighted that, in contrast to probability samples, there is no single method that can be used for all non-probability samples. Based on the methods available in the literature several statistical software solutions have been developed. These pieces of software will be presented and compared with **nonprobsvy** [@nonprobsvy] in the next section. ## Software for non-probability samples Table 2 presents a comparison of selected packages in terms of the availability of different inference methods. The reason for this is the availability of software designed for inference based on non-probability samples, not a general way to estimate propensity scores or model $\mathbb{E}(Y \mid X)$. For such reviews, one can refer to @valliant2018survey [Chapter 6] who discussed weighting and matching in this context in Stata [@Stata2025] using `logit`, `svy` and `teffects nnmatch` commands or @lohr2021sas [e.g., Chapter 8] who covers SAS software [@sas94; e.g., `PROC LOGISTIC` procedure]. We focus on packages available as open source or via CRAN or PyPI (for non-CRAN or non-PyPI software see @cobo2024software). Our comparison includes five packages that focus specifically on non-probability sampling: in R -- **NonProbEst** [@NonProbEst] and our **nonprobsvy**, and in Python -- **balance** [@sarig2023balancepythonpackage], **inps** [@castro2024inps], and in Stan [@carpenter2017stan] e.g., through the R package **rstanarm**, which allows you to easily apply the **rstanarm**'s MRP method [@rstanarm]. Table: **Table 2.** A comparison of packages and implemented methods. | | Stan | balance (Python) | inps (Python) | NonProbEst (R) | nonprobsvy (R) | |:--------------------------------------|:----:|:----------------:|:-------------:|:--------------:|:--------------:| | IPW | -- | ✓ | ✓ | ✓ | ✓ | | Calibrated IPW | -- | -- | -- | -- | ✓ | | MI | -- | -- | -- | ✓ | ✓ | | DR | -- | -- | ✓ | -- | ✓ | | MRP | ✓ | -- | -- | -- | -- | | Variable selection | ✓ | ✓ | ✓ | ✓ | ✓ | | Analytical variance | -- | -- | -- | -- | ✓ | | Bootstrap variance | -- | -- | -- | ✓ | ✓ | | Integration with **survey**/**samplics** | -- | -- | -- | -- | ✓ | The **NonProbEst** package, while more limited than **nonprobsvy**, offers various techniques, such as PS or prediction approaches (e.g., model-calibrated). Users can choose several different settings for PS weights, variable selection and can estimate variance of the mean using the leave-one-out Jackknife procedure. Unfortunately, has had no new CRAN release since version 0.2.4. While the package contains functions designed for specific methods, it does not allow users to leverage the **survey** package [@survey-pkg] for estimation. The **balance** package is solely dedicated to the PS approach (as of version 0.20.0). It assumes that the reference probability sample is available and the authors have implemented the variance estimator of the weighted mean as a measure of uncertainty for the IPW estimator. The weights for the IPW estimator are constructed using the approach proposed by @Schonlau2017. The **inps** package supports the use of unit-level data from a probability sample or the population, implements IPW, MI and DR estimators, and offers users the possibility of selecting variables but is still at a very early stage of development. It also implements kernel weighting and a simple bootstrap approach via the **scipy.stats** package [@scipy2020]. Neither **balance** nor **inps** supports the use of the **samplics** package [@Diallo2021]. Finally, we note that the MRP approach is implemented solely in Stan with variable selection specified by an appropriate prior, but it can be easily implemented in other probabilistic programming languages, e.g., using the **Turing.jl** package [@turing] in the Julia language [@julia]. The **nonprobsvy** package has several advantages over those presented above. Firstly, it implements state-of-the-art methods recently proposed in the literature, along with valid statistical inference procedures. Secondly, it offers other approaches, such as calibrated IPW (where PS weights match either the population or estimated totals), NN and PMM matching, various IPW and DR estimators using generalized linear models (GLM) with different link functions for probability estimation. Thirdly, it supports the functions included in the **survey** package to account for the design of the probability sample. Finally, we provide a user-friendly API that mimics `glm()`, `survey::calibrate()` and other functions known in R, together with the main function to specify the approach and estimators. Among the packages reviewed here, is the only one that offers such functionalities. The remaining part of the paper is structured as follows. The next section is dedicated to the theory of statistical inference based on non-probability samples. We provide the basic set-up and introduce specific methods in separate sections. We follow the notation used by @wu2022statistical throughout the paper. We then describe the main function and the package functionalities, and present an empirical study showcasing the process of integrating data from the Polish Job Vacancy Survey with voluntary administrative data from the Central Job Offers Database in order to estimate the share of companies with at least one vacancy offered on a single shift. We then present classes and S3 methods implemented in the package. The paper ends with a summary and plans for future development. The Appendix presents algorithms for selected MI estimators. The package **nonprobsvy** [@nonprobsvy] is available from the Comprehensive R Archive Network (CRAN) at . # Methods for non-probability samples ## Basic setup Let $U = \{1, \ldots, N\}$ denote the target population consisting of $N$ labeled units. Each unit $i$ has an associated $L$-dimensional vector of auxiliary variables $\boldsymbol{x}_{i}$ and the study (target) variable $y_{i}$. Let $\{ (y_i, \boldsymbol{x}_i), i \in \mathcal{I}_{\text{NP}}\}$ be a non-probability sample $S_{\text{NP}}$ of size $n_{\text{NP}}$, where $\mathcal{I}_{\text{NP}}$ is an index set for $S_{\text{NP}}$. Let $\left\{\left(\boldsymbol{x}_i, \pi_{\text{P},i}\right), i \in \mathcal{I}_{\text{P}}\right\}$ be a probability sample $S_{\text{P}}$ of size $n_{\text{P}}$, where $\mathcal{I}_{\text{P}}$ is an index set for $S_{\text{P}}$, and the only information known for all units in the population refers to auxiliary variables $\boldsymbol{x}$ and inclusion probabilities $\pi_{\text{P}}$. Each unit in the $S_{\text{P}}$ sample has been assigned a design-based weight given by $d_{\text{P},i} = 1/\pi_{\text{P},i}$. Let $I_{\text{NP},i}=I(i \in S_{\text{NP}})$ and $I_{\text{P},i}=I(i \in S_{\text{P}})$ be indicators of inclusion in the non-probability sample $S_{\text{NP}}$ and the probability sample $S_{\text{P}}$, respectively, which are defined for all units in the target population. Let $\pi_{\text{NP},i}=\mathbb{P}(I_{\text{NP},i}=1 \mid \boldsymbol{x}_i)$ be propensity scores, which characterize the $S_{\text{NP}}$ sample's inclusion and participation mechanisms. Unlike $\pi_{\text{P},i}$, the $\pi_{\text{NP},i}$ and $d_{\text{NP},i}=1/\pi_{\text{NP},i}$ are unknown. The description of the data is presented in a more concise form in Table 3. Table: **Table 3.** Two-sample setting with no overlap between samples. | Sample | ID | Inclusion ($I_{\text{NP}}$) | Design weight ($d$) | Covariates ($\boldsymbol{x}$) | Study variable ($y$) | |:----------------------|:------------------|:---------------------------:|:-------------------:|:-----------------------------:|:--------------------:| | Non-probability | 1 | 1 | ? | ✓ | ✓ | | $S_{\text{NP}}$ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | | | $n_{\text{NP}}$ | 1 | ? | ✓ | ✓ | | Probability | $n_{\text{NP}}+1$ | 0 | ✓ | ✓ | ? | | $S_{\text{P}}$ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | | | $n_{\text{NP}}+n_{\text{P}}$ | 0 | ✓ | ✓ | ? | The goal is to estimate the finite population mean $\mu_{y} = N^{-1}\sum_{i=1}^{N} y_{i}$ of the target variable $y$. As values of $y_{i}$ are not observed in the probability sample, they cannot be used to estimate the target quantity. Instead, one could try combining the non-probability and probability samples (or known population totals) to estimate $\mu_{y}$. Given the absence of a universally accepted method for achieving this objective, assumptions vary considerably, as outlined by @wu2022statistical. However, the following are the main assumptions that apply to the majority of methods presented in this section: * **A1** $I_{\text{NP}, i}$ and the study variable $y_{i}$ are independent given the set of covariates $\boldsymbol{x}_{i}$ (the missing at random mechanism). * **A2** All the units in the target population have non-zero PS, i.e., $\pi_{\text{NP}, i}>0$, $i=1, 2, \ldots, N$ (i.e., no coverage error). * **A3** The indicator variables $I_{\text{NP}, 1}, I_{\text{NP}, 2}, \ldots, I_{\text{NP}, N}$ are independent given the set of auxiliary variables $\left(\boldsymbol{x}_1, \boldsymbol{x}_2, \ldots, \boldsymbol{x}_N\right)$ (i.e., no clustering). Currently, we ignore overlap between $S_{\text{NP}}$ and $S_{\text{P}}$, and assume no measurement error in $y_i$ or $\boldsymbol{x}_i$ and the fact that values of $\boldsymbol{x}_i$ are known. The setting presented in Table 3 can also be extended to calibrated $d_{\text{P}, i}$ weights (i.e., $d_{\text{P}, i}$ adjusted for under-coverage, non-contact or non-response; cf. @sarndal2005estimation) but this requires additional developments in the theory about the consistency of the MI, IPW and DR estimators. In the next sections we briefly present the methods implemented in the package. ## Prediction-based approach ### Prediction estimators In the model-based prediction approach the following semi-parametric model for the finite population is assumed: $$ \mathbb{E}_{\xi}\left(y_i \mid \boldsymbol{x}_i\right)=m\left(\boldsymbol{x}_i, \boldsymbol{\beta}\right), \text { and } \quad V_{\xi}\left(y_i \mid \boldsymbol{x}_i\right)=v\left(\boldsymbol{x}_i\right) \sigma^2, \quad i=1, 2, \ldots, N, \tag{1} $$ where the mean function $m(\cdot, \cdot)$ and the variance $v(\cdot)$ have known forms, $\boldsymbol{\beta}$ and $\sigma^2$ are model parameters and $y_i$ are also assumed to be conditionally independent given the $\boldsymbol{x}_i$. The model in Equation 1 is assumed to hold for all units in the non-probability sample $S_{\text{NP}}$. The parameters of the model in Equation 1 can be estimated using the quasi maximum likelihood estimation method, which includes linear and non-linear models, such as GLM. Let $\boldsymbol{\beta}_0$ and $\sigma^2_0$ be the true values of the model parameters $\boldsymbol{\beta}$ and $\sigma^2$ under the adopted model and $\hat{\boldsymbol{\beta}}$ be the quasi maximum likelihood estimator of $\boldsymbol{\beta}_0$. Let $m_i=m(\boldsymbol{x}_i, \boldsymbol{\beta}_0)$ and $\hat{m}_i=m(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}})$ be calculated for all units $i=1, \ldots, N$. Under this setting, as @wu2022statistical notes, there are two commonly used prediction estimators: $$ \hat{\mu}_{y, \text{PR1}}=\frac{1}{N} \sum_{i=1}^N \hat{m}_i \quad \text { and } \quad \hat{\mu}_{y, \text{PR2}}=\frac{1}{N}\left\{\sum_{i \in S_{\text{NP}}} y_i-\sum_{i \in S_{\text{NP}}} \hat{m}_i+\sum_{i=1}^N \hat{m}_i\right\}. \tag{2} $$ Under linear models, where $m(\boldsymbol{x}_i, \boldsymbol{\beta})=\boldsymbol{x}_i^{\top}\boldsymbol{\beta}$, the two estimators in Equation 2 reduce to: $$ \hat{\mu}_{y, \text{PR1}}=\mu_{\boldsymbol{x}}^{\top} \hat{\boldsymbol{\beta}} \quad \text { and } \quad \hat{\mu}_{y, \text{PR2}}=\frac{n_{\text{NP}}}{N}\left(\overline{y}_{\text{NP}}-\overline{\boldsymbol{x}}_{\text{NP}}^{\top} \hat{\boldsymbol{\beta}}\right)+\mu_{\boldsymbol{x}}^{\top} \hat{\boldsymbol{\beta}}, \tag{3} $$ where $\mu_{\boldsymbol{x}} = N^{-1}\sum_{i=1}^N\boldsymbol{x}_i$ is the vector of the population means of the $\boldsymbol{x}$ variables and $\overline{\boldsymbol{x}}_{\text{NP}}=n_{\text{NP}}^{-1}\sum_{i \in S_{\text{NP}}}\boldsymbol{x}_i$ is the vector of the simple means of $\boldsymbol{x}$ from the non-probability $S_{\text{NP}}$ sample. If the linear model contains an intercept and $\hat{\boldsymbol{\beta}}$ is the ordinary least squares estimator, then $\hat{\mu}_{y, \text{PR1}}=\hat{\mu}_{y, \text{PR2}}$. The estimator in Equation 3 is appealing as it only requires the non-probability sample $S_{\text{NP}}$ and reference population means (or totals and population size $N$). If the population means are unknown, they can be replaced by estimates provided by the reference probability sample $S_{\text{P}}$, i.e., $\sum_{i=1}^N \hat{m}_i$ is replaced with $\sum_{i \in S_{\text{P}}} d_{\text{P}, i}\hat{m}_i$ for Equation 2 and $\mu_{\boldsymbol{x}}$ is replaced by $\hat{\mu}_{\boldsymbol{x}}=\hat{N}_{\text{P}}^{-1}\sum_{i \in S_{\text{P}}}d_{\text{P}, i}\boldsymbol{x}_i$ for Equation 3 where $\hat{N}_{\text{P}}=\sum_{i \in S_{\text{P}}}d_{\text{P}, i}$. ### Mass imputation estimators Model-based prediction estimators of $\mu$ can be treated as mass imputation estimators, since the information on $y_i$ is missing entirely in the reference probability sample $S_{\text{P}}$ (but $\boldsymbol{x}_i$ is available) and $y_i$ can be imputed based on the non-probability sample as values of $\{ (y_i, \boldsymbol{x}_i), i \in \mathcal{I}_{\text{NP}}\}$ are known. The general form of the MI estimator is given by: $$ \hat{\mu}_{y, \text{MI}}=\frac{1}{\hat{N}_{\text{P}}} \sum_{i \in S_{\text{P}}} d_{\text{P}, i} y_i^*, $$ where $y_i^*$ is the imputed value of $y_i$ and $\hat{N}_{\text{P}}$ is defined as previously. Under deterministic regression imputation, the $\hat{\mu}_{y, \text{MI}}$ estimator reduces to the estimators in Equation 2. There are several ways of imputing $y_i^*$ and in the **nonprobsvy** package we have implemented the following MI estimators: the semi-parametric approach based on GLM (MI-GLM), kernel smoothing (MI-NPAR), nearest neighbor matching (MI-NN) and predictive mean matching (MI-PMM). The properties of the MI-GLM estimator, where $y_i^*$ are imputed with $\hat{m}_i$ from the semi-parametric model, were studied by @kim_combining_2021. In the **nonprobsvy** package, we account for the following GLM families: `gaussian()`, `binomial()` and `poisson()`. The MI-NPAR estimator was proposed by @chen_nonparametric_2022 who studied its finite population properties. Currently, we implement this approach using local polynomial regression using `stats::loess()` function. In the next releases we will allow users to use the **np** [@nppkg] or **KernSmooth** [@KernSmooth] packages. The MI-NN estimator was initially proposed by @rivers2007sampling under the name sample matching and theoretical properties of the MI-NN estimator for large non-probability samples (big data, i.e., covering a significant part of the target population) were studied by @yang2021integration. The basic idea of NN matching is as follows: 1. for each unit $i$ in the probability sample $S_{\text{P}}$ find a donor $j$ (or multiple donors) in the $S_{\text{NP}}$ sample based on some distance between $\boldsymbol{x}_i$ and $\boldsymbol{x}_j$, 2. use the matched values $y_j$ from $S_{\text{NP}}$ to impute missing $y_i$ values in the probability sample $S_{\text{P}}$. Imputed values of $y_i^*$ depend on the number of selected $k$ neighbors: for $k=1$, the closest unit is selected, and for $k>1$, one can calculate a simple average over a vector of selected $y$ values. A detailed description of the procedure is presented in Algorithm 1 in the Appendix. The MI-NN estimator suffers from the curse of dimensionality, i.e., asymptotic bias of the MI estimator increases as the number of covariates $\boldsymbol{x}$ increases with a fixed $k$ [@abadie2006large; @yang_asymptotic_2020]. The single PMM imputation approach has been proposed to overcome this issue [@chlebicki2025]. In the PMM approach, matching is done using predicted values of $\hat{m}_i=m\left(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}}\right)$ instead of $\boldsymbol{x}_i$, thus the NN algorithm is modified as follows: 1. fit the model $m\left(\boldsymbol{x}_i, \boldsymbol{\beta}\right)$ to non-probability $S_{\text{NP}}$ sample, 2. assign predicted values $\hat{m}_i$ to all units in $S_{\text{NP}}$ and $S_{\text{P}}$, 3. match all units from the $S_{\text{P}}$ sample to donor units from the $S_{\text{NP}}$ sample based on $\hat{m}$ values. The MI-PMM estimator is the same as in the NN approach but differs from the multiple imputation PMM implemented in the **mice** package [@mice]. In **nonprobsvy** we have implemented a single, deterministic imputation method (i.e., one dataset is created). @chlebicki2025 studied properties of two variants of the MI-PMM estimator for non-probability samples: matching predicted to predicted ($\hat{m}-\hat{m}$ matching; denoted as MI-PMM-A) and matching predicted to observed ($\hat{m}-y$ matching; denoted as MI-PMM-B). Details of the procedure can be found in Algorithm 2 and Algorithm 3 in the Appendix. @chlebicki2025 also prove the consistency of the MI-PMM-A estimator under model misspecification, i.e., the assumed model may differ from the true one. ### Variance estimators for the prediction approach Variance of the MI estimators can be estimated analytically or using the bootstrap approach. The analytical estimator of the variance of the MI-GLM estimator proposed by @kim_combining_2021 [p. 950] contains two components: $\hat{V}_1$ (based on the information from both samples $S_{\text{NP}}$ and $S_{\text{P}}$) and $\hat{V}_2$ (based exclusively on the probability sample $S_{\text{P}}$). For the MI-NN estimator @yang2021integration proposed a variance estimator for large $S_{\text{NP}}$ samples, which reduces to the part for the probability sample $S_{\text{P}}$ (i.e., the design-based variance estimator of the mean, which can easily be obtained from the **survey** package) and a version for smaller samples can be found in @chlebicki2025. With respect to the variance of MI-PMM estimators, @chlebicki2025 propose formulas which are the same as those for the MI-NN estimators. The analytical variance estimator for the MI-NPAR was proposed by @chen_nonparametric_2022. In the bootstrap approach each bootstrap replication $b=1, \ldots, B$ consists of the following steps. 1. Independently: * Draw a simple random sampling with replacement from the non-probability sample $S_{\text{NP}}$ (using `base::sample()` function). * Draw a sample according to the declared sampling design from the probability sample $S_{\text{P}}$ (e.g., one can use the `as.svrepdesign()` function from the **survey** package). 2. Estimate $\mu_{y, \text{MI}}^b$ using an appropriate approach (e.g., MI-GLM, MI-NN or MI-PMM). After obtaining $B$ bootstrap replicates, estimate variance using the following equation: $$ \hat{V}_{\text{boot}} = \frac{1}{B-1}\sum_{b=1}^B\left(\hat{\mu}^b_{y, \text{MI}} - \hat{\mu}_{y, \text{MI}}\right)^2. \tag{4} $$ The above approaches are applied when unit-level data from the probability sample $S_{\text{P}}$ are available. If only population means or totals and the population size are available, the current implementation supports MI-GLM, using the first component $\hat V_1$ of the @kim_combining_2021 variance estimator, with survey-based population quantities substituted when available. MI-NPAR, MI-NN, and MI-PMM require unit-level auxiliary data. To estimate the variance of the MI-NN and MI-PMM estimators we only allow the bootstrap approach with known population means. For the probability sample, bootstrap replicates can be generated internally through from the package. The current version does not, however, accept user-supplied replicate-weight designs. ## Inverse probability weighting Inverse probability weighting, another popular estimation approach, involves estimating PS given by $\pi_{\text{NP}, i}=\mathbb{P}\left(i \in S_{\text{NP}} \right)$. As in the case of the prediction-based approach, there are two variants of the IPW estimator, given by: $$ \hat{\mu}_{y, \text{IPW-HT}}=\frac{1}{N} \sum_{i \in S_{\text{NP}}} \frac{y_i}{\hat{\pi}_{\text{NP}, i}} \quad \text{and} \quad \hat{\mu}_{y, \text{IPW-Hájek}}=\frac{1}{\hat{N}_{\text{NP}}} \sum_{i \in S_{\text{NP}}} \frac{y_i}{\hat{\pi}_{\text{NP}, i}}, \tag{5} $$ where the $\hat{\mu}_{y, \text{IPW-HT}}$ is an adjustment of the Horvitz-Thompson estimator, and the $\hat{\mu}_{y, \text{IPW-Hájek}}$ is an adjustment of the Hájek estimator, where the estimated population size is given by $\hat{N}_{\text{NP}} = \sum_{i \in S_{\text{NP}}} \hat{\pi}_{\text{NP}, i}^{-1}$. The use of this estimator with respect to non-probability samples is discussed by @lee2006propensity and @biffignandi2021handbook [Chapter 13] and there are several approaches to using propensity scores along with alternative versions of the weights [cf. @elliott_inference_2017, Section 3]. In an article by @chen2020doubly dedicated to the properties of the estimators in Equation 5, the authors proved their consistency and derived their closed form versions. @wu2022statistical [Section 4.2] argues that the $\hat{\mu}_{y, \text{IPW-Hájek}}$ estimator performs better than $\hat{\mu}_{y, \text{IPW-HT}}$ even if the population size is known. The construction of the IPW estimator involves two steps: 1. estimating the PS, 2. deriving $d_{\text{NP}, i}$ weights (which, in our case, are equal to $\pi_{\text{NP}, i}^{-1}$) To estimate the propensity scores $\pi_{\text{NP}, i}=\pi(\boldsymbol{x}_i, \boldsymbol{\gamma})$ one can use the likelihood approach assuming that the information about $\boldsymbol{x}_i$ is available for each unit in the population given by $$ \begin{aligned} \ell(\boldsymbol{\gamma}) &= \log\left\{\prod_{i=1}^N \left(\pi_{\text{NP}, i}\right)^{I_{\text{NP}, i}}\left(1-\pi_{\text{NP}, i}\right)^{1-I_{\text{NP}, i}}\right\} \\ &= \sum_{i \in S_{\text{NP}}} \log \left\{\frac{\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)}{1-\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)}\right\}+\sum_{i=1}^N \log \left\{1-\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)\right\}. \end{aligned} \tag{6} $$ In practice, a function of this form cannot be used because not all units from the population are observed. A more realistic approach consists in using a reference probability sample $S_{\text{P}}$, which means that the second component of the equation in Equation 6 is replaced, yielding a pseudo log-likelihood function given by $$ \ell^*(\boldsymbol{\gamma}) = \sum_{i \in S_{\text{NP}}} \log \left\{\frac{\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)}{1-\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)}\right\}+ \sum_{i \in S_{\text{P}}} d_{\text{P}, i} \log \left\{1-\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)\right\}. $$ The maximum pseudo-likelihood estimator $\hat{\boldsymbol{\gamma}}$ can be obtained as the solution to the pseudo score equation, which, under the logit function assumed for $\pi_{\text{NP}, i}$, is given by $$ \boldsymbol{U}(\boldsymbol{\gamma}) = \sum_{i \in S_{\text{NP}}} \boldsymbol{x}_i - \sum_{i \in S_{\text{P}}} d_{\text{P}, i} \pi(\boldsymbol{x}_i, \boldsymbol{\gamma}) \boldsymbol{x}_i. \tag{7} $$ In general, pseudo score functions $\boldsymbol{U}(\boldsymbol{\gamma})$ for true values of model parameters $\boldsymbol{\gamma}_0$ are unbiased under the joint $q p$ randomization in the sense that $\mathbb{E}_{q p}\left\{\boldsymbol{U}\left(\boldsymbol{\gamma}_0\right)\right\}=\boldsymbol{0}$, which implies that the estimator $\hat{\boldsymbol{\gamma}}$ is consistent for $\boldsymbol{\gamma}_0$ [@wu2022statistical]. The terms in Equation 7 can be replaced by general estimation equations. Let $\boldsymbol{h}(\boldsymbol{x}, \boldsymbol{\gamma})$ be a user-specified vector of functions with the same dimension of $\boldsymbol{\gamma}$ and $\boldsymbol{G}(\boldsymbol{\gamma})$ is defined as $$ \boldsymbol{G}(\boldsymbol{\gamma})=\sum_{i \in S_{\text{NP}}} \boldsymbol{h}\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)-\sum_{i \in S_{\text{P}}} d_{\text{P}, i} \pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right) \boldsymbol{h}\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right), $$ then solving for $\boldsymbol{G}(\boldsymbol{\gamma})=\boldsymbol{0}$ with the chosen parametric form of $\pi_{\text{NP}, i}$ and the chosen $\boldsymbol{h}(\boldsymbol{x}, \boldsymbol{\gamma})$ produces the consistent estimator of $\hat{\boldsymbol{\gamma}}$. In the literature, the most commonly considered functions are $\boldsymbol{h}\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right) = \boldsymbol{x}_i$ and $\boldsymbol{h}\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right) = \boldsymbol{x}_i \pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)^{-1}$. Note that if the function $\boldsymbol{h}\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right)=\boldsymbol{x}_i$, then $\boldsymbol{G}$ reduces to $\boldsymbol{U}$ and for the second variant of the $\boldsymbol{h}$ function we get $$ \boldsymbol{G}(\boldsymbol{\gamma}) = \sum_{i \in S_{\text{NP}}} \frac{\boldsymbol{x}_i}{\pi\left(\boldsymbol{x}_i, \boldsymbol{\gamma}\right) }-\sum_{i \in S_{\text{P}}} d_{\text{P}, i} \boldsymbol{x}_{i}, \tag{8} $$ which can be viewed as calibrated IPW. Equation 8 only requires the knowledge of population totals for auxiliary variables $\boldsymbol{x}$. Moreover, the use of Equation 8 yields a doubly robust estimator that is consistent if either the propensity model is correct or the outcome is linear in the calibration variables [cf. @kim_theory_2012]. ### Variance estimators for the inverse probability weighting approach @chen2020doubly [Section 3.2] derived asymptotic variance estimators for both IPW estimators in Equation 5 and presented the plug-in variance estimator for the $\hat{\mu}_{y, \text{IPW-Hájek}}$ estimator assuming logistic regression. In the package we have implemented this approach for `"logit"`, `"probit"` and `"cloglog"` link functions of the `binomial()` family. We refer the reader to @wu2022statistical [Section 6.2] and @chrostowski2024statistical [Chapter 3] for more details on how these estimators are derived based on the general estimating equations approach. Another approach is to use a bootstrap, which is essentially the same as the one presented in Equation 4, where $\hat{\mu}_{y}$ is replaced by one of the estimators in Equation 5. ## Doubly robust approach The IPW and MI estimators are suited to correctly specified models for PS and outcome regression models, respectively. The DR approach was proposed to improve robustness and efficiency [cf. @robins1994estimation]. It incorporates a prediction model for the response variable $y_i$ and a PS model for participation $I_{\text{NP}, i}$. This approach is doubly robust in the sense that the DR estimator remains consistent even if one of the models is misspecified. We need to consider a joint randomization approach involving a non-probability sample $S_{\text{NP}}$ and a probability sample $S_{\text{P}}$ and DR inference is conducted within the $qp$ or $\xi p$ framework without specifying which one is correct. The general formula for the DR estimator is given by $$ \tilde{\mu}_{y, \text{DR}} = \frac{1}{N}\sum_{i \in S_{\text{NP}}}\frac{y_i - m_i}{\pi_{\text{NP}, i}} + \frac{1}{N}\sum_{i =1}^N m_i, $$ where $\pi_{\text{NP}, i}$ and $m_i$ are defined as previously. In the next sections we discuss two approaches to the DR estimation. ### Parameters estimated separately @chen2020doubly proposed two DR estimators given in Equation 9 and Equation 10 assuming that the population size is either known or estimated: $$ \hat{\mu}_{y, \text{DR-HT}}=\frac{1}{N} \sum_{i \in S_{\text{NP}}} d_{\text{NP}, i}\left\{y_i-m\left(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}}\right)\right\}+\frac{1}{N} \sum_{i \in S_{\text{P}}} d_{\text{P}, i} m\left(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}}\right), \tag{9} $$ and $$ \hat{\mu}_{y, \text{DR-Hájek}}=\frac{1}{\hat{N}_{\text{NP}}} \sum_{i \in S_{\text{NP}}} d_{\text{NP}, i}\left\{y_i-m\left(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}}\right)\right\}+\frac{1}{\hat{N}_{\text{P}}} \sum_{i \in S_{\text{P}}} d_{\text{P}, i} m\left(\boldsymbol{x}_i, \hat{\boldsymbol{\beta}}\right), \tag{10} $$ where $d_{\text{NP}, i}=\pi\left(\boldsymbol{x}_i, \hat{\boldsymbol{\gamma}}\right)^{-1}$, $\hat{N}_{\text{NP}}=\sum_{i \in S_{\text{NP}}} d_{\text{NP}, i}$ and $\hat{N}_{\text{P}}=\sum_{i \in S_{\text{P}}} d_{\text{P}, i}$. The estimator in Equation 10, including the estimated population size, has a better performance in terms of bias and the mean squared error and should be used in practice. However, the main limitation is associated with variance estimation, which is discussed at the end of this section. @chen2020doubly suggested constructing the estimators in Equation 9 or 10 based on the two models estimated separately, which is an attractive option, since one can specify a different number of variables for the propensity and outcome model. An alternative approach, proposed by @yang_doubly_2020, and similar to the one described by @kim2014doubly, is discussed in the next section. ### Minimization of the bias for doubly robust methods @yang_doubly_2020 discussed variable selection for a high-dimensional setting and noted that post-selection asymptotic bias may be non-negligible. Therefore, according to @yang_doubly_2020, the idea is to develop estimating equations for the $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ parameters based on the bias of the population mean estimator. In this way, the parameters can be estimated in a single step, rather than in two separate steps. First, the authors derived the asymptotic bias of $\hat{\mu}_{y,\text{DR}}$, which is given by $$ \begin{aligned} a\!\operatorname{bias}(\boldsymbol\gamma,\boldsymbol\beta) &=\E\{\hat\mu_{y, \text{DR}}(\boldsymbol\gamma,\boldsymbol\beta)-\mu\} \\ &=\E\left[ \frac{1}{N}\sum_{i=1}^{N} \left\{\frac{I_{\text{NP},i}}{\pi(\boldsymbol x_i,\boldsymbol\gamma)}-1\right\} \{y_i-m(\boldsymbol x_i,\boldsymbol\beta)\} \right] \\ &\quad+ \E\left[ \frac{1}{N}\sum_{i=1}^{N} \{I_{\text{P},i}d_{\text{P},i}-1\}m(\boldsymbol x_i,\boldsymbol\beta) \right]. \end{aligned}$$ The goal of this approach is to minimize $\{a\!\operatorname{bias}(\boldsymbol\gamma,\boldsymbol\beta)\}^2$, which leads to the following system of empirical estimating equations: $$ \left(\begin{array}{c} \sum_{i=1}^N I_{\text{NP}, i} \frac{ \dot{\boldsymbol\pi}\left(\boldsymbol{x}_i,\boldsymbol{\gamma}\right) }{ \pi\left(\boldsymbol{x}_i,\boldsymbol{\gamma}\right)^2 } \left\{ y_i-m\left(\boldsymbol{x}_i,\boldsymbol{\beta}\right) \right\} \\[1.2em] \sum_{i=1}^N \frac{I_{\text{NP}, i}} {\pi\left(\boldsymbol{x}_i,\boldsymbol{\gamma}\right)} \dot{m}\left(\boldsymbol{x}_i,\boldsymbol{\beta}\right) - \sum_{i \in S_{\text{P}}} d_{\text{P}, i} \dot{m}\left(\boldsymbol{x}_i,\boldsymbol{\beta}\right) \end{array}\right) = \boldsymbol{0} \tag{11} $$ where $\dot{m}\left(\boldsymbol{x}_i,\boldsymbol{\beta}\right)= \frac{\partial m\left(\boldsymbol{x}_i,\boldsymbol{\beta}\right)} {\partial \boldsymbol{\beta}}$ and $\dot{\boldsymbol\pi}\left(\boldsymbol{x}_i,\boldsymbol{\gamma}\right)=\frac{\partial \pi\left(\boldsymbol{x}_i,\boldsymbol{\gamma}\right)} {\partial \boldsymbol{\gamma}}$. This formulation is valid for any differentiable link function used in the sampling-score model. For the logit link, Equation 11 reduces to the estimating-equation form used by @yang_doubly_2020. The system in Equation 11 can be solved using the Newton--Raphson method. Without variable selection and under the logit sampling-score model, this approach reduces to the @kim2014doubly estimating-equation system used by @chen2020doubly and @wu2022statistical for doubly robust inference and doubly robust variance estimation in non-probability samples. The main limitation of this approach is that the joint estimating-equation blocks must be formulated on compatible covariate vectors and have the same dimension. Following @yang_doubly_2020, this is achieved after variable selection by re-estimating both nuisance models on the union of selected variables, $C=\widehat M_\gamma\cup\widehat M_\beta$, where $\widehat M_\gamma$ and $\widehat M_\beta$ denote the sets of variables selected for the propensity score and outcome models, respectively. This construction does not imply unbiasedness when both nuisance models are misspecified; its role is to reduce the first-order post-selection effect and to preserve root-$n_{\min}$ inference, with $n_{\min}=\min(|S_{\text{P}}|,|S_{\text{NP}}|)$, under correct specification of either the sampling-score model or the outcome model. In the package, the bias-minimization approach is implemented using the general derivative form in Equation 11, which allows different link functions for the propensity-score model. For the logit link, this formulation coincides with the estimating equation of @yang_doubly_2020. Note that we use a different convention for $\boldsymbol h$: our choice $\boldsymbol h(\boldsymbol{x};\boldsymbol\gamma)=\boldsymbol{x}/\pi(\boldsymbol{x};\boldsymbol\gamma)$ corresponds to the IPW estimating equation of @yang_doubly_2020 written with $\boldsymbol h(\boldsymbol X;\boldsymbol\alpha)=\boldsymbol X$. As mentioned above, the choice of either Equation 9 or 10 results in a different approach to variance estimation, which is discussed in the next section. ### Variance estimators for the doubly robust approach @yang_doubly_2020 derived a closed-form variance estimator for the known-$N$ DR estimator in Equation 9. This result requires knowledge of the population size $N$ and the corresponding bias-correction term to obtain a valid estimator of $V_{\xi p}\left(\hat{\mu}_{y,\text{DR}}-\mu_y\right)$ under the outcome regression model $\xi$. @chen2020doubly also consider Hájek-type DR estimators with estimated denominators, which motivates the version in Equation 10. However, an analytical variance estimator for the exact Hájek-type DR estimator implemented in the package has not, to our knowledge, been formally established. Therefore, in the package we provide a plug-in approximation obtained by replacing $N$ in the known-$N$ variance formula by the estimated denominators $\hat{N}_{\text{NP}}$ and $\hat{N}_{\text{P}}$. This option should be interpreted as heuristic and used with caution. Alternatively, one can use the bootstrap approach. @chen2020doubly demonstrated that the bootstrap approach presented in Section [2.2](#sec-prediction){reference-type="ref" reference="sec-prediction"} performs well in terms of the coverage rate when one of the working models is correctly specified. This is why this approach is recommended for all users. ## Variable selection algorithms @yang_asymptotic_2020 point out that it is crucial to use variable selection techniques during estimation, especially when dealing with high-dimensional non-probability samples. Variable selection not only improves model stability and computational feasibility, but also reduces variance, which can increase when irrelevant auxiliary variables are included. The most popular approaches described in the literature are penalization methods, such as least absolute shrinkage and selection operator (LASSO), smoothly clipped absolute deviation (SCAD) or minimax concave penalty (MCP), which, thanks to the appropriate loss functions, shrink the coefficients in variables that have no significant effect on the dependent variable [cf. @tibshirani1996regression; @ncvreg]. The selection procedure for non-probability methods works in a similar way, with loss functions modified to account for external data sources, such as sample or population totals or averages. In particular, the technique consists of two steps: 1. select the relevant variables using an appropriately constructed loss function (and possibly using the approach shown in Equation 11 to obtain the final estimates of the model parameters), 2. construct the selected estimator using variables selected from the first step. For instance, @yang_doubly_2020 used Equation 12 as a loss function for estimating outcome equation parameters: $$ \operatorname{Loss}\left(\lambda_{\boldsymbol{\beta}}\right)=\sum_{i=1}^N I_{\text{NP}, i}\left[y_i-m\left\{\boldsymbol{x}_i, \boldsymbol{\beta}(\lambda_{\boldsymbol{\beta}})\right\}\right]^2, \tag{12} $$ where $m\left\{\boldsymbol{x}_i, \boldsymbol{\beta}(\lambda_{\boldsymbol{\beta}})\right\}$ is the penalized function for the $\boldsymbol{\beta}$ parameters with a tuning parameter $\lambda_{\boldsymbol{\beta}}$ and loss functions for the PS function presented in Table 4, where $\lambda_{\boldsymbol{\gamma}}$ is the tuning parameter and $\boldsymbol{x}_{i, l}$ denotes $l$-th component of $\boldsymbol{x}_{i}$. Table: **Table 4.** Loss functions for the PS function depending on the $\boldsymbol{h}(\cdot, \cdot)$ function. | $\boldsymbol{h}(\boldsymbol{x}_i, \boldsymbol{\gamma})$ function | Loss function | |:---|:---| | $\boldsymbol{x}_i$ | $\sum_{l=1}^L\left(\sum_{i=1}^N\left[I_{\text{NP}, i} - \frac{I_{\text{P}, i}\pi\left\{\boldsymbol{x}_i, \boldsymbol{\gamma}(\lambda_{\boldsymbol{\gamma}})\right\}}{\pi_{\text{P}, i}}\right] \boldsymbol{x}_{i, l}\right)^2$ | | $\boldsymbol{x}_i \pi(\boldsymbol{x}_i, \boldsymbol{\gamma})^{-1}$ | $\sum_{l=1}^L\left(\sum_{i=1}^N\left[\frac{I_{\text{NP}, i}}{\pi\left\{\boldsymbol{x}_i, \boldsymbol{\gamma}(\lambda_{\boldsymbol{\gamma}})\right\}}-\frac{I_{\text{P}, i}}{\pi_{\text{P}, i}}\right] \boldsymbol{x}_{i, l}\right)^2$ | @yang_doubly_2020 only discussed the SCAD penalty and the $\boldsymbol{h}(\boldsymbol{x}_i, \boldsymbol{\gamma})=\boldsymbol{x}_i$ function for the DR estimator. In the **nonprobsvy** package we have extended this approach to the first variant of $\boldsymbol{h}(\boldsymbol{x}_i, \boldsymbol{\gamma})$, shown in the first row of Table 4, and have allowed the user to select other link functions for the $\pi_{\text{NP}, i}$, implemented other penalty functions and extended the possibility of selecting variables to MI and IPW estimators. In the next section we discuss how to define the approaches presented above. # The main function and the package functionalities ## The `nonprob` function The **nonprobsvy** package is built around the `nonprob()` function. The main design objective was to make the use of `nonprob` as similar as possible to standard R functions for fitting statistical models, such as `stats::glm()`, while incorporating survey design features from the **survey** package. The most important arguments are given in Table 5 and the obligatory ones include `data` as well as one of the following three -- `selection`, `outcome`, or `target` -- depending on which method has been selected. In the case of `outcome` and `target` multiple $y$ variables can be specified. Table: **Table 5.** A description of the `nonprob()` function arguments. | Argument | Description | |:---|:---| | `data` | A `data.frame` with data from the non-probability sample of size $n_{\text{NP}}$. The data frame must contain all variables referenced in the `selection`, `outcome`, and `target` formulas. | | `selection` | A `formula` for the selection equation (e.g., `~ x1 + x2`). This formula defines the PS model (for the IPW estimator). | | `outcome` | A `formula` for the outcome equation (e.g., `y1 + y2 ~ x1 + x2`). We allow for multiple target variables under the same prediction model (for the MI and DR estimator). | | `target` | A `formula` with target variables (e.g., `~ y1 + y2 + y3`). | | `svydesign` | An optional `svydesign2` object containing the probability sample of size $n_{\text{P}}$. Must include all variables specified in `selection`, `outcome`, and `target` formulas. Used for calibration and validation of the non-probability sample estimates. | | `pop_totals`, `pop_means`, `pop_size` | An optional named `vector` with population totals or means of the covariates and population size. For `pop_totals` and `pop_means`, length must equal the number of covariates in the selection model, with names matching variable names in the `data` frame. | | `method_selection` | A link function for the PS (`c("logit", "probit", "cloglog")`). | | `method_outcome` | Specification of the MI approach (`c("glm", "nn", "pmm", "npar")`). | | `family_outcome` | A GLM family for the MI approach (`c("gaussian", "binomial", "poisson")`). | | `subset` | An optional `vector` specifying a subset of observations to be used in the fitting process. | | `strata` | An optional parameter for estimation for sub-populations (currently not supported). | | `case_weights` | An optional `vector` of prior case (frequency) weights to be used in the fitting process. | | `na_action` | A `function` indicating what to do with `NA`'s (default `na.omit`). | | `control_selection`, `control_outcome`, `control_inference` | Control functions with parameters for PS, the outcome model and variance estimation, respectively. | | `start_selection`, `start_outcome` | An optional `vector` with starting values for the parameters of the PS and the outcome equation. | | `verbose` | A `logical` value indicating if information should be printed. | | `se` | A `logical` value indicating whether to calculate and return the standard error of the estimated mean. | | `...` | Additional optional argument for further development (currently not supported). | The package allows the user to provide either reference population data (via the `pop_totals`, or `pop_means` and `pop_size`) or a probability sample declared by the `svydesign` argument (`svydesign2` class from the **survey** package). The `nonprob()` function is used to specify inference methods through the `selection` and `outcome` arguments. If a `svydesign2` object is provided and the `selection` argument is specified, then the IPW estimators are used (by default parameters of the PS model employ Equation 7), if the `outcome` argument is specified, then the MI approach is used (the default option is the MI-GLM with the `gaussian` family) and if both are specified, then the DR approach is applied (parameters $(\boldsymbol{\beta}, \boldsymbol{\gamma})$ are estimated separately and the Equation 10 is used). A particular inference method is selected through the `method_selection`, `method_outcome`, `family_outcome`, `control_selection` and `control_outcome` arguments. The variance estimation method is selected through the `control_inference` argument. Resulting object of class `nonprob` is a list that contains the following (most important) elements: * `data` -- a `data.frame` containing the non-probability sample. * `X` -- a `matrix` containing both samples. * `y` -- a `list` containing all variables declared in either the `target` or `outcome` arguments. * `R` -- a numeric `vector` informing about inclusion in the non-probability sample. * `ps_scores` -- propensity scores or `NULL` (for the MI estimators). * `ipw_weights` -- inverse probability weights or `NULL` (for the MI estimators). * `output` -- a `data.frame` containing point and standard error estimates. * `confidence_interval` -- a `data.frame` containing confidence intervals for the mean. * `outcome` -- a `list` of results for each `outcome` model. * `selection` -- a `list` of results for the `selection` model. * `svydesign` -- a `svydesign2` object passed by the `svydesign` argument. * `ys_rand_pred, ys_nons_pred` -- a `list` of predicted values from the MI model. ## Controlling the type of estimators The `control_out()` function can be used to specify various aspects of the estimation process, including the variable selection methods (through different `penalty` options like SCAD, LASSO, and MCP with their respective tuning parameters defined in the same way as in the `control_sel()` function; cf. @yang_doubly_2020), and detailed configuration for NN, PMM and NPAR approaches (using parameters like `pmm_match_type`, `pmm_weights`, `pmm_k_choice` or `npar_loess`). For NN and PMM approaches we use the **RANN** package [@rann-pkg] with the kd-tree algorithm and the Euclidean distance as default. We currently do not support other distances (e.g., Gower). Table 6 presents example usage of the `control_out()` function for three types of MI estimators. Table: **Table 6.** Example declarations of the MI estimators. | Estimator | Declaration with `control_out()` | |:---|:---| | MI-GLM with the LASSO penalty and 5 folds | `nonprob(outcome = y1 ~ x1 + x2, data = df, svydesign = prob, control_outcome = control_out(penalty = "lasso", nfolds = 5))` | | MI-NN with the `bd` algorithm | `nonprob(outcome = y1 ~ x1 + x2, data = df, svydesign = prob, control_outcome = control_out(treetype = "bd"))` | | MI-PMM-B with $k=3$ | `nonprob(outcome = y1 ~ x1 + x2, data = df, svydesign = prob, control_outcome = control_out(k = 3, pmm_match_type = 2))` | The `control_sel()` function provides essential control parameters for fitting the selection model in the `nonprob()` function. It allows users to select between MLE given by Equation 7 or GEE defined in Equation 8 through the `est_method` argument, specify the $\boldsymbol{h}(\cdot, \cdot)$ function through the `gee_h_fun` argument, specify the optimizer (`optimizer` argument) and which variable selection method should be applied (using different penalty functions like SCAD, lasso, and MCP by specifying the `penalty` argument) along with parameters (e.g., the number of folds through the `nfolds` argument). The parameters of the PS for the calibrated IPW is estimated by using the **nleqslv** package and fitting parameters (arguments starting with the `nleqslv_*`). Table 7 presents example usage of the `control_sel()` function for two types of IPW estimators. Table: **Table 7.** Example declarations of the IPW estimators. | Estimator | Declaration with `control_sel()` | |:---|:---| | Calibrated IPW | `nonprob(selection = ~ x1 + x2, target = ~ y1, data = df, svydesign = prob, control_selection = control_sel(est_method = "gee"))` | | IPW with the MCP penalty and 5 folds | `nonprob(selection = ~ x1 + x2, target = ~ y1, data = df, svydesign = prob, control_selection = control_sel(penalty = "MCP", nfolds = 5))` | ## Controlling variance estimation Finally, the `control_inf()` function configures the parameters for variance estimation in the `nonprob()` function. It allows users to specify whether the analytical or bootstrap approach should be used (the `var_method` argument), whether the variable selection method should be applied (the `vars_selection` argument) and what type of bootstrap should be applied for the probability sample (the `rep_type` argument). This function is also used to specify the inference procedure for the DR approach: if a union of variables should be applied (the `vars_combine` argument) and if the bias correction (the `bias_correction` argument) should be applied after variable selection (the `vars_selection`) and variable combination. Table 8 presents example usage of the `control_inf()` function for the IPW and DR estimators. Table: **Table 8.** Example declarations of the IPW and DR estimators. | Estimator | Declaration with `control_inf()` | |:---|:---| | Calibrated IPW with variable selection, bootstrap and $B=50$ | `nonprob(selection = ~ x1 + x2, target = ~ y1, data = df, svydesign = prob, control_selection = control_sel(est_method = "gee"), control_inference = control_inf(vars_selection = TRUE, var_method = "bootstrap", rep_type = "subbootstrap", B = 50))` | | The DR with the SCAD penalty, 5 folds and bias correction | `nonprob(selection = ~ x1 + x2, outcome = y1 ~ x1 + x2, data = df, svydesign = prob, control_selection = control_sel(penalty = "SCAD", nfolds = 5), control_inference = control_inf(vars_selection = TRUE, bias_correction = TRUE, vars_combine = TRUE))` | In the next sections we present a case study illustrating the process of integrating a non-probability sample with a reference probability sample. We present various estimators and compare them. Finally, we describe more advanced options available in the package. # Data analysis example ## Description of the data The package can be installed in the standard manner using: ```{r install, eval = FALSE} install.packages("nonprobsvy") ``` Before we explain the case study let us first load the necessary packages for this paper. ```{r libs, message = FALSE, warning = FALSE} library("nonprobsvy") library("ggplot2") ``` The goal of the case study was to integrate survey (`jvs`) and administrative (`admin`) data about job vacancies in Poland. The first source, the job vacancy survey (JVS), contains 6,523 units. The JVS provides a probability sample drawn according to a stratified sampling design. More details can be found in @jvs2022. The dataset contains information about the industry code (14 levels, the `nace` column), `region` (16 levels), sector (2 levels, the `private` column), company size (3 levels: Small up to 9, Medium 10-49 and Large 50+) and the final weight (i.e., the design weight corrected for non-contact and non-response), which is treated as the $d$ weight. ```{r jvs-data} data("jvs", package = "nonprobsvy") head(jvs) ``` Since the **nonprobsvy** package relies on the functionalities of the **survey** package, we need to define the `svydesign2` object via the `svydesign` function, as shown below. The dataset does not contain the true stratification variable, so we use a simplified version by specifying `~ size + nace + region`; similarly, since we do not have information regarding non-response and its correction, we simply assume that the `weight` sums up to the population size. ```{r jvs-design} jvs_svy <- svydesign(ids = ~ 1, weights = ~ weight, strata = ~ size + nace + region, data = jvs) ``` Our second source (`admin`), the Central Job Offers Database, is a register containing all vacancies submitted to Public Employment Offices (see ). We treat this register as a non-probability sample since it contains administrative data provided on a voluntary basis, so the inclusion mechanism is unknown. This dataset was prepared in such a way that records deemed to be out of scope (either in terms of the definition of vacancy or the population of entities) were excluded. In addition to the same variables found in the JVS, the dataset contains one called `single_shift`, which is our target variable, defined as: whether a company seeks at least one employee for a single-shift job. The goal of this case study is to estimate the share of companies that seek employees for a single-shift job in Poland in a given quarter. ```{r admin-data} data("admin", package = "nonprobsvy") head(admin) ``` One should keep in mind that this paper does not aim to provide a complete tutorial on how to use non-probability samples for statistical inference. We therefore do not include the stage of aligning variables to meet the same definitions, assessing the strength of the relation between auxiliary variables and the target variable, the selection mechanism and the distribution mis-matches between both samples. In the examples below we assume that there is no overlap between both sources and the naive, reference estimate, given by the mean of the `single_shift` column of `admin`, which is equal to 66.1%. ## Estimation ### Propensity score approach First, we start with the IPW approach, which offers the choice between two estimation methods: MLE (default) and GEE (calibrated to the estimated survey totals). We start by calling the `nonprob` function, where we define the `selection` argument indicating which variables are to be included, the `target` argument, which specifies the variable of interest, i.e., `single_shift`. The remaining arguments specify the `svydesign` object, the dataset and the link function (`method_selection`; `"logit"` is the default value). As the `pop_size` argument was not specified, the $\hat{\mu}_{y, \text{IPW-Hájek}}$ is returned. ```{r ipw1} ipw_est1 <- nonprob( selection = ~ region + private + nace + size, target = ~ single_shift, svydesign = jvs_svy, data = admin, method_selection = "logit" ) ``` In order to get the basic information about the estimated target quantity we can use the `print` method to display the object. It provides information about the methods, source of auxiliary variable, variable selection, the naive (uncorrected estimator) and the corrected along with the standard error and confidence interval. The last line, i.e., `selected estimator`, is limited to the maximum of 5 variables. ```{r ipw1-print} ipw_est1 ``` If we want to see detailed information about the approach, we can use the `summary` method. This function provides information on the sample sizes, sum of the IPW weights along with the distribution of IPW weights (for non-probability sample) and propensity scores (denoted as IPW probabilities) for both samples. A more advanced user may be interested in inspecting the details of the fit by accessing the `selection` element (i.e., `ipw_est1$selection`). ```{r ipw1-summary} summary(ipw_est1) ``` To extract the point estimate along with its standard error and 95% confidence interval in a form of `data.frame` we can use the `extract` method as shown below. ```{r ipw1-extract} extract(ipw_est1) ``` If we want to use the calibrated IPW approach, it is necessary to define the `control_sel()` function in the `control_selection` argument by setting the `est_method` argument equal to `"gee"` (the default is `"mle"`) and the value of `gee_h_fun`. ```{r ipw2} ipw_est2 <- nonprob( selection = ~ region + private + nace + size, target = ~ single_shift, svydesign = jvs_svy, data = admin, method_selection = "logit", control_selection = control_sel(gee_h_fun = 1, est_method = "gee") ) ``` Results are comparable to the standard IPW point estimate (70.4 vs 72.2) while the standard error is slightly higher. ```{r ipw2-print} ipw_est2 ``` The calibrated IPW significantly improves the balance, which can be accessed via the `check_balance()` function: ```{r ipw-balance} data.frame(ipw_mle = check_balance(~ size - 1, ipw_est1, 1)$balance, ipw_gee = check_balance(~ size - 1, ipw_est2, 1)$balance) ``` Notice that neither in the package nor in this paper do we provide a detailed description of the post-hoc results, such as the covariate balance. This can be done using existing CRAN packages, e.g., through the `bal.tab()` function from the **cobalt** package [@cobalt]. ### Prediction-based approach If the user is interested in the prediction-based approach, in particular involving MI estimators, then, they should specify the argument `outcome` as a formula (as in the case of the `glm()` function). We allow a single outcome (specified as `y ~ x1 + x2 + ... + xk`) and multiple outcomes (as `y1 + y2 + y3 ~ x1 + x2 + ... + xk`). Note that if the `outcome` argument is specified, then there is no need to specify the `target` argument. By default, the GLM type of an MI estimator is used (i.e., `method_outcome = "glm"`). In the code below we show how this type of an MI estimator can be declared. ```{r mi1} mi_est1 <- nonprob( outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_outcome = "glm", family_outcome = "binomial" ) mi_est1 ``` In order to employ an MI estimator based on NN matching, one can specify `method_outcome = "nn"` for the nearest neighbors search using all variables specified in the `outcome` argument, `method_outcome = "pmm"` to use PMM or `method_outcome = "npar"` to use non-parametric approach. For the NN and PMM estimators, we employ $k=5$ nearest neighbors (i.e., `control_out(k=5)`) but the variance of the PMM estimator should be estimated using bootstrap approach suggested by @chlebicki2025. However, in this example we have decided not to use bootstrap as all variables are categorical and this results in a large number of NNs having the same distance and the numerical approximation defined in `control_out(eps)` may give slightly different results between platforms. Therefore the reported variance is based on the probability sample $S_{\text{P}}$ only. In the case of the MI-NN estimator there is no need to specify the `family_outcome` argument as no model is estimated underneath. ```{r mi23} mi_est2 <- nonprob( outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_outcome = "nn", control_outcome = control_out(k = 5) ) mi_est3 <- nonprob( outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_outcome = "pmm", family_outcome = "binomial", control_outcome = control_out(k = 5) ) ``` Results of both estimators are more or less similar, but it should be noted that the NN version suffers from the curse of dimensionality, so the PMM version seems to be more reliable. ```{r mi23-extract} rbind("NN" = extract(mi_est2)[, 2:3], "PMM" = extract(mi_est3)[, 2:3]) ``` As discussed in the methods section, IPW and MI estimators are consistent only when the model and auxiliary variables are correctly specified. To overcome this problem, the user can turn to doubly robust estimators. ### The doubly robust approach In order to choose doubly robust estimation the user needs to specify both the `selection` and `outcome` arguments. These formulas can be specified with the same or varying number of auxiliary variables. As in the MI approach, we also allow multiple outcomes. In the following example code we have specified the non-calibrated IPW and the MI-GLM estimator. ```{r dr1} dr_est1 <- nonprob( selection = ~ region + private + nace + size, outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_selection = "logit", method_outcome = "glm", family_outcome = "binomial" ) dr_est1 ``` Detailed results can be displayed by using the `summary` method, which prints both information about the distribution of the IPW weights as well as the outcome residuals and predictions. For more details one can access the `"outcome"` and `"selection"` fields in the `nonprob` object. ```{r dr1-summary} summary(dr_est1) ``` Finally, we can use the bias minimization approach, as proposed by @yang_doubly_2020, by specifying the `control_inference = control_inf(bias_correction = TRUE)` argument together with the `vars_combine = TRUE` and `vars_selection = TRUE` as this approach requires variable selection followed by variable union from both equations. To keep the example fast we use 3 folds and 10 values of $\lambda$ for both models. ```{r dr2} set.seed(2026) dr_est2 <- nonprob( selection = ~ region + private + nace + size, outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_selection = "logit", method_outcome = "glm", family_outcome = "binomial", control_selection = control_sel(nfolds = 3, nlambda = 10), control_outcome = control_out(nfolds = 3, nlambda = 10), control_inference = control_inf(bias_correction = TRUE, vars_combine = TRUE, vars_selection = TRUE) ) dr_est2 ``` ## Comparison of estimates Finally, as there is no single method for non-probability samples, we suggest comparing results in a single table or a plot. Figure [1](#fig:comparison-of-est){reference-type="ref" reference="fig:comparison-of-est"} presents point estimates along with 95% confidence intervals. The various estimators show interesting patterns compared to the naive estimate (red dotted line). Both IPW estimators result in the point estimates over 70% but IPW GEE is characterized by large standard error. The MI estimators demonstrate similar pattern but MI-PMM and MI-NN produces the widest confidence interval and MI-NN is characterized with the lowest point estimate. The estimates from MI-GLM and the DR estimators, with and without bias minimization, are more closely clustered and have similar confidence interval widths. All adjusted methods indicate a population parameter higher than the naive estimate, although the magnitude of the adjustment differs across methods. This agreement provides useful sensitivity evidence that the substantive conclusion is not highly dependent on the particular adjustment method, but it should not be interpreted as formal validation that selection bias has been fully removed. ```{r comparison-of-est, fig.cap = "Figure 1: Comparison of estimates of the share of job vacancies offered on a single-shift."} df_s <- rbind(extract(ipw_est1), extract(ipw_est2), extract(mi_est1), extract(mi_est2), extract(mi_est3), extract(dr_est1), extract(dr_est2)) df_s$est <- c("IPW (MLE)", "IPW (GEE)", "MI (GLM)", "MI (NN)", "MI (PMM)", "DR", "DR (BM)") ggplot(data = df_s, aes(y = est, x = mean, xmin = lower_bound, xmax = upper_bound)) + geom_point() + geom_vline(xintercept = mean(admin$single_shift), linetype = "dotted", color = "red") + geom_errorbar() + labs(x = "Point estimator and confidence interval", y = "Estimators") + theme_bw() ``` ## Advanced usage ### Bootstrap approach for variance estimation In the package we allow the user to estimate the variance of the mean analytically (by default) or using the bootstrap approach, as described in the prediction approach section. We use analytical variance estimators proposed in the papers referenced in the methods section. The calculation of the standard error can be disabled using `nonprob(se = FALSE)`. The bootstrap approach implemented in the package refers to: * the non-probability sample -- currently only simple random sampling with replacement is available via the `base::sample()`, * the probability sample -- all the approaches implemented in the `as.svrepdesign()` function of the **survey** package are supported and we refer the reader to the relevant help file. The bootstrap approach is specified via the `control_inf()` function with `var_method = "bootstrap"`. The bootstrap method for the probability sample is controlled via the `rep_type` argument, which passes the method to the `as.svrepdesign()` function. The number of iterations is set in the `num_boot` argument (100 by default). If the samples are large or the estimation method is complicated (e.g., involves variable selection) one can set `verbose = TRUE` to track progress. By default bootstrap results are stored in the `boot_sample` element of the resulting list (to disable this option, `keep_boot` should be set to `FALSE`). The following code is an example of applying the IPW approach with the bootstrap approach specified by the argument `control_inference` of the `nonprob()` function. ```{r ipw-boot} set.seed(2024) ipw_est1_boot <- nonprob( selection = ~ region + private + nace + size, target = ~ single_shift, svydesign = jvs_svy, data = admin, method_selection = "logit", control_inference = control_inf(var_method = "bootstrap", num_boot = 50), verbose = FALSE ) ``` Next, we compare the estimated standard error of variance estimation for the analytical and the bootstrap approach. ```{r ipw-boot-compare} rbind("IPW analytic variance" = extract(ipw_est1)[, 2:3], "IPW bootstrap variance" = extract(ipw_est1_boot)[, 2:3]) ``` The boot samples can be accessed via the `boot_sample` element of the output list of the `nonprob()` function. Note that the output is returned as a `matrix` because we allow multiple `target` variables. ```{r ipw-boot-sample} head(ipw_est1_boot$boot_sample, n = 3) ``` ### Variable selection algorithms In this section we briefly show how to use variable selection algorithms. In order to indicate that a variable selection algorithm should be used one should specify the `control_inference = control_inf(vars_selection = TRUE)` argument. Then, the user should either leave the default setting or specify the outcome parameters via the `control_out()` function or the `control_sel()` function. Both functions have the same parameters: * `penalty` -- The penalization function used during variables selection (possible values: `c("SCAD", "lasso", "MCP")`). * `nlambda` -- The number of $\lambda$ values; by default set to 50 (grid search). * `lambda_min` -- The smallest value for $\lambda$, as a fraction of `lambda.max`; 0.001 by default. * `lambda` -- A user specified vector of lambdas (only for the `control_sel()` function). * `nfolds` -- The number of folds for cross validation; by default set to 10. * `a_SCAD, a_MCP` -- The tuning parameter of the SCAD and MCP penalty for the selection model; by default set to 3.7 and 3, respectively. In the case of the MI approach we rely on the **ncvreg** package [@ncvreg], which is the only R package that employs the SCAD method. For the IPW and DR approaches, we have developed our own codes in C++ via the **Rcpp** and **RcppArmadillo** packages. In the code below we apply variable selection for the MI-GLM estimator using only 3 folds, 10 possible values of $\lambda$ parameters and the LASSO penalty. ```{r mi-sel} set.seed(2024) mi_est1_sel <- nonprob( outcome = single_shift ~ region + private + nace + size, svydesign = jvs_svy, data = admin, method_outcome = "glm", family_outcome = "binomial", control_outcome = control_out(nfolds = 3, nlambda = 10, penalty = "lasso"), control_inference = control_inf(vars_selection = TRUE), verbose = TRUE ) ``` In this case study, the MI-GLM estimator with variable selection yields almost the same results as the approach without it. Point estimates and standard errors differ at the fourth and third digit, respectively. ```{r mi-sel-compare} rbind("MI without var sel" = extract(mi_est1)[, 2:3], "MI with var sel" = extract(mi_est1_sel)[, 2:3]) ``` The result object of the `cv.ncvreg` class is stored in the `"outcome"` element of the result, and you can access the selected variables using the `coef` generic method as follows. ```{r mi-sel-coef} round(coef(mi_est1_sel)$coef_out[, 1], 4) ``` If a user is interested in viewing the PS estimation results, then access to the `"selection"` element is required. ```{r ipw-coef} round(coef(ipw_est1)$coef_sel[, 1], 4) ``` # Classes and S3 methods The package contains the main class `nonprob` and the supplementary class `nonprob_method` and the related `nonprob_summary` class. All available S3 methods can be obtained by calling `methods(class = "nonprob")`. For instance, the `check_balance()` function, already mentioned in the case study, is used to view the balance by checking how the PS weights reproduce known or estimated population totals; the `nobs()` function returns the sample size of the probability and non-probability samples. ```{r nobs} nobs(dr_est1) ``` Table 9 presents S3 methods implemented for objects of class `nonprob`. We have intentionally limited the number of implemented methods as the goal of the package is to provide point and interval estimates. If users are interested in assessing the quality of the models or the covariate balance, they should use existing R packages. Table: **Table 9.** S3 methods implemented for objects of class `nonprob` in **nonprobsvy**. | Function | Description | |:---|:---| | `check_balance` | Returns a `list` of totals for IPW and DR estimators. | | `coef` | Returns a list of coefficients for outcome or selection model. | | `confint` | Returns the confidence interval for the target variable(s). | | `extract` | Returns a `data.frame` with estimation results. | | `nobs` | Returns the number of observations for both samples. | | `plot` | Plots a comparison of naive (uncorrected) and estimated mean along with CI. | | `pop_size` | Returns the population size (fixed or not). | | `summary` | Returns a `nonprob_summary` class with additional information. | | `update` | Returns an `updated` object. | | `weights` | Returns IPW weights for the IPW and DR estimator. | The confidence interval for the target variable(s) can be estimated using the generic `confint()` function. ```{r confint} confint(dr_est1, level = 0.99) ``` If a user is interested in assessing the distribution of the IPW weights, we suggest using the `weights()` function. ```{r weights} summary(weights(dr_est1)) ``` In the case of mass imputation estimators it is possible to use special functions `method_glm()` that are of class `nonprob_method`. An example is given below ```{r method-glm} res_glm <- method_glm( y_nons = admin$single_shift, X_nons = model.matrix(~ region + private + nace + size, admin), X_rand = model.matrix(~ region + private + nace + size, jvs), svydesign = jvs_svy) res_glm ``` For the IPW we have created the function `method_ps()`, which returns requested functions for a specific link. This approach is motivated by the following reasons: 1) we allow different estimation techniques, such as maximum likelihood and general estimation equations; 2) the IPW estimator is also used for the DR estimator; and 3) variance estimators depend on whether the IPW or the DR estimator is applied (with combination with the MI estimator). ```{r method-ps} method_ps() ``` Details can be found in help page. ```{r method-ps-help, eval = FALSE} ?method_ps() ``` # Summary and future work The **nonprobsvy** package provides a comprehensive R software solution that addresses inference challenges connected with non-probability samples by integrating them with probability samples or known population totals/means. As non-probability data sources, like administrative registers, voluntary online panels, and social media data become increasingly available, statisticians need robust methods to produce reliable population estimates. The package implements state-of-the-art approaches including mass imputation, inverse probability weighting, and doubly robust methods, each designed to correct selection bias by leveraging auxiliary data. By providing a unified framework and its integration with the **survey** package, the **nonprobsvy** makes complex statistical methods for non-probability samples more accessible, enabling researchers to produce robust estimates even when working with non-representative data. There are several avenues for future development of the **nonprobsvy** package. One key priority is to implement model-based calibration and additional methods for estimating propensity scores and weights. The package currently assumes no overlap between probability and non-probability samples, so accounting for potential overlap (e.g., in big data sources and registers) is another important extension. Additional planned developments include handling non-ignorable sample selection mechanisms, developing a theory for maintaining consistency with calibration weights, and supporting multiple non-probability samples from various sources for the purpose of data integration. Further methodological extensions under consideration include empirical likelihood approaches for doubly/multiply robust estimation, integration of machine learning methods like debiased/double machine learning from causal inference, handling measurement errors in big data variables, and expanding the bootstrap approach beyond simple random sampling with replacement. The package will also be extended to handle the `svyrep.design` class from the **survey** package and the **svrep** package. These developments will enhance its capabilities for handling complex survey data structures and modern estimation challenges. # Acknowledgements The authors' work has been financed by the National Science Centre in Poland, OPUS 20, grant no. 2020/39/B/HS4/00941. Łukasz Chrostowski was the main developer and maintainer of the package up to version 0.1.0. Parts of this paper are based on Łukasz's Master's thesis (available at ). Piotr Chlebicki has contributed to the package and has implemented MI-PMM estimators. Maciej Beręsewicz was responsible for the design of the package and for testing, reviewing and contributing to the source code. He also prepared the manuscript. He was also responsible for a significant restructuring of the package between versions 0.1.0 and 0.2.0+. We also thank the authors of StatsClaw [@qinxu2026statsclaw], which helped us identify bugs and inconsistencies in our code. The authors thank the editors and reviewers for their constructive suggestions. # Appendix: Algorithms for the MI-NN and MI-PMM estimators **Algorithm 1.** Mass imputation using the $k$ nearest neighbor algorithm. 1. If $k=1$, then for each $i \in S_{\text{P}}$ match $\hat{\nu}(i,1)=\hat{\nu}(i)$ such that \[ \hat{\nu}(i)= \operatornamewithlimits{arg\, min}_{j\in S_{\text{NP}}} d\left(\boldsymbol{x}_i,\boldsymbol{x}_j\right) \] (e.g., Euclidean distance). 2. If $k>1$, then \[ \hat{\nu}(i, z) = \operatornamewithlimits{arg\, min}_{j\in S_{\text{NP}}\setminus\bigcup_{t=1}^{z-1} \{\hat{\nu}(i, t)\}} d\left(\boldsymbol{x}_i, \boldsymbol{x}_j\right), \] i.e., $\hat{\nu}(i, z)$ is the \(z\)th nearest neighbor from the sample $S_{\text{NP}}$. 3. For each $i \in S_{\text{P}}$, calculate the imputed value as \[ y_i^* = \frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i, t)}. \] **Algorithm 2.** Mass imputation using predictive mean matching variant: $\hat{y}-\hat{y}$ matching. 1. Estimate regression model $m(\boldsymbol{x}, \boldsymbol{\beta})$ parameters. 2. Predict \[ \hat{y}_{i}=m\left(\boldsymbol{x}_{i}, \hat{\boldsymbol{\beta}}\right), \qquad \hat{y}_{j}=m\left(\boldsymbol{x}_{j}, \hat{\boldsymbol{\beta}}\right) \] for $i\in S_{\text{P}}$ and $j\in S_{\text{NP}}$. 3. If $k=1$, then for each $i\in S_{\text{P}}$ match $\hat{\nu}(i,1)=\hat{\nu}(i)$ such that \[ \hat{\nu}(i)= \operatornamewithlimits{arg\, min}_{j\in S_{\text{NP}}} d\left(\hat{y}_{i}, \hat{y}_{j}\right). \] 4. If $k>1$, then \[ \hat{\nu}(i, z) = \operatornamewithlimits{arg\, min}_{j\in S_{\text{NP}}\setminus\bigcup_{t=1}^{z-1} \{\hat{\nu}(i, t)\}} d\left(\hat{y}_{i}, \hat{y}_{j}\right), \] i.e., $\hat{\nu}(i, z)$ is the \(z\)th nearest neighbor from the sample $S_{\text{NP}}$. 5. For each $i \in S_{\text{P}}$, calculate the imputed value as \[ y_i^* = \frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i, t)}. \] **Algorithm 3.** Mass imputation using predictive mean matching variant: $\hat{y}-y$ matching. 1. Estimate regression model $m(\boldsymbol{x}, \boldsymbol{\beta})$ parameters. 2. Predict \[ \hat{y}_{i}=m\left(\boldsymbol{x}_{i}, \hat{\boldsymbol{\beta}}\right) \] for $i \in S_{\text{P}}$. 3. If $k=1$, then for each $i\in S_{\text{P}}$ match $\hat{\nu}(i,1)=\hat{\nu}(i)$ such that \[ \hat{\nu}(i)= \operatornamewithlimits{arg\, min}_{j \in S_{\text{NP}}} d\left(\hat{y}_{i}, y_{j}\right). \] 4. If $k>1$, then \[ \hat{\nu}(i, z) = \operatornamewithlimits{arg\, min}_{j \in S_{\text{NP}} \setminus \bigcup_{t=1}^{z-1} \{\hat{\nu}(i, t)\}} d\left(\hat{y}_{i}, y_{j}\right), \] i.e., $\hat{\nu}(i, z)$ is the \(z\)th nearest neighbor from the sample $S_{\text{NP}}$. 5. For each $i \in S_{\text{P}}$, calculate the imputed value as \[ y_i^* = \frac{1}{k}\sum_{t=1}^{k}y_{\hat{\nu}(i, t)}. \] # References