Package {weightflow}

Title:	Declarative Recipes for Staged Survey Weighting with Recipe-Aware Replicate Variances
Version:	1.3.0
Description:	Builds survey analysis weights by declaring the whole weighting process as an ordered recipe of explicit adjustments, estimated in a single call. Steps cover within-cluster selection, subsampling for two-phase designs, nonresponse by weighting classes or response-propensity models (optionally machine-learning, with cross-fitting), calibration to known totals following Deville and Sarndal (1992) <doi:10.2307/2290268>, optionally model-assisted, non-probability samples by pseudo-weighting, mass imputation and doubly robust estimators, and range-restricted trimming. Rotating and pure panels add panel-selection probabilities, attrition, longitudinal weights, gross flows and composite estimation. Variances come from a recipe-aware bootstrap and jackknife that resample or delete primary sampling units and re-apply the entire cascade on each replicate, following Rao and Wu (1988) <doi:10.1080/01621459.1988.10478591>; panel replicates are coordinated across waves, so the sample overlap enters the variance of a net change as covariance, and two-phase variances split into first- and second-phase components (V = V1 + V2). A self-contained HTML report documents each step, and the weights bridge to the 'survey' and 'srvyr' packages. The methods, and the simulation evidence behind the variance estimators, are described in Ferreira (2026) <doi:10.1177/18747655261484262>.
License:	MIT + file LICENSE
Encoding:	UTF-8
Language:	en-US
Depends:	R (≥ 4.1.0)
Imports:	stats, utils, graphics, parallel
Suggests:	MASS, rpart, ranger, testthat (≥ 3.0.0), survey, srvyr, dplyr, tidyr, ggplot2, haven, archive, knitr, rmarkdown, spelling, xgboost, yaml
Config/roxygen2/version:	8.0.0
URL:	https://github.com/jpferreira33/weightflow, https://jpferreira33.github.io/weightflow/
BugReports:	https://github.com/jpferreira33/weightflow/issues
Config/testthat/edition:	3
LazyData:	true
VignetteBuilder:	knitr
NeedsCompilation:	no
Packaged:	2026-09-11 17:41:02 UTC; jp
Author:	Juan Pablo Ferreira [aut, cre, cph], Andrés Gutiérrez [aut]
Maintainer:	Juan Pablo Ferreira <juanpablo.ferreira@fcea.edu.uy>
Repository:	CRAN
Date/Publication:	2026-09-11 20:10:02 UTC

weightflow: declarative survey weighting

Description

Builds analysis weights from design base weights by declaring the weighting process as an ordered recipe of explicit adjustments – unknown eligibility, within-cluster selection (e.g. within household), nonresponse, calibration, trimming, rounding, rescaling, assertions – and then estimating that recipe in one call. The package also produces replicate weights and design-based standard errors that carry the variability of the whole cascade, so a weighting project no longer has to end at the weights.

Details

Start with weighting_spec(), add ⁠step_*()⁠ adjustments, estimate the cascade with prep(), and extract the weights with collect_weights(). Inspect with summary(), plot() and report_weighting().

Author(s)

Maintainer: Juan Pablo Ferreira juanpablo.ferreira@fcea.edu.uy (ORCID) [copyright holder]

Authors:

• Juan Pablo Ferreira juanpablo.ferreira@fcea.edu.uy (ORCID) [copyright holder]

• Andrés Gutiérrez (ORCID)

See Also

Useful links:

• https://github.com/jpferreira33/weightflow

• https://jpferreira33.github.io/weightflow/

• Report bugs at https://github.com/jpferreira33/weightflow/issues

Direct estimates and design SEs per domain, ready for small-area estimation

Description

Small-area estimation (SAE) area-level models (Fay-Herriot) need, for each domain, the direct estimate, its design-based variance and the effective sample size. This function computes exactly those from a recipe-aware replicate object, so the design-based ingredients flow into emdi, sae or hbsae without leaving the weightflow variance machinery. It does not fit any SAE model itself.

Usage

as_sae_input(
  object,
  variable,
  by,
  type = c("mean", "total"),
  level = 0.95,
  cv_breaks = c(0.165, 0.33)
)

Arguments

Details

The domain standard error is the recipe-aware replicate SE (it re-runs the whole recipe per replicate), so it already reflects nonresponse and calibration, not just the final weights.

Value

A data frame with one row per domain: domain, n (active units), n_eff (Kish effective sample size), estimate, se, cv, ci_lower, ci_upper and rating. Pass estimate and se^2 (the sampling variance) to a Fay-Herriot model.

Note for very small domains: the domain estimates share one bootstrap, and a replicate in which any domain has no active unit is dropped for all domains (a conservative choice). With many tiny domains this wastes replicates; use more replicates in the bootstrap, or estimate sparse domains in a separate call.

See Also

domain_summary(), bootstrap_estimate(), design_effect()

Other cascade audit: collect_propensities(), collect_step_detail(), collect_weights(), domain_summary(), weight_factors(), weighting_alerts()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
boot <- bootstrap_weights(spec, replicates = 50, strata = "region",
                          psu = "psu", seed = 1)
as_sae_input(boot, "responded", by = "region")

Export weightflow weights to a survey design

Description

as_svydesign() builds a linearization (ultimate-cluster) survey.design from a prepped recipe, treating the final weights as fixed constants. as_svrepdesign() builds a replicate-weights svyrep.design from a bootstrap_weights() or jackknife_weights() object. Both are the bridge to the survey package, and therefore to svytotal(), svymean(), svyratio(), svyby(), svyglm() and domain estimation generally.

Usage

as_svydesign(object, ids, strata = NULL, weight_name = ".weight", ...)

as_svrepdesign(object, ...)

Arguments

Details

Only as_svrepdesign() propagates the variability of the weighting adjustments (nonresponse, calibration, ...), because each replicate re-runs the whole recipe. as_svydesign() is design-based linearization on the fixed final weights: its standard errors reflect the sampling design but treat the adjustments as known without error, so they are usually smaller. Use as_svrepdesign() (with bootstrap_weights() / jackknife_weights()) when the adjustment variability should be included.

Value

A survey.design / svyrep.design object.

See Also

Other variance estimation: bootstrap_estimate(), bootstrap_weights(), collect_replicate_weights(), jackknife_estimate(), jackknife_weights()

Gross-flow TOTALS with standard errors, plus net flows and margins

Description

The gross change is about the number of people who move between states, i.e. totals. From a longitudinal-weight bootstrap_weights() object this returns, all with a bootstrap standard error computed within each replicate: the from x to count matrix (the gross flows, in population totals), the net flow matrix ⁠i->j minus j->i⁠, and the margins – the origin totals (how many started in each state), the destination totals (how many ended in each), the stayers (the diagonal) and the movers (off-diagonal).

Usage

boot_flows(boot, from, to, states = NULL)

Arguments

Value

a weightflow_flows object with counts/counts_se, net/net_se, origin/origin_se, dest/dest_se, and stayers/movers (each with its SE).

See Also

boot_transition(), transition_matrix()

Transition matrix with per-cell bootstrap standard errors

Description

Like transition_matrix(), but every cell also gets a standard error and confidence interval from the bootstrap replicate weights of a longitudinal-weight bootstrap_weights() object, by re-tabulating the flow within each replicate. The recipe is re-run per replicate, so the uncertainty of the longitudinal weighting propagates into the flow SEs.

Usage

boot_transition(
  boot,
  from,
  to,
  states = NULL,
  format = c("row", "col", "joint", "counts")
)

Arguments

Value

a weightflow_transition_boot with estimate, se, ci_lower, ci_upper matrices.

See Also

transition_matrix(), bootstrap_weights()

Bootstrap estimate, standard error and confidence interval

Description

Applies a statistic to the point weights and to every bootstrap replicate, and returns the estimate with its bootstrap standard error and a normal confidence interval. boot_total() and boot_mean() are the two shortcuts you will use most: a weighted total and a weighted mean of one column.

Usage

bootstrap_estimate(
  boot,
  statistic,
  level = 0.95,
  ci_type = c("normal", "t", "percentile"),
  df = NULL
)

boot_total(boot, variable)

boot_mean(boot, variable)

Arguments

Details

The bootstrap variance takes the replicate estimates \hat\theta^{*}_b around the point estimate \hat\theta (the mse = TRUE convention of survey), over the R valid replicates (a failed replicate is dropped, not counted),

\widehat V_{\mathrm{boot}}(\hat\theta) = \frac{1}{R}\sum_{b=1}^{R}\big(\hat\theta^{*}_b - \hat\theta\big)^2.

Value

A data frame with estimate, se, ci_lower, ci_upper.

See Also

Other variance estimation: as_svydesign(), bootstrap_weights(), collect_replicate_weights(), jackknife_estimate(), jackknife_weights()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
boot <- bootstrap_weights(spec, replicates = 50, strata = "region",
                          psu = "psu", seed = 1)
# a t interval with the design degrees of freedom (safer with few PSUs)
bootstrap_estimate(boot, function(w, d) sum(w * d$responded), ci_type = "t")

Recipe-aware bootstrap replicate weights

Description

Builds bootstrap replicate weights by resampling primary sampling units (PSUs) with replacement within strata and re-running the entire weighting recipe on each replicate – every estimated stage (nonresponse, calibration, model calibration, trimming), not just one. Reach for this when those stages are estimated from the sample and you want their uncertainty inside the standard error, instead of conditioning on them as if they were known.

Usage

bootstrap_weights(
  object,
  replicates = 200L,
  strata = NULL,
  psu = NULL,
  m = NULL,
  fpc = NULL,
  lonely_psu = c("certainty", "collapse"),
  seed = NULL,
  cores = 1L,
  progress = TRUE,
  .tp_component = c("full", "phase1", "phase2")
)

Arguments

Details

The multiplier is the Rao-Wu rescaling bootstrap: within a stratum with n PSUs, m PSUs are drawn with replacement (default m = n - 1) and unit i in PSU k gets \lambda = 1 - \sqrt{m/(n-1)} + \sqrt{m/(n-1)}\,(n/m)\,t_k, with t_k the number of times its PSU was drawn.

Value

An object of class weightflow_boot with the replicates matrix (units x replicates), the point weights, and the design metadata.

See Also

Other variance estimation: as_svydesign(), bootstrap_estimate(), collect_replicate_weights(), jackknife_estimate(), jackknife_weights()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
boot <- bootstrap_weights(spec, replicates = 50, strata = "region",
                          psu = "psu", seed = 1)
boot_total(boot, "responded")
# with a first-stage finite-population correction (per-stratum sampling fraction)
d <- sample_survey; d$f <- 0.1
spec_f <- weighting_spec(d, base_weights = pw)
bootstrap_weights(spec_f, replicates = 50, strata = "region", psu = "psu",
                  fpc = "f", seed = 1)

Net change between two panel waves, with honest variance

Description

Estimates a net change (a statistic in wave 2 minus the same statistic in wave 1) from a wave_bootstrap(), decomposing its variance into V1, V2 and the sampling covariance the overlap induces: V = V1 + V2 - 2*Cov. Because the bootstrap replicates are coordinated by PSU, Cov emerges from them – no analytic covariance formula. Reports the correlation rho and deff_change = V / (V1 + V2), i.e. how much the overlap lowered the variance of the change relative to treating the waves as independent.

Usage

change_estimate(
  wb,
  statistic,
  waves = NULL,
  level = 0.95,
  type = c("absolute", "relative"),
  by = NULL,
  ci_type = c("normal", "t", "percentile"),
  df = NULL
)

change_mean(
  wb,
  variable,
  waves = NULL,
  level = 0.95,
  type = c("absolute", "relative"),
  by = NULL
)

change_total(
  wb,
  variable,
  waves = NULL,
  level = 0.95,
  type = c("absolute", "relative"),
  by = NULL
)

Arguments

Value

an object of class weightflow_change with the estimate, se, V, V1, V2, Vind = V1 + V2, cov, rho, deff_change and the confidence interval; or, with ⁠by=⁠, a weightflow_change_by holding one row per domain.

See Also

wave_bootstrap(), change_mean(), change_total(), level_estimate()

Evaluate an estimation pipeline

Description

Runs every estimand of a weightflow_estimation across the cross of its step_domain() columns, using the coordinated replicates for the variance, and returns a tidy table.

Usage

collect_estimates(est)

Arguments

Value

a weightflow_estimation_result: one row per (estimand x domain cell) with estimate, se, ci_lower, ci_upper and (for changes) rho. The object also carries a row-aligned detail frame (confidence level, effective replicates, the two wave levels behind a change and the overlap design effect) and the pipeline's filters, so report_panel() can document the estimates without re-running them.

See Also

step_estimate(), report_panel()

Recover the fitted response propensities of a nonresponse step

Description

A step_nonresponse(method = "propensity") step fits a response-propensity model and adjusts the weights by 1/\hat p. prep() keeps the full per-unit propensity vector \hat p (out-of-fold when cross-fitting is used) on the step, but it is not returned by collect_weights(). This accessor extracts it aligned to the sample, so you can inspect its distribution and confirm the nonresponse model is well fitted before trusting the adjusted weights. It works the same way whether the adjustment was made at the unit level or, through cluster, at the household level (there the household propensity is broadcast to its members).

Usage

collect_propensities(object, step = NULL)

Arguments

Value

