"Horvitz-Thompson: Same results from conversion vs. creating from scratch"
"Yates-Grundy: Same results from conversion vs. creating from scratch"
"Ultimate Cluster: Same results from conversion vs. creating from scratch"
"Multistage Stratified SRS: Same results from conversion vs. creating from scratch"
"SD1: Same results from conversion vs. creating from scratch"
"SD2: Same results from conversion vs. creating from scratch"
"Two-phase Design: Same results from conversion vs. creating from scratch"
"The sample quadratic form matrix for this design and variance estimator is not positive semidefinite."
"`exact_vcov = TRUE` gives exact variance estimate for totals"
"`exact_vcov = TRUE` throws informative error message if there too few replicates"
"Bootstrap estimate with many replicates close to expected value"
"Returns `tbl_svy` if the input is a `tbl_svy` and 'srvyr' is loaded"

`exact_vcov = TRUE` gives exact variance estimate for totals

Bootstrap estimate with many replicates close to expected value

# Check basic successive-difference quadratic forms ----
successive_diffs <- function(y, type = "SD2") {
result <- sum_of_squares * (n/(2*(n-1)))
sum_of_squares <- sum_of_squares + (tail(y, 1) - head(y, 1))^2
"Helper function for successive difference quadratic forms works correctly", {
object = as.matrix(t(wtd_y) %*% svrep:::make_sd_matrix(n = length(wtd_y), type = "SD2") %*% wtd_y),
expected = successive_diffs(wtd_y, type = "SD2") |> as.matrix()
object = as.matrix(t(wtd_y) %*% svrep:::make_sd_matrix(n = length(wtd_y), type = "SD1") %*% wtd_y),
expected = successive_diffs(wtd_y, type = "SD1") |> as.matrix()
object = as.matrix(t(wtd_y) %*% svrep:::make_sd_matrix(n = length(wtd_y), type = "SD1",
expected = (1-0.9) * successive_diffs(wtd_y, type = "SD1") |> as.matrix()
object = as.matrix(t(wtd_y) %*% svrep:::make_sd_matrix(n = length(wtd_y), type = "SD2",
expected = (1-0.9) * successive_diffs(wtd_y, type = "SD2") |> as.matrix()
object = svrep:::make_sd_matrix(n = 1) |> as.matrix(),
expected = matrix(0, nrow = 1, ncol = 1)
object = as.matrix(c(99, 16) %*% svrep:::make_sd_matrix(n = 2, type = "SD1") %*% c(99,16)),
expected = as.matrix(successive_diffs(c(99,16), type = "SD1"))
object = as.matrix(c(99, 16) %*% svrep:::make_sd_matrix(n = 2, type = "SD2") %*% c(99,16)),
expected = as.matrix(successive_diffs(c(99,16), type = "SD2"))
regexp = "must be an integer greater than"
object = svrep:::make_sd_matrix(n = 2, f = 1.1),
regexp = "must be a single number between"
# Check basic SRSWOR quadratic form ----
"Helper function for SRSWOR quadratic forms works correctly", {
object = as.matrix(t(wtd_y) %*% svrep:::make_srswor_matrix(n = length(wtd_y)) %*% wtd_y),
object = as.matrix(t(wtd_y_matrix) %*% svrep:::make_srswor_matrix(n = length(wtd_y)) %*% wtd_y_matrix),
fpcs = list('sampsize' = as.matrix(rep(n, times = n)),
'popsize' = as.matrix(rep(Inf, times = n))),
object = as.matrix(t(wtd_y) %*% svrep:::make_srswor_matrix(n = length(wtd_y),
expected = (1-0.9) * n *  cov(as.matrix(wtd_y)) |> as.matrix()
object = svrep:::make_srswor_matrix(n = 1) |> as.matrix(),
object = svrep:::make_srswor_matrix(n = 0),
object = svrep:::make_srswor_matrix(n = 2, f = 1.1),
## Load an example dataset that uses unequal probability sampling
## Create matrix to represent the Horvitz-Thompson estimator as a quadratic form
## This is the correct result that we expect the 'svrep' package to match
horvitz_thompson_matrix <- matrix(nrow = n, ncol = n)
horvitz_thompson_matrix[i,j] <- 1 - (pi[i,i] * pi[j,j])/pi[i,j]
## Get the quadratic form using the 'svrep' pacakge
## Compare quadratic form matrix from 'svrep' to the correct result
"Generates correct quadratic form for Horvitz-Thompson", {
expect_equal(object = as.matrix(ht_quad_form), expected = horvitz_thompson_matrix)
## Use the 'survey' package to obtain the Horvitz-Thompson estimate
## Use the 'svrep' package to obtain the Horvitz-Thompson estimate
quad_form_result <- t(y_wtd) %*% ht_quad_form %*% y_wtd
"Horvitz-Thompson estimate matches the 'survey' package", {
expect_equal(object = as.matrix(quad_form_result), expected = svy_pkg_result)
contribution <- -0.5*(1 - (pi_i*pi_j)/pi_ij) * (y_wtd[i] - y_wtd[j])^2
yg_quad_form_matrix <- make_quad_form_matrix(variance_estimator = "Yates-Grundy",
quad_form_result <-  as.matrix(t(y_wtd) %*% yg_quad_form_matrix %*% y_wtd)
"Generates correct quadratic form for Yates-Grundy", {
expect_equal(object = quad_form_result, expected = syg_result)
expect_equal(object = quad_form_result, expected = svy_pkg_result)
v_h_sd1 = (1 - unique(SAMPLING_PROB)) * successive_diffs(y = TOTCIR/SAMPLING_PROB, type = "SD1"),
v_h_sd2 = (1 - unique(SAMPLING_PROB)) * successive_diffs(y = TOTCIR/SAMPLING_PROB, type = "SD2")
v_sd1 = sum(v_h_sd1), v_sd2 = sum(v_h_sd2)
"`make_quad_form_matrix()` works correctly for 'SD1' and 'SD2' with stratified single-stage samples", {
object = as.numeric(t(wtd_y) %*% sd1_quad_form %*% wtd_y),
object = as.numeric(t(wtd_y) %*% sd2_quad_form %*% wtd_y),
##_ Check that works if rows are not in sampling sort order
"`make_quad_form_matrix()` with 'SD1' and 'SD2' works correctly when data frame isn't sorted", {
sorted_result <- t(y_wtd_sorted) %*% sorted_quad_form_matrix %*% y_wtd_sorted
unsorted_result <- t(y_wtd_unsorted) %*% unsorted_quad_form_matrix %*% y_wtd_unsorted
"successive-differences estimate less than ST-SRSWOR approximation, for sys sample from sorted frame", {
"Correct results for ultimate cluster", {
# Correct result for a single-stage stratified sample
object = as.matrix(t(wtd_y_matrix) %*% quad_UC %*% wtd_y_matrix),
) |> svytotal(x = ~ TOTCIR + TOTSTAFF, na.rm = TRUE) |>
# Correct result for a multistage cluster sample
wtd_y_matrix[rowSums(is.na(wtd_y_matrix)) > 0] <- 0
strata_ids = data.frame('PSU_STRATUM' = rep(1, times = nrow(library_multistage_sample)),
'SSU_STRATUM' = rep(1, times = nrow(library_multistage_sample))),
"Fills in implicit cluster or strata IDs", {
strata_pop_sizes = library_stsys_sample |> transmute(N = 50000),
# Check "Stratified Multistage SRS" results ----
"Correct results for stratified multistage SRS", {
variance_estimator = "Stratified Multistage SRS",
# Use function to produce quadratic form matrix
# Determine expected correct quadratic form matrix
expected_quad_form <- quad_form_UC + (0.1 * later_stage_quad_form)
object = as.matrix(t(wtd_y_matrix) %*% qf_matrix %*% wtd_y_matrix),
# Check "Deville-1" and "Deville-2" results ----
"Correct results for Deville-1 and Deville-2", {
diag(exp_dev1_quad_form) <- c_i * (1 - c_i/sum(c_i))
expect_equal(object = as.matrix(dev1_quad_form), expected = exp_dev1_quad_form)
diag(exp_dev2_quad_form) <- c_i * (1 - c_i/sum(c_i))
expect_equal(object = as.matrix(dev2_quad_form), expected = exp_dev2_quad_form)
test_that("Deville-Tille quadratic form gives expected estimate", {
aux_vars <- matrix(rnorm(n = 15), nrow = 5, ncol = 3)
c_values <- ((length(probs))/(length(probs) - ncol(aux_vars))) * (1-probs)
denom_matrix <- lapply(seq_along(probs), function(k) {
c_values[k] * (aux_vars[k,] %*% t(aux_vars[k,])) / (probs[k]^2)
right_vector <- lapply(seq_along(probs), function(l) {
hat_matrix <- inv_denom_matrix %*% right_vector
y_star <- as.vector(aux_vars %*% hat_matrix)
var_est <- sapply(seq_along(probs), function(k) {
## Compare the "by-hand" result to the quadratic form result
Q_var_est <- as.vector(t(y/probs) %*% Q %*% (y/probs))
## If every unit selected with certainty, variance estimate should be 0
object = make_deville_tille_matrix(probs = c(1,1),
expected = matrix(c(0,0,0,0), nrow = 2, ncol = 2)
# Ensure function checks inputs for issues ----
regexp = "`made-up` is not a supported variance estimator"
object = make_quad_form_matrix(variance_estimator = "Yates-Grundy"),
regexp = "must supply a matrix.+joint_probs"
object = make_quad_form_matrix(variance_estimator = "Yates-Grundy",
regexp = "must be a matrix.+no missing values"
object = make_quad_form_matrix(variance_estimator = "SD1"),
regexp = "must supply a matrix or data frame to `cluster_ids`"
object = make_quad_form_matrix(variance_estimator = "SD1",
regexp = "must supply a vector to `sort_order`"
object = make_quad_form_matrix(variance_estimator = "Stratified Multistage SRS",
strata_ids = data.frame(STRATUM = c(1,1))),
regexp = "matrix or data frame to `strata_pop_sizes`"
object = make_quad_form_matrix(variance_estimator = "Deville-1",
regexp = "must supply a matrix or data frame to `probs`"
regexp = "Must set `rows=TRUE` and/or `cols=TRUE`"