Bootstrap Methods for Surveys

Quadratic Forms

Let \(v( \hat{T_y})\) be the textbook variance estimator for an estimated population total \(\hat{T}_y\) of some variable \(y\). The base weight for case \(i\) in our sample is \(w_i\), and we let \(\breve{y}_i\) denote the weighted value \(w_iy_i\).

Suppose we can represent our textbook variance estimator as a quadratic form: \(v(\hat{T}_y) = \breve{y}\Sigma\breve{y}^T\), for some \(n \times n\) matrix \(\Sigma\). Our only constraint on \(\Sigma\) is that, for our sample, it must be symmetric and positive semi-definite (in other words, it should never lead to a negative variance estimate, no matter what the value of \(\breve{y}\) is).

For example, the popular Horvitz-Thompson estimator based on first-order inclusion probabilities \(\pi_k\) and second-order inclusion probabilities \(\pi_{kl}\) can be represented as a positive semi-definite matrix with entries \((1-\pi_k)\) along the main diagonal and entries \((1 - \frac{\pi_k \pi_l}{\pi_{kl}})\) everywhere else. An illustration for a sample with \(n=3\) is shown below:

\[ \Sigma_{HT} = \begin{bmatrix} (1-\pi_1) & (1 - \frac{\pi_1 \pi_2}{\pi_{12}}) & (1 - \frac{\pi_1 \pi_3}{\pi_{13}}) \\ (1 - \frac{\pi_2 \pi_1}{\pi_{21}}) & (1 - \pi_2) & (1 - \frac{\pi_2 \pi_3}{\pi_{23}}) \\ (1 - \frac{\pi_3 \pi_1}{\pi_{31}}) & (1 - \frac{\pi_3 \pi_2}{\pi_{32}}) & (1 - \pi_3) \end{bmatrix} \]

As another example, the successive-difference variance estimator for a systematic sample can be represented as a positive semi-definite matrix whose diagonal entries are all \(1\), whose superdiagonal and subdiagonal entries are all \(-1/2\), and whose top right and bottom left entries are \(-1/2\) (Ash 2014). An illustration for a sample with \(n=5\) is shown below:

\[ \Sigma_{SD2} = \begin{bmatrix} 1 & -1/2 & 0 & 0 & -1/2\\ -1/2 & 1 & -1/2 & 0 & 0 \\ 0 & -1/2 & 1 & -1/2 & 0 \\ 0 & 0 & -1/2 & 1& -1/2 \\ -1/2 & 0 & 0 & -1/2 & 1 \end{bmatrix} \]

To obtain a quadratic form matrix for a variance estimator, we can use the function make_quad_form_matrix(), which takes as inputs the name of a variance estimator and relevant survey design information. For example, the following code produces the quadratic form matrix of the “SD2” variance estimator we saw earlier.

make_quad_form_matrix(
    variance_estimator = "SD2",
    cluster_ids = c(1,2,3,4,5) |> data.frame(),
    strata_ids = c(1,1,1,1,1) |> data.frame(),
    sort_order = c(1,2,3,4,5)
)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]  1.0 -0.5  0.0  0.0 -0.5
#> [2,] -0.5  1.0 -0.5  0.0  0.0
#> [3,]  0.0 -0.5  1.0 -0.5  0.0
#> [4,]  0.0  0.0 -0.5  1.0 -0.5
#> [5,] -0.5  0.0  0.0 -0.5  1.0

In the following example, we use this method to estimate the variance for a stratified systematic sample of U.S. public libraries. First, we create a quadratic form matrix to represent the SD2 successive-difference estimator.

# Load an example dataset of a stratified systematic sample
data('library_stsys_sample', package = 'svrep')

# Represent the SD2 successive-difference estimator as a quadratic form,
# and obtain the matrix of that quadratic form
sd2_quad_form <- make_quad_form_matrix(
  variance_estimator = 'SD2',
  cluster_ids = library_stsys_sample |> select(FSCSKEY),
  strata_ids = library_stsys_sample |> select(SAMPLING_STRATUM),
  strata_pop_sizes = library_stsys_sample |> select(STRATUM_POP_SIZE),
  sort_order = library_stsys_sample |> pull("SAMPLING_SORT_ORDER")
)

class(sd2_quad_form)
#> [1] "matrix" "array"
dim(sd2_quad_form)
#> [1] 219 219