The sample data.frame with columns appended: .propensity (the fitted response propensity \hat p, NA for units outside the model, i.e. ineligible / already dropped), .responded (the response indicator the model used), .weight_in (the weight reaching the step, see below), .factor (the multiplier the step actually applied to the unit), .status (a factor that labels each unit as "eligible respondent", "eligible nonrespondent" or "not in propensity model"), and, when the step uses propensity classes (num_classes), .class (the assigned class). Units not in the propensity model carry NA in the per-unit columns. .weight_in is the weight reaching the nonresponse step – it already carries any earlier adjustment (unknown-eligibility redistribution, within-cluster selection), not the raw base weight. At the unit level it is also the weight the propensity model is fitted with (unless weight_model = FALSE); with cluster, the model is fitted at the household level with the household weight, which equals .weight_in only when weights are uniform within the household. .factor equals 1/\hat p only when num_classes = NULL; with propensity classes it is the class-level adjustment, so 1/.propensity does not reconstruct the applied factor – use .factor. The stage-by-stage weights are available through weight_factors() and domain_summary().

See Also

step_nonresponse(), collect_weights()

Other cascade audit: as_sae_input(), collect_step_detail(), collect_weights(), domain_summary(), weight_factors(), weighting_alerts()

Examples

fit <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "propensity",
                   formula = ~ sex + region, engine = "logit") |>
  prep()
p <- collect_propensities(fit)
summary(p$.propensity)

Collect replicate weights into a data frame ready for srvyr

Description

Returns the data with the point weight and every replicate weight as ordinary columns, plus the replication design as attributes. This is the form srvyr::as_survey_rep() and survey::svrepdesign() expect, and the form to write out when the analysis continues in another session, another script or another language.

Usage

collect_replicate_weights(
  object,
  weight_name = ".weight",
  prefix = "rep_",
  drop_zero = TRUE,
  scramble = FALSE
)

Arguments

Value

A data frame: the original columns, weight_name, and one column per replicate. The number of replicates is in attribute "R", and the replication design in attributes "type", "scale" and "rscales".

See Also

Other variance estimation: as_svydesign(), bootstrap_estimate(), bootstrap_weights(), jackknife_estimate(), jackknife_weights()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
boot <- bootstrap_weights(spec, replicates = 30, strata = "region",
                          psu = "psu", seed = 1, progress = FALSE)
df <- collect_replicate_weights(boot)   # or a weightflow_jack object
if (requireNamespace("srvyr", quietly = TRUE) &&
    requireNamespace("dplyr", quietly = TRUE)) {
  srvyr::as_survey_rep(df, weights = .weight,
                       repweights = dplyr::starts_with("rep_"),
                       type = attr(df, "type"), combined.weights = TRUE,
                       scale = attr(df, "scale"), rscales = attr(df, "rscales"),
                       mse = TRUE)
}

Per-unit detail of one step of the cascade

Description

A generic companion to collect_weights() that returns, for a single step, the weight it received and the multiplier it applied to every unit, plus any quantities that step computed internally (for a propensity step, the fitted propensity and its class). .weight_in and .factor are read from the stage-by-stage weights that prep() already stores, so .weight_in * .factor equals the weight leaving the step by construction, for any step. Native columns (those a step exposes on its own) are NA for units the step did not touch.

Usage

collect_step_detail(object, step = NULL)

Arguments

Value

The sample data.frame with .weight_in (the weight reaching the step, carrying every earlier adjustment) and .factor (the multiplier the step applied to each unit, NA where the incoming weight is zero) appended, plus any native columns of the chosen step (for a propensity step: .propensity, .responded, and .class when propensity classes are used), which are NA outside the units the step covers. Here .factor is defined for every unit with a nonzero incoming weight, so an active unit the step did not touch reports .factor = 1; this differs from collect_propensities(), where .factor is NA outside the propensity model (see its .status).

See Also

collect_weights(), collect_propensities(), weight_factors()

Other cascade audit: as_sae_input(), collect_propensities(), collect_weights(), domain_summary(), weight_factors(), weighting_alerts()

Examples

fit <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "propensity",
                   formula = ~ sex + region, engine = "logit") |>
  prep()
d <- collect_step_detail(fit, step = 1)
head(d[!is.na(d$.factor), c(".weight_in", ".factor", ".propensity")])

Extract the data with the computed weights

Description

Returns the sample as a data.frame with the final analysis weight attached as a column, ready to hand to an estimation routine. By default the units that left the cascade (weight 0) are dropped, so what comes back is the responding, in-scope sample and its weights.

Usage

collect_weights(
  object,
  drop_zero = TRUE,
  keep_intermediate = FALSE,
  weight_name = ".weight"
)

Arguments

Value

data.frame.

See Also

Other cascade audit: as_sae_input(), collect_propensities(), collect_step_detail(), domain_summary(), weight_factors(), weighting_alerts()

Examples

fitted <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
head(collect_weights(fitted))

Data-defect diagnostics for a non-probability sample

Description

For a non-probability sample, Meng (2018) decomposes the error of the sample mean into the data-defect correlation (rho, the population correlation between the target variable and the participation indicator), the sampling fraction and the problem difficulty. The effective sample size a probability sample would need to match the same mean-squared error is

n_{\mathrm{eff}} = \frac{f/(1-f)}{\rho^2}, \qquad f = n/N,

which does not depend on the outcome except through rho. Because rho on the target variable is not observable from the sample alone, data_defect() returns the effective size across a grid of plausible residual rho (read it as an ignorance range, not a single number), plus the measurable selection strength on the covariates used for pseudo-weighting (the largest correlation between an auxiliary and participation, which pseudo-weighting neutralises).

Usage

data_defect(object, ddc_grid = c(0.001, 0.005, 0.01, 0.05, 0.1))

Arguments

Value

a list (class weightflow_data_defect) with the sample size n, the estimated population size N (the sum of the final weights), the fraction f, the sensitivity grid (ddc, n_eff), and aux, a data frame of the covariate-participation correlations (or NULL when the recipe has no pseudo-weighting step).

References

Meng, X.-L. (2018). Statistical paradises and paradoxes in big data (I). Annals of Applied Statistics 12(2), 685-726.

See Also

step_pseudoweight(), report_weighting()

Kish design effect from unequal weighting

Description

Computes Kish's design effect due to unequal weighting, deff = 1 + CV^2(w) = m \sum w^2 / (\sum w)^2, and the effective sample size n_\mathrm{eff} = m / deff it implies. It is the standard one-number summary of what a weighting cascade cost in precision, and it is what the summary() and plot() methods of a prepped recipe report step by step.

Usage

design_effect(w)

Arguments

Details

Zero weights are dropped (they are the "dropped-unit" marker); negative weights – a valid but unusual output of unbounded linear/GREG calibration – are kept active, so the count n matches collect_weights(). Be aware, however, that the Kish formula assumes non-negative weights: a negative weight shrinks \sum w and enlarges \sum w^2 at once, so with negatives present deff is inflated and no longer interpretable as an effective-sample summary. prep() raises an alert when a calibration produces negative weights; prefer bounds to keep the factor positive if you need the design effect to be meaningful.

Value

list with deff, n_eff, cv and n.

Examples

design_effect(sample_survey$pw)

Flag re-identification risk from outlier weights within a publication cell

Description

A unit whose final weight is far larger than the rest of its publication cell is a disclosure risk in a public-use file: an extreme weight makes a rare unit stand out. disclosure_risk() flags, within each cell defined by by, the units whose final weight exceeds ratio times the cell's median weight, and reports the unit's share of the cell's total weight. Trimming (step_trim_weights() or the totals-preserving step_trim_calibrated()) is the usual remedy.

Usage

disclosure_risk(object, by, ratio = 10)

Arguments

Value

A data.frame, one row per flagged unit, with the row index (.row), the cell, the unit weight, the cell_median, the cell_n (active units in the cell) and cell_share (the unit's fraction of the cell's total weight), ordered from the largest weight down. Zero rows when nothing is flagged.

See Also

step_trim_weights(), step_trim_calibrated(), collect_replicate_weights()

Examples

fit <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
disclosure_risk(fit, by = "region")

Per-domain weight summary at every stage of the cascade

Description

For quality control by study domain (for example a department / DAM), this summarises how the weights move within each domain at every stage of the recipe: the base weights, then the weights after each step. It reads the stage-by-stage weights that prep() already stores, so it adds no computation to the cascade and never changes a weight.

Usage

domain_summary(object, by, min_n_eff = NULL)

Arguments

Value

A data.frame with one row per stage x domain and the columns stage (an ordered factor: base weights, then ⁠1. <step>⁠, ⁠2. <step>⁠, ...), domain (an ordered factor), n_active (active units in the domain at that stage), sum_w (sum of the active weights), mean_w, deff (the Kish design effect within the domain) and n_eff; and, when min_n_eff is given, publishable. Reading down a domain shows how its weight total and dispersion evolve step by step.

See Also

design_effect(), weight_factors(), summary.prepped_weighting_spec()

Other cascade audit: as_sae_input(), collect_propensities(), collect_step_detail(), collect_weights(), weight_factors(), weighting_alerts()

Examples

fit <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_calibrate(method = "raking", margins = list(region = c(table(population$region)))) |>
  prep()
domain_summary(fit, by = "region")

Jackknife estimate, standard error and confidence interval

Description

Applies a statistic to the point weights and to every delete-a-PSU replicate, and returns the estimate with its stratified jackknife (JKn) standard error and a normal confidence interval. jack_total() and jack_mean() are the shortcuts for a weighted total and a weighted mean of one column.

Usage

jackknife_estimate(
  jack,
  statistic,
  level = 0.95,
  ci_type = c("normal", "t"),
  df = NULL
)

jack_total(jack, variable)

jack_mean(jack, variable)

Arguments

Details

The stratified (JKn) variance sums each stratum's delete-a-PSU spread,

\widehat V_{JK} = \sum_h \frac{n_h - 1}{n_h}\sum_{i \in h}\big(\hat\theta_{(hi)} - \hat\theta_h\big)^2,

with \hat\theta_{(hi)} the estimate with PSU i of stratum h deleted and \hat\theta_h their within-stratum mean; the unstratified JK1 uses a single stratum. No finite population correction is applied.

Value

A data frame with estimate, se, ci_lower, ci_upper.

Note

jack_total() / jack_mean() center the replicate deviations on the per-stratum mean of the deleted-PSU estimates (the standard JKn). The survey design built by as_svrepdesign() instead uses mse = TRUE, which centers on the point estimate. Both are legitimate, so the standard errors from jack_total() and from svytotal() on the same object can differ slightly.

See Also

Other variance estimation: as_svydesign(), bootstrap_estimate(), bootstrap_weights(), collect_replicate_weights(), jackknife_weights()

Examples

spec <- weighting_spec(sample_one, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
jk <- jackknife_weights(spec, strata = "region", psu = "psu", progress = FALSE)
jackknife_estimate(jk, function(w, d) sum(w * d$employed, na.rm = TRUE))

Recipe-aware delete-a-PSU jackknife replicate weights

Description

Builds jackknife replicate weights by deleting one primary sampling unit (PSU) at a time and re-running the entire weighting recipe on each replicate. This is the deterministic sibling of bootstrap_weights(): same recipe-aware variance, no random number generation, and a replicate count fixed by the design rather than chosen by the analyst.

Usage

jackknife_weights(
  object,
  strata = NULL,
  psu = NULL,
  lonely_psu = c("certainty", "collapse"),
  cores = 1L,
  progress = TRUE
)

Arguments

Details

For a stratum h with n_h PSUs, the replicate that deletes PSU i zeros the base weight of that PSU and inflates the remaining PSUs of the stratum by n_h/(n_h-1); other strata are unchanged. There is one replicate per PSU. Strata with a single PSU contribute no variance and are skipped. This is the stratified jackknife (JKn); with strata = NULL it is the unstratified jackknife (JK1), and with psu = NULL each unit is its own PSU (delete-one-unit jackknife).

Value

An object of class weightflow_jack with the replicates matrix (units x replicates), the point weights, the per-replicate stratum and stratum size (used by jackknife_estimate()), and the design metadata.

See Also

Other variance estimation: as_svydesign(), bootstrap_estimate(), bootstrap_weights(), collect_replicate_weights(), jackknife_estimate()

Examples

spec <- weighting_spec(sample_one, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region))))
jk <- jackknife_weights(spec, strata = "region", psu = "psu", progress = FALSE)
jack_total(jk, "employed")

Level estimate for a single panel wave, with its replicate variance

Description

Estimates a statistic in ONE wave of a wave_bootstrap() or wave_jackknife() object, with the variance from that wave's replicates. Lets you report a level and a change from the same coordinated object.

Usage

level_estimate(
  wb,
  statistic,
  wave = NULL,
  level = 0.95,
  ci_type = c("normal", "t"),
  df = NULL
)

level_mean(wb, variable, wave = NULL, level = 0.95)

level_total(wb, variable, wave = NULL, level = 0.95)

Arguments

Value

an object of class weightflow_level with estimate, se, V and the interval.

See Also

change_estimate(), level_mean(), level_total()

Read the nonresponse-sensitivity analysis from a prepped recipe

Description

Returns the proxy pattern-mixture ignorance analysis stored by step_nr_sensitivity(): the adjusted mean of the study variable at each phi, with the ignorance interval and the proxy strength.

Usage

nr_sensitivity(object, step = NULL)

Arguments

Value

a list (class weightflow_nr_sensitivity) with table (phi, mu), rho, ybar_r, ignorance (the min-max interval over phi), mu_mar (the phi = 0 estimate) and the respondent/nonrespondent counts.

See Also

step_nr_sensitivity()

Synthetic rotating- and pure-panel datasets

Description

Four small, reproducible household-panel datasets that ship with weightflow to test and illustrate the panel tools. They share one structure and differ only in the rotation system, so the same code runs on a pure panel and on each rotating design. All are in long format: one row per person and per wave the person is in sample. Continuing units keep the same household_id / person_id across waves, so they link the panel; the household is the natural cluster and c(household_id, person_no) the person-level key.

