knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(weightflow)
Calibration adjusts the weights so that the weighted sample reproduces known
population totals of auxiliary variables. In practice those totals almost always
arrive as a table: a census cross-tabulation, a projection, a spreadsheet
read into R. step_calibrate() accepts totals in two shapes:
margins as a named
list, or a totals vector aligned with the model matrix); andcount. Any column that is not the
counts column is treated as a post-stratification variable.This vignette shows both, side by side, for the three calibration methods, so you can pick whichever matches the data you already have.
We use the bundled example data throughout.
data(population) data(sample_survey)
Post-stratification calibrates to the joint distribution of one or more categorical variables (the cells of their cross-classification).
The tidy way: a data frame with the category column and a counts column. Here we
build the population counts with table() (its data-frame form has a Freq
column), then pass count = "Freq".
region_totals <- as.data.frame(table(region = population$region)) region_totals ps <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "poststratify", totals = region_totals, count = "Freq") |> prep() sum(collect_weights(ps)$.weight) # sums to the population size N
The classic equivalent uses margins, a named list of named vectors:
ps_classic <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "poststratify", margins = list(region = c(table(population$region)))) |> prep() sum(collect_weights(ps_classic)$.weight)
Both produce the same weights.
This is where the tidy format helps most. If your table has several category columns, weightflow crosses them automatically to form the post-strata. You do not need to build a single collapsed cell variable by hand.
rs_totals <- as.data.frame(table(region = population$region, sex = population$sex)) head(rs_totals) ps_cross <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "poststratify", totals = rs_totals, count = "Freq") |> prep() sum(collect_weights(ps_cross)$.weight)
If you had pre-combined the categories into a single column (for example a
"North-F" cell label), that is simply the one-column case above, so it is
handled by the same rule with no extra work.
Raking calibrates to several independent margins (each variable separately, iterated), which is what you use when you do not have the full cross-tabulation, only the marginal totals. In the tidy format you pass 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)) rk <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "raking", totals = list(m_region, m_sex), count = "Freq") |> prep() sum(collect_weights(rk)$.weight)
The classic equivalent uses margins:
rk_classic <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "raking", margins = list(region = c(table(population$region)), sex = c(table(population$sex)))) |> prep() sum(collect_weights(rk_classic)$.weight)
Because raking is iterative, weightflow warns you if the margins are mutually inconsistent (they do not all sum to the same population size) or if the iteration does not converge, rather than silently returning weights that do not satisfy the margins.
Linear (GREG) calibration handles categorical and continuous auxiliaries together, through a model formula.
Linear calibration is built on model.matrix(), and traditionally the totals
must be supplied in exactly that internal shape: a vector that includes the
intercept (the population size $N$) and, for each factor, drops one reference
category absorbed into the intercept, using treatment-contrast column names.
Knowing that a category is silently omitted, and reproducing the exact model-matrix names, is a common source of mistakes. Providing totals this way is the norm in established survey-calibration tools, and it is precisely the part that trips people up.
# The classic model-matrix vector: intercept = N, and region *without* its # reference level (the first, "North"), with model.matrix column names. 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")) lin_classic <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "linear", formula = ~ region + sex, totals = pop_tot) |> prep() sum(collect_weights(lin_classic)$.weight)
In the tidy format you give the complete categories (no reference dropped) and a single number for each continuous total, as a named list matching the formula terms. weightflow builds the model-matrix totals internally, including the intercept and the omitted reference category, so you never handle them.
lin_tidy <- weighting_spec(sample_survey, base_weights = pw) |> step_calibrate(method = "linear", formula = ~ region + sex, totals = list(region = m_region, sex = m_sex), count = "Freq") |> prep() sum(collect_weights(lin_tidy)$.weight)
This yields the same weights as the model-matrix vector above, but you supplied the totals in the natural, complete form.
For a continuous auxiliary, give its population total as a single number. Here
we add income; because income is observed only for respondents in this
example, we calibrate the respondent subsample.
resp <- subset(sample_survey, responded == 1) lin_mixed <- 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() # the calibrated weights reproduce every target X <- model.matrix(~ region + sex + income, data = resp) colSums(collect_weights(lin_mixed)$.weight * X)
Ridge (penalized) calibration works with the tidy format too: add penalty as
usual. Under ridge the achieved totals are deliberately not exact, so weightflow
reports the deviation instead of warning about it.
Sometimes the benchmarks are known by domain and you want the weights to
reproduce them within each domain, not only overall. Pass by = with the
domain (partition) column: step_calibrate() then calibrates independently
inside each domain, each to its own totals. The tidy totals carry the domain as
an extra column, and the domain variable does not go in the formula or the
margins (it is the partition). This is also called partitioned calibration.
Domain calibration earns its keep when, within each domain, you calibrate to
more than one target: two or more raked margins, or a continuous total. (With
a single margin per domain it would collapse to an ordinary global
cross-tabulation, so by = would add nothing.) Here we rake to two margins, sex
and age group, within each region: the region-specific sex × age structure
is not assumed to be the same across regions, which is exactly what a global cross
could not express with these marginal-only benchmarks.
# benchmarks known BY REGION: a sex margin and an age-group margin per region. 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)) dom <- weighting_spec(samp, base_weights = pw) |> step_calibrate(method = "raking", totals = list(sex_by_region, age_by_region), count = "Freq", by = "region") |> prep() # within each region the weights reproduce BOTH margins w <- dom$final_weight round(xtabs(w ~ region + sex, data = cbind(samp, w = w))) round(xtabs(w ~ region + age_grp, data = cbind(samp, w = w)))
Mixing a categorical and a continuous auxiliary by domain works too. The
continuous total is a data frame domain, value (one total per domain). Here we
use the exponential (raking) distance so the weights stay positive:
inc_by_region <- aggregate(income ~ region, population, sum) # region, income resp <- subset(sample_survey, responded == 1) lin_dom <- weighting_spec(resp, base_weights = pw) |> step_calibrate(method = "linear", formula = ~ sex + income, totals = list(sex = sex_by_region, income = inc_by_region), count = "Freq", calfun = "raking", by = "region") |> prep() # the income total is reproduced within each region w <- lin_dom$final_weight got <- tapply(w * resp$income, resp$region, sum) cbind(calibrated = round(got), benchmark = inc_by_region$income[match(names(got), inc_by_region$region)])
by = composes with calfun, bounds, penalty and the integrative option,
all applied within each domain. With by = NULL (the default) calibration is
global, as in the sections above.
step_model_calibration() (model-assisted, Wu & Sitter 2001) fits a working
model for each study variable, predicts it over the population, and calibrates so
the weighted sample reproduces two kinds of target at once: the population total
of each prediction (the model-assisted part) and the totals of a set of
consistency auxiliaries given by x_formula (exactly as in linear
calibration).
By default those consistency totals are read from the population frame. But the
control totals of the auxiliaries often come from another source: an official
published total, or a variable that is not even in the frame. For that,
x_totals accepts the same two shapes as linear calibration (tidy named
list, or classic model-matrix vector), paired with count for the tidy form. The
model predictors and the consistency auxiliaries are independent: a variable can
drive the model without being a control total, and a control total need not enter
the model.
resp <- subset(sample_survey, responded == 1) mc <- weighting_spec(resp, base_weights = pw) |> step_model_calibration( x_formula = ~ region + age, # consistency block models = list(income = y_model(income ~ age + sex, # model block engine = "glm")), population = population, # used for prediction x_totals = list(region = m_region, age = sum(population$age)), count = "Freq") |> prep() # both blocks are reproduced: the X totals and the model prediction total mc$steps[[1]]$diagnostics
Here region is a tidy data frame of category counts and age is a single
external number, just like the mixed linear example. population is still
required because the model must predict over every population unit, but the
consistency auxiliaries only need to exist in the sample: when x_totals is
given, weightflow does not read them from the frame, so age (or a variable
absent from population altogether) can be controlled from an outside total.
Leaving x_totals = NULL keeps the earlier behaviour, reading the X totals from
population.
The tidy format is not only more convenient; it checks the totals against the sample and explains problems in survey terms.
NA) is an error. Calibration
requires every unit to have a value for each auxiliary; impute first, or use a
complete frame variable.margins/totals inputs remain fully supported; existing code
keeps working unchanged.In all cases the calibration itself is identical: the two shapes are just different ways to hand weightflow the same population totals.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.