Next, we estimate the sampling variance of the estimated total for the TOTCIR variable using the quadratic form.

# Obtain weighted values
wtd_y <- as.matrix(library_stsys_sample[['TOTCIR']] /
                     library_stsys_sample[['SAMPLING_PROB']])
wtd_y[is.na(wtd_y)] <- 0

# Obtain point estimate for a population total
point_estimate <- sum(wtd_y)

# Obtain the variance estimate using the quadratic form
variance_estimate <- t(wtd_y) %*% sd2_quad_form %*% wtd_y

# Summarize results
sprintf("Estimate: %s", point_estimate)
#> [1] "Estimate: 2195846830.39647"
sprintf("Standard Error: %s", sqrt(variance_estimate[1,1]))
#> [1] "Standard Error: 564706867.832222"
sprintf("Coefficient of Variation: %s",
        round(sqrt(variance_estimate[1,1])/point_estimate, 2))
#> [1] "Coefficient of Variation: 0.26"

Forming Adjustment Factors

Our goal is to form \(B\) sets of bootstrap weights, where the \(b\)-th set of bootstrap weights is a vector of length \(n\) denoted \(\mathbf{a}^{(b)}\), whose \(k\)-th value is denoted \(a_k^{(b)}\). This gives us \(B\) replicate estimates of the population total, \(\hat{T}_y^{*(b)}=\sum_{k \in s} a_k^{(b)} \breve{y}_k\), for \(b=1, \ldots B\), from which we can easily calculate an estimate of the sampling variance.

\[ v_B\left(\hat{T}_y\right)=\frac{\sum_{b=1}^B\left(\hat{T}_y^{*(b)}-\hat{T}_y\right)^2}{B} \]

We can write this bootstrap variance estimator as a quadratic form:

\[ \begin{aligned} v_B\left(\hat{T}_y\right) &=\mathbf{\breve{y}}^{\prime}\Sigma_B \mathbf{\breve{y}} \\ \textit{where}& \\ \boldsymbol{\Sigma}_B &= \frac{\sum_{b=1}^B\left(\mathbf{a}^{(b)}-\mathbf{1}_n\right)\left(\mathbf{a}^{(b)}-\mathbf{1}_n\right)^{\prime}}{B} \end{aligned} \]

Note that if the vector of adjustment factors \(\mathbf{a}^{(b)}\) has expectation \(\mathbf{1}_n\) and variance-covariance matrix \(\boldsymbol{\Sigma}\), then we have the bootstrap expectation \(E_{*}\left( \boldsymbol{\Sigma}_B \right) = \boldsymbol{\Sigma}\). Since the bootstrap process takes the sample values \(\breve{y}\) as fixed, the bootstrap expectation of the variance estimator is \(E_{*} \left( \mathbf{\breve{y}}^{\prime}\Sigma_B \mathbf{\breve{y}}\right)= \mathbf{\breve{y}}^{\prime}\Sigma \mathbf{\breve{y}}\). Thus, we can produce a bootstrap variance estimator with the same expectation as the textbook variance estimator simply by randomly generating \(\mathbf{a}^{(b)}\) from a distribution with the following two conditions:

Condition 1: \(\quad \mathbf{E}_*(\mathbf{a})=\mathbf{1}_n\)

Condition 2: \(\quad \mathbf{E}_*\left(\mathbf{a}-\mathbf{1}_n\right)\left(\mathbf{a}-\mathbf{1}_n\right)^{\prime}=\mathbf{\Sigma}\)

The simplest, general way to generate the adjustment factors is to simulate from a multivariate normal distribution \(\mathbf{a} \sim MVN(\mathbf{1}_n, \boldsymbol{\Sigma})\), which is the method used in this package. However, this method can lead to negative adjustment factors and hence negative bootstrap weights, which–while perfectly valid for variance estimation–may be undesirable from a practical point of view. Thus, in the following subsection, we describe one method for adjusting the replicate factors so that they are nonnegative and \(E_{*} \left( \mathbf{\breve{y}}^{\prime}\Sigma_B \mathbf{\breve{y}}\right) =\mathbf{\breve{y}}^{\prime}\Sigma \mathbf{\breve{y}}\).

Adjusting Generalized Survey Bootstrap Replicates to Avoid Negative Weights