Usage

panel_puro

panel_cl

panel_ine

panel_us

Format

A data.frame with one row per person-wave and the columns:

household_id

household id, persistent across waves (the panel link / cluster).

person_id, person_no

person id and person-number within household; c(household_id, person_no) is the person key.

stratum, psu

design stratum and primary sampling unit (PSU nested in stratum, at least two PSUs per stratum), for the coordinated bootstrap / jackknife.

region, sex, age

covariates usable as estimation domains.

wave

wave (month) index.

rotation_group, month_in_sample

rotation group and order-in-sample.

pw

design (base) weight.

disposition

between-wave disposition: "R", "NR", "OS", "UNK".

lf_status

labour status within the working-age population, a factor with levels "emp" / "unemp" / "inact"; NA for non-respondents. This is the status argument of step_cre() (the previous-wave composite auxiliary).

employed, unemployed

employed / unemployed indicators (labour force; NA if not "R" or not in the labour force).

income

labour income (NA if not "R").

Details

The between-wave disposition disposition follows the four-state taxonomy the longitudinal cascade needs: "R" responded, "NR" eligible nonresponse (reweight), "OS" out of scope – left the target population between waves, so the household exits permanently and is not reweighted – and "UNK" unknown eligibility. The variables of interest (employed, unemployed, income) are observed only when disposition == "R" (and, for the labour-force items, when the person is in the labour force), otherwise NA; they repeat across waves for continuing persons, with within-person persistence, so net change and gross flows are meaningful.

The datasets differ only in the rotation calendar (which waves each rotation group is in sample), giving different overlap structures:

panel_puro

Pure panel, no rotation: every unit is followed across all 4 waves; the sample shrinks only through attrition. Style of EU-SILC / SLID pure panels.

panel_cl

Chile ENE, 2-2-2 (in-out-in): 3 waves, consecutive overlap ~1/2, and units that return (in sample in waves 1 and 3 but not 2). The real ENE has no public rotation group; rotation_group is included for teaching, but the panel also links through the persistent ids alone.

panel_ine

INE Uruguay ECH / StatCan LFS, 6-month rotation: 6 groups in sample each wave, one sixth rotating out per wave, so consecutive overlap ~5/6.

panel_us

US CPS, 4-8-4: 3 waves, consecutive overlap ~3/4, with a cohort that leaves and returns after the 8-month gap.

Examples

# rotation structure of the 6-month panel
panel_design(panel_ine, unit = c("household_id", "person_no"), wave = "wave",
             rotation_group = "rotation_group", pattern = "6")
# coordinated change of the unemployment rate between two waves
t1 <- subset(panel_ine, wave == 1 & disposition == "R")
t2 <- subset(panel_ine, wave == 2 & disposition == "R")
wb <- wave_bootstrap(
  list(T1 = weighting_spec(t1, base_weights = pw),
       T2 = weighting_spec(t2, base_weights = pw)),
  replicates = 100, strata = "stratum", psu = "psu", seed = 1, progress = FALSE)
change_mean(wb, "unemployed")

Describe the rotating-panel structure of a survey

Description

Tags stacked per-wave survey data with its rotation structure so the panel steps and the report can read it, and so the overlap and linkage quality can be inspected before any weighting. panel_design() computes only; it does not create or change weights. It is the panel analogue of reference_sample(): the descriptor lives in attr(data, "wf_panel") and the returned object is still an ordinary data.frame.

Usage

panel_design(
  data,
  unit,
  wave,
  rotation_group = NULL,
  cluster = NULL,
  pattern = NULL,
  waves = NULL,
  reference_wave = NULL
)

Arguments

Details

