survey.R
In senseweight: Sensitivity Analysis for Weighted Estimators

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)

## ----setup, message = FALSE, warning=FALSE------------------------------------
library(senseweight)
library(survey)

## -----------------------------------------------------------------------------
data(poll.data)
poll.data |> head()

## ----warning=FALSE, message=FALSE---------------------------------------------
poll_srs <- svydesign(ids = ~ 1, data = poll.data)

## -----------------------------------------------------------------------------
pop_targets = c(1, 0.212, 0.264, 0.236, 0.310, 
                0.114, 0.360, 0.528, 0.114, 
                0.021, 0.034, 0.805, 
                0.266, 0.075, 0.312, 0.349)
#Match covariate names in polling data 
names(pop_targets) = model.matrix(~.-Y, data = poll.data) |> colnames()
print(pop_targets)

## -----------------------------------------------------------------------------
#Set up raking formula:
formula_rake <- ~ age_buckets + educ + gender + race + pid + bornagain

#PERFORM RAKING:
model_rake <- calibrate(
  design = poll_srs,
  formula = formula_rake,
  population = pop_targets,
  calfun = "raking",
  force = TRUE
)


rake_results <- svydesign( ~ 1, data = poll.data, weights = stats::weights(model_rake))
#Estimate from raking results:
weights = stats::weights(rake_results) * nrow(model_rake)

unweighted_estimate = svymean(~ Y, poll_srs, na.rm = TRUE)
weighted_estimate = svymean(~ Y, model_rake, na.rm = TRUE)

## -----------------------------------------------------------------------------
print(unweighted_estimate)

## -----------------------------------------------------------------------------
print(weighted_estimate)

## -----------------------------------------------------------------------------
summarize_sensitivity(estimand = 'Survey',
                      Y = poll.data$Y,
                      weights = weights,
                      svy_srs = unweighted_estimate, 
                      svy_wt = weighted_estimate,
                      b_star = 0.5)


## -----------------------------------------------------------------------------
robustness_value(estimate = as.numeric(weighted_estimate[1]),
                 b_star = 0.5,
                 sigma2 = var(poll.data$Y), 
                 weights = weights)

## -----------------------------------------------------------------------------
benchmark_survey('educ', 
                 formula = formula_rake,
                 weights = weights,
                 population_targets = pop_targets,
                 sample_svy = poll_srs,
                 Y = poll.data$Y)

## -----------------------------------------------------------------------------
covariates = c("age_buckets", "educ", "gender", "race",
               "educ", "pid", "bornagain")

benchmark_results = run_benchmarking(estimate = as.numeric(weighted_estimate[1]),
                 RV = 0.05,
                 formula = formula_rake,
                 weights = weights,
                 Y = poll.data$Y,
                 sample_svy = poll_srs,
                 population_targets = pop_targets,
                 estimand= "Survey")
print(benchmark_results)


## ----fig.width=6.5, fig.height=5----------------------------------------------
contour_plot(varW = var(weights), 
             sigma2 = var(poll.data$Y),
             killer_confounder = 0.5, 
             df_benchmark = benchmark_results,
             shade = TRUE, 
             label_size = 4)