Quadratic Forms

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

Suppose we can represent our textbook variance estimator as a quadratic form: v(T̂_y) = y̆Σy̆^T, for some n × n matrix Σ. Our only constraint on Σ 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 y̆ is).

For example, the popular Horvitz-Thompson estimator based on first-order inclusion probabilities π_k and second-order inclusion probabilities π_kl can be represented as a positive semi-definite matrix with entries (1 − π_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)
)
#> 5 x 5 sparse Matrix of class "dsCMatrix"
#>                              
#> [1,]  1.0 -0.5  .    .   -0.5
#> [2,] -0.5  1.0 -0.5  .    .  
#> [3,]  .   -0.5  1.0 -0.5  .  
#> [4,]  .    .   -0.5  1.0 -0.5
#> [5,] -0.5  .    .   -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. This can be done by using the svydesign() function to describe the survey design and then using get_design_quad_form() to obtain the quadratic form for a specified variance estimator.

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

# First, sort the rows in the order used in sampling
library_stsys_sample <- library_stsys_sample |>
  dplyr::arrange(SAMPLING_SORT_ORDER)

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

# Obtain the quadratic form for the target estimator
sd2_quad_form <- get_design_quad_form(
  design = survey_design,
  variance_estimator = "SD2"
)
#> For `variance_estimator='SD2', assumes rows of data are sorted in the same order used in sampling.

class(sd2_quad_form)
#> [1] "dsCMatrix"
#> attr(,"package")
#> [1] "Matrix"
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[['LIBRARIA']] /
                     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
std_error <- sqrt(variance_estimate[1,1])

# Summarize results
sprintf("Estimate: %s", round(point_estimate))
#> [1] "Estimate: 65642"
sprintf("Standard Error: %s", round(std_error))
#> [1] "Standard Error: 13972"

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 a^(b), whose k-th value is denoted a_k^(b). This gives us B replicate estimates of the population total, T̂_y^*(b) = ∑_k ∈ sa_k^(b)y̆_k, for b = 1, …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 a^(b) has expectation 1_n and variance-covariance matrix Σ, then we have the bootstrap expectation E_*(Σ_B) = Σ. Since the bootstrap process takes the sample values y̆ as fixed, the bootstrap expectation of the variance estimator is E_*(y̆^′Σ_By̆) = y̆^′Σy̆. Thus, we can produce a bootstrap variance estimator with the same expectation as the textbook variance estimator simply by randomly generating a^(b) from a distribution with the following two conditions:

Condition 1:   E_*(a) = 1_n

Condition 2:   E_*(a − 1_n)(a − 1_n)^′ = Σ

The simplest, general way to generate the adjustment factors is to simulate from a multivariate normal distribution a ∼ MVN(1_n, Σ), 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_*(y̆^′Σ_By̆) = y̆^′Σy̆.

Adjusting Generalized Survey Bootstrap Replicates to Avoid Negative Weights

Let A = [a⁽¹⁾⋯a^(b)⋯a^(B)] denote the (n × B) matrix of bootstrap adjustment factors. To eliminate negative adjustment factors, Beaumont and Patak (2012) propose forming a rescaled matrix of nonnegative replicate factors 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 τ can be set based on the realized adjustment factor matrix A or by choosing τ prior to generating the adjustment factor matrix A so that τ 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.238

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 generalized 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)
svytotal(x = ~ Bush + Kerry, design = gen_boot_design_ht)

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 250754549

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

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

# Check that the scale factor was calculated correctly
tau <- multistage_boot_design$tau
print(tau)
#> [1] 4.79
B <- ncol(multistage_boot_design$repweights)
print(B)
#> [1] 500

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

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.

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

The actual value of tau used can be extracted from the function’s output using the attr() function.

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

For convenience, the values to use for the scale and rscales arguments of svrepdesign() are included as attributes of the adjustment factors created by make_gen_boot_factors().

# Retrieve value of 'scale'
rep_adj_factors |>
  attr('scale') 
#> [1] 0.041405

# Compare to manually-calculated value
  (tau^2) / B
#> [1] 0.041405

# Retrieve value of 'rscales'
rep_adj_factors |>
  attr('rscales') |> 
  head() # Only show first 5 values
#> [1] 1 1 1 1 1 1

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 τ²/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 = attr(rep_adj_factors, 'scale'),
  rscales = attr(rep_adj_factors, 'rscales')
)

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.756  4.149 0.9455