data must be in stacked (long) form: one row per unit per wave, with a unit column that is stable across waves and a wave column identifying the period. From the unit x wave membership it derives the observed overlap matrix (the fraction of each wave retained in every other wave, i.e. CEPAL's traslape), and, when rotation_group is given, the panel-selection probability ⁠Pr(panel selection)⁠ for each adjacent pair and for the full combination – the reciprocal of which is the CEPAL panel base-weight factor.

The linkage quality is measured, not assumed: if a pattern is supplied the observed adjacent overlap is compared against the one it implies, and a large gap (alert PN-01) is the early warning that the linkage key is unstable (relabelled ids, a redesigned frame, duplicated panels). Uneven rotation-group sizes, which break the scalar reciprocal of ⁠Pr(panel selection)⁠, raise PN-02.

Value

data, unchanged as a data frame, with the descriptor in attr(data, "wf_panel") and class "wf_panel_design" prepended.

See Also

panel_merge(), reference_sample()

Examples

# person-level tracking: the key is household id + person line number
pd <- panel_design(panel_ine, unit = c("household_id", "person_no"), wave = "wave",
                   rotation_group = "rotation_group", cluster = "household_id",
                   pattern = "6")
pd            # overlap matrix, Pr(panel selection), linkage rate, alerts
summary(pd)

Linear combination of panel waves, with honest between-wave variance

Description

Estimates any linear combination psi = sum(contrast * theta_wave) of a statistic across panel waves – an annual average (contrast = rep(1/W, W)), a net change (contrast = c(-1, 1), i.e. what change_estimate() does), a semester contrast, etc. – and its variance ⁠a' Sigma a⁠, where Sigma is the between-wave covariance matrix that the overlap of a rotating panel induces. That covariance is taken straight from the coordinated replicates of wave_bootstrap() or wave_jackknife(), so it is correct without any analytic covariance formula. Treating the waves as independent (the diagonal of Sigma only) understates the variance of an average and overstates that of a change.

Usage

panel_estimate(
  wb,
  statistic,
  contrast = NULL,
  waves = NULL,
  level = 0.95,
  ci_type = c("normal", "t"),
  df = NULL
)

panel_mean(wb, variable, contrast = NULL, waves = NULL, level = 0.95)

panel_total(wb, variable, contrast = NULL, waves = NULL, level = 0.95)

Arguments

Value

an object of class weightflow_panel_estimate with estimate, se, V, Vind (the independence-assuming variance), deff = V / Vind, the covariance matrix Sigma, the per-wave point estimates, the contrast, and the confidence interval.

See Also

change_estimate(), wave_bootstrap(), wave_jackknife()

Examples

waves <- lapply(1:3, function(t)
  weighting_spec(subset(panel_ine, wave == t & disposition == "R"), base_weights = pw))
names(waves) <- c("T1", "T2", "T3")
wb <- wave_bootstrap(waves, replicates = 100, strata = "stratum", psu = "psu",
                     seed = 1, progress = FALSE)
rate <- function(w, d) weighted.mean(d$unemployed, w, na.rm = TRUE)
panel_estimate(wb, rate)                       # average unemployment level over the waves
panel_estimate(wb, rate, contrast = c(-1, 0, 1))  # T3 - T1 change

Build the wide longitudinal file from per-wave surveys

Description

Joins a named list of per-wave data.frames on a stable unit key into the wide, one-row-per-unit file that the longitudinal recipe consumes. Each wave's non-key columns are suffixed with the wave name, and per-wave presence indicators ⁠.wf_in_<wave>⁠ (and, if responded is given, response indicators ⁠.wf_resp_<wave>⁠) are added so a later step_attrition() can model the response pattern.

Usage

panel_merge(
  waves,
  by,
  responded = NULL,
  require = c("any", "all"),
  suffix = "_"
)

Arguments

Details

require = "any" (the default) keeps every unit seen in at least one wave – which the attrition model needs, because it has to be fitted over responders and non-responders. require = "all" keeps only the intersection (CEPAL's ⁠s(2) = s1 cap s2 cap ...⁠); prefer to reach it with step_attrition() after estimating, not by dropping units before the model can see them.

Value

a wide data.frame, one row per unit, ready for panel_design() / weighting_spec().

See Also

panel_design()

Examples

# reshape two waves of the long panel into one wide, one-row-per-unit file
wide <- panel_merge(
  list(T1 = subset(panel_ine, wave == 1), T2 = subset(panel_ine, wave == 2)),
  by = c("household_id", "person_no"), require = "any")
names(wide)                       # per-wave columns are suffixed _T1 / _T2

Panel-selection probability for a set of combined waves

Description

Returns ⁠Pr(panel selection)⁠ for combining waves in a panel_design() that carries a rotation_group: the fraction of rotation-group cohorts present in all the combined waves. Its reciprocal is the CEPAL panel base-weight factor used by step_panel_overlap() (prob = panel_pr(pd, c("T1","T2"))). Surveys without a public rotation-group variable (e.g. the Chilean ENE) cannot use this – derive the probability another way and pass it to step_panel_overlap().

Usage

panel_pr(object, waves = NULL)

Arguments

Value

a single probability in ⁠(0, 1]⁠.

See Also

panel_design(), step_panel_overlap()

Examples

pd <- panel_design(panel_ine, unit = c("household_id", "person_no"), wave = "wave",
                   rotation_group = "rotation_group")
panel_pr(pd)                    # over all waves
panel_pr(pd, c("1", "2"))       # combining waves 1 and 2

Diagnostic plots for the weights

Description

Draws the weighting cascade: one histogram of the adjustment factor per step, plus a four-panel summary of the final weights, the cumulative factor, base against final weight, and the design effect by stage. Base graphics only, no dependencies. Use it after prep() to see how the weights moved, where summary() tells you by how much.

Usage

## S3 method for class 'prepped_weighting_spec'
plot(x, type = c("all", "factors", "summary"), ...)

Arguments

Value

Invisibly, the prepped object x. Called for its side effect of drawing the diagnostic plots described above.

Examples

fitted <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
plot(fitted)

Synthetic target population (sampling frame)

Description

The simulated frame every weightflow example draws from: 4,495 persons nested in 1,882 households, in 120 primary sampling units, in 4 regions used as strata. Because the whole population is observed, it supplies the known population totals a calibration step needs, and the true values against which a weighted estimate can be checked.

Usage

population

Format

A data frame with one row per person:

person_id

individual identifier

household_id

household identifier (cluster)

psu

primary sampling unit (segment) within the stratum

region

stratum: North, South, East or West

sex

F or M

age

age in years (18-95)

income

annual income

employed

employment indicator (0/1)

Estimate the weighting cascade

Description

Runs an inert weighting_spec() recipe. Starting from the design base weights, prep() applies each step in the order it was piped, multiplying the current weight by that step's adjustment factor, and returns an object holding the weight at every stage, the per-step diagnostics and the quality alerts. This is the only function that computes weights.

Usage

prep(spec, min_cell_n = 30, max_factor = 2.5, warn = FALSE)

Arguments

Value

a "prepped_weighting_spec" object. Every quality incident is recorded in ⁠$alerts⁠ (readable with weighting_alerts() / has_alerts()), regardless of warn. This includes warnings a step raises internally, such as a calibration that could not meet its constraints: they are captured into ⁠$alerts⁠ even when the surrounding warnings are suppressed, so ⁠$alerts⁠ is the single reliable channel for programmatic quality control.

See Also

weightflow-alerts for the catalogue of quality alerts prep() can raise, weighting_alerts() / has_alerts() to read them, and vignette("inspecting-auditing") for the full quality-control workflow.

Examples

rec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region")
prep(rec)

Print a bootstrap replicate-weight object

Description

Compact one-screen summary of a weightflow_boot object: how many replicates were requested, how many units are in the data, how many of them are still active (final weight above zero), and which columns defined the resampling design.

Usage

## S3 method for class 'weightflow_boot'
print(x, ...)

Arguments

Value

(invisibly) the object.

Print a jackknife replicate-weight object

Description

Compact one-screen summary of a weightflow_jack object: how many delete-a-PSU replicates were built, how many units are in the data, how many of them are still active (final weight above zero), and which columns defined the deletion design.

Usage

## S3 method for class 'weightflow_jack'
print(x, ...)

Arguments

Value

(invisibly) the object.

Read a weighting recipe from a YAML file

Description

Reads a recipe written by write_recipe(). With data = NULL (the default) it returns an inspectable recipe manifest (for review or archival); pass data to reconstruct an executable weighting_spec bound to that data, ready for prep().

Usage

read_recipe(file, data = NULL, references = NULL, allow_code = FALSE)

Arguments

Details

A reconstructed recipe is validated only when you call prep(): a hand-edited recipe with an out-of-range or mistyped value surfaces its error there, not at read_recipe() time. And because reading a recipe evaluates the stored expressions, only read recipes you trust, as you would source() an R script.

Value

With data = NULL, a weightflow_recipe manifest (a list with a print method). With data, a weighting_spec.

See Also

write_recipe()

Other recipe serialization: write_recipe()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region")
f <- tempfile(fileext = ".yml"); write_recipe(spec, f)
read_recipe(f)                       # inspect the manifest
spec2 <- read_recipe(f, data = sample_survey)   # rebuild an executable recipe

Use a weighted survey as the calibration reference instead of a frame

Description

Wraps a reference-survey microdata data.frame together with its design weights so it can be passed as the population argument of step_model_calibration() (and any step that takes a population frame). The calibration totals are then the weighted sums over the reference survey – an estimate of the population totals – instead of unweighted sums over a full frame. This is the model-assisted / two-survey setup: fit the model on your sample, project it onto a larger reference survey, and calibrate to the weighted totals of the projection (Wu and Sitter 2001; Kim and Rao 2012).

Usage

reference_sample(data, weights, replicates = NULL)

Arguments

Details

A reference survey with all weights equal to 1 reproduces the plain-frame behaviour exactly. To propagate the reference survey's own sampling variance into the recipe-aware bootstrap, pass its replicate weights through replicates: each bootstrap replicate then re-estimates the totals from the paired reference replicate (Opsomer and Erciulescu 2021), so the extra variance from estimating the totals is captured. Without replicates the totals are treated as fixed (a reasonable approximation when the reference is much larger than the sample, and the same assumption made when calibrating to another survey's published totals).

Value

data tagged so that step_model_calibration() weights its totals by weights. It is still an ordinary data.frame.

See Also

step_model_calibration()

Examples

ref <- reference_sample(population, weights = rep(1, nrow(population)))

Panel / longitudinal HTML report

Description

Builds an HTML report by ADDING the panel cards to the standard report. When a longitudinal weight is given, it renders the full report_weighting() page for that weight (cascade, weight distribution, deff, ...) and injects a "Panel / longitudinal" section with the rotation structure (panel_design()), the coordinated net-change variance (change_estimate()), the attrition/retention summary, and the gross flows (transition_matrix() / boot_transition()). Without a longitudinal weight it writes a small standalone page with those cards. Any argument may be NULL; its card is skipped.

Usage

report_panel(
  design = NULL,
  change = NULL,
  longitudinal = NULL,
  transition = NULL,
  coordinated = NULL,
  estimates = NULL,
  variance = NULL,
  file = NULL,
  open = TRUE,
  lang = c("en", "es")
)

Arguments

Value

the path to the written HTML file, invisibly.

See Also

report_weighting(), panel_design(), change_estimate(), boot_transition()

Self-contained HTML quality report for a weighting recipe

Description

Writes a single, self-contained HTML file documenting how a prepped recipe turned the design base weights into the final weights: the cascade, the parameters requested at each step, the per-stage weight summary, the step-specific diagnostics (calibration, nonresponse, trimming), the fieldwork outcome rates and an audit trail. It is the deliverable to attach to a methodological report or a weighting annex: everything is inline, so the file opens offline and can be archived or emailed as one artifact.

Usage

report_weighting(
  object,
  file = NULL,
  open = TRUE,
  plots = TRUE,
  narrative = TRUE,
  lang = c("en", "es"),
  metadata = NULL,
  replicates = NULL,
  domains = NULL,
  y_vars = NULL
)

Arguments

Value

(invisibly) the path to the HTML file.

Examples

fitted <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
# writes a self-contained HTML report to a temporary file (open = FALSE so
# nothing is launched); use open = TRUE to view it in the browser.
path <- report_weighting(fitted, open = FALSE)

Synthetic address sample with one selected person per household

Description

A multistage sample of 417 addresses (stratum, then PSU, then household, then one person inside the reached household) carrying every complication a household survey meets: unresolved eligibility, out-of-scope addresses, household nonresponse, unequal within-household selection and person nonresponse. It is the dataset that exercises the complete weighting cascade.

Usage

sample_one

Format

A data frame with one row per sampled household (the selected person, or a single placeholder row for non-roster cases):

person_id, household_id, psu

identifiers

region

stratum

sex, age

selected person's attributes (NA on non-roster rows)

pw

design base weight (product of the stage selection probabilities)

status

"eligible", "ineligible" or "unknown"

disposition

full field disposition as a single factor (a recode of the indicator columns): "eligible respondent", "eligible nonrespondent", "household nonresponse", "ineligible" or "unknown eligibility"

unknown_elig

1 if eligibility is unknown (no roster)

ineligible

1 if the address is out of scope (no roster)

hh_responded

1 reached, 0 household nonresponse, NA for non-eligible

responded

1 if the selected person responded (NA on non-roster rows)

n_elig

number of eligible persons in the household (NA on non-roster rows)

p_within

within-household selection probability of the selected person

income, employed

survey outcomes; NA unless the person responded

Synthetic person sample with a take-all household roster

Description

A stratified two-stage sample of 467 persons drawn from population: PSUs within region, then households within PSU, then every adult of the selected household (take-all roster). It carries unequal design base weights, an unknown-eligibility flag and a person-level response indicator, and it is the dataset the short examples in this package use.

Usage

sample_survey

Format

A data frame with one row per sampled person:

person_id, household_id, psu

identifiers

region, sex, age

frame auxiliaries, known for all units

pw

design base weight (inverse sampling fraction)

unknown_elig

1 if eligibility is unknown

responded

1 if the person responded

income, employed

survey outcomes; NA for nonrespondents

Assert quality conditions on the weights

Description

A checkpoint that leaves the weights untouched and instead verifies that they meet quality thresholds at this point of the cascade, raising an error or a warning when they do not. Use it to stop a production pipeline before bad weights are published, in the spirit of a validation step inside a recipe.

Usage

step_assert(
  spec,
  max_deff = NULL,
  max_weight_ratio = NULL,
  min_n_eff = NULL,
  on_fail = c("error", "warning"),
  id = NULL
)

Arguments

Value

The input weighting_spec with this checkpoint appended to its recipe. The check is recorded only; it is evaluated when prep() is called and does not modify the weights.

See Also

Other weighting steps: step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_assert(max_deff = 5, on_fail = "warning") |> prep()

Attrition adjustment for panel waves

Description

The panel-facing nonresponse step: it adjusts for attrition (nonresponse between waves) when building a longitudinal weight, reweighting the units that stayed in to also represent those that dropped out, among the units that remain eligible. It is a thin wrapper over step_nonresponse() – so the theory stays visible in the recipe – that uses the same estimators but under a name that reads correctly in a panel cascade and with the panel conventions: the covariates should come from a wave where the unit was observed (e.g. the first period), and it does NOT absorb eligibility (out-of-scope and unknown-eligibility go in step_drop_ineligible() / step_unknown_eligibility()).

Usage

step_attrition(
  spec,
  respondent,
  method = c("propensity", "rhg", "weighting_class", "calibration"),
  formula = NULL,
  by = NULL,
  engine = c("logit", "tree", "forest", "boost"),
  num_classes = 5L,
  weight_model = TRUE,
  cluster = NULL,
  crossfit = NULL,
  crossfit_seed = NULL,
  totals = NULL,
  count = NULL,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  penalty = NULL,
  equal_within_cluster = FALSE,
  maxit = 50L,
  tol = 1e-06,
  id = NULL
)

Arguments

Details

The attrition adjustment prescribed by ECLAC's household-survey manual (ch. XVI) and used by Statistics Canada's SLID (Naud 2002; LaRoche 2003) is response-propensity weighting (method = "propensity", i.e. inverse of the estimated response probability; Little 1986; Rosenbaum 1987) – not an ECLAC invention. method = "rhg" is the response-homogeneity-group variant (propensity stratified into num_classes groups, then the class-mean rate), which stabilises the weights.

Two ECLAC prescriptions are not implemented, and the difference matters when they apply. The manual (ch. XVI, sec. B.1.b) allows two fallback imputations for units with no auxiliary information at all: a nonrespondent whose rotation panel does not overlap (impute the overall effective response rate as its phi), and a newly incorporated nonrespondent (impute the adjusted expansion factor of its household). Here a covariate that is NA for an eligible unit is an error, not an imputation – an NA propensity would let that nonrespondent survive the adjustment silently, which is the failure the error exists to prevent, and the package will not silently substitute a value of its own for a missing input.

So when the error fires, the fix belongs in the data, and either ECLAC route is available to you there: impute the missing covariates before the step (the manual's own fallback is the overall effective response rate, i.e. a constant, which is what a model with no covariates for those units amounts to), give a newly incorporated nonrespondent its household's adjusted factor, or restrict formula to a covariate set observed for every eligible unit. What the package will not do is choose one of those for you and leave no trace in the recipe. Multi-wave retention chaining (decomposing Pr(in at T3) into Pr(reach T2) x Pr(T2 to T3)) is likewise absent: this step adjusts one transition at a time, which is what ch. XVI specifies for two consecutive periods.

Value

the weighting_spec with the attrition step appended.

The arguments of step_nonresponse()

Attrition is nonresponse over waves, so this step delegates to step_nonresponse() and takes its arguments too. That was not true before: the twelve remaining ones were missing, including crossfit – which the step's own quality alert recommends when a flexible engine overfits the propensity, a remedy that could not be applied – and totals / calfun / bounds, i.e. the whole Sarndal-Lundstrom calibration route with nothing to configure.

See Also

step_nonresponse(), step_panel_overlap(), step_drop_ineligible()

Examples

wide <- panel_merge(
  list(T1 = subset(panel_ine, wave == 1), T2 = subset(panel_ine, wave == 2)),
  by = c("household_id", "person_no"), require = "all")
weighting_spec(wide, base_weights = pw_T1) |>
  step_panel_overlap(prob = 5 / 6) |>
  step_drop_ineligible(disposition_T2 == "OS", reason = "left the target population") |>
  step_attrition(respondent = disposition_T2 == "R", method = "propensity",
                 formula = ~ age_T1 + sex_T1 + region_T1) |>
  prep()

Calibration to population totals

Description

Adjusts the weights so that the weighted sample reproduces known population totals of auxiliary variables, while moving the incoming weights as little as possible. This is the calibration estimator of Deville and Sarndal (1992), with raking, post-stratification and linear (GREG) calibration, optional bounds on the adjustment factor, ridge relaxation of the targets, domain partitions and one-weight-per-cluster (integrative) calibration.

Usage

step_calibrate(
  spec,
  margins = NULL,
  method = c("raking", "poststratify", "linear"),
  formula = NULL,
  totals = NULL,
  count = NULL,
  by = NULL,
  cluster = NULL,
  equal_within_cluster = FALSE,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  maxit = 50L,
  tol = 1e-06,
  penalty = NULL,
  population = NULL,
  id = NULL
)

Arguments

Details

The calibrated weight is w_i = d_i g_i, close to the incoming weights d_i and reproducing the totals \sum_{i \in s} w_i \mathbf{x}_i = \mathbf{X}, with the factor g_i minimizing a distance to d_i: linear g_i = 1 + \mathbf{x}_i'\boldsymbol\lambda or exponential ("raking") g_i = \exp(\mathbf{x}_i'\boldsymbol\lambda). The penalty (ridge) relaxes the exact constraint to steady extreme weights when the auxiliaries are many or collinear, minimizing

\sum_i \frac{(w_i - d_i)^2}{d_i q_i} + \frac{1}{s}\sum_j c_j (\hat X_j - X_j)^2,

where c_j is the cost of missing constraint j: as c_j \to \infty the constraint is met exactly, as c_j \to 0 the weights return to d_i.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

# Raking to population margins
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_calibrate(method = "raking", id = "calib_main",
                 margins = list(sex    = c(table(population$sex)),
                                region = c(table(population$region)))) |>
  prep()
# the id ("calib_main") labels the step in the print-out and selects it in
# collect_step_detail(fit, "calib_main")

# ridge (penalized) calibration: relaxes the targets to control extreme
# weights; a smaller penalty relaxes more. Uses only base R.
pop_tot <- c("(Intercept)" = nrow(population),
             regionSouth = sum(population$region == "South"),
             regionEast  = sum(population$region == "East"),
             regionWest  = sum(population$region == "West"),
             sexM        = sum(population$sex == "M"))
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_calibrate(method = "linear", formula = ~ region + sex,
                 totals = pop_tot, penalty = 1) |>
  prep()

# --- Tidy `totals` format (recommended) ---------------------------------
# Post-stratification: give the population counts as a data frame with one or
# more category columns plus a counts column named by `count`.
ps_totals <- as.data.frame(table(region = population$region, sex = population$sex))
weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "poststratify", totals = ps_totals, count = "Freq") |>
  prep()

# Raking: a list of data frames, one per margin.
m_region <- as.data.frame(table(region = population$region))
m_sex    <- as.data.frame(table(sex = population$sex))
weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 totals = list(m_region, m_sex), count = "Freq") |>
  prep()

# Linear/GREG with mixed auxiliaries: data frames for categoricals (all
# categories) and a single number for a continuous total. weightflow builds
# the model.matrix totals internally, so you never drop a reference category.
resp <- subset(sample_survey, responded == 1)
weighting_spec(resp, base_weights = pw) |>
  step_calibrate(method = "linear", formula = ~ region + sex + income,
                 totals = list(region = m_region, sex = m_sex,
                               income = sum(population$income)),
                 count = "Freq") |>
  prep()

# Domain (partitioned) calibration: `by` calibrates independently within each
# domain, each to its own totals. The domain is an extra column in the tidy
# totals, not a term in the formula/margins. Here we rake two margins (sex and
# age group) within each region, which a single global cross could not express.
pop  <- transform(population,
  age_grp = cut(age, c(0, 30, 45, 60, Inf), labels = c("18-30","31-45","46-60","60+")))
samp <- transform(sample_survey,
  age_grp = cut(age, c(0, 30, 45, 60, Inf), labels = c("18-30","31-45","46-60","60+")))
sex_by_region <- as.data.frame(table(region = pop$region, sex     = pop$sex))
age_by_region <- as.data.frame(table(region = pop$region, age_grp = pop$age_grp))
weighting_spec(samp, base_weights = pw) |>
  step_calibrate(method = "raking",
                 totals = list(sex_by_region, age_by_region),
                 count = "Freq", by = "region") |>
  prep()

Composite regression estimator (CRE / regression composite estimation)

Description

Final-weight step for rotating panels that exploits the sample overlap in the estimation, not only in the variance. It is the composite calibration of the Canadian LFS (Gambino, Kennedy and Singh 2001; Fuller and Rao 2001) and the composite calibration of the Uruguayan ECH (INE, methodology sec. 8.4): the same estimator. On top of the ordinary calibration to known population totals X, it adds a second set of constraints to control totals Zhat estimated from the previous wave (the composite auxiliaries: previous-month labour status, optionally crossed by domains). Because Zhat uses the previous wave's composite weights, the estimator is recursive.

Usage

