Package 'weightflow'

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")

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()

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()

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()

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,
  level = 0.95,
  ci_type = c("normal", "t", "percentile"),
  df = NULL
)

boot_mean(
  boot,
  variable,
  level = 0.95,
  ci_type = c("normal", "t", "percentile"),
  df = NULL
)

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")

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)

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,
  ci_type = c("normal", "t", "percentile"),
  df = NULL,
  type = c("absolute", "relative"),
  by = NULL
)

change_total(
  wb,
  variable,
  waves = NULL,
  level = 0.95,
  ci_type = c("normal", "t", "percentile"),
  df = NULL,
  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()

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()

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)

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)
}

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")])

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))

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()

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)

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")

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")

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, level = 0.95, ci_type = c("normal", "t"), df = NULL)

jack_mean(jack, variable, level = 0.95, ci_type = c("normal", "t"), df = NULL)

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))

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,
  groups = NULL,
  seed = NULL,
  lonely_psu = c("certainty", "collapse"),
  max_replicates = 5000L,
  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).

With a direct (element) sample the delete-one-unit jackknife is the right estimator, but it needs one replicate per unit: the whole recipe is re-prepped n times and the replicate matrix is ⁠n x n⁠, which is 80 GB at n = 100,000. groups gives the standard device for that case, the delete-a-group jackknife (Kott 2001; Rust and Rao 1996): the sample is partitioned at random into G groups within each stratum and one whole group is deleted at a time, so there are G replicates instead of n. It is the same estimator the engine already computes – the group simply plays the part of the PSU – with $(G-1)/G$ in place of $(n_h-1)/n_h$, and df = G - strata, which is the honest precision of the variance. G between 30 and 100 is the usual range; it is not defaulted, because the choice is the analyst's and the delete-one and delete-a-group variances are not the same number (the second is noisier).

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")

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,
  ci_type = c("normal", "t"),
  df = NULL
)

level_total(
  wb,
  variable,
  wave = NULL,
  level = 0.95,
  ci_type = c("normal", "t"),
  df = NULL
)

Arguments

Value

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

See Also

change_estimate(), level_mean(), level_total()

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()

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")

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)

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,
  ci_type = c("normal", "t"),
  df = NULL
)

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

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

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

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

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)

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)

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)

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.

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.

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

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 – the same assumption made when calibrating to another survey's published totals.

How much that assumption costs depends on what the reference is used for, and the two uses are not equivalent. In step_model_calibration() the reference supplies totals, a smooth functional of it. In step_pseudoweight() every reference unit is a row of the pooled participation fit, so the fitted propensities – and therefore every pseudo-weight – depend on the reference unit by unit. Treating it as fixed there understates the variance even when the reference is several times the size of the sample, and the recipe cannot repair this on its own: building the reference's replicates needs its strata, PSUs and finite-population correction (and its control totals, if it was calibrated), which only the reference's own producer has. So pass them, from the survey's published replicate file or built with bootstrap_weights() on the reference's design; without them bootstrap_weights() and jackknife_weights() say which assumption they fell back on.

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)))

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()

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)

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

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

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()

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()

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()

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",
  rescale_previous = FALSE,
  on_missing_prev = c("carry_backward", "zero"),
  cluster = NULL,
  equal_within_cluster = FALSE,
  calfun = c("linear", "logit", "raking"),
  bounds = NULL,
  rotation_group = NULL,
  n_groups = 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)

# ONE population vector for both waves, as a series of projections would be.
# `Zhat` is a total on the previous wave's scale, so calibrating each wave to its
# own design-weighted total puts two population scales in one system and the
# difference lands on the status estimate (see `rescale_previous`).
Xpop <- colSums(t1$pw * stats::model.matrix(~ sex, data = t1))

# 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 = Xpop, 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 = Xpop,
           alpha = 2/3, status_ref = "inact") |>
  prep()
fit2

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()

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"))

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()

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",
  by = NULL,
  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

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()