Let \(\mathbf{A} = \left[ \mathbf{a}^{(1)} \cdots \mathbf{a}^{(b)} \cdots \mathbf{a}^{(B)} \right]\) denote the \((n \times B)\) matrix of bootstrap adjustment factors. To eliminate negative adjustment factors, Beaumont and Patak (2012) propose forming a rescaled matrix of nonnegative replicate factors \(\mathbf{A}^S\) by rescaling each adjustment factor \(a_k^{(b)}\) as follows:

\[ \begin{aligned} a_k^{S,(b)} &= \frac{a_k^{(b)} + \tau - 1}{\tau} \\ \textit{where } \tau &\geq 1 - a_k^{(b)} \geq 1 \\ &\textit{for all }k \textit{ in } \left\{ 1,\ldots,n \right\} \\ &\textit{and all }b \textit{ in } \left\{1, \ldots, B\right\} \\ \end{aligned} \]

The value of \(\tau\) can be set based on the realized adjustment factor matrix \(\mathbf{A}\) or by choosing \(\tau\) prior to generating the adjustment factor matrix \(\mathbf{A}\) so that \(\tau\) is likely to be large enough to prevent negative bootstrap weights.

If the adjustment factors are rescaled in this manner, it is important to adjust the scale factor used in estimating the variance with the bootstrap replicates, which becomes \(\frac{\tau^2}{B}\) instead of \(\frac{1}{B}\).

\[ \begin{aligned} \textbf{Prior to rescaling: } v_B\left(\hat{T}_y\right) &= \frac{1}{B}\sum_{b=1}^B\left(\hat{T}_y^{*(b)}-\hat{T}_y\right)^2 \\ \textbf{After rescaling: } v_B\left(\hat{T}_y\right) &= \frac{\tau^2}{B}\sum_{b=1}^B\left(\hat{T}_y^{S*(b)}-\hat{T}_y\right)^2 \\ \end{aligned} \]

When sharing a dataset that uses rescaled weights from a generalized survey bootstrap, the documentation for the dataset should instruct the user to use replication scale factor \(\frac{\tau^2}{B}\) rather than \(\frac{1}{B}\) when estimating sampling variances.

Option 1: Convert an existing design to a generalized bootstrap design

The simplest method is to convert an existing survey design object to a generalized bootstrap design.

With this approach, we create a survey design object using the svydesign() function. This allows us to represent information about stratification and clustering (potentially at multiple stages), as well as information about finite population corrections.

# Load example data from stratified systematic sample
data('library_stsys_sample', package = 'svrep')

# First, ensure data are sorted in same order as was used in sampling
library_stsys_sample <- library_stsys_sample[
  order(library_stsys_sample$SAMPLING_SORT_ORDER),
]

# Create a survey design object
design_obj <- svydesign(
  data = library_stsys_sample,
  strata = ~ SAMPLING_STRATUM,
  ids = ~ 1,
  fpc = ~ STRATUM_POP_SIZE
)

Next, we convert the survey design object to a replicate design using the function as_gen_boot_design(). The function argument variance_estimator allows us to specify the name of a variance estimator to use as the basis for creating replicate weights.

# Convert to generalized bootstrap replicate design
gen_boot_design_sd2 <- as_gen_boot_design(
  design = design_obj,
  variance_estimator = "SD2",
  replicates = 2000
)
#> For `variance_estimator='SD2', assumes rows of data are sorted in the same order used in sampling.

# Estimate sampling variances
svymean(x = ~ TOTSTAFF, na.rm = TRUE, design = gen_boot_design_sd2)
#>            mean     SE
#> TOTSTAFF 19.756 4.2644

For a PPS design that uses the Horvitz-Thompson or Yates-Grundy estimator, we can create a generalized bootstrap estimator with the same expectation. In the example below, we create a PPS design with the ‘survey’ package and convert it to a generalied bootstrap design.

# Load example data of a PPS survey of counties and states
   data('election', package = 'survey')

# Create survey design object
   pps_design_ht <- svydesign(
     data = election_pps,
     id = ~1, fpc = ~p,
     pps = ppsmat(election_jointprob),
     variance = "HT"
   )

  
# Convert to generalized bootstrap replicate design
gen_boot_design_ht <- pps_design_ht |>
  as_gen_boot_design(variance_estimator = "Horvitz-Thompson",
                     replicates = 5000, tau = "auto")