step_cre(
  spec,
  previous = NULL,
  status,
  composite = list(NULL),
  id_unit,
  formula,
  totals = NULL,
  count = NULL,
  birth = NULL,
  alpha = 2/3,
  overlap = "auto",
  on_missing_prev = c("carry_backward", "zero"),
  cluster = NULL,
  equal_within_cluster = FALSE,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  rotation_group = NULL,
  status_ref = NULL,
  id = NULL
)

Arguments

Details

The composite auxiliary z is known for the overlap sample (units present in t-1) and imputed for the birth rotation group (units entering at t, with no t-1 value) by combining two imputations: MR1 (mean imputation, aimed at the level) and MR2 (carry-backward, aimed at the change), mixed as ⁠z = (1 - alpha) z1 + alpha z2⁠ with alpha = 2/3.

The composite weights are w_cre = w_nr * g, with g the GREG factor of the augmented calibration ⁠[x | z]⁠ to ⁠[X | Zhat]⁠; the same solver as step_calibrate(). With equal_within_cluster = TRUE the auxiliaries are averaged to the household (integrated method of weighting, Lemaitre-Dufour 1987) so every household member shares one final weight.

Value

the input weighting_spec with the CRE step appended (evaluated at prep()). Chain waves by passing each prepped wave as the previous of the next.

References

Gambino J, Kennedy B, Singh MP (2001). Regression composite estimation for the Canadian Labour Force Survey. Survey Methodology 27(1):65-74. Fuller WA, Rao JNK (2001). A regression composite estimator with application to the Canadian Labour Force Survey. Survey Methodology 27(1):45-51. Instituto Nacional de Estadistica (Uruguay). Metodologia de la Encuesta Continua de Hogares, seccion 8.4 (calibracion compuesta).

See Also

Other weighting steps: step_assert(), step_calibrate(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

# Composite (CRE) estimation on the 6-month rotating panel `panel_ine`.
# `lf_status` is the previous-wave labour status (emp / unemp / inact).
t1 <- subset(panel_ine, wave == 1 & disposition == "R")
t2 <- subset(panel_ine, wave == 2 & disposition == "R")
t1$sex <- factor(t1$sex); t2$sex <- factor(t2$sex)
Xtot <- function(d) colSums(d$pw * stats::model.matrix(~ sex, data = d))

# seed wave: no previous month, so step_cre() reduces to a linear calibration to X
seed <- weighting_spec(t1, base_weights = pw) |>
  step_cre(previous = NULL, status = lf_status, formula = ~ sex,
           totals = Xtot(t1), status_ref = "inact") |>
  prep()

# composite wave: augment X with the previous-wave status, country-level and by sex
fit2 <- weighting_spec(t2, base_weights = pw) |>
  step_cre(previous = seed, status = lf_status, composite = list(NULL, "sex"),
           id_unit = c("household_id", "person_no"), formula = ~ sex, totals = Xtot(t2),
           alpha = 2/3, status_ref = "inact") |>
  prep()
fit2

Declare the recipe's scope: cross-sectional or longitudinal weights

Description

Two argument-free declarative steps that say what the recipe is building. They do not touch the weights (like step_assert()); they record the intent so the rest of the recipe, the report and the estimators behave accordingly, and so a wrong-purpose estimate can be flagged. All the panel detail – waves, rotation group, reference wave – already lives in panel_design(), so these steps need no arguments.

Usage

step_cross_sectional(spec, id = NULL)

step_longitudinal(spec, id = NULL)

Arguments

Details

Put the scope step first. step_cross_sectional() is the default (a plain recipe with no scope step behaves as cross-sectional). step_longitudinal() requires the data to be tagged by panel_design() and marks the recipe as building the panel (longitudinal) weight, so steps like step_panel_overlap() apply and the calibration targets the reference wave.

Value

The input weighting_spec with the scope declared.

ECLAC conventions for the longitudinal weight

The three points below are conventions of ECLAC's household-survey manual (ch. XVI). The package cannot enforce them – a vector of control totals carries no label saying which period it belongs to – so they are the reader's to apply.

Who is in the longitudinal population. It is made up of the units that were in the target population at the first period and remained in it through the last. Units that left (death, migration, institutionalization) are out, and so are units that entered after the first period: a panel adds no elements over time. Drop the first group with step_drop_ineligible() before the attrition step – leaving the target population is not nonresponse and must not be compensated by redistributing weight.

Calibrate to the FIRST period's totals. ECLAC is explicit that the auxiliary totals used in the calibration must represent the population of the first period of interest, because a panel that adds no elements over the measurement periods is representative only of the period in which it was selected (ch. XVI, sec. B.2.c). In the original:

"los totales auxiliares utilizados en la calibracion deben representar la
 poblacion del primer periodo de interes, puesto que, al conformar un panel
 que no suma elementos a lo largo de los periodos de medicion, la muestra
 sera representativa unicamente del periodo en el que fue seleccionada"

Passing the later period's projections is a silent error: the recipe converges, the totals close, and the weights represent a population the sample never had a chance to cover.

Cross-sectional estimates from a longitudinal file are reference only. They will not match the published cross-sectional figures and should not: the target population of the panel combination is not the target population of the cross-section (ch. XVI, sec. C). Use the longitudinal weight for gross flows, transitions and durations – what it exists for – and the cross-sectional weight for levels.

See Also

panel_design(), step_panel_overlap()

Examples

wide <- panel_merge(
  list(T1 = subset(panel_ine, wave == 1), T2 = subset(panel_ine, wave == 2)),
  by = c("household_id", "person_no"), require = "all")
weighting_spec(wide, base_weights = pw_T1) |>
  step_longitudinal() |>
  step_panel_overlap(prob = 5 / 6) |> prep()

Declarative estimation over a coordinated panel object

Description

step_domain() and step_estimate() build an estimation pipeline (weightflow_estimation) on top of a saved wave_bootstrap() / wave_jackknife() object, so that changes, levels and linear contrasts – with the honest overlap covariance – can be requested declaratively and disaggregated, the way the weighting recipe is built with ⁠step_*⁠. It is a thin front end over change_estimate(), panel_estimate() and level_estimate(); the coordinated replicates do the variance. Evaluate with collect_estimates() (printing does it too).

Usage

step_domain(x, ...)

step_filter(x, condition)

step_estimate(
  x,
  statistic,
  over = NULL,
  type = c("absolute", "relative"),
  waves = NULL,
  contrast = NULL,
  level = 0.95,
  label = NULL
)

step_transition(x, from, to, format = c("row", "col", "joint", "counts"))

Arguments

Value

a weightflow_estimation.

See Also

collect_estimates(), change_estimate(), panel_estimate(), level_estimate()

Examples

t1 <- subset(panel_ine, wave == 1 & disposition == "R")
t2 <- subset(panel_ine, wave == 2 & disposition == "R")
wb <- wave_bootstrap(
  list(T1 = weighting_spec(t1, base_weights = pw),
       T2 = weighting_spec(t2, base_weights = pw)),
  replicates = 100, strata = "stratum", psu = "psu", seed = 1, progress = FALSE)

# net change of the unemployment rate, by region
collect_estimates(wb |> step_domain(region) |>
  step_estimate(mean(unemployed), over = "change"))

# subpopulation: same change among the working-age population (rows masked, not dropped)
collect_estimates(wb |> step_filter(age >= 25 & age <= 54) |>
  step_estimate(mean(unemployed), over = "change"))

Drop ineligible (out-of-scope) units

Description

Sets the weight of the units known to be outside the target population to zero, so they leave the cascade and take no part in any later step or in collect_weights(). Their weight is discarded, not redistributed: the weight total is meant to fall by exactly the mass they carried. Use it once eligibility has been resolved, immediately after step_unknown_eligibility().

Usage

step_drop_ineligible(spec, ineligible, reason = NULL, id = NULL)

Arguments

Details

Apply it AFTER step_unknown_eligibility: ineligibles must be present and NOT flagged as unknown during that step, so they take part in the known-eligibility group and receive their share of the redistributed unknown weight. Their weight is then correctly discarded here (it represents the ineligible share of the unknown units, which are out of scope).

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

df <- transform(sample_survey,
                ineligible = as.integer(region == "West" & age > 90))
weighting_spec(df, base_weights = pw) |>
  step_drop_ineligible(ineligible = ineligible) |>
  prep()

Model-assisted calibration (Wu and Sitter 2001)

Description

Fits a working model for each study variable, predicts it over the whole population, and calibrates the weights so that the sample total of every prediction matches its population total, on top of the usual auxiliary totals – which may come from the population frame itself or from an external source (a census table, an administrative register). Reach for it when you hold, or can supply, those control totals and the outcome is well predicted by the auxiliaries: the predictions act as extra, highly relevant controls and buy precision that calibrating on x alone cannot.

Usage

step_model_calibration(
  spec,
  x_formula,
  models,
  population,
  x_totals = NULL,
  count = "Freq",
  cluster = NULL,
  equal_within_cluster = FALSE,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  maxit = 100L,
  tol = 1e-07,
  crossfit = NULL,
  crossfit_seed = NULL,
  id = NULL
)

Arguments

Details

Requires COMPLETE auxiliary information: a data.frame population with the x_formula columns and the model predictors for the whole population (or a reference frame/census).

The predictions \hat y_i enter as extra constraints, \sum_{i \in s} w_i \hat y_i = \sum_{i \in U} \hat y_i, solved together with the benchmark auxiliary totals \mathbf{X}. When the working model is linear this reduces to GREG; a nonlinear learner adds efficiency through the prediction constraint while the totals \mathbf{X} preserve design consistency even if the model is misspecified.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

References

Wu, C. and Sitter, R. R. (2001). A model-calibration approach to using complete auxiliary information from survey data. Journal of the American Statistical Association, 96(453), 185-193. doi:10.1198/016214501750333054.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_model_calibration(
    x_formula  = ~ sex + region,
    models     = list(income = y_model(income ~ age + sex, engine = "glm")),
    population = population) |>
  prep()

# with cross-fitting (out-of-fold predictions, avoids overfitting)
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_model_calibration(
    x_formula  = ~ sex + region,
    models     = list(income = y_model(income ~ age + sex, engine = "glm")),
    population = population, crossfit = 5, crossfit_seed = 1) |>
  prep()

# consistency totals from an external source (tidy format): a data frame per
# factor and a single number per continuous total. `population` is still used
# for the model predictions. Adjust for nonresponse first, since the outcome
# is only observed for respondents.
m_region <- as.data.frame(table(region = population$region))
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  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 = sum(population$age)),
    count      = "Freq") |>
  prep()

# equal weights within a household (integrative, Lemaitre-Dufour): one weight
# per cluster, so person and household estimates stay coherent. The final
# weights are constant within each cluster among its active members.
fit_hh <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_model_calibration(
    x_formula  = ~ sex + region,
    models     = list(income = y_model(income ~ age + sex, engine = "glm")),
    population  = population,
    cluster = "household_id", equal_within_cluster = TRUE) |>
  prep()
w <- fit_hh$final_weight
max(tapply(w[w > 0], sample_survey$household_id[w > 0],
           function(x) diff(range(x))))    # 0: one weight per household

Nonresponse adjustment

Description

Inflates the weights of the eligible respondents so they also represent the eligible nonrespondents, under the assumption that response is ignorable given the information used. Three estimators are available – weighting classes, a response-propensity model (four engines, with optional cross-fitting), and calibration of the respondents to auxiliary totals – at the unit level or, through cluster, at a coarser level (e.g. the household).

Usage

step_nonresponse(
  spec,
  respondent,
  method = c("weighting_class", "propensity", "calibration"),
  by = NULL,
  formula = NULL,
  engine = c("logit", "tree", "forest", "boost"),
  weight_model = TRUE,
  num_classes = 5L,
  cluster = NULL,
  crossfit = NULL,
  crossfit_seed = NULL,
  totals = NULL,
  count = NULL,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  penalty = NULL,
  equal_within_cluster = FALSE,
  maxit = 50L,
  tol = 1e-06,
  id = NULL
)

Arguments

Details

The three estimators are the same operation with a different inflation factor. Weighting classes inflate the responding weight to the cell total, f_c = \sum_{i \in c} w_i / \sum_{i \in c} r_i w_i (with r_i = 1 if i responds), applied as w_i \leftarrow f_c w_i. A response-propensity model instead adjusts unit by unit, w_i \leftarrow w_i / \hat\phi_i, with \hat\phi_i the estimated response propensity. Calibration of the respondents solves for a factor v_i that makes the respondents reproduce a reference total (the two-phase / Sarndal-Lundstrom approach).

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class",
                   by = "region")

# household-level nonresponse (whole household responds or not)
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class",
                   by = "region", cluster = "household_id") |>
  prep()
# propensity with cross-fitting (out-of-sample, avoids overfitting)
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "propensity",
                   formula = ~ region + sex, engine = "logit",
                   num_classes = 5, crossfit = 5, crossfit_seed = 1) |>
  prep()

# gradient boosting engine (requires the 'xgboost' package)
if (requireNamespace("xgboost", quietly = TRUE)) {
  weighting_spec(sample_survey, base_weights = pw) |>
    step_nonresponse(respondent = responded, method = "propensity",
                     formula = ~ region + sex + age, engine = "boost",
                     num_classes = 5, crossfit = 5) |>
    prep()
}

# nonresponse by calibration (two-phase): calibrate the respondents to the
# R+NR design-weighted totals of the auxiliaries at that stage, so their
# estimates reproduce the pre-nonresponse ones (Sarndal & Lundstrom 2005)
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "calibration",
                   formula = ~ region + sex) |>
  prep()

Sensitivity of a mean to nonignorable nonresponse or selection

Description

