The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
nonprobsvy
:
an R package for modern statistical inference methods based on
non-probability samplesThe goal of this package is to provide R users access to modern methods for non-probability samples when auxiliary information from the population or probability sample is available:
The package allows for:
ncvreg
, Rcpp
,
RcppArmadillo
packages),survey
package when probability
sample is available Lumley (2023),logit
,
probit
and cloglog
) and outcome
(gaussian
, binomial
and poisson
)
variables.Details on use of the package be found:
You can install the recent version of nonprobsvy
package
from main branch Github with:
::install_github("ncn-foreigners/nonprobsvy") remotes
or install the stable version from CRAN
install.packages("nonprobsvy")
or development version from the dev
branch
::install_github("ncn-foreigners/nonprobsvy@dev") remotes
Consider the following setting where two samples are available: non-probability (denoted as \(S_A\) ) and probability (denoted as \(S_B\)) where set of auxiliary variables (denoted as \(\boldsymbol{X}\)) is available for both sources while \(Y\) and \(\boldsymbol{d}\) (or \(\boldsymbol{w}\)) is present only in probability sample.
Sample | Auxiliary variables \(\boldsymbol{X}\) | Target variable \(Y\) | Design (\(\boldsymbol{d}\)) or calibrated (\(\boldsymbol{w}\)) weights | |
---|---|---|---|---|
\(S_A\) (non-probability) | 1 | \(\checkmark\) | \(\checkmark\) | ? |
… | \(\checkmark\) | \(\checkmark\) | ? | |
\(n_A\) | \(\checkmark\) | \(\checkmark\) | ? | |
\(S_B\) (probability) | \(n_A+1\) | \(\checkmark\) | ? | \(\checkmark\) |
… | \(\checkmark\) | ? | \(\checkmark\) | |
\(n_A+n_B\) | \(\checkmark\) | ? | \(\checkmark\) |
Suppose \(Y\) is the target variable, \(\boldsymbol{X}\) is a matrix of auxiliary variables, \(R\) is the inclusion indicator. Then, if we are interested in estimating the mean \(\bar{\tau}_Y\) or the sum \(\tau_Y\) of the of the target variable given the observed data set \((y_k, \boldsymbol{x}_k, R_k)\), we can approach this problem with the possible scenarios:
Estimator | Example code |
---|---|
Mass imputation based on regression imputation |
|
Inverse probability weighting |
|
Inverse probability weighting with calibration constraint |
|
Doubly robust estimator |
|
Estimator | Example code |
---|---|
Mass imputation based on regression imputation |
|
Mass imputation based on nearest neighbour imputation |
|
Mass imputation based on predictive mean matching |
|
Mass imputation based on regression imputation with variable selection (LASSO) |
|
Inverse probability weighting |
|
Inverse probability weighting with calibration constraint |
|
Inverse probability weighting with calibration constraint with variable selection (SCAD) |
|
Doubly robust estimator |
|
Doubly robust estimator with variable selection (SCAD) and bias minimization |
|
Simulate example data from the following paper: Kim, Jae Kwang, and Zhonglei Wang. “Sampling techniques for big data analysis.” International Statistical Review 87 (2019): S177-S191 [section 5.2]
library(survey)
library(nonprobsvy)
set.seed(1234567890)
<- 1e6 ## 1000000
N <- 1000
n <- rnorm(n = N, mean = 1, sd = 1)
x1 <- rexp(n = N, rate = 1)
x2 <- rnorm(n = N) # rnorm(N)
epsilon <- 1 + x1 + x2 + epsilon
y1 <- 0.5*(x1 - 0.5)^2 + x2 + epsilon
y2 <- exp(x2)/(1+exp(x2))
p1 <- exp(-0.5+0.5*(x2-2)^2)/(1+exp(-0.5+0.5*(x2-2)^2))
p2 <- rbinom(n = N, size = 1, prob = p1)
flag_bd1 <- as.numeric(1:N %in% sample(1:N, size = n))
flag_srs <- N/n
base_w_srs <- data.frame(x1,x2,y1,y2,p1,p2,base_w_srs, flag_bd1, flag_srs)
population <- N/sum(population$flag_bd1) base_w_bd
Declare svydesign
object with survey
package
<- svydesign(ids= ~1, weights = ~ base_w_srs,
sample_prob data = subset(population, flag_srs == 1))
Estimate population mean of y1
based on doubly robust
estimator using IPW with calibration constraints.
<- nonprob(
result_dr selection = ~ x2,
outcome = y1 ~ x1 + x2,
data = subset(population, flag_bd1 == 1),
svydesign = sample_prob
)
Results
summary(result_dr)
#>
#> Call:
#> nonprob(data = subset(population, flag_bd1 == 1), selection = ~x2,
#> outcome = y1 ~ x1 + x2, svydesign = sample_prob)
#>
#> -------------------------
#> Estimated population mean: 2.95 with overall std.err of: 0.04195
#> And std.err for nonprobability and probability samples being respectively:
#> 0.000783 and 0.04195
#>
#> 95% Confidence inverval for popualtion mean:
#> lower_bound upper_bound
#> y1 2.867789 3.03224
#>
#>
#> Based on: Doubly-Robust method
#> For a population of estimate size: 1025063
#> Obtained on a nonprobability sample of size: 693011
#> With an auxiliary probability sample of size: 1000
#> -------------------------
#>
#> Regression coefficients:
#> -----------------------
#> For glm regression on outcome variable:
#> Estimate Std. Error z value P(>|z|)
#> (Intercept) 0.996282 0.002139 465.8 <2e-16 ***
#> x1 1.001931 0.001200 835.3 <2e-16 ***
#> x2 0.999125 0.001098 910.2 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -----------------------
#> For glm regression on selection variable:
#> Estimate Std. Error z value P(>|z|)
#> (Intercept) -0.498997 0.003702 -134.8 <2e-16 ***
#> x2 1.885629 0.005303 355.6 <2e-16 ***
#> -------------------------
#>
#> Weights:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.000 1.071 1.313 1.479 1.798 2.647
#> -------------------------
#>
#> Covariate balance:
#> (Intercept) x2
#> 25062.8473 -517.5862
#> -------------------------
#>
#> Residuals:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.99999 0.06603 0.23778 0.26046 0.44358 0.62222
#>
#> AIC: 1010622
#> BIC: 1010645
#> Log-Likelihood: -505309 on 694009 Degrees of freedom
Mass imputation estimator
<- nonprob(
result_mi outcome = y1 ~ x1 + x2,
data = subset(population, flag_bd1 == 1),
svydesign = sample_prob
)
Results
summary(result_mi)
#>
#> Call:
#> nonprob(data = subset(population, flag_bd1 == 1), outcome = y1 ~
#> x1 + x2, svydesign = sample_prob)
#>
#> -------------------------
#> Estimated population mean: 2.95 with overall std.err of: 0.04203
#> And std.err for nonprobability and probability samples being respectively:
#> 0.001227 and 0.04201
#>
#> 95% Confidence inverval for popualtion mean:
#> lower_bound upper_bound
#> y1 2.867433 3.032186
#>
#>
#> Based on: Mass Imputation method
#> For a population of estimate size: 1e+06
#> Obtained on a nonprobability sample of size: 693011
#> With an auxiliary probability sample of size: 1000
#> -------------------------
#>
#> Regression coefficients:
#> -----------------------
#> For glm regression on outcome variable:
#> Estimate Std. Error z value P(>|z|)
#> (Intercept) 0.996282 0.002139 465.8 <2e-16 ***
#> x1 1.001931 0.001200 835.3 <2e-16 ***
#> x2 0.999125 0.001098 910.2 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> -------------------------
Inverse probability weighting estimator
<- nonprob(
result_ipw selection = ~ x2,
target = ~y1,
data = subset(population, flag_bd1 == 1),
svydesign = sample_prob)
Results
summary(result_ipw)
#>
#> Call:
#> nonprob(data = subset(population, flag_bd1 == 1), selection = ~x2,
#> target = ~y1, svydesign = sample_prob)
#>
#> -------------------------
#> Estimated population mean: 2.925 with overall std.err of: 0.05
#> And std.err for nonprobability and probability samples being respectively:
#> 0.001586 and 0.04997
#>
#> 95% Confidence inverval for popualtion mean:
#> lower_bound upper_bound
#> y1 2.82679 3.022776
#>
#>
#> Based on: Inverse probability weighted method
#> For a population of estimate size: 1025063
#> Obtained on a nonprobability sample of size: 693011
#> With an auxiliary probability sample of size: 1000
#> -------------------------
#>
#> Regression coefficients:
#> -----------------------
#> For glm regression on selection variable:
#> Estimate Std. Error z value P(>|z|)
#> (Intercept) -0.498997 0.003702 -134.8 <2e-16 ***
#> x2 1.885629 0.005303 355.6 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> -------------------------
#>
#> Weights:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.000 1.071 1.313 1.479 1.798 2.647
#> -------------------------
#>
#> Covariate balance:
#> (Intercept) x2
#> 25062.8473 -517.5862
#> -------------------------
#>
#> Residuals:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.99999 0.06603 0.23778 0.26046 0.44358 0.62222
#>
#> AIC: 1010622
#> BIC: 1010645
#> Log-Likelihood: -505309 on 694009 Degrees of freedom
Work on this package is supported by the National Science Centre, OPUS 20 grant no. 2020/39/B/HS4/00941.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.