The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.
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.
Health stats visible at Monitor.