A diagnostic step (it does not change any weight) that gauges how much the weighted mean of a study variable could move if response, or participation in a non-probability sample, depended on the outcome itself beyond the observed auxiliaries. It implements the proxy pattern-mixture model of Andridge and Little (2011): the auxiliaries are reduced to a single proxy (the respondent regression prediction of y), and a single sensitivity parameter phi in ⁠[0, 1]⁠ moves the mechanism from ignorable given the proxy (phi = 0, MAR) to depending only on the outcome (phi = 1). Evaluated over a grid of phi, the adjusted means form an ignorance interval to read alongside the sampling confidence interval; see nr_sensitivity() and the report block.

Usage

step_nr_sensitivity(
  spec,
  y,
  formula,
  respondent = NULL,
  eligible = NULL,
  phi = c(0, 0.25, 0.5, 0.75, 1),
  id = NULL
)

Arguments

Details

The adjusted mean at sensitivity \phi is

\mu(\phi) = \bar{y}_r + (1 - \pi)\,\frac{s_{yr}}{s_{xr}}\,m(\phi)\,(\bar{x}_{nr} - \bar{x}_r),

with slope m(\phi) = \frac{(1-\phi)\rho + \phi}{(1-\phi) + \phi\rho}, so that m(0) = \rho (ignorable given the proxy) and m(1) = 1/\rho.

The proxy correlation rho (the multiple correlation of y on the auxiliaries among respondents) sets how informative the auxiliaries are: a weak proxy widens the ignorance interval (at phi = 1 the slope is 1/rho). The step reads the base design weights, so place it anywhere in the recipe; it needs the nonrespondents still present (a study variable that is NA for them, or an explicit respondent indicator).

Value

the input weighting_spec with this diagnostic step appended.

References

Andridge, R. R. and Little, R. J. A. (2011). Proxy pattern-mixture analysis for survey nonresponse. Journal of Official Statistics 27(2), 153-180.

See Also

nr_sensitivity(), step_assert()

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Adjust base weights by the panel-selection probability (CEPAL ch. XVI)

Description

Divides each incoming weight by ⁠Pr(panel selection)⁠ – the probability that a rotation panel belongs to the combined-wave sample – turning the design weight of the reference wave into the panel base weight of the longitudinal sample (d_base = d1 / Pr). This is the first step of the Verma-Betti-Ghellini longitudinal-weight sequence: combining k waves in a design that keeps g of its G rotation groups in all of them gives Pr = g / G, so the factor is its reciprocal (e.g. 4/3 when 3 of 4 groups persist, 4 when only 1 does).

Usage

step_panel_overlap(spec, prob, id = NULL)

Arguments

Details

prob is a probability in ⁠(0, 1]⁠: either a per-unit column (or expression) or a single constant applied to every active unit. When the data has a rotation_group, get the constant from panel_pr(): prob = panel_pr(pd, c("T1", "T2")). When it does not (e.g. the Chilean ENE), derive the probability yourself (from the design fraction or the cluster) and pass it here. Run this on the longitudinal sample (the intersection of the combined waves), before the attrition and calibration steps.

Value

The input weighting_spec with this step appended to its recipe.

See Also

panel_pr(), panel_design(), step_nonresponse(), step_calibrate()

Examples

# wide longitudinal file of the units in sample in both waves
wide <- panel_merge(
  list(T1 = subset(panel_ine, wave == 1), T2 = subset(panel_ine, wave == 2)),
  by = c("household_id", "person_no"), require = "all")

# rotation group known -> derive Pr(panel selection) from the panel design
pd <- panel_design(panel_ine, unit = c("household_id", "person_no"), wave = "wave",
                   rotation_group = "rotation_group")
weighting_spec(wide, base_weights = pw_T1) |>
  step_panel_overlap(prob = panel_pr(pd, c("1", "2"))) |> prep()

# no rotation group (e.g. Chile ENE) -> pass the probability directly
weighting_spec(wide, base_weights = pw_T1) |>
  step_panel_overlap(prob = 5 / 6) |> prep()

Pseudo-weights for a non-probability sample against a reference

Description