# Compare sampling variances from bootstrap vs. Horvitz-Thompson estimator
svytotal(x = ~ Bush + Kerry, design = pps_design_ht)
#>          total      SE
#> Bush  64518472 2604404
#> Kerry 51202102 2523712
svytotal(x = ~ Bush + Kerry, design = gen_boot_design_ht)
#>          total      SE
#> Bush  64518472 2627361
#> Kerry 51202102 2499754

We can also use the generalized bootstrap for designs that use multistage, stratified simple random sampling without replacement.

library(dplyr) # For data manipulation

# Create a multistage survey design
  multistage_design <- svydesign(
    data = library_multistage_sample |>
      mutate(Weight = 1/SAMPLING_PROB),
    ids = ~ PSU_ID + SSU_ID,
    fpc = ~ PSU_POP_SIZE + SSU_POP_SIZE,
    weights = ~ Weight
  )

# Convert to a generalized bootstrap design
  multistage_boot_design <- as_gen_boot_design(
    design = multistage_design,
    variance_estimator = "Stratified Multistage SRS"
  )

# Compare variance estimates
  svytotal(x = ~ TOTCIR, na.rm = TRUE, design = multistage_design)
#>             total        SE
#> TOTCIR 1634739229 251589313
  svytotal(x = ~ TOTCIR, na.rm = TRUE, design = multistage_boot_design)
#>             total        SE
#> TOTCIR 1634739229 236364980

Unless specified otherwise, as_gen_boot_design() automatically selects a rescaling value \(\tau\) to use for eliminating negative adjustment factors. The scale attribute of the resulting replicate survey design object is thus set to equal \(\tau^2/B\). The specific value of \(\tau\) can be retrieved from the replicate design object, as follows.

# View overall scale factor
overall_scale_factor <- gen_boot_design_ht$scale
print(overall_scale_factor)
#> [1] 0.00464648

# Check that the scale factor was calculated correctly
tau <- gen_boot_design_ht$tau
print(tau)
#> [1] 4.82
B <- ncol(gen_boot_design_ht$repweights)
print(B)
#> [1] 5000

print( (tau^2) / B )
#> [1] 0.00464648

Option 2: Create the quadratic form matrix and then use it to create bootstrap weights

The generalized survey bootstrap can be implemented with a two-step process:

• Step 1: Use make_quad_form_matrix() to represent the variance estimator as a quadratic form’s matrix.

# Load an example dataset of a stratified systematic sample
data('library_stsys_sample', package = 'svrep')

# Represent the SD2 successive-difference estimator as a quadratic form,
# and obtain the matrix of that quadratic form
sd2_quad_form <- make_quad_form_matrix(
  variance_estimator = 'SD2',
  cluster_ids = library_stsys_sample |> select(FSCSKEY),
  strata_ids = library_stsys_sample |> select(SAMPLING_STRATUM),
  strata_pop_sizes = library_stsys_sample |> select(STRATUM_POP_SIZE),
  sort_order = library_stsys_sample |> pull("SAMPLING_SORT_ORDER")
)

• Step 2: Use make_gen_boot_factors() to generate replicate factors based on the target quadratic form. The function argument tau can be used to avoid negative adjustment factors using the previously-described method. The actual value of tau used can be extracted from the function’s output using the attr() function.

rep_adj_factors <- make_gen_boot_factors(
  Sigma = sd2_quad_form,
  num_replicates = 500,
  tau = "auto"
)

tau <- attr(rep_adj_factors, 'tau')
B <- ncol(rep_adj_factors)

Using the adjustment factors thus created, we can create a replicate survey design object by using the function svrepdesign(), with arguments type = "other" and specifying the scale argument to use the factor \(\tau^2/B\).

gen_boot_design <- svrepdesign(
  data = library_stsys_sample |>
    mutate(SAMPLING_WEIGHT = 1/SAMPLING_PROB),
  repweights = rep_adj_factors,
  weights = ~ SAMPLING_WEIGHT,
  combined.weights = FALSE,
  type = "other",
  scale = tau^2/B,
  rscales = rep(1, times = B)
)

This allows us to estimate sampling variances, even for quite complex sampling designs.

gen_boot_design |>
  svymean(x = ~ TOTSTAFF,
          na.rm = TRUE, deff = TRUE)
#>             mean      SE DEff
#> TOTSTAFF 19.7565  4.2024 0.97