Example of FPBB with the Louisville Vaccination Survey

finite-population-bayesian-bootstrap.R

suppressPackageStartupMessages({
  library(survey)
  library(svrep)
  library(polyapost)
})

set.seed(1999)

# Load example survey data ----
  data("lou_vax_survey", package = 'svrep')
  lou_vax_survey_design <- svydesign(ids = ~ 1, weights = ~ SAMPLING_WEIGHT,
                                     data = lou_vax_survey)

# Create bootstrap replicate weights, represented with a survey design object ----
  M <- sum(lou_vax_survey$variables$SAMPLING_WEIGHT) # Population size
  L <- 30 # Number of usual bootstrap replicates
  B <- 5 # Number of synthetic populations per bootstrap replicate
  
  ## Create Rao-Wu bootstrap replicates
  lou_vax_boot <- lou_vax_survey_design |>
    as.svrepdesign(type = "subbootstrap",
                   replicates = L,
                   mse = TRUE)
  
  boot_rep_weights <- weights(lou_vax_boot, type = "analysis")
  
  ## For each replicate:
  
    ### Identify number of observed population elements in the replicate
    nonzero_indices <- apply(X = boot_rep_weights, MARGIN = 2,
                             FUN = function(wts) {
                               wts > 0
                             })
    boot_n_obs <- colSums(nonzero_indices)
    
    ### Identify number of population elements *NOT* in the replicate
    boot_n_unobserved <- M - boot_n_obs
    
    ### Rescale each replicate's weights to sum to number of ultimate sampling units
    boot_weight_sums <- apply(boot_rep_weights, MARGIN = 2,
                              FUN = sum)
  
  
    rescaled_rep_weights <- sapply(X = seq_len(ncol(boot_rep_weights)),
                                   function(rep_index) {
                                     wts <- as.vector(boot_rep_weights[,rep_index])
                                     rescaling_factor <- boot_n_obs[rep_index]/boot_weight_sums[rep_index]
                                     wts <- rescaling_factor * wts
                                     return(wts)
                                   })

  ## Create Finite-Population Bayesian Bootstrap (FPBB) replicate weights
  ## For the l-th Rao-Wu bootstrap sample,
  ## create b bootstrap population using resampling from the posterior predictive distribution
  ## of elements in the nonsampled population, given the elements in the l-th bootstrap sample.
  ## Operationally, we draw a Polya sample from an urn with unequal selection probabilities
  ## and add that Polya sample to the original bootstrap sample.
  fpbb_weights <- matrix(data = 0, nrow = nrow(lou_vax_boot), ncol = L * B)
  lb <- 1L
  for (l in seq(L)) {
    for (b in seq(B)) {
      sampled_nonzero_indices <- polyapost::wtpolyap(
        ysamp = seq(boot_n_obs[l]),
        wts = rescaled_rep_weights[nonzero_indices[,l],l],
        k = boot_n_unobserved[l]
      )
      fpbb_weights[nonzero_indices[,l],lb] <- as.vector(table(sampled_nonzero_indices))
      lb <- lb + 1L
    }
  }
  
  lou_vax_fpbb_design <- survey::svrepdesign(
    data = lou_vax_survey,
    weights = lou_vax_survey[['SAMPLING_WEIGHT']],
    repweights = fpbb_weights,
    type = 'other',
    scale = 1/(L * B), # Need to check that this scale factor is correct
    rscales = 1
  )
  
  lou_vax_fpbb_design |> svymean(x = ~ VAX_STATUS, na.rm = TRUE)
  lou_vax_boot |> svymean(x = ~ VAX_STATUS, na.rm = TRUE)