if (FALSE) {
# Load example data for primary survey ----
suppressPackageStartupMessages(library(survey))
data(api)
primary_survey <- svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc) |>
as.svrepdesign(type = "JK1")
# Load example data for control survey ----
control_survey <- svydesign(id = ~ 1, fpc = ~fpc, data = apisrs) |>
as.svrepdesign(type = "JK1")
# Estimate control totals ----
estimated_controls <- svytotal(x = ~ stype + enroll,
design = control_survey)
control_point_estimates <- coef(estimated_controls)
control_vcov_estimate <- vcov(estimated_controls)
# Calibrate totals for one categorical variable and one numeric ----
calibrated_rep_design <- calibrate_to_estimate(
rep_design = primary_survey,
estimate = control_point_estimates,
vcov_estimate = control_vcov_estimate,
cal_formula = ~ stype + enroll
)
# Inspect estimates before and after calibration ----
##_ For the calibration variables, estimates and standard errors
##_ from calibrated design will match those of the control survey
svytotal(x = ~ stype + enroll, design = primary_survey)
svytotal(x = ~ stype + enroll, design = control_survey)
svytotal(x = ~ stype + enroll, design = calibrated_rep_design)
##_ Estimates from other variables will be changed as well
svymean(x = ~ api00 + api99, design = primary_survey)
svymean(x = ~ api00 + api99, design = control_survey)
svymean(x = ~ api00 + api99, design = calibrated_rep_design)
# Inspect weights before and after calibration ----
summarize_rep_weights(primary_survey, type = 'overall')
summarize_rep_weights(calibrated_rep_design, type = 'overall')
# For reproducibility, specify which columns are randomly selected for Fuller method ----
column_selection <- calibrated_rep_design$col_selection
print(column_selection)
calibrated_rep_design <- calibrate_to_estimate(
rep_design = primary_survey,
estimate = control_point_estimates,
vcov_estimate = control_vcov_estimate,
cal_formula = ~ stype + enroll,
col_selection = column_selection
)
}
Run the code above in your browser using DataLab