Declarative, pipeable survey weighting in base R: from design weights to calibrated, model-assisted, variance-ready weights.
weightflow builds survey weights by chaining
hierarchical adjustments with a tidymodels-style API, and
estimates their variances with a bootstrap that re-applies the whole
recipe on each replicate. It has no hard dependencies
(base R, R >= 4.1) and bridges to
survey/srvyr for design-based inference.
Where does it fit? survey and srvyr are the
standard tools for analysing data once you already have
weights. weightflow sits one step earlier: it builds those
weights from the design base weights, making every adjustment
(eligibility, nonresponse, calibration, trimming) an explicit, auditable
step, and then hands the result to
survey/srvyr for inference.
weightflow expresses the whole weighting process as a sequence of explicit steps. The diagram below summarizes the flow and the choices that depend on the design and on the available auxiliary information.
# From CRAN
install.packages("weightflow")
# Development version (latest changes)
# install.packages("remotes")
remotes::install_github("jpferreira33/weightflow")A recipe is inert: building it computes nothing.
prep() walks the steps in order and estimates the
cascade of factors; collect_weights() extracts the final
weights. Separating define from apply makes the whole
process reproducible and auditable, and it is exactly what lets the
bootstrap re-run the entire cascade per replicate.
library(weightflow)
recipe <- weighting_spec(sample_one, base_weights = pw) |>
step_unknown_eligibility(unknown = unknown_elig, by = "region") |>
step_drop_ineligible(ineligible = ineligible) |>
# household nonresponse: the whole dwelling is lost (no roster), so the
# adjustment is at the household level and uses only frame information
step_nonresponse(respondent = hh_responded, method = "weighting_class",
by = "region", cluster = "household_id") |>
step_select_within(prob = p_within) |>
# person nonresponse: among the selected persons, the roster gives sex and age
# even for those who did not respond, so a propensity model can use them
step_nonresponse(respondent = responded, method = "propensity",
formula = ~ region + sex + age, engine = "logit",
num_classes = 10) |>
step_calibrate(method = "raking",
margins = list(region = c(table(population$region)),
sex = c(table(population$sex)))) |>
step_trim_weights() |>
step_assert(max_deff = 3)
fitted <- prep(recipe) # estimate the cascade
summary(fitted) # per-stage diagnostics + Kish deff
wts <- collect_weights(fitted) # data.frame with .weightThe article A full weighting pipeline on a real household survey (ECH 2019) runs the whole workflow on open microdata from Uruguay’s continuous household survey: it induces realistic eligibility and nonresponse, weights the survivors back with integrated household calibration, validates the poverty-rate estimate against a known truth, and attaches design-based confidence intervals with the bootstrap.
The methods below are what set weightflow apart. Each is opt-in: the defaults reproduce classic survey weighting, and one argument switches the method on.
Estimate the response propensity with gradient boosting instead of
logistic regression, useful when nonresponse depends on the covariates
in nonlinear or interacting ways. The engine also drives the outcome
models in step_model_calibration().
step_nonresponse(respondent = responded, method = "propensity",
formula = ~ region + sex + age, engine = "boost")A flexible learner that predicts the same units it trained on
overfits the propensity, which inflates the weights and the variance.
Cross-fitting estimates each unit from a model trained on the
other folds; folds are formed by cluster when a
cluster is set, so household members never leak across
folds.
step_nonresponse(respondent = responded, method = "propensity",
formula = ~ region + sex + age, engine = "boost",
crossfit = 5, crossfit_seed = 1)In practice this is the difference between a stable adjustment and one dominated by a few extreme weights: on the bundled data, boosting without cross-fitting inflates the design effect, while cross-fitting brings it back down (the Machine learning, cross-fitting and robust calibration article shows the two side by side).
When you calibrate to many margins, forcing every constraint exactly
can produce extreme weights. Ridge calibration relaxes the targets in a
controlled way: a single, scale-free penalty trades a
little accuracy on the totals for much steadier weights.
step_calibrate(method = "linear", formula = ~ region + sex,
totals = pop_totals, penalty = 1) # smaller = more relaxationInstead of a hand-picked cutoff, choose the trimming threshold that minimizes an estimate of bias^2 + variance (Potter 1990), balancing the bias of trimming against the variance from extreme weights.
step_trim_weights(method = "potter")Hand weightflow the population totals the way they actually arrive,
as a data frame (a census cross-tab, a projection, a spreadsheet),
instead of a fiddly model-matrix vector. Name the counts column with
count; several category columns are crossed automatically,
and weightflow builds the intercept and the dropped reference levels for
you.
region_sex <- as.data.frame(table(region = population$region, sex = population$sex))
step_calibrate(method = "poststratify", totals = region_sex, count = "Freq")Calibrate independently within each domain, each to its own
totals, with one argument (by). The domain is just a column
in the tidy totals, not a term in the formula, and it composes with
calfun, bounds, penalty and the
integrative cluster option.
It earns its keep with a quantitative control total that differs by domain, awkward to express by hand, since it needs domain-by-covariate interactions. Here each region is calibrated to its sex counts and to its own income total:
sex_by_region <- as.data.frame(table(region = population$region, sex = population$sex))
income_by_region <- aggregate(income ~ region, population, sum) # region -> income total
step_calibrate(method = "linear", formula = ~ sex + income,
totals = list(sex = sex_by_region, income = income_by_region),
count = "Freq", by = "region", calfun = "raking")Raking fits the case where, within each region, you know the margins separately (each region’s sex totals and its age-band totals, not their cross):
sex_by_region <- as.data.frame(table(region = population$region, sex = population$sex))
age_by_region <- as.data.frame(table(region = population$region, age_grp = population$age_grp))
step_calibrate(method = "raking",
totals = list(sex_by_region, age_by_region),
count = "Freq", by = "region")A calfun = "raking" distance (g = exp(u)) keeps the
calibrated weights positive without explicit bounds while still hitting
the targets exactly, on categorical and continuous auxiliaries alike,
and with the integrative option.
step_calibrate(method = "linear", formula = ~ region + income,
totals = list(region = m_region, income = 1.2e6),
count = "Freq", calfun = "raking")The control totals of the model-calibration auxiliaries often come
from an outside source (an official figure, a variable not in the
frame). Pass them through x_totals, in the same tidy shape
as linear calibration; population is then used only for the
model predictions.
step_model_calibration(
x_formula = ~ region + age,
models = list(income = y_model(income ~ age + sex, engine = "glm")),
population = population,
x_totals = list(region = m_region, age = 5.1e5), count = "Freq")Alongside the bootstrap, a delete-a-PSU jackknife re-runs the whole recipe on each replicate, so the replicate weights carry every adjustment. Stratified (JKn) or unstratified (JK1), and it bridges to survey/srvyr for any estimand or domain.
jk <- jackknife_weights(spec, strata = "region", psu = "psu")
jack_total(jk, "employed")After a nonresponse adjustment, summary() and
report_weighting() automatically report the R-indicator
(Schouten, Cobben & Bethlehem) plus the partial R-indicators: how
representative the response is, and which variable drives the gap. No
new function to call.
# printed by summary() when the recipe adjusts for nonresponse:
# R-indicator (representativity of response): 0.890 (on region, sex)Adjustment steps, applied in the order you pipe them:
| Step | What it does |
|---|---|
step_unknown_eligibility() |
Redistribute unknown-eligibility cases among the known ones (person-
or household-level via cluster). |
step_drop_ineligible() |
Zero out out-of-scope units. |
step_select_within() |
Within-household selection (unequal prob or equal
n_eligible). |
step_nonresponse() |
Weighting classes or propensity (logit / CART / random forest / xgboost), with optional k-fold cross-fitting, person- or household-level. |
step_calibrate() |
Raking, post-stratification, linear/GREG; bounded (Deville-Särndal), integrative (one weight per household), and ridge (penalized) options. |
step_model_calibration() |
Wu-Sitter model calibration with working models for the outcomes (any engine, with cross-fitting). |
step_trim(), step_trim_weights() |
Manual or automatic trimming (Tukey fence or Potter MSE-optimal), insertable anywhere. |
step_round(), step_rescale() |
Integer rounding and rescaling to a size or total. |
step_assert() |
Quality checkpoint on deff, weight ratio or effective n. |
Eligibility and response accept 0/1 dummy columns or any logical condition.
Diagnostics and reporting: summary()
and plot() show the per-stage cascade with the Kish
design effect (deff = 1 + CV^2) and effective sample size;
weight_factors() returns the per-unit, per-step factors;
report_weighting() writes a self-contained HTML report
(pipeline diagram, variables used, per-stage summaries and per-step
visuals) with no graphics device or server required.
Variance estimation (see the Variance estimation article). Once the weights are built, get design-based standard errors with a bootstrap that re-runs the whole recipe on each replicate:
boot <- bootstrap_weights(recipe, replicates = 500, strata = "region", psu = "psu")
boot_mean(boot, "income") # estimate, SE and 95% CI
# hand the replicate weights to survey / srvyr for the rest of the analysis
rep_design <- as_svrepdesign(boot) # a svyrep.design object
collect_replicate_weights(boot) # replicate weights as a data.frameThe bootstrap resamples PSUs within strata (Rao-Wu rescaling) and
then re-applies the entire cascade (eligibility, nonresponse,
calibration, trimming) on each replicate. So the replicate weights carry
two sources of variability at once: the sampling design
(the resampling of PSUs within strata) and every weighting adjustment
(each one is re-estimated on each replicate). Re-running the full recipe
per replicate is automatic here, rather than something you
re-orchestrate by hand on top of the replicate weights, and the result
plugs straight into survey/srvyr through
as_svrepdesign() for any downstream estimator.
Three bundled datasets: population (the frame),
sample_survey (take-all roster) and sample_one
(multistage select-one design), all with stratum, PSU and design weight,
so the full pipeline and the variance methods run natively.
apply_step() is the internal S3 generic behind each
step. To add an adjustment, define a step_*() constructor
(inert) and its apply_step.<class>() method; nothing
else changes.
General framework
Nonresponse and machine-learning propensities
Calibration
Design effect and trimming
Variance estimation
MIT © Juan Pablo Ferreira