Estimate the variance of non-response bias estimate, using the survey package

variance-of-nr-bias.R

# Generate example population and sample ----
  population <- data.frame(
    vax_status = sample(x = c(0,1), prob = c(0.25, 0.75), size = 1000, replace = TRUE),
    response_status = sample(x = c("Respondent", "Nonrespondent"),
                             size = 1000, replace = TRUE, prob = c(0.8, 0.2))
  )
  
  sample_data <- population[sample(x = 1000,size=150),]

# Create a survey design object ----

  library(survey)
  sample_data$pop_size <- 1000
  svy_design <- survey::svydesign(
    data = sample_data,
    ids = ~ 1, fpc = ~ pop_size
  )

# Estimate means for respondents and for overall population ----
  svy_design <- transform(
    svy_design,
    is_respondent = ifelse(response_status == "Respondent", 1, 0),
    any_case = 1,
    respondent_vax_status = vax_status * is_respondent
  )
  
  estimated_means <- svyratio(
    numerator = ~ respondent_vax_status + vax_status,
    denominator = ~ is_respondent + any_case,
    design = svy_design, covmat = TRUE
  )
  
  domain_means <- coef(estimated_means)[c(1,4)]
  print(domain_means)
  
# Estimate non-response bias ----
  
  svycontrast(
    stat = estimated_means,
    contrast = list('bias_estimate' = c(1,0,0,-1))
  )