Introduction to the nonprobsampling Package
Source:vignettes/nonprobsampling.Rmd
nonprobsampling.RmdOverview
The nonprobsampling package provides methods for
correcting participation bias in nonprobability samples using one or
multiple external probability reference surveys. The package is based on
the generalized estimating equation framework proposed by Landsman et
al. (2026), in which pseudo-weights are constructed by aligning
auxiliary information from the nonprobability sample with corresponding
information from probability reference surveys.
The package implements three one-reference methods within the framework: the raking ratio calibration method of Landsman et al. (2026), the adjusted logistic propensity (ALP) method of Wang et al. (2021), and the Chen–Li–Wu (CLW) method of Chen et al. (2020). It also implements the multi-reference calibration method of Landsman et al. (2026), which integrates auxiliary information from multiple reference surveys when no single reference survey contains all variables relevant to participation.
Datasets
The package includes three datasets derived from real studies and used throughout this vignette to illustrate estimation of the mean serum PSA (prostate-specific antigen) level among U.S. men aged 55 years and older.
| Dataset | Type | Source | n |
|---|---|---|---|
sc |
Nonprobability sample | A synthetic sample drawn from a finite population, which was constructed based on NHANES survey data from the 1999–2010 cycles | 2,404 |
sp1 |
First probability reference survey | National Health and Nutrition Examination Survey (NHANES) 1999–2010 | 3,494 |
sp2 |
Second probability reference survey | National Health Interview Survey (NHIS) 1997–2008 | 35,525 |
The dataset sc represents the nonprobability sample and
contains the outcome variable psa_level, defined as serum
PSA level (ng/mL). The datasets sp1 and sp2
are probability reference surveys used to estimate pseudo-weights.
For illustration purposes, all shared variables have already been
harmonized across the example datasets. Variable names, factor levels,
and category definitions are consistent. The variables shared between
sc and sp1 are agecat,
race, education, BMI,
comorbidity, diabetes and
pros_enlarged. The variables shared between sc
and sp2 are agecat, race,
BMI, comorbidity, and
diabetes.
For user-supplied data, shared variables must be harmonized before estimation. Variable names must match exactly, and continuous variables included in the participation model should be measured on comparable scales across datasets.
library(nonprobsampling)
data("sc")
data("sp1")
data("sp1_bootstrap")
data("sp2")
str(sc)
#> 'data.frame': 2404 obs. of 8 variables:
#> $ psa_level : num 0.3 1 1 1.2 0.3 1 0.3 1 1 0.3 ...
#> $ BMI : Factor w/ 4 levels "Morbidly Obese",..: 1 1 4 3 2 1 2 4 4 1 ...
#> $ race : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
#> $ agecat : Factor w/ 4 levels "1","2","3","4": 1 2 2 4 2 3 2 2 2 1 ...
#> $ education : Factor w/ 5 levels "1","2","3","4",..: 4 3 4 4 1 3 1 4 4 4 ...
#> $ pros_enlarged: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
#> $ comorbidity : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 1 1 2 ...
#> $ diabetes : Factor w/ 2 levels "0","1": 2 2 1 1 1 2 1 1 1 2 ...
str(sp1)
#> 'data.frame': 3494 obs. of 14 variables:
#> $ agecat : Factor w/ 4 levels "1","2","3","4": 4 2 4 3 2 2 1 1 4 3 ...
#> $ marital : Factor w/ 4 levels "1","2","3","4": 1 3 1 1 1 1 1 1 1 1 ...
#> $ race : Factor w/ 4 levels "1","2","3","4": 3 1 3 3 1 1 1 1 1 1 ...
#> $ education : Factor w/ 5 levels "1","2","3","4",..: 1 3 2 3 3 2 2 4 5 2 ...
#> $ employment : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 2 1 1 ...
#> $ smoking : Factor w/ 3 levels "1","2","3": 2 2 2 2 2 2 2 1 3 2 ...
#> $ comorbidity : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 2 1 2 1 ...
#> $ psa_level : num NA NA NA NA NA NA NA NA NA NA ...
#> $ BMI : Factor w/ 4 levels "Morbidly Obese",..: 2 3 4 3 4 3 4 4 2 4 ...
#> $ diabetes : Factor w/ 2 levels "0","1": 2 2 1 2 1 1 2 1 1 1 ...
#> $ pros_enlarged: Factor w/ 2 levels "0","1": 2 1 1 2 1 2 1 1 2 1 ...
#> $ strata_sp1 : int 13 12 6 13 5 12 10 11 1 7 ...
#> $ psu_sp1 : int 2 2 1 2 2 1 1 1 2 1 ...
#> $ wts_sp1 : num 401 8923 481 498 9314 ...
str(sp2)
#> 'data.frame': 35525 obs. of 11 variables:
#> $ agecat : Factor w/ 4 levels "1","2","3","4": 3 1 3 1 4 4 1 1 1 1 ...
#> $ marital : Factor w/ 4 levels "1","2","3","4": 2 1 4 1 1 2 2 3 1 1 ...
#> $ race : Factor w/ 4 levels "1","2","3","4": 3 3 3 3 3 3 2 2 2 3 ...
#> $ employment : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 1 1 1 ...
#> $ diabetes : Factor w/ 2 levels "0","1": 2 2 1 1 1 1 1 2 2 1 ...
#> $ BMI : Factor w/ 4 levels "Morbidly Obese",..: 2 2 4 2 3 4 2 4 1 3 ...
#> $ smoking : Factor w/ 3 levels "1","2","3": 1 1 2 3 2 2 3 3 2 1 ...
#> $ comorbidity: Factor w/ 2 levels "0","1": 2 2 2 1 1 1 2 2 2 2 ...
#> $ wts_sp2 : int 2466 8271 2913 9368 3401 1701 2253 2255 5829 1452 ...
#> $ strata_sp2 : num 5095 5152 5152 5325 5171 ...
#> $ psu_sp2 : int 2 1 1 2 1 1 2 2 2 1 ...
head(sc)
#> psa_level BMI race agecat education pros_enlarged comorbidity
#> 1 0.3 Morbidly Obese 1 1 4 0 1
#> 2 1.0 Morbidly Obese 1 2 3 0 1
#> 3 1.0 Overweight 1 2 4 0 0
#> 4 1.2 Obese 1 4 4 0 1
#> 5 0.3 Normal 1 2 1 0 1
#> 6 1.0 Morbidly Obese 1 3 3 0 1
#> diabetes
#> 1 1
#> 2 1
#> 3 0
#> 4 0
#> 5 0
#> 6 1Step 1 — Create Survey Design Objects
Before calling est_pw(), each probability reference
survey must be converted to a survey design object using either
survey::svydesign() or survey::svrepdesign().
These objects retain the sampling information needed for design-based
variance estimation.
ref1_design <- survey::svydesign(
ids = ~psu_sp1,
strata = ~strata_sp1,
weights = ~wts_sp1,
data = sp1,
nest = TRUE
)
ref2_design <- survey::svydesign(
ids = ~psu_sp2,
strata = ~strata_sp2,
weights = ~wts_sp2,
data = sp2,
nest = TRUE
)The replicate weights in the example reference survey
sp1_bootstrap were generated using the Rao-Wu rescaling
bootstrap method (Rust and Rao, 1996). They are provided to illustrate
variance estimation with a replicate-weight survey design object. In
real applications, users should use the replicate weights supplied by
the survey data provider when available.
ref1_design_rep_wts <- survey::svrepdesign(
data = sp1_bootstrap,
weights = ~wts_sp1,
repweights = "bw[0-9]+",
type = "bootstrap"
)Step 2.1 — Estimate Pseudo-Weights (One Reference Survey)
est_pw() is the main function of the package. It takes a
list in which the first element is the nonprobability sample data frame
and the remaining elements are probability reference surveys represented
as survey design objects.
The method argument specifies the pseudo-weighting
method to use, and p_formula defines the participation
model by selecting the shared auxiliary variables used to estimate the
pseudo-weights.
Calibration Estimator
The raking ratio calibration estimator is a one-reference pseudo-weighting method within the generalized estimating equation framework.
fit_cali <- est_pw(
data = list(sc, ref1_design),
p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
method = "calibration",
control = pw_solver_control(ftol = 1e-6)
)
print(fit_cali)
#> Pseudo-weight fit ("pw_fit")
#>
#> Call:
#> est_pw(data = list(sc, ref1_design), p_formula = ~agecat + race +
#> education + comorbidity + BMI + diabetes, method = "calibration",
#> control = pw_solver_control(ftol = 1e-06))
#>
#> Method: One reference calibration
#> Participation model: 16 parameters (incl. intercept)
#> Convergence: converged (nleqslv, 5 iter, max|EE| = 1.30e-08)
#>
#> Pseudo-weights (n = 2404):
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 4099 7244 8819 8651 10139 13958
#> Sum: 20,797,012
#>
#> Use summary() for coefficients and diagnostics; pwmean() to estimate means.
summary(fit_cali)
#> Call:
#> est_pw(data = list(sc, ref1_design), p_formula = ~agecat + race +
#> education + comorbidity + BMI + diabetes, method = "calibration",
#> control = pw_solver_control(ftol = 1e-06))
#>
#> Method: One reference calibration
#>
#> Participation model involves the following variables:
#> agecat2 agecat3 agecat4 race2 race3 race4 education2 education3 education4 education5 comorbidity1 BMINormal BMIObese BMIOverweight diabetes1
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 5
#> Max |estimating equation|: 1.304e-08
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) agecat2 agecat3 agecat4 race2 race3 race4
#> -9.1480 -0.0344 -0.1093 -0.2328 0.3702 0.6331 0.4028
#> education2 education3 education4 education5 comorbidity1 BMINormal BMIObese
#> -0.1751 -0.1227 -0.0394 0.0951 0.0121 0.0484 0.0858
#> BMIOverweight diabetes1
#> 0.0345 0.0546When no p_formula is provided, it is constructed
automatically from all variables shared between the nonprobability
sample and the reference survey(s). The generated formula is shown in
the model summary.
fit_cali_auto <- est_pw(
data = list(sc, ref1_design),
method = "calibration",
control = pw_solver_control(ftol = 1e-6)
)
summary(fit_cali_auto)
#> Call:
#> est_pw(data = list(sc, ref1_design), method = "calibration",
#> control = pw_solver_control(ftol = 1e-06))
#>
#> Method: One reference calibration
#> Generated default p_formula: ~psa_level + BMI + race + agecat + education + pros_enlarged + comorbidity + diabetes
#>
#> Participation model involves the following variables:
#> psa_level BMINormal BMIObese BMIOverweight race2 race3 race4 agecat2 agecat3 agecat4 education2 education3 education4 education5 pros_enlarged1 comorbidity1 diabetes1
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 7
#> Max |estimating equation|: 1.863e-08
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) psa_level BMINormal BMIObese BMIOverweight race2
#> -8.7089 -0.0800 0.1392 0.0121 0.0087 0.4721
#> race3 race4 agecat2 agecat3 agecat4 education2 education3
#> 0.6846 0.3769 0.0457 0.0226 -0.0967 -0.0844 -0.1381
#> education4 education5 pros_enlarged1 comorbidity1 diabetes1
#> -0.0412 -0.0172 0.1179 0.0135 -0.0013
#>
#> (1190 observations deleted due to missingness in sp)Interaction and higher-order terms can be specified through
p_formula.
fit_cali_interaction <- est_pw(
data = list(sc, ref1_design),
p_formula = ~ BMI * psa_level + diabetes + pros_enlarged + I(psa_level^2),
method = "calibration",
control = pw_solver_control(ftol = 1e-6)
)
summary(fit_cali_interaction)
#> Call:
#> est_pw(data = list(sc, ref1_design), p_formula = ~BMI * psa_level +
#> diabetes + pros_enlarged + I(psa_level^2), method = "calibration",
#> control = pw_solver_control(ftol = 1e-06))
#>
#> Method: One reference calibration
#>
#> Participation model involves the following variables:
#> BMINormal BMIObese BMIOverweight psa_level diabetes1 pros_enlarged1 I(psa_level^2) BMINormal:psa_level BMIObese:psa_level BMIOverweight:psa_level
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 12
#> Max |estimating equation|: 6.519e-09
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) BMINormal BMIObese BMIOverweight psa_level diabetes1
#> -8.7225 0.1445 0.1695 0.1073 -0.0369 0.0861
#> pros_enlarged1 I(psa_level^2) BMINormal:psa_level BMIObese:psa_level
#> 0.0330 -0.0019 0.0360 -0.0839
#> BMIOverweight:psa_level
#> -0.0291
#>
#> (1190 observations deleted due to missingness in sp)Setting verbose = TRUE prints a brief processing log,
including the number of observations used after missing-data filtering,
the estimation method, participation model dimension, and solver
convergence diagnostics such as iteration count and the final maximum
absolute estimating equation value.
fit_cali <- est_pw(
data = list(sc, ref1_design),
p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
method = "calibration",
verbose = TRUE,
control = pw_solver_control(ftol = 1e-6)
)
#> sc: 2404 rows, no exclusions.
#> sp: 3494 rows, no exclusions.
#> Fitting Calibration participation model with 16 variables.
#> Converged after 5 iterations (max |eq| = 1.304e-08).The updated sc data frame is returned with a
pseudo_wts column added.
head(fit_cali$sc_updated)
#> psa_level BMI race agecat education pros_enlarged comorbidity
#> 1 0.3 Morbidly Obese 1 1 4 0 1
#> 2 1.0 Morbidly Obese 1 2 3 0 1
#> 3 1.0 Overweight 1 2 4 0 0
#> 4 1.2 Obese 1 4 4 0 1
#> 5 0.3 Normal 1 2 1 0 1
#> 6 1.0 Morbidly Obese 1 3 3 0 1
#> diabetes pseudo_wts
#> 1 1 9142.734
#> 2 1 10285.206
#> 3 0 9772.742
#> 4 0 11185.397
#> 5 0 9154.969
#> 6 1 11084.897The estimated pseudo-weights are also returned as a numeric vector.
head(fit_cali$pseudo_weights)
#> [1] 9142.734 10285.206 9772.742 11185.397 9154.969 11084.897ALP Estimator
The adjusted logistic propensity (ALP) estimator is another one-reference pseudo-weighting method within the generalized estimating equation framework.
fit_alp <- est_pw(
data = list(sc, ref1_design_rep_wts),
p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
method = "alp"
)
print(fit_alp)
#> Pseudo-weight fit ("pw_fit")
#>
#> Call:
#> est_pw(data = list(sc, ref1_design_rep_wts), p_formula = ~agecat +
#> race + education + comorbidity + BMI + diabetes, method = "alp")
#>
#> Method: One reference ALP
#> Participation model: 16 parameters (incl. intercept)
#> Convergence: converged (nleqslv, 5 iter, max|EE| = 1.59e-12)
#>
#> Pseudo-weights (n = 2404):
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 4056 7079 8772 8659 10229 14067
#> Sum: 20,815,867
#>
#> Use summary() for coefficients and diagnostics; pwmean() to estimate means.CLW Estimator
The Chen–Li–Wu (CLW) estimator is another one-reference pseudo-weighting method within the generalized estimating equation framework.
fit_clw <- est_pw(
data = list(sc, ref1_design),
p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
method = "clw"
)
print(fit_clw)
#> Pseudo-weight fit ("pw_fit")
#>
#> Call:
#> est_pw(data = list(sc, ref1_design), p_formula = ~agecat + race +
#> education + comorbidity + BMI + diabetes, method = "clw")
#>
#> Method: One reference CLW
#> Participation model: 16 parameters (incl. intercept)
#> Convergence: converged (nleqslv, 5 iter, max|EE| = 1.53e-10)
#>
#> Pseudo-weights (n = 2404):
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 4057 7079 8772 8659 10229 14067
#> Sum: 20,815,875
#>
#> Use summary() for coefficients and diagnostics; pwmean() to estimate means.Step 2.2 — Estimate Pseudo-Weights (Multiple Reference Surveys)
When two or more probability reference surveys are available, use
method = "multi" (or leave method = NULL for
automatic selection). The argument p_formula should be
supplied as a list containing one participation model formula per
reference survey.
By default, precali = TRUE applies cumulative
precalibration before pseudo-weight estimation. Under this procedure,
the reference surveys are sequentially calibrated to previously
processed surveys using their overlapping variables in the
p_formula. This helps preserve information on the
relationships between overlapping and survey-specific auxiliary
variables within each reference survey.
The argument sp_order controls the order in which
reference surveys are processed during cumulative precalibration and
pseudo-weight estimation. With sp_order = "size" (default),
surveys are ordered from largest to smallest sample size, so larger
surveys are used as reference and overlapping variables are
preferentially taken from them. With sp_order = "given",
surveys are processed in the user-supplied order, which may be
preferable when a particular reference survey is considered more
reliable or important.
Order Surveys by Sample Size (Default)
fit_multi_order_by_size <- est_pw(
data = list(sc, ref1_design, ref2_design),
p_formula = list(
~ agecat + race + education + psa_level + pros_enlarged + comorbidity + BMI + diabetes,
~ agecat + race + diabetes + BMI + comorbidity
),
sp_order = "size",
precali = TRUE,
control = pw_solver_control(ftol=1e-6)
)
summary(fit_multi_order_by_size)
#> Call:
#> est_pw(data = list(sc, ref1_design, ref2_design), sp_order = "size",
#> precali = TRUE, p_formula = list(~agecat + race + education +
#> psa_level + pros_enlarged + comorbidity + BMI + diabetes,
#> ~agecat + race + diabetes + BMI + comorbidity), control = pw_solver_control(ftol = 1e-06))
#>
#> Method: Multi-reference calibration
#>
#> Precalibration summary:
#> Non-calibrated sample (largest): sp[[2]]
#> Calibrated sample: sp[[1]]
#> Calibration variables: survey weights total, agecat2, agecat3, agecat4, BMINormal, BMIObese, BMIOverweight, comorbidity1, diabetes1, race2, race3, race4
#>
#> Reference samples summary:
#> Order of samples by size (largest to smallest):
#> sp[[2]] (n = 35525)
#> sp[[1]] (n = 2304)
#> Shared variables in sp[[2]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, diabetes1, BMINormal, BMIObese, BMIOverweight, comorbidity1
#> Shared variables in sp[[1]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, education2, education3, education4, education5, psa_level, pros_enlarged1, comorbidity1, BMINormal, BMIObese, BMIOverweight, diabetes1
#> Variables used for calculation in sp[[2]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, diabetes1, BMINormal, BMIObese, BMIOverweight, comorbidity1
#> Variables used for calculation in sp[[1]]:
#> education2, education3, education4, education5, psa_level, pros_enlarged1
#>
#> Participation model involves the following variables:
#> agecat2 agecat3 agecat4 race2 race3 race4 diabetes1 BMINormal BMIObese BMIOverweight comorbidity1 education2 education3 education4 education5 psa_level pros_enlarged1
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 7
#> Max |estimating equation|: 4.768e-07
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) agecat2 agecat3 agecat4 race2 race3 race4
#> -11.1368 0.0493 0.0508 -0.1023 0.4634 0.7063 0.6381
#> diabetes1 BMINormal BMIObese BMIOverweight comorbidity1 education2
#> -0.0666 -0.3844 -0.1567 -0.6146 0.0394 -0.0858
#> education3 education4 education5 psa_level pros_enlarged1
#> -0.1519 -0.0425 -0.0212 -0.0849 0.0883
#>
#> (1190 observations deleted due to missingness in sp[[1]])Use the User-Specified Survey Order
fit_multi_as_given <- est_pw(
data = list(sc, ref1_design, ref2_design),
p_formula = list(
~ agecat + race + education + psa_level + pros_enlarged + comorbidity + diabetes,
~ agecat + race + BMI + diabetes + comorbidity
),
sp_order = "given",
precali = TRUE
)
summary(fit_multi_as_given)
#> Call:
#> est_pw(data = list(sc, ref1_design, ref2_design), sp_order = "given",
#> precali = TRUE, p_formula = list(~agecat + race + education +
#> psa_level + pros_enlarged + comorbidity + diabetes, ~agecat +
#> race + BMI + diabetes + comorbidity))
#>
#> Method: Multi-reference calibration
#>
#> Precalibration summary:
#> Non-calibrated sample (first): sp[[1]]
#> Calibrated sample: sp[[2]]
#> Calibration variables: survey weights total, agecat2, agecat3, agecat4, comorbidity1, diabetes1, race2, race3, race4
#>
#> Reference samples summary:
#> Order of samples kept as provided:
#> sp[[1]] (n = 2304)
#> sp[[2]] (n = 35525)
#> Shared variables in sp[[1]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, education2, education3, education4, education5, psa_level, pros_enlarged1, comorbidity1, diabetes1
#> Shared variables in sp[[2]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, BMINormal, BMIObese, BMIOverweight, diabetes1, comorbidity1
#> Variables used for calculation in sp[[1]]:
#> agecat2, agecat3, agecat4, race2, race3, race4, education2, education3, education4, education5, psa_level, pros_enlarged1, comorbidity1, diabetes1
#> Variables used for calculation in sp[[2]]:
#> BMINormal, BMIObese, BMIOverweight
#>
#> Participation model involves the following variables:
#> agecat2 agecat3 agecat4 race2 race3 race4 education2 education3 education4 education5 psa_level pros_enlarged1 comorbidity1 diabetes1 BMINormal BMIObese BMIOverweight
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 7
#> Max |estimating equation|: 3.725e-09
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) agecat2 agecat3 agecat4 race2 race3 race4
#> -8.2265 0.0475 0.0593 -0.0428 0.5086 0.7317 0.4247
#> education2 education3 education4 education5 psa_level pros_enlarged1
#> -0.0908 -0.1519 -0.0670 -0.0045 -0.0799 0.0985
#> comorbidity1 diabetes1 BMINormal BMIObese BMIOverweight
#> -0.0263 -0.1061 -0.3840 -0.1522 -0.6174
#>
#> (1190 observations deleted due to missingness in sp[[1]])Step 3 — Estimate Population Means and Prevalences
Once pseudo-weights have been estimated, pwmean() uses
them to estimate population means for numeric outcomes and population
prevalences of each observed category for categorical outcomes. In both
cases, the returned object includes unweighted and pseudo-weighted
estimates, standard errors, and confidence intervals.
For a numeric outcome, pwmean() estimates the overall
population mean. Since psa_level is numeric, the estimate
is reported as a mean. Binary 0/1 outcomes are also treated as numeric
and reported as means.
est_cali_numeric <- pwmean(fit_cali, y = "psa_level")
print(est_cali_numeric)
#>
#> Pseudo-weighted (calibration) Estimators:
#> Domain: Overall
#> Mean: 1.541785
#> Std. Error: 0.035178
#> 95% CI: [1.472837, 1.610733]
summary(est_cali_numeric)
#> Call:
#> pwmean(object = fit_cali, y = "psa_level")
#>
#> Method: One reference calibration
#>
#> Domain: Overall
#>
#> Unweighted estimators:
#> Mean: 1.536194
#> Std. Error: 0.034852
#> 95% CI: [1.467886, 1.604502]
#>
#> Pseudo-weighted (calibration) estimators:
#> Mean: 1.541785
#> Std. Error: 0.035178
#> 95% CI: [1.472837, 1.610733]For a categorical outcome, pwmean() estimates the
population prevalence of each observed category. Here, BMI
is a categorical outcome, so the output is printed with category and
prevalence labels.
est_cali_factor <- pwmean(fit_cali, y = "BMI")
print(est_cali_factor)
#>
#> Pseudo-weighted (calibration) Estimators:
#> Category: BMI = Morbidly Obese
#> Domain: Overall
#> Prevalence: 0.087922
#> Std. Error: 0.006216
#> 95% CI: [0.075739, 0.100104]
#>
#> Category: BMI = Normal
#> Domain: Overall
#> Prevalence: 0.248901
#> Std. Error: 0.009127
#> 95% CI: [0.231012, 0.266790]
#>
#> Category: BMI = Obese
#> Domain: Overall
#> Prevalence: 0.231551
#> Std. Error: 0.010053
#> 95% CI: [0.211848, 0.251255]
#>
#> Category: BMI = Overweight
#> Domain: Overall
#> Prevalence: 0.431626
#> Std. Error: 0.010305
#> 95% CI: [0.411429, 0.451823]
summary(est_cali_factor)
#> Call:
#> pwmean(object = fit_cali, y = "BMI")
#>
#> Method: One reference calibration
#>
#> Unweighted estimators:
#> category domain prevalence se CI_lower CI_upper
#> BMI = Morbidly Obese Overall 0.08111481 0.005569341 0.0701991 0.09203052
#> BMI = Normal Overall 0.25707155 0.008915047 0.2395984 0.27454472
#> BMI = Obese Overall 0.23835275 0.008691808 0.2213171 0.25538838
#> BMI = Overweight Overall 0.42346090 0.010079620 0.4037052 0.44321659
#>
#> Pseudo-weighted (calibration) estimators:
#> category domain prevalence se CI_lower CI_upper
#> BMI = Morbidly Obese Overall 0.0879218 0.006215612 0.07573942 0.1001042
#> BMI = Normal Overall 0.2489010 0.009127032 0.23101234 0.2667896
#> BMI = Obese Overall 0.2315515 0.010053047 0.21184786 0.2512551
#> BMI = Overweight Overall 0.4316257 0.010304846 0.41142861 0.4518229Domain (Subgroup) Estimates
When a categorical domain variable is passed to zcol,
the function returns separate estimates for each domain.
With a numeric outcome and a domain variable, pwmean()
estimates a separate population mean within each domain. In this
example, the mean of psa_level is estimated within each
BMI domain.
est_by_bmi_numeric <- pwmean(fit_cali, y = "psa_level", zcol = "BMI")
print(est_by_bmi_numeric)
#>
#> Pseudo-weighted (calibration) Estimators:
#> Domain: BMI = Morbidly Obese
#> Mean: 1.732436
#> Std. Error: 0.140720
#> 95% CI: [1.456629, 2.008243]
#>
#> Domain: BMI = Normal
#> Mean: 1.728680
#> Std. Error: 0.083032
#> 95% CI: [1.565940, 1.891420]
#>
#> Domain: BMI = Obese
#> Mean: 1.362181
#> Std. Error: 0.057429
#> 95% CI: [1.249623, 1.474740]
#>
#> Domain: BMI = Overweight
#> Mean: 1.491525
#> Std. Error: 0.049834
#> 95% CI: [1.393851, 1.589199]
summary(est_by_bmi_numeric)
#> Call:
#> pwmean(object = fit_cali, y = "psa_level", zcol = "BMI")
#>
#> Method: One reference calibration
#>
#> Unweighted estimators:
#> domain mean se CI_lower CI_upper
#> BMI = Morbidly Obese 1.673128 0.12901739 1.420259 1.925998
#> BMI = Normal 1.692864 0.08216978 1.531814 1.853914
#> BMI = Obese 1.360960 0.05707777 1.249089 1.472830
#> BMI = Overweight 1.513487 0.05114169 1.413251 1.613723
#>
#> Pseudo-weighted (calibration) estimators:
#> domain mean se CI_lower CI_upper
#> BMI = Morbidly Obese 1.732436 0.14072032 1.456629 2.008243
#> BMI = Normal 1.728680 0.08303227 1.565940 1.891420
#> BMI = Obese 1.362181 0.05742881 1.249623 1.474740
#> BMI = Overweight 1.491525 0.04983443 1.393851 1.589199With a categorical outcome and a domain variable,
pwmean() estimates a prevalence of each outcome category
within each domain. In this example, the prevalence of each
agecat category is estimated within each BMI
domain.
est_by_bmi_factor <- pwmean(fit_cali, y = "agecat", zcol = "BMI")
print(est_by_bmi_factor)
#>
#> Pseudo-weighted (calibration) Estimators:
#> Category: agecat = 1
#> Domain: BMI = Morbidly Obese
#> Prevalence: 0.282210
#> Std. Error: 0.032150
#> 95% CI: [0.219196, 0.345223]
#>
#> Category: agecat = 1
#> Domain: BMI = Normal
#> Prevalence: 0.317870
#> Std. Error: 0.019683
#> 95% CI: [0.279293, 0.356447]
#>
#> Category: agecat = 1
#> Domain: BMI = Obese
#> Prevalence: 0.350691
#> Std. Error: 0.020744
#> 95% CI: [0.310033, 0.391349]
#>
#> Category: agecat = 1
#> Domain: BMI = Overweight
#> Prevalence: 0.314797
#> Std. Error: 0.015112
#> 95% CI: [0.285177, 0.344417]
#>
#> Category: agecat = 2
#> Domain: BMI = Morbidly Obese
#> Prevalence: 0.325504
#> Std. Error: 0.034144
#> 95% CI: [0.258582, 0.392425]
#>
#> Category: agecat = 2
#> Domain: BMI = Normal
#> Prevalence: 0.256639
#> Std. Error: 0.018350
#> 95% CI: [0.220675, 0.292604]
#>
#> Category: agecat = 2
#> Domain: BMI = Obese
#> Prevalence: 0.295654
#> Std. Error: 0.019828
#> 95% CI: [0.256792, 0.334515]
#>
#> Category: agecat = 2
#> Domain: BMI = Overweight
#> Prevalence: 0.262580
#> Std. Error: 0.014320
#> 95% CI: [0.234513, 0.290646]
#>
#> Category: agecat = 3
#> Domain: BMI = Morbidly Obese
#> Prevalence: 0.212538
#> Std. Error: 0.030268
#> 95% CI: [0.153214, 0.271862]
#>
#> Category: agecat = 3
#> Domain: BMI = Normal
#> Prevalence: 0.244468
#> Std. Error: 0.018005
#> 95% CI: [0.209178, 0.279757]
#>
#> Category: agecat = 3
#> Domain: BMI = Obese
#> Prevalence: 0.211416
#> Std. Error: 0.017904
#> 95% CI: [0.176325, 0.246508]
#>
#> Category: agecat = 3
#> Domain: BMI = Overweight
#> Prevalence: 0.229366
#> Std. Error: 0.013103
#> 95% CI: [0.203685, 0.255047]
#>
#> Category: agecat = 4
#> Domain: BMI = Morbidly Obese
#> Prevalence: 0.179749
#> Std. Error: 0.030531
#> 95% CI: [0.119909, 0.239589]
#>
#> Category: agecat = 4
#> Domain: BMI = Normal
#> Prevalence: 0.181023
#> Std. Error: 0.016755
#> 95% CI: [0.148183, 0.213863]
#>
#> Category: agecat = 4
#> Domain: BMI = Obese
#> Prevalence: 0.142239
#> Std. Error: 0.015901
#> 95% CI: [0.111073, 0.173405]
#>
#> Category: agecat = 4
#> Domain: BMI = Overweight
#> Prevalence: 0.193257
#> Std. Error: 0.012796
#> 95% CI: [0.168177, 0.218337]
summary(est_by_bmi_factor)
#> Call:
#> pwmean(object = fit_cali, y = "agecat", zcol = "BMI")
#>
#> Method: One reference calibration
#>
#> Unweighted estimators:
#> category domain prevalence se CI_lower CI_upper
#> agecat = 1 BMI = Morbidly Obese 0.3282051 0.03371244 0.26212996 0.3942803
#> agecat = 1 BMI = Normal 0.3446602 0.01913314 0.30715993 0.3821605
#> agecat = 1 BMI = Obese 0.3891798 0.02038608 0.34922376 0.4291357
#> agecat = 1 BMI = Overweight 0.3467583 0.01492416 0.31750753 0.3760092
#> agecat = 2 BMI = Morbidly Obese 0.3384615 0.03397280 0.27187607 0.4050470
#> agecat = 2 BMI = Normal 0.2734628 0.01794467 0.23829188 0.3086337
#> agecat = 2 BMI = Obese 0.2966841 0.01909960 0.25924958 0.3341187
#> agecat = 2 BMI = Overweight 0.2681729 0.01389157 0.24094592 0.2953999
#> agecat = 3 BMI = Morbidly Obese 0.2000000 0.02871833 0.14371311 0.2562869
#> agecat = 3 BMI = Normal 0.2346278 0.01706018 0.20119050 0.2680652
#> agecat = 3 BMI = Obese 0.1954625 0.01658085 0.16296461 0.2279604
#> agecat = 3 BMI = Overweight 0.2229862 0.01305248 0.19740386 0.2485686
#> agecat = 4 BMI = Morbidly Obese 0.1333333 0.02440588 0.08549868 0.1811680
#> agecat = 4 BMI = Normal 0.1472492 0.01426576 0.11928881 0.1752096
#> agecat = 4 BMI = Obese 0.1186736 0.01352220 0.09217061 0.1451767
#> agecat = 4 BMI = Overweight 0.1620825 0.01155602 0.13943313 0.1847319
#>
#> Pseudo-weighted (calibration) estimators:
#> category domain prevalence se CI_lower CI_upper
#> agecat = 1 BMI = Morbidly Obese 0.2822096 0.03215020 0.2191964 0.3452229
#> agecat = 1 BMI = Normal 0.3178697 0.01968255 0.2792926 0.3564468
#> agecat = 1 BMI = Obese 0.3506910 0.02074425 0.3100330 0.3913490
#> agecat = 1 BMI = Overweight 0.3147969 0.01511246 0.2851770 0.3444168
#> agecat = 2 BMI = Morbidly Obese 0.3255038 0.03414428 0.2585822 0.3924253
#> agecat = 2 BMI = Normal 0.2566394 0.01834973 0.2206746 0.2926042
#> agecat = 2 BMI = Obese 0.2956536 0.01982781 0.2567918 0.3345154
#> agecat = 2 BMI = Overweight 0.2625798 0.01431989 0.2345133 0.2906462
#> agecat = 3 BMI = Morbidly Obese 0.2125376 0.03026785 0.1532137 0.2718615
#> agecat = 3 BMI = Normal 0.2444675 0.01800519 0.2091780 0.2797571
#> agecat = 3 BMI = Obese 0.2114163 0.01790418 0.1763248 0.2465079
#> agecat = 3 BMI = Overweight 0.2293664 0.01310283 0.2036853 0.2550475
#> agecat = 4 BMI = Morbidly Obese 0.1797490 0.03053137 0.1199086 0.2395894
#> agecat = 4 BMI = Normal 0.1810233 0.01675539 0.1481834 0.2138633
#> agecat = 4 BMI = Obese 0.1422391 0.01590142 0.1110728 0.1734053
#> agecat = 4 BMI = Overweight 0.1932569 0.01279605 0.1681771 0.2183367Handling Missing Data
nonprobsampling handles missing data in two separate
stages.
First, est_pw() handles missing values in the
participation model variables specified by p_formula. By
default, est_pw() uses
na.action = stats::na.omit, which removes observations with
missing values before estimating pseudo-weights.
Alternatively, users may specify
na.action = stats::na.exclude. In this case, excluded
observations from the nonprobability sample are retained in the returned
objects, but their pseudo-weights are set to NA. This is
useful when users want to keep the original row alignment of the
nonprobability sample for later merging or imputation.
Second, pwmean() handles missing values in the outcome
variable and, if specified, the categorical domain variable. This second
layer of missing-data handling is applied after pseudo-weights have
already been estimated.
For illustration, we create temporary versions of the example
datasets. In the nonprobability sample, we insert three missing values
into the participation variable agecat and four missing
values into the outcome variable psa_level. In the
reference survey sp1, we insert two missing values into the
participation variable education.
sc_with_na <- sc
sc_with_na$agecat[c(1, 2, 4)] <- NA
sc_with_na$psa_level[c(6, 7, 8, 9)] <- NA
sp1_with_na <- sp1
sp1_with_na$education[c(3, 4)] <- NA
ref1_design_with_na <- survey::svydesign(
ids = ~psu_sp1,
strata = ~strata_sp1,
weights = ~wts_sp1,
data = sp1_with_na,
nest = TRUE
)Missing-Data Handling in est_pw() with na.exclude
The following example uses the calibration method with
na.action = stats::na.exclude. A brief summary of
missing-data handling is displayed at the bottom of the model
summary.
fit_na_exclude <- est_pw(
data = list(sc_with_na, ref1_design_with_na),
method = "calibration",
p_formula = ~ agecat + race + education + comorbidity + BMI + diabetes,
na.action = stats::na.exclude,
control = pw_solver_control(ftol = 1e-6)
)
summary(fit_na_exclude)
#> Call:
#> est_pw(data = list(sc_with_na, ref1_design_with_na), p_formula = ~agecat +
#> race + education + comorbidity + BMI + diabetes, method = "calibration",
#> na.action = stats::na.exclude, control = pw_solver_control(ftol = 1e-06))
#>
#> Method: One reference calibration
#>
#> Participation model involves the following variables:
#> agecat2 agecat3 agecat4 race2 race3 race4 education2 education3 education4 education5 comorbidity1 BMINormal BMIObese BMIOverweight diabetes1
#>
#> Solver diagnostics:
#> Solver: nleqslv
#> Method: Newton
#> Termination code: 1
#> Iterations: 5
#> Max |estimating equation|: 1.490e-08
#> Message: Function criterion near zero
#>
#> Participation model coefficients:
#> (Intercept) agecat2 agecat3 agecat4 race2 race3 race4
#> -9.1557 -0.0345 -0.1075 -0.2338 0.3725 0.6348 0.4038
#> education2 education3 education4 education5 comorbidity1 BMINormal BMIObese
#> -0.1751 -0.1244 -0.0437 0.0944 0.0112 0.0576 0.0940
#> BMIOverweight diabetes1
#> 0.0442 0.0513
#>
#> (3 observations deleted due to missingness in sc)
#> (2 observations deleted due to missingness in sp)A summary is available in the na_summary component.
fit_na_exclude$na_summary
#> NA processing summary:
#> n_orig n_used n_excluded
#> sc 2404 2401 3
#> sp 3494 3492 2Information about observations excluded because of missing values in
the participation model variables can be obtained using the S3
na.action() method.
na.action(fit_na_exclude)
#> [1] 1 2 4
#> attr(,"class")
#> [1] "exclude"With na.exclude, the returned pseudo-weight vector and
sc_updated data frame preserve the original row structure.
Observations excluded from pseudo-weight estimation remain in the
output, with their pseudo-weights recorded as NA.
head(fit_na_exclude$pseudo_weights)
#> [1] NA NA 9796.551 NA 9149.513 11215.230
head(fit_na_exclude$sc_updated)
#> psa_level BMI race agecat education pros_enlarged comorbidity
#> 1 0.3 Morbidly Obese 1 <NA> 4 0 1
#> 2 1.0 Morbidly Obese 1 <NA> 3 0 1
#> 3 1.0 Overweight 1 2 4 0 0
#> 4 1.2 Obese 1 <NA> 4 0 1
#> 5 0.3 Normal 1 2 1 0 1
#> 6 NA Morbidly Obese 1 3 3 0 1
#> diabetes pseudo_wts
#> 1 1 NA
#> 2 1 NA
#> 3 0 9796.551
#> 4 0 NA
#> 5 0 9149.513
#> 6 1 11215.230Missing-Data Handling in pwmean()
Missing values in the outcome variable are handled when
pwmean() is called. This is separate from the missing-value
handling performed by est_pw().
result <- pwmean(fit_na_exclude, "psa_level")
summary(result)
#> Call:
#> pwmean(object = fit_na_exclude, y = "psa_level")
#>
#> Method: One reference calibration
#>
#> Domain: Overall
#>
#> Unweighted estimators:
#> Mean: 1.537379
#> Std. Error: 0.034904
#> 95% CI: [1.468969, 1.605790]
#>
#> Pseudo-weighted (calibration) estimators:
#> Mean: 1.544108
#> Std. Error: 0.035273
#> 95% CI: [1.474975, 1.613241]
#>
#> (4 observations deleted due to missingness)A summary of outcome-level missing-value handling is stored in the
returned object. The corresponding na.action information
can also be obtained using the S3 na.action() method.
result$na.action
#> [1] 6 7 8 9
#> attr(,"class")
#> [1] "omit"
na.action(result)
#> [1] 6 7 8 9
#> attr(,"class")
#> [1] "omit"Missing-Data Handling in est_pw() with na.omit (Default)
Using na.action = stats::na.omit removes observations
with missing participation model variables from the returned objects.
This is the default behavior.
fit_na_omit <- est_pw(
data = list(sc_with_na, ref1_design_with_na),
method = "calibration",
na.action = stats::na.omit,
control = pw_solver_control(ftol = 1e-6)
)Information about the excluded observations can be obtained using the
S3 na.action() method.
na.action(fit_na_omit)
#> [1] 1 2 4 6 7 8 9
#> attr(,"class")
#> [1] "omit"With na.omit, the returned pseudo-weight vector and
sc_updated contain only the observations used for
pseudo-weight estimation.
head(fit_na_omit$pseudo_weights)
#> [1] 6555.138 5119.046 6552.331 7517.718 5119.046 7878.904
head(fit_na_omit$sc_updated)
#> psa_level BMI race agecat education pros_enlarged comorbidity
#> 3 1.0 Overweight 1 2 4 0 0
#> 5 0.3 Normal 1 2 1 0 1
#> 10 0.3 Morbidly Obese 1 1 4 0 1
#> 11 1.2 Obese 1 4 4 0 1
#> 12 0.3 Normal 1 2 1 0 1
#> 13 3.2 Overweight 1 1 5 0 0
#> diabetes pseudo_wts
#> 3 0 6555.138
#> 5 0 5119.046
#> 10 1 6552.331
#> 11 0 7517.718
#> 12 0 5119.046
#> 13 0 7878.904In summary, est_pw() handles missing values in
participation model variables, while pwmean() handles
missing values in the outcome and domain variables. With
na.omit, excluded observations are removed from the
returned objects; with na.exclude, the original row
structure is preserved and excluded observations receive NA
pseudo-weights.
Controlling the Numerical Solver
The control argument accepts an object created by
pw_solver_control(), which is passed to the underlying
nleqslv nonlinear equation solver. It controls the
numerical algorithm used to solve the estimating equations, including
convergence tolerances, iteration limits, and diagnostic output. The
default settings are suitable for most applications, but advanced users
may adjust them when troubleshooting convergence issues.
ctrl <- pw_solver_control(
method = "Broyden",
xtol = 1e-10,
ftol = 1e-4,
maxit = 30,
trace = TRUE
)
fit_ctrl <- est_pw(
data = list(sc, ref1_design),
p_formula = ~ agecat + race + education + comorbidity + BMI,
method = "calibration",
control = ctrl
)
#> Algorithm parameters
#> --------------------
#> Method: Broyden Global strategy: double dogleg (initial trust region = -2)
#> Maximum stepsize = 1.79769e+308
#> Scaling: fixed
#> ftol = 0.0001 xtol = 1e-10 btol = 0.001 cndtol = 1e-12
#>
#> Iteration report
#> ----------------
#> Iter Jac Lambda Eta Dlt0 Dltn Fnorm Largest |f|
#> 0 1.344466e+12 1.015181e+06
#> 1 N(4.6e-03) N 0.3165 0.8407 0.8407 3.697989e+11 5.861171e+05
#> 2 B(4.6e-03) N 0.9663 0.0769 0.1538 3.161761e+09 5.858978e+04
#> 3 B(4.6e-03) N 0.2224 0.0508 0.0508 1.669546e+09 3.711881e+04
#> 4 B(4.6e-03) N 0.8740 0.0072 0.0144 2.958202e+06 1.158335e+03
#> 5 B(4.5e-03) N 0.6003 0.0008 0.0016 1.232564e+05 3.730882e+02
#> 6 B(4.5e-03) N 0.2706 0.0003 0.0007 7.010684e+02 2.222404e+01
#> 7 B(4.5e-03) N 0.2724 0.0000 0.0000 2.802099e+01 3.941451e+00
#> 8 B(4.6e-03) N 0.4647 0.0000 0.0000 8.489106e-01 1.015415e+00
#> 9 B(4.5e-03) N 0.2540 0.0000 0.0000 1.891429e-03 4.723867e-02
#> 10 B(4.5e-03) N 0.2325 0.0000 0.0000 1.414529e-06 1.055873e-03
#> 11 B(4.5e-03) N 0.2497 0.0000 0.0000 4.135542e-09 5.835295e-05Quick Reference
| Function | Purpose |
|---|---|
est_pw() |
Estimate pseudo-weights from a nonprobability sample and one or multiple reference surveys |
pwmean() |
Estimate pseudo-weighted population means for numeric outcomes and population prevalences for categorical outcomes, overall or by domain |
pw_solver_control() |
Configure numerical solver settings |
| Argument | Values | Notes |
|---|---|---|
method |
"alp", "clw", "calibration",
"multi", NULL
|
NULL auto-selects based on number of reference
surveys |
sp_order |
"size", "given"
|
Only relevant for multi-reference calibration method;
"size" sorts largest first |
precali |
TRUE, FALSE
|
Precalibration before multi-reference calibration step |
na.action |
na.omit, na.exclude,
na.fail
|
Controls missing-data handling |
References
Chen, Y., Li, P., and Wu, C. (2020). Doubly robust inference with nonprobability survey samples. Journal of the American Statistical Association, 115(532), 2011–2021. doi:10.1080/01621459.2019.1677241
Landsman, V., Wang, L., Carrillo-Garcia, I., Mitani, A. A., Smith, P. M., Graubard, B. I., Bui, T., and Carnide, N. (2026). Correction for Participation Bias in Nonprobability Samples Using Multiple Reference Surveys. Statistics in Medicine, 45(3–5). doi:10.1002/sim.70403
Rust, K. F., and Rao, J. N. K. (1996). Variance estimation for complex surveys using replication techniques. Statistical Methods in Medical Research, 5(3), 283–310. doi:10.1177/096228029600500305
Wang, L., Valliant, R., and Li, Y. (2021). Adjusted logistic propensity weighting methods for population inference using nonprobability volunteer-based epidemiologic cohorts. Statistics in Medicine, 40(24), 5237–5250. doi:10.1002/sim.9122