For a non-probability sample (opt-in panel, volunteer or river sample) with no design weights, step_pseudoweight() estimates each unit's participation propensity \hat p against a probability reference_sample() and assigns the pseudo-weight (1 - \hat{p})/\hat{p} (the participation odds; Elliott and Valliant 2017), which inflates each unit to the population so the weights sum to the reference's estimated population size. It stacks the non-probability sample and the reference internally (the participation indicator and the two samples' weights are built for you), fits the propensity, and returns the pseudo-weight on the non-probability units only; the reference is used to train the model and then dropped.

Usage

step_pseudoweight(
  spec,
  reference,
  formula,
  engine = c("logit", "tree", "forest", "boost"),
  num_classes = NULL,
  crossfit = NULL,
  crossfit_seed = NULL,
  id = NULL
)

Arguments

Details

The recipe must be a non-probability spec: weighting_spec(..., nonprob = TRUE). This step is the inverse-propensity (IPW) route; you can instead, or additionally, calibrate to a reference_sample() with step_calibrate() / step_model_calibration() (mass imputation / model-based), and combining both gives the doubly robust estimator.

Value

the input weighting_spec with this step appended.

References

Elliott, M. R. and Valliant, R. (2017). Inference for non-probability samples. Statistical Science 32(2), 249-264.

See Also

reference_sample(), step_calibrate(), step_model_calibration()

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

set.seed(1)
N   <- nrow(population)
# a biased volunteer sample (men over-participate) and a probability reference
vol <- population[rbinom(N, 1, plogis(-2 + 0.9 * (population$sex == "M"))) == 1,
                  c("region", "sex", "income")]
ref <- population[sample(N, 600), c("region", "sex")]
ref$d <- N / 600                                   # its design weights
fit <- weighting_spec(vol, base_weights = NULL, nonprob = TRUE) |>
  step_pseudoweight(reference = reference_sample(ref, "d"),
                    formula = ~ region + sex, engine = "logit") |>
  prep()
# the pseudo-weighted mean corrects the volunteer bias
c(naive = mean(vol$income),
  pseudo = weighted.mean(vol$income, fit$final_weight),
  truth = mean(population$income))

Rescale the weights to a fixed sum

Description

Multiplies the active weights by a single constant so that they add up to a chosen total: either the number of active units (mean weight 1) or an arbitrary number. Use it as a presentation step, when the analysis wants normalized weights rather than population-scale ones.

Usage

step_rescale(spec, to = c("n", "total"), total = NULL, by = NULL, id = NULL)

Arguments

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_rescale(to = "n") |> prep()

Round the final weights

Description

Rounds the weights to a given number of decimals, either unit by unit ("nearest"), with the largest-remainder method ("preserve_total"), which keeps the weighted total exactly, or with the cube method ("balanced"), which keeps the calibrated totals – by domain, not only the grand total – as close as the integer grid allows. Typically the last step of a recipe, after calibration, when the weights have to be delivered as integers or with a fixed number of decimals.

Usage

step_round(
  spec,
  digits = 0L,
  method = c("nearest", "preserve_total", "balanced"),
  by = NULL,
  id = NULL
)

Arguments

Details

The "balanced" method implements balanced rounding by the cube method (Deville and Tille 2004; ECLAC/CEPAL household-survey methodology, chapter 9, section F.2) natively, with no external sampling dependency. It is randomized: call set.seed() before prep() for a reproducible result.

"preserve_total" is randomized too, but only where it has to be. Ties in the fractional part are the rule rather than the exception – a self-weighting design, weights that land on .5, calibrated weights on a grid – and they are broken at random, so the extra unit is allocated without regard to the order of the file. Deciding ties by row order instead moves mass systematically towards whatever the file is sorted by, usually region, while the grand total (the one thing the method promises) stays exactly right and hides it. Weights whose fractional parts are distinct are rounded exactly as before.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

References

Deville J-C, Tille Y (2004). Efficient balanced sampling: the cube method. Biometrika 91(4):893-912.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_round(digits = 0) |> prep()

Within-cluster selection adjustment

Description

Undoes one stage of subsampling inside a cluster: when only some of the eligible units of a household (or dwelling, or area segment) were selected, the selected ones must represent the whole cluster, so their weight is multiplied by the inverse of the within-cluster selection probability. Apply it after the cluster-level eligibility and nonresponse steps and before the person-level nonresponse step.

Usage

step_select_within(
  spec,
  prob = NULL,
  n_eligible = NULL,
  n_selected = NULL,
  id = NULL
)

Arguments

Details

Despite the name, the cluster need not be a household and the unit need not be a person: the step is the generic within-cluster subsampling adjustment. In a multi-stage design it can appear more than once – e.g. dwellings selected within sampled area segments, then persons selected within dwellings – each occurrence undoing one stage of subsampling.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

# simple random selection of one eligible person per household
df <- transform(sample_survey,
                n_elig = ave(person_id, household_id, FUN = length))
weighting_spec(df, base_weights = pw) |>
  step_select_within(n_eligible = n_elig)

# simple random selection of two eligible persons per household
weighting_spec(df, base_weights = pw) |>
  step_select_within(n_eligible = n_elig, n_selected = 2)

Second-phase subsampling (two-phase sampling)

Description

Undoes a second phase of sampling: when a subsample of the first-phase units was drawn for a more expensive follow-up (measuring an outcome, a longer questionnaire), the subsampled units must represent the whole first-phase sample. Their weight is multiplied by the inverse of the phase-2 selection probability, and the not-subsampled units leave the cascade (weight 0).

Usage

step_subsample(spec, selected, prob, psu, id = NULL)

Arguments

Details

The step also records the phase-2 design (the selection probability and the phase-2 sampling unit) so that bootstrap_weights() can reproduce the two-phase variance V = V_1 + V_2: the phase-1 sampling variance plus the expected conditional variance of the phase-2 subsample. A single-phase bootstrap would only capture V_1 and undercover.

The coupling is additive, not multiplicative: the per-unit resampling factor has variance (1-f_1)\pi_2 + (1-\pi_2), the sum of the phase-1 component (1-f_1)\pi_2 (seen through the subsample) and the phase-2 conditional component 1-\pi_2. A naive product of two factors adds a spurious interaction term and is too wide. In practice the factor is drawn once per phase-2 sampling unit from a strictly positive Gamma of that mean and variance, so every replicate weight stays positive (a downstream propensity/GLM step re-runs cleanly). See vignette("two-phase-sampling") for the methodology and its Monte Carlo validation.

The second phase is modelled as Poisson (independent / Bernoulli) selection of the sampling unit nested in the first phase (e.g. households subsampled from a first-phase household sample). This is the general-purpose choice: a Poisson second phase is conservative for, and closely approximates, the without- replacement and stratified subsampling schemes used in practice when the phase-2 sampling fraction is small – which is the usual case, since a costly follow-up subsamples only a fraction of the first phase. The phase-1 sampling fraction f_1 is taken from the fpc argument of bootstrap_weights() and defaults to 0 (negligible, the usual case in household surveys), which reduces the coupling to a single independent per-unit factor of variance 1 – a Gamma of variance 1, i.e. an Exponential(1) (the Bayesian-bootstrap multiplier): valid and strictly positive, but right-skewed, which is part of why replicate-to-replicate variance estimates have heavier tails.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

bootstrap_weights() for the two-phase variance; step_select_within() for within-cluster subsampling that is not a separate sampling phase.

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_trim(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

# households subsampled for a follow-up module, selected with prob p2
df <- transform(sample_survey,
                in_phase2 = as.integer(runif(nrow(sample_survey)) < 0.3),
                p2 = 0.3)
weighting_spec(df, base_weights = pw) |>
  step_subsample(selected = in_phase2, prob = p2, psu = "household_id")

Trim extreme weights against a ratio

Description

Caps the current weights at a multiple of a reference (each unit's base weight, the group median, or an absolute value) and, by default, redistributes the removed mass among the untrimmed units so the weighted total survives the trim. With by, the reference and the cap are computed separately within each subgroup. An optional step that can appear anywhere in the recipe, more than once.

Usage

step_trim(
  spec,
  max_ratio,
  min_ratio = NULL,
  reference = c("base", "median", "value"),
  redistribute = TRUE,
  by = NULL,
  maxit = 50L,
  id = NULL
)

Arguments

Details

There is no standard threshold: max_ratio is an analyst decision, a bias-variance trade-off. Use Kish's design effect (see summary) to judge whether trimming is worth it.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim_calibrated(), step_trim_weights(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_trim(max_ratio = 3, reference = "base")

Trimmed calibration (range-restricted, totals-preserving)

Description

Pulls already-calibrated weights into an absolute interval ⁠[lower, upper]⁠ without breaking the calibration: instead of capping and redistributing, it re-solves a bounded calibration whose targets are the totals the incoming weights already reproduce, optionally with its own band per subgroup through by. It is the only one of the three trimming steps that leaves the calibration totals intact.

Usage

step_trim_calibrated(
  spec,
  formula,
  lower = NULL,
  upper = NULL,
  calfun = c("linear", "raking"),
  by = NULL,
  cluster = NULL,
  equal_within_cluster = FALSE,
  maxit = 100L,
  tol = 1e-07,
  id = NULL
)

Arguments

Details

The absolute-weight bound is imposed as a per-unit factor bound ⁠w_new / w in [lower/w, upper/w]⁠ on top of the incoming weights, using a bounded (range-restricted) calibration with the truncated Deville-Sarndal distances: the range-restricted Euclidean distance (calfun = "linear", the default) or the multiplicative one (calfun = "raking"). Weights inside the range that are not needed to restore the totals stay put; the out-of-range ones saturate at their bound and the rest move as little as possible. If the range is too tight to preserve every total, the totals that cannot be met are relaxed and a warning is raised.

This step is meant to run after a step_calibrate() or a step_model_calibration(): it acts on the active incoming weights (including any negative weights an unbounded linear calibration produced, which it can bring back into ⁠[lower, upper]⁠) and leaves dropped units (weight 0) alone. After a step_model_calibration() it preserves both the known-margin totals and the model-prediction totals: that step saves its prediction columns, which this step appends to formula's design so the trimmed weights keep every total the model calibration reproduced (pass the same x_formula as formula).

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

References

Deville, J.-C. and Sarndal, C.-E. (1992). Calibration estimators in survey sampling. Journal of the American Statistical Association, 87, 376-382. doi:10.2307/2290268. The totals-preserving trimming solves a bounded (range-restricted) calibration with the truncated distances introduced there. Folsom, R. E. and Singh, A. C. (2000). The generalized exponential model for sampling weight calibration for extreme values, nonresponse and poststratification. Proceedings of the ASA Survey Research Methods Section, 598-603, formalises the same range-restricted (generalized exponential) family.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_weights(), step_unknown_eligibility()

Examples

# calibrate, then trim the calibrated weights into [5.5, 13.5] without breaking
# the region/sex totals (the calibrated weights of sample_survey live in ~[5.4, 14])
weighting_spec(sample_survey, base_weights = pw) |>
  step_calibrate(method = "raking",
                 margins = list(region = c(table(population$region)),
                                sex    = c(table(population$sex)))) |>
  step_trim_calibrated(~ region + sex, lower = 5.5, upper = 13.5) |>
  prep()

Automatic weight trimming to an absolute band

Description

Caps the weights into an absolute interval ⁠[lower, upper]⁠ and hands the removed mass back to the units that were not capped, so the weighted total is preserved. This is the step to use when you have not calibrated yet (or will calibrate afterwards) and you want the cutoff chosen from the data rather than argued for: with upper = NULL it picks one by the Tukey far-out fence or by Potter's MSE rule.

Usage

step_trim_weights(
  spec,
  lower = 1,
  upper = NULL,
  method = c("tukey", "potter"),
  redistribute = c("proportional", "uniform"),
  strict = TRUE,
  maxit = 50L,
  kappa = 1,
  id = NULL
)

Arguments

Details

Two things are worth knowing about method = "potter" before reading its cutoff as optimal, because both are assumptions rather than results.

The criterion is the mean squared error of a total, so its two terms are deliberately on different orders: the bias of capping is the trimmed mass, which grows like the sample size, and the variance term is the sum of squared remaining weights, which also grows like the sample size, so bias squared grows like its square. The cutoff therefore rises with the sample size for the same weight distribution – replicating a 300-unit sample to 30,000 moved it from 126 to 183 in one test, capping a third as many units. That is the correct behaviour for a total (the relative bias stays put while the relative variance shrinks, so trimming buys less), not a defect, but it does mean the rule is not a property of the weight distribution alone. The criterion is invariant to the scale of the weights.

The bias it charges is the bias of capping without redistribution, while this step always redistributes: the weighted total is preserved exactly, so the real bias is only the difference between the study variable's mean among the capped units and among the units receiving their mass. The criterion therefore overstates the bias and, other things equal, caps less than the MSE it is named after would. kappa is the handle: set it below 1 to buy back that conservatism.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_unknown_eligibility()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_trim_weights(lower = 1, strict = TRUE) |> prep()

# Potter MSE-optimal cutoff chosen from the data
weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_trim_weights(method = "potter") |> prep()

Unknown-eligibility adjustment

Description

Redistributes the weight of the cases whose eligibility was never resolved onto the resolved cases of the same adjustment cell, so the resolved units stand in for the unresolved share of the frame. Reach for it as the first step of the cascade, while the known-ineligible units are still in the data.

Usage

step_unknown_eligibility(spec, unknown, by = NULL, cluster = NULL, id = NULL)

Arguments

Details

Within each cell c the resolved cases are scaled up to also carry the unresolved ones,

w_i^{\mathrm{out}} = w_i \, \frac{\sum_{j \in c} w_j}{\sum_{j \in c,\, \mathrm{resolved}} w_j},

with every weight on the right the weight entering the step, so the ratio is one number per cell, the cell total is conserved exactly, and the result does not depend on the order in which units are updated.

Value

The input weighting_spec with this step appended to its recipe. The step is recorded only; it is evaluated when prep() is called.

See Also

Other weighting steps: step_assert(), step_calibrate(), step_cre(), step_drop_ineligible(), step_model_calibration(), step_nonresponse(), step_nr_sensitivity(), step_pseudoweight(), step_rescale(), step_round(), step_select_within(), step_subsample(), step_trim(), step_trim_calibrated(), step_trim_weights()

Examples

weighting_spec(sample_survey, base_weights = pw) |>
  step_unknown_eligibility(unknown = unknown_elig, by = "region")

# household-level redistribution (unknown units without roster)
weighting_spec(sample_survey, base_weights = pw) |>
  step_unknown_eligibility(unknown = unknown_elig, by = "region",
                           cluster = "household_id")

Detailed per-step diagnostics

Description

Prints the full audit of an estimated recipe: the stage-by-stage evolution of the weights, then one block per step with that step's own diagnostics table and the design effect before and after it, and finally the R-indicator when the recipe adjusted for nonresponse. Where print() answers "what is in this recipe?", summary() answers "what did each step do, and what did it cost?".

Usage

## S3 method for class 'prepped_weighting_spec'
summary(object, ...)

Arguments

Value

(invisibly) the prepped object.

Examples

fitted <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
summary(fitted)

Gross-flow transition matrix between two panel waves

Description

Weighted cross-tabulation of a categorical state at an earlier wave (from) against a later wave (to) – e.g. the labour-force status of the same people across two periods – computed on the longitudinal weight. This is the gross flow (CEPAL ch. XVII): who moved between which states. boot_transition() adds a per-cell standard error from the ordinary bootstrap.

Usage

transition_matrix(
  data,
  from,
  to,
  weights = NULL,
  states = NULL,
  format = c("row", "col", "joint", "counts")
)

Arguments

Value

a weightflow_transition object holding the matrix.

Why there is a single weight, and not two

Feinberg and Stasny (1983) describe a gross-change table built from the two cross-sectional weights: when w_{k,t-1} \neq w_{k,t}, the smaller weight goes to the (i, j) cell and the difference goes to an "out of population" cell – (Outside, j) or (i, Outside) depending on the sign – on the assumption that the weights differ only because of natural entries to and exits from the target population. ECLAC cites it (ch. XVI, sec. B) as the state of things before a longitudinal weight is built.

This function takes one weight because the package builds that longitudinal weight instead, by the sequence the same chapter prescribes: panel base weight, adjustment for nonresponse in the first period, an explicit definition of the longitudinal population, the attrition adjustment and the final calibration. Once that weight exists the discrepancy Feinberg-Stasny reconstructs no longer arises – who left the population was decided explicitly at step_drop_ineligible() rather than inferred from a difference between two weights. The two-weight construction is the alternative to the longitudinal weight, not a complement to it.

See Also

boot_transition()

Decompose a two-phase variance into V = V1 + V2

Description

For a recipe containing step_subsample(), two_phase_variance() splits the recipe-aware bootstrap variance of an estimate into its first-phase and second-phase components, V = V_1 + V_2. The per-unit coupling factor has variance d = (1-f_1)\pi_2 + (1-\pi_2), the sum of the phase-1 term (1-f_1)\pi_2 and the phase-2 conditional term 1-\pi_2; running the same bootstrap with each term in turn yields the two components.

Usage

two_phase_variance(
  object,
  variable,
  estimator = c("mean", "total"),
  replicates = 500L,
  seed = NULL,
  fpc = NULL
)

Arguments

Details

The share prop_phase2 = V2 / V is an operational read: a large share means the second-phase subsampling drives the uncertainty, so subsampling more would pay off; a small share means a denser (more expensive) subsample would buy little, and the first phase is the binding constraint.

Value

An object of class weightflow_tp_variance: V1, V2, V, the matching standard errors se1, se2, se, and prop_phase2 = V2 / V.

See Also

bootstrap_weights(), step_subsample().

Examples

df <- transform(sample_survey,
                in2 = as.integer(runif(nrow(sample_survey)) < 0.3), p2 = 0.3)
spec <- weighting_spec(df, base_weights = pw) |>
  step_subsample(selected = in2, prob = p2, psu = "household_id")
two_phase_variance(spec, "income", replicates = 100)

Coordinated bootstrap across panel waves

Description

Builds recipe-aware bootstrap replicate weights for two or more waves that are coordinated by PSU: a PSU present in several waves gets the same resampling in all of them, so the sampling covariance between waves – which the overlap of a rotating panel induces – is captured. Feed the result to change_estimate() for the honest variance of a net change. Ignoring the coordination (a separate per-wave bootstrap) over-estimates the variance of change whenever the waves overlap and are correlated.

Usage

wave_bootstrap(
  specs,
  replicates = 500L,
  strata = NULL,
  psu = NULL,
  m = NULL,
  seed = NULL,
  refit_steps = "all",
  resample = c("multinom", "binom"),
  lonely_psu = c("certainty", "collapse"),
  progress = TRUE
)

Arguments

Details

Self-contained: it does not modify bootstrap_weights(); each wave's recipe is re-run through prep() on perturbed base weights, exactly as the single-sample bootstrap does, but with a PSU-coordinated Rao-Wu draw shared across waves.

Composite (CRE) recursion is handled automatically: when a wave's recipe holds a step_cre() whose previous is the preceding wave, that wave's estimated control totals ⁠Zhat*⁠ are re-estimated per replicate from the previous wave's replicate-b weights (not the frozen point weights), so the reported variance includes the sampling variance of Zhat itself. Waves are processed in order and each wave's replicate matrix is carried into the next, giving the full coordinated variance of the recursive estimator. (Assumes each step_cre's previous is the immediately preceding wave in specs; otherwise that step keeps its fixed point previous.)

Value

an object of class weightflow_wave_boot: per-wave point weights, replicate matrices (aligned by replicate index across waves), the per-wave data, and the design.

See Also

change_estimate(), bootstrap_weights(), panel_design()

Examples

t1 <- subset(panel_ine, wave == 1 & disposition == "R")
t2 <- subset(panel_ine, wave == 2 & disposition == "R")
wb <- wave_bootstrap(
  list(T1 = weighting_spec(t1, base_weights = pw),
       T2 = weighting_spec(t2, base_weights = pw)),
  replicates = 200, strata = "stratum", psu = "psu", seed = 1, progress = FALSE)
change_estimate(wb, function(w, d) weighted.mean(d$unemployed, w, na.rm = TRUE))

Extract the carry artifact of a period

Description

The object wave_step() leaves for the next run: the PSU multiplicities that coordinate the next period's bootstrap, the replicate values of the declared estimands, and – only when the recipe needs them – the replicate weights.

Usage

wave_carry(x)

Arguments

Value

a wf_wave_carry, to be saved (saveRDS()) and passed as previous next period.

See Also

wave_step()

Linear combination of an estimand across a chain of periods

Description

wave_step() reports the level of its own period and the net change against each carry it was given – that is, pairwise contrasts. Many published series are not pairwise: a rolling quarter is the average of three consecutive months, a semester-on-semester contrast is c(-1/2, -1/2, 1/2, 1/2), an annual average is rep(1/W, W). wave_contrast() estimates any such combination \psi = a'\theta directly from the carries the chain already wrote to disk, with V(\psi) = a' \Sigma a.

Usage

wave_contrast(carries, estimand, contrast = NULL, level = 0.95)

Arguments

Details

The carries make this possible without keeping the waves in memory: each one stores the R replicate values of every declared estimand, and those replicates are paired across periods because the coordination transferred the PSU multiplicities. Stacking them into an ⁠R x W⁠ matrix and taking its (uncentred-at-R) covariance recovers the full between-period covariance matrix, from which any linear combination follows. With contrast = c(-1, 1) the result reproduces the ⁠$change⁠ row of wave_step() exactly.

Value

a one-row data.frame with estimate, se, V, lo, hi, R_used, the periods and the contrast used, plus the covariance matrix Sigma as an attribute.

See Also

wave_step(), wave_carry(); panel_estimate() does the same on a wave_bootstrap() object, when all waves are held together.

Coordinated delete-one jackknife across panel waves

Description

Deterministic counterpart of wave_bootstrap(): builds delete-one-PSU replicate weights that are coordinated across waves (the same PSU is removed from every wave simultaneously), so the sampling covariance the overlap induces is captured without any randomness. Feed the result to change_estimate() for an honest, reproducible variance of a net change. Meant as the cheap exact oracle to validate wave_bootstrap(); it is also a valid variance estimator in its own right.

Usage

wave_jackknife(
  specs,
  strata = NULL,
  psu = NULL,
  refit_steps = "all",
  lonely_psu = c("certainty", "collapse"),
  progress = TRUE
)

Arguments

Details

Self-contained: each wave's recipe is re-run through prep() on base weights with the deleted PSU zeroed and its stratum mates reweighted by n_h/(n_h-1); it does not modify bootstrap_weights() or touch R/variance.R.

Value

an object of class weightflow_wave_jack accepted by change_estimate(); replicate columns are aligned across waves by the union of PSU ids (a PSU absent from a wave, or alone in its stratum, yields that wave's point weights, i.e. a zero jackknife contribution).

See Also

wave_bootstrap(), change_estimate()

Examples

t1 <- subset(panel_ine, wave == 1 & disposition == "R")
t2 <- subset(panel_ine, wave == 2 & disposition == "R")
wj <- wave_jackknife(
  list(T1 = weighting_spec(t1, base_weights = pw),
       T2 = weighting_spec(t2, base_weights = pw)),
  strata = "stratum", psu = "psu", progress = FALSE)
change_mean(wj, "unemployed")

One period of a coordinated panel bootstrap, chained from the previous ones

Description

Runs the whole weighting recipe of a single period, coordinates its bootstrap with the periods that share sample with it, and returns the cross-sectional weights, the level and net-change estimates, and a compact carry artifact for the next run. It is the production counterpart of wave_bootstrap(): that one needs every wave in memory at once, this one needs only what the earlier runs left behind.

Usage

wave_step(
  spec,
  previous = NULL,
  estimands = NULL,
  by = NULL,
  replicates = 500L,
  strata = NULL,
  psu = NULL,
  period = NULL,
  rotation_map = NULL,
  m = NULL,
  seed = NULL,
  refit_steps = "all",
  resample = c("multinom", "binom"),
  carry = c("auto", "thin", "fat"),
  progress = TRUE
)

Arguments

Details

Coordination transfers the multiplicity of each primary sampling unit, following the Statistics Canada LFS (cat. 71-526-X, section 7.2.2, after Roberts, Kovacevic, Mantel and Phillips 2001): current-period PSUs are paired with earlier ones, common PSUs with themselves, and each inherits its partner's multiplicity. When a stratum's PSU count is unchanged this is a permutation, so the stratum total is preserved exactly and every shared PSU keeps exactly its resampling; the coordination is then exact. When the count changes, multiplicity is added or dropped at random until the stratum total closes, and the coordination is approximate – ⁠$strata⁠ reports, per stratum, the share of replicates that needed no adjustment.

previous is a list of carries, not a single one, because which earlier periods share sample is a property of the rotation design: a 6-consecutive design overlaps at lags 1 to 5 and nowhere else, a 4-(8)-4 design overlaps at lags 1-3 and again at 9-15, and a 1-(3)-1-(3)-1-(3)-1 design has no overlap at lag 1 at all. Each PSU inherits from the most recent carry that holds it, so gaps and returning cohorts need no window parameter. Supply every carry whose period shares sample with this one.

The cross-sectional weights are not touched by any of this: ⁠$weights⁠ is exactly prep(spec)$final_weight.

Value

an object of class wf_wave_step with ⁠$weights⁠ (the publishable cross-sectional weights), ⁠$level⁠, ⁠$change⁠, ⁠$strata⁠ (the per-stratum coordination report), and ⁠$carry⁠ – the artifact to save for the next period.

See Also

wave_carry(), wave_bootstrap(), change_estimate()

Per-unit adjustment factors table

Description

Unrolls the cascade into a data.frame: the weight of every unit at every stage, plus the factor each step applied to it. This is the tidy form of prep()'s ⁠$history⁠, and the starting point for any diagnostic that plot() does not already draw.

Usage

weight_factors(object)

Arguments

Value

data.frame with one weight column per stage and one factor per step.

See Also

Other cascade audit: as_sae_input(), collect_propensities(), collect_step_detail(), collect_weights(), domain_summary(), weighting_alerts()

Examples

fitted <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  prep()
head(weight_factors(fitted))

Quality alerts raised while preparing a recipe

Description

prep() computes non-fatal quality alerts at every stage of the cascade, not only at the end. Each alert names a specific risk, its trigger and a remedy. Alerts are always stored on the prepped object (read them with weighting_alerts() / has_alerts()) and shown in the HTML report; with prep(warn = TRUE) they are also raised as R warnings. Every alert is tagged with the step that produced it, e.g. "[step_calibrate] ...".

Details

This page catalogues the alerts by theme. Thresholds are the prep() defaults unless noted; min_cell_n and max_factor are arguments of prep().

Weight distribution

• Negative calibration weights: an unbounded linear/GREG calibration produced weight(s) below zero. They stay active (the totals are honest), but a negative survey weight is rarely wanted and makes the estimator erratic. Use a bounded distance (calfun = "logit" with bounds) or step_trim_calibrated().

• Sub-one weights: unit(s) with a final weight below 1, so they represent less than themselves. Usually a sign of an over-aggressive adjustment.

• Extreme g-factors: a calibration adjustment factor outside the Deville-Sarndal range ⁠[0.1, 10]⁠. Consider bounds, a ridge penalty, or coarser auxiliaries.

Adjustment cells

• Small cells: an adjustment cell (weighting class or post-stratum) with fewer than min_cell_n (default 30, Kalton and Flores-Cervantes 2003) cases. Collapse cells or switch to raking.

• Empty cells: a cell with no unit to adjust to (no respondents, or all of unknown eligibility); the affected units are set to weight 0.

• Large adjustment factors: a per-cell factor above max_factor (default 2.5), which inflates the variance.

Nonresponse and response propensities

• Tiny propensities: a minimum fitted response propensity below 0.01 blows up the 1/p weights; check the model or trim.

• Collapsed propensity classes: the fitted propensities were nearly constant, so the requested num_classes could not be formed and the class correction did nothing.

• Miscalibrated propensities: a calibration slope far from 1; honest probabilities matter for 1/p. Use num_classes (rank-robust) or revise the model.

• Non-positive nonresponse g-weights: an implied response probability of zero or below; use a bounded distance or other auxiliaries.

Calibration

• Ill-conditioned system: near-collinear auxiliaries (large condition number) make the weights unstable; drop a redundant auxiliary or set a ridge penalty.

• Did not converge: raking, linear or bounded calibration stopped without meeting every target; increase maxit or check the margins are mutually consistent.

• Reconciled control totals: the tidy totals did not all sum to the same population size and were rescaled to a common N (informative, not an error).

Machine-learning adjustments

• Learner without cross-fitting: a flexible learner (tree / forest / boost) run without crossfit can understate the variance through same-sample prediction. Add crossfit to estimate each unit out of sample.

Two-phase (double) sampling

• Few phase-2 sampling units: fewer than 30 units were subsampled at the second phase. The recipe-aware bootstrap resamples the phase-2 variance component (V2) at this level, so few units leave V2 with few degrees of freedom and an unstable phase-2 standard error. Inspect the split with two_phase_variance().

• Tiny phase-2 probability: a minimum phase-2 selection probability below 0.02 expands the subsampled weights sharply (up to 1/pi2), inflating V2. Check the phase-2 design or trim the expanded weights.

Finalising

• Rounded to zero: rounding pushed a small or negative weight to exactly 0, so the unit left the active set and collect_weights(); round to more decimals or resolve those weights first.

See Also

prep(), weighting_alerts(), has_alerts(), step_assert() for a hard quality gate, and vignette("inspecting-auditing") for the full programmatic quality-control workflow.

Conventions shared by every weightflow step

Description

Five things behave the same way in all twelve ⁠step_*()⁠ functions and are easier to learn once than twelve times: what the active set is and how a unit leaves it, why the order of the cascade is not arbitrary, the three different scales on which weight bounds are expressed, which arguments take a bare column name and which take a string, and how to read the diagnostics that prep() stores.

Details

The active set. A unit is active when its weight is finite and non-zero (is.finite(w) & w != 0). A weight of exactly 0 is the "dropped" marker: step_drop_ineligible() sets out-of-scope units to zero, an empty adjustment cell collapses to zero, and a unit that leaves the active set takes no part in any later step and is not returned by collect_weights(). A negative weight, by contrast, is a valid (if unusual) output of unbounded linear/GREG calibration and stays active: it is counted by collect_weights(), the stage funnel and design_effect(), so the reported totals match the weights actually returned. (One caveat: the Kish design effect assumes non-negative weights, so with negatives present its value is inflated – see design_effect().)

Why the order is not arbitrary. Each step multiplies the current weight, i.e. the weight leaving the previous step, so the cascade is read top to bottom. The methodological order of a household survey is: resolve unknown eligibility (step_unknown_eligibility()) while the ineligible units are still present, then drop the ineligible (step_drop_ineligible()), undo any within-cluster subsampling (step_select_within()), adjust for nonresponse (step_nonresponse()), calibrate to population totals (step_calibrate()), and finally trim, round or rescale for delivery. Putting calibration before a trim, or a trim before the nonresponse adjustment, changes the estimator, which is why the steps are explicit rather than inferred.

Three scales for weight bounds. Bounds are expressed on three different scales, and mixing them up is the most common mistake for users coming from survey:

• step_trim() max_ratio / min_ratio are a ratio to a reference (each unit's base weight, the group median, or an absolute value): max_ratio = 3 caps a weight at three times its reference.

• step_trim_weights() and step_trim_calibrated() lower / upper are absolute weights: upper = 400 caps the weight itself at 400.

• step_calibrate() bounds = c(L, U) bound the calibration g-factor w_new / d (the multiplicative adjustment), not the weight.

Bare column name vs string. Arguments that name a single variable to be evaluated on the data take an unquoted (bare) column name or condition: base_weights in weighting_spec(), and respondent, unknown, ineligible, prob, n_eligible, n_selected in the steps. Arguments that name grouping or design structure take character strings: by, cluster, and strata / psu in the variance functions. So it is step_nonresponse(respondent = responded, by = "region") – responded bare, "region" quoted.

The same concept, different argument names. A few concepts are named differently across functions for historical reasons. The correspondence, so you do not have to guess:

• Weight bounds: step_trim() uses max_ratio / min_ratio (a ratio), while step_trim_weights() and step_trim_calibrated() use upper / lower (absolute weights). See "Three scales for weight bounds" above.

• Primary sampling unit: as_svydesign() names it ids, while bootstrap_weights() and jackknife_weights() name it psu. Both mean the cluster identifier.

• Calibration form vs distance: in step_calibrate(), method selects the family ("raking", "poststratify", "linear"); calfun selects the distance function ("linear", "logit", "raking") and applies only when method = "linear". So method = "linear" and calfun = "linear" are not the same knob.

Reading the diagnostics. prep() returns an object that carries the weight at every stage (⁠$history⁠, a named list of vectors), one entry per step (⁠$steps⁠, each with its own ⁠$diagnostics⁠ table and ⁠$alerts⁠), and the recipe-level ⁠$alerts⁠. Rather than read those directly, use summary() for the stage-by-stage audit, weighting_alerts() / has_alerts() for the quality incidents, plot() for the visual cascade, weight_factors() for the per-unit factor table, design_effect() for the Kish design effect, and report_weighting() for the full self-contained HTML report.

Quality alerts recorded while preparing a recipe

Description

weighting_alerts() returns the character vector of quality incidents recorded by prep(), each tagged with the step that produced it. has_alerts() is a convenience predicate. Every incident is captured here, including warnings a step raises internally (for example a calibration that could not meet its constraints), so this is the reliable channel for programmatic quality control even when warnings were suppressed.

Usage

weighting_alerts(object)

has_alerts(object)

Arguments

Value

weighting_alerts(): a character vector (empty if the recipe ran clean). has_alerts(): a single logical.

See Also

Other cascade audit: as_sae_input(), collect_propensities(), collect_step_detail(), collect_weights(), domain_summary(), weight_factors()

Examples

fit <- weighting_spec(sample_survey, base_weights = pw) |>
  step_trim(max_ratio = 3) |>
  prep()
weighting_alerts(fit)
has_alerts(fit)

Start a weighting specification

Description

Opens a weighting recipe on a sample and its design base weights. The object it returns is inert: it holds the data, the name of the base-weight column and an empty list of steps, and computes nothing. Every ⁠step_*()⁠ function takes such an object and returns it with one more step appended; prep() estimates the result.

Usage

weighting_spec(data, base_weights = NULL, nonprob = FALSE)

Arguments

Value

an object of class "weighting_spec".

Examples

rec <- weighting_spec(sample_survey, base_weights = pw)
rec
# a non-probability sample: no design weights, base weight 1
np <- weighting_spec(sample_survey, base_weights = NULL, nonprob = TRUE)

Write a weighting recipe to a YAML file

Description

Serializes the recipe (the weighting method, not the data) to a human-readable YAML file: the base-weight column, the non-probability flag, and every step's id, type and parameters. The file is a versionable metadata artifact you can review in a pull request, archive next to the quality report, or reference from a metadata system. Read it back with read_recipe().

Usage

write_recipe(object, file, timestamp = TRUE)

Arguments

Details

Formulas and captured column expressions are stored as text. A reference_sample() is stored as a descriptor only (its microdata is not metadata), so a step that calibrates or pseudo-weights against a reference must be given that reference again when the recipe is reconstructed. Small control-totals tables (for example a tidy poststratification total) are serialized in full; a data frame larger than 10,000 rows is treated as microdata and rejected (route it through reference_sample()).

A recipe is portable only if its captured expressions reference columns of the data (for example respondent = responded). An expression that referenced objects from your R session (say respondent = id %in% ids_resp) is stored as text but cannot be reconstructed elsewhere, and will error at prep().

Value

the file path, invisibly.

See Also

read_recipe()

Other recipe serialization: read_recipe()

Examples

spec <- weighting_spec(sample_survey, base_weights = pw) |>
  step_nonresponse(respondent = responded, method = "weighting_class", by = "region") |>
  step_calibrate(method = "raking", margins = list(region = c(table(population$region))))
f <- tempfile(fileext = ".yml")
write_recipe(spec, f)

Specify a working model for a study variable y

Description

Declares the working model that step_model_calibration() fits for one study variable: a formula, a learner (a linear/glm model or a machine-learning method such as a regression tree, a random forest or gradient boosting) and, for glm, a family. It builds no model and touches no data – it records the specification that the calibration step will fit on the sample and predict over the population.

Usage

y_model(formula, engine = c("glm", "tree", "forest", "boost"), family = NULL)

Arguments

Value

a model specification list.

Examples

y_model(income ~ age + sex, engine = "glm")