Unequal Probability Sampling (UPS) improves estimation precision by sampling units with probabilities proportional to their size (PPS). If a cluster or county is larger, it has a higher chance of being included in the sample. This approach prevents extreme variance when cluster sizes vary drastically.
Notation. The population has \(N\) units (clusters, counties, …). Unit \(i\) has
We write \(t_i\) rather than \(y_i\) so that \(\bar{t}_U\), the mean per unit, is kept apart from \(\bar{y}_U\), the mean per unit of size that ratio estimators target.
1. General Formulas for UPS With Replacement (UPSWR)
We make \(n\) independent draws. On each draw, unit \(i\) is selected with probability \(\psi_i\), where \(\sum_{i=1}^N \psi_i = 1\). A unit may be selected more than once, and it then appears in the sample once per draw.
Hansen–Hurwitz estimator of the total. Each draw gives an unbiased estimate \(u_i = t_i/\psi_i\) of \(t\). The HH estimator is their average, and the SE comes from their sample standard deviation \(s_u\): \[ \hat{t}_{\mathrm{HH}} = \frac{1}{n}\sum_{i=1}^n \frac{t_i}{\psi_i}, \qquad \widehat{\mathrm{SE}}(\hat{t}_{\mathrm{HH}}) = \frac{s_u}{\sqrt{n}}, \qquad s_u^2 = \frac{1}{n-1}\sum_{i=1}^n \left(\frac{t_i}{\psi_i} - \hat{t}_{\mathrm{HH}}\right)^2 . \] Dividing by \(N\) gives the HH estimator of the mean per unit, and its SE: \[ \hat{\bar{t}}_{\mathrm{HH}} = \frac{\hat{t}_{\mathrm{HH}}}{N} = \frac{1}{n}\sum_{i=1}^n \frac{t_i}{N\psi_i}, \qquad \widehat{\mathrm{SE}}(\hat{\bar{t}}_{\mathrm{HH}}) = \frac{\widehat{\mathrm{SE}}(\hat{t}_{\mathrm{HH}})}{N}. \]
Hansen–Hurwitz ratio estimator. Dividing the HH estimator of \(t\) by the HH estimator of the total of an auxiliary variable \(x\) gives \[ \hat{\bar{y}}_{\mathrm{HH},r} = \frac{\hat{t}_{\mathrm{HH}}}{\hat{t}_{x,\mathrm{HH}}} = \frac{\sum_{i=1}^n t_i/\psi_i}{\sum_{i=1}^n x_i/\psi_i}, \] which estimates \(t/t_x\). Multiplying by \(\bar{x}_U\) or \(t_x\) gives estimators of \(\bar{t}_U\) or \(t\). Its SE is that of an ordinary ratio estimator, applied to the pairs \((t_i/\psi_i,\; x_i/\psi_i)\). Note that \(\hat{\bar{y}}_{\mathrm{HH},r}\) itself never uses \(t_x\) — only the \(\psi_i\); the total \(t_x\) enters only when converting the ratio to \(\bar{t}_U\) or \(t\).
2. The PPS Special Case
PPS is the general design with the auxiliary variable taken to be the size, and the selection probabilities taken proportional to it: \[ x_i = M_i, \qquad t_x = M, \qquad \bar{x}_U = \bar{M}, \qquad \psi_i = \frac{x_i}{t_x} = \frac{M_i}{M} . \]
The ratio estimator collapses to an average. Every draw now has \(x_i/\psi_i = t_x\), a constant, so the denominator of the ratio estimator is exactly \(n\,t_x = nM\) and \[ \hat{\bar{y}}_{\mathrm{HH},r} = \frac{\sum_{i=1}^n t_i/\psi_i}{\sum_{i=1}^n x_i/\psi_i} = \frac{1}{n}\sum_{i=1}^n \frac{t_i}{x_i} = \frac{1}{n}\sum_{i=1}^n \frac{t_i}{M_i} = \frac{1}{n}\sum_{i=1}^n \bar{y}_i , \] the plain average of the sampled means per element \(\bar{y}_i = t_i/M_i\). Its SE comes from the sample variance \(s_{\bar{y}}^2\) of the \(\bar{y}_i\): \[ \widehat{\mathrm{SE}}(\hat{\bar{y}}_{\mathrm{HH},r}) = \frac{s_{\bar{y}}}{\sqrt{n}} , \qquad s_{\bar{y}}^2 = \frac{1}{n-1}\sum_{i=1}^n \left(\bar{y}_i - \hat{\bar{y}}_{\mathrm{HH},r}\right)^2 . \]
Converting to a total or a per-unit mean uses the general factors \(t_x\) and \(\bar{x}_U\), which here are \(M\) and \(\bar{M}\): \[ \hat{t}_{\mathrm{HH}} = t_x\,\hat{\bar{y}}_{\mathrm{HH},r} = M\,\hat{\bar{y}}_{\mathrm{HH},r}, \qquad \hat{\bar{t}}_{\mathrm{HH}} = \bar{x}_U\,\hat{\bar{y}}_{\mathrm{HH},r} = \bar{M}\,\hat{\bar{y}}_{\mathrm{HH},r} = \frac{M}{N}\,\hat{\bar{y}}_{\mathrm{HH},r} . \]
So under PPS the ratio step adds nothing new: the design already uses \(M_i\) fully, and HH is most precise when \(t_i/\psi_i\) is nearly constant — that is, when \(t_i\) is nearly proportional to \(M_i\).
Comparison with the SRS Ratio Estimator of \(\bar{t}_U\)
Under a simple random sample of \(n\) units drawn without replacement, the usual estimators of \(\bar{t}_U\) are \[ \hat{\bar{t}} = \frac{1}{n}\sum_{i=1}^n t_i, \qquad \hat{\bar{t}}_r = \hat{B}\,\bar{x}_U, \quad \hat{B} = \frac{\sum_{i=1}^n t_i}{\sum_{i=1}^n x_i}, \] with SEs \(\sqrt{(1 - n/N)\,s_t^2/n}\) and \((\bar{x}_U/\bar{x})\sqrt{(1 - n/N)\,s_e^2/n}\), where \(e_i = t_i - \hat{B}x_i\). UPSWR and the SRS ratio estimator both exploit a size variable, but in different ways: UPSWR uses it at the design stage, while the ratio estimator uses it at the estimation stage.
A direct comparison of the SRS ratio estimator with UPSWR, both using the same size variable, is given with the app in Section 6.7.
3. UPS Without Replacement (UPSWO)
When sampling without replacement, the inclusion probability of unit \(i\) across all draws is \(\pi_i\).
The estimating functions used throughout this book (including ratio_est, reg_est, srs_est, upswr_total, upswr_ratio, cluster_upswr_ratio, cluster_upswo_ratio, and the format_est_gt table-formatting helper used below) live in a single shared file, samplingestimate.r, which we source below.
samplingestimate.r)
# =============================================================
# Shared estimation functions for the Elements of Sampling Survey
# book. Source this file from any chapter that needs point/SE/CI
# estimators or their gt-table formatting helpers, e.g.:
#
# source("samplingestimate.r")
#
# =============================================================
suppressPackageStartupMessages({
library(gt)
})
# ---------------------------------------------------------
# Shared gt-table helpers
# ---------------------------------------------------------
## Formats a gt table column as ordinary fixed-point, UNLESS its values are
## too small (< `small` in absolute value) or too large (>= `large`), in
## which case it switches that column to scientific notation (rendered as
## "m x 10^n" via gt's exp_style = "x10n") instead.
fmt_auto <- function (gt_tbl, data, columns, decimals = 2, small = 1e-3, large = 1e5)
{
for (col in columns)
{
x <- data[[col]]
x <- x[is.finite (x) & x != 0]
if (length (x) == 0) next
## judge by the column's TYPICAL magnitude, not by its extremes: the
## Sum row is n times the size of a data row, and one residual that
## lands near zero should not push a whole column into scientific
typical <- stats::median (abs (x))
if (typical < small || typical >= large)
gt_tbl <- fmt_scientific (gt_tbl, columns = tidyselect::all_of (col), decimals = decimals, exp_style = "x10n")
else
gt_tbl <- fmt_number (gt_tbl, columns = tidyselect::all_of (col), decimals = decimals)
}
gt_tbl
}
## Formats a single scalar for inline formula display (fixed-point with
## thousands separators, switching to "m \times 10^{n}" LaTeX scientific
## notation for very small or very large magnitudes, mirroring fmt_auto()'s
## use of exp_style = "x10n") --- used to show worked-example steps such as
## \bar y = sum/n with the actual numbers plugged in.
fmt_num <- function (x, decimals = 2, small = 1e-3, large = 1e5)
{
use_sci <- is.finite (x) & x != 0 & (abs (x) < small | abs (x) >= large)
exps <- ifelse (x == 0, 0, floor (log10 (abs (x))))
mant <- ifelse (x == 0, 0, x / 10^exps)
ifelse (use_sci,
sprintf ("%s \\times 10^{%d}", formatC (mant, format = "f", digits = decimals), exps),
formatC (x, format = "f", digits = decimals, big.mark = ","))
}
## The agriculture files (agpop.csv and the samples drawn from it: agsrs.csv,
## agstrat.csv, ...) record county acreages in acres, which run into the
## millions and make every working table below hard to read. This converts the
## acreage columns to THOUSANDS of acres, so a mean reads as 297.9 rather than
## 297,897, and a variance as 2.0e5 rather than 2.0e11. The file's missing-value
## code, -99, becomes NA first, so it cannot be silently rescaled to -0.099.
## d --- data frame just read from one of those files
## cols --- acreage columns to convert (those absent from d are ignored)
acres_in_thousands <- function (d, cols = c ("acres92", "acres87", "acres82"))
{
cols <- intersect (cols, names (d))
d[cols] <- lapply (d[cols], function (x) { x[x == -99] <- NA; x / 1000 })
d
}
## Truncates a row-level working data.frame whose LAST row is a summary row
## (e.g. "Sum"/"Total"): if there are more than head_n data rows (excluding
## that summary row), keeps only the first head_n data rows, inserts one
## "..." placeholder row, then the summary row --- so long working tables
## stay readable. id_col --- name of the row-label column.
truncate_working <- function (working, id_col, head_n = 6)
{
n <- nrow (working) - 1
if (n <= head_n + 1) return (working)
dots <- working[1, ]
dots[1, ] <- NA
dots[1, id_col] <- "..."
rbind (working[seq_len (head_n), ], dots, working[nrow (working), ])
}
# Helper function to format estimation outputs with appropriate mathematical headers
format_est_gt <- function(est_data, est_type = c("mean", "total", "ratio", "reg_mean", "reg_total", "ratio_mean", "domain_mean",
"hh_total", "hh_mean", "hh_ratio_mean")) {
est_type <- match.arg(est_type)
sym_map <- list(
"mean" = c(est = "$\\bar{y}$", se = "$\\mathrm{SE}(\\bar{y})$"),
"total" = c(est = "$\\hat{t}$", se = "$\\mathrm{SE}(\\hat{t})$"),
"ratio" = c(est = "$\\hat{B}$", se = "$\\mathrm{SE}(\\hat{B})$"),
"reg_mean" = c(est = "$\\bar{y}_{reg}$", se = "$\\mathrm{SE}(\\bar{y}_{reg})$"),
"reg_total" = c(est = "$\\hat{t}_{reg}$", se = "$\\mathrm{SE}(\\hat{t}_{reg})$"),
"ratio_mean" = c(est = "$\\bar{y}_{r}$", se = "$\\mathrm{SE}(\\bar{y}_{r})$"),
"domain_mean" = c(est = "$\\bar{y}_{d}$", se = "$\\mathrm{SE}(\\bar{y}_{d})$"),
## UPSWR (Hansen-Hurwitz): total, mean per unit, and ratio (mean per element)
"hh_total" = c(est = "$\\hat{t}_{\\mathrm{HH}}$",
se = "$\\mathrm{SE}(\\hat{t}_{\\mathrm{HH}})$"),
"hh_mean" = c(est = "$\\hat{\\bar{t}}_{\\mathrm{HH}}$",
se = "$\\mathrm{SE}(\\hat{\\bar{t}}_{\\mathrm{HH}})$"),
"hh_ratio_mean" = c(est = "$\\hat{\\bar{y}}_{\\mathrm{HH},r}$",
se = "$\\mathrm{SE}(\\hat{\\bar{y}}_{\\mathrm{HH},r})$")
)
est_sym <- sym_map[[est_type]]["est"]
se_sym <- sym_map[[est_type]]["se"]
# Convert to data frame. Use unname() to strip hidden lm() coefficient names
if (is.vector(est_data) && !is.list(est_data)) {
df <- as.data.frame(t(unname(est_data)))
} else {
df <- as.data.frame(est_data)
}
# Force standard column names so gt() never gets confused
colnames(df) <- c("Est.", "S.E.", "ci.low", "ci.upp")
# Check if we should use rownames as a stub
use_stub <- !is.null(rownames(df)) && !all(rownames(df) == as.character(1:nrow(df)))
df |>
gt(rownames_to_stub = use_stub) |>
cols_label(
Est. = md(est_sym),
S.E. = md(se_sym),
ci.low = md("95% CI Lower"),
ci.upp = md("95% CI Upper")
) |>
fmt_number(
columns = everything(),
decimals = 2
) |>
tab_options(table.width = pct(70))
}
## Formats an lm object as two gt tables, replacing the plain-text
## summary(lmfit) output:
## $coef --- coefficient table (Estimate, S.E., t value, p value)
## $fit --- one-row table of model fit statistics, including R^2
## returns list($coef, $fit)
format_lm_gt <- function (lmfit)
{
s <- summary (lmfit)
coef_df <- as.data.frame (s$coefficients)
colnames (coef_df) <- c ("Estimate", "S.E.", "t value", "p value")
coef_gt <- coef_df |>
gt (rownames_to_stub = TRUE) |>
fmt_number (columns = c ("Estimate", "S.E.", "t value"), decimals = 4) |>
fmt_number (columns = "p value", decimals = 4) |>
tab_header (title = "Coefficient Estimates") |>
tab_options (table.width = pct(70))
fstat <- s$fstatistic
fit_df <- data.frame (
R.squared = s$r.squared,
Adj.R.squared = s$adj.r.squared,
Sigma = s$sigma,
F.statistic = unname (fstat["value"]),
p.value = unname (pf (fstat["value"], fstat["numdf"], fstat["dendf"], lower.tail = FALSE))
)
fit_gt <- fit_df |>
gt () |>
cols_label (
R.squared = md ("$R^2$"),
Adj.R.squared = md ("Adj. $R^2$"),
Sigma = md ("$\\hat\\sigma$"),
F.statistic = "F-statistic",
p.value = "p-value"
) |>
fmt_number (columns = everything (), decimals = 4) |>
tab_header (title = "Model Fit Summary") |>
tab_options (table.width = pct(70))
list (coef = coef_gt, fit = fit_gt)
}
# ---------------------------------------------------------
# Base Estimators
# ---------------------------------------------------------
## sdata --- a vector of original survey data
## N --- population size
## estimate --- "mean" (default) returns the population MEAN estimate;
## "total" returns the population TOTAL estimate (N * mean,
## with SE/CI scaled accordingly); requires finite N
## show.details --- if TRUE (default), also return a gt table of the
## row-level working values y_i, the fitted value \hat y_i
## (constant, and equal to \bar y under the SRS model),
## e_i = y_i-\hat y_i, and e_i^2, with a Sum row and a
## footnote showing s_y^2; NULL otherwise
## col_labels --- named list (y, yhat, dev, ybar, s) of RAW latex: the
## data, fitted and deviation columns, the symbol for the
## sample mean, and for the sample SD; a partial list fills in
## the remaining defaults, so callers can override just one label
## title --- the working table's title
## returns list ($estimate = c(Est., S.E., ci.low, ci.upp), $table)
srs_est <- function (sdata, N = Inf, estimate = c ("mean", "total"),
show.details = TRUE, col_labels = list (),
title = "SRS Mean Estimation: Working Table")
{
estimate <- match.arg (estimate)
col_labels <- modifyList (list (y = "y_i", yhat = "\\hat y_i", dev = "e_i",
ybar = "\\bar y", s = "s_y"), col_labels)
ybar_lab <- col_labels$ybar
s_lab <- col_labels$s
n <- length (sdata)
ybar <- mean (sdata)
dev <- sdata - ybar
s2y <- sum (dev^2) / (n - 1)
se.ybar <- sqrt((1 - n / N)) * sd (sdata) / sqrt(n)
if (estimate == "total")
{
if (!is.finite (N))
stop ("N must be finite to compute estimate = \"total\"")
est <- N * ybar
se.est <- N * se.ybar
}
else
{
est <- ybar
se.est <- se.ybar
}
mem <- qt (0.975, df = n - 1) * se.est
estimate_vec <- c (Est. = est, S.E. = se.est, ci.low = est - mem, ci.upp = est + mem)
if (!show.details)
return (list (estimate = estimate_vec, table = NULL))
working <- data.frame (
i = c (seq_len (n), "Sum"),
y = c (sdata, sum (sdata)),
yhat = rep (ybar, n + 1),
dev = c (dev, sum (dev)),
dev2 = c (dev^2, sum (dev^2))
)
working <- truncate_working (working, id_col = "i")
table <- gt (working) |>
tab_header (title = title, subtitle = md (sprintf ("$n = %d$", n))) |>
cols_label (
i = md ("$i$"),
y = md (sprintf ("$%s$", col_labels$y)),
yhat = md (sprintf ("$%s$", col_labels$yhat)),
dev = md (sprintf ("$%s$", col_labels$dev)),
dev2 = md (sprintf ("$(%s-%s)^2$", col_labels$y, col_labels$yhat))
) |>
cols_width (i ~ px (40), everything () ~ px (110)) |>
fmt_auto (data = working, columns = c ("y", "yhat", "dev", "dev2"), decimals = 2) |>
sub_missing (missing_text = "...") |>
tab_style (
style = cell_text (weight = "bold"),
locations = cells_body (rows = i == "Sum")
) |>
tab_footnote (
footnote = md (paste0 (
if (col_labels$yhat != ybar_lab)
sprintf ("$%s = %s$ for every unit under the SRS model; ", col_labels$yhat, ybar_lab),
sprintf ("$%s^2 = \\sum_i (%s-%s)^2/(n-1) = %s$.",
s_lab, col_labels$y, col_labels$yhat, fmt_num (s2y)))),
locations = cells_column_labels (columns = dev2)
) |>
tab_source_note (
source_note = md (if (estimate == "total")
sprintf ("$\\hat t = N%s = %s \\times %s = %s$",
ybar_lab, fmt_num (N, 0), fmt_num (ybar), fmt_num (est))
else
sprintf ("$%s = \\dfrac{\\sum_i %s}{n} = \\dfrac{%s}{%d} = %s$",
ybar_lab, col_labels$y, fmt_num (sum (sdata)), n, fmt_num (ybar)))
) |>
tab_source_note (
source_note = md (if (estimate == "total")
sprintf ("$\\mathrm{SE}(\\hat t) = N \\times \\mathrm{SE}(%s) = %s \\times %s = %s$",
ybar_lab, fmt_num (N, 0), fmt_num (se.ybar), fmt_num (se.est))
else if (is.finite (N))
sprintf ("$\\mathrm{SE}(%s) = \\sqrt{1-\\dfrac{n}{N}}\\,\\dfrac{%s}{\\sqrt n} = \\sqrt{1-\\dfrac{%d}{%d}}\\,\\dfrac{%s}{\\sqrt{%d}} = %s$",
ybar_lab, s_lab, n, N, fmt_num (sqrt (s2y)), n, fmt_num (se.ybar))
else
sprintf ("$\\mathrm{SE}(%s) = \\dfrac{%s}{\\sqrt n} = \\dfrac{%s}{\\sqrt{%d}} = %s$",
ybar_lab, s_lab, fmt_num (sqrt (s2y)), n, fmt_num (se.ybar)))
)
list (estimate = estimate_vec, table = table)
}
## ydata --- observations of the variable of interest
## xdata --- observations of the auxiliary variable
## xbarU --- population mean of the auxiliary variable; required when
## estimate is "mean" or "total"; ignored when estimate = "model"
## N --- population size
## estimate --- "mean" (default) returns the ratio estimate of the
## population MEAN of y, Bhat*xbarU; "total" returns the ratio
## estimate of the population TOTAL of y (requires finite N);
## "model" returns the raw ratio coefficient Bhat = ybar/xbar
## itself (with its own SE)
## show.details --- if TRUE (default), also return a gt table of the
## row-level working values y_i, x_i, fitted yhat_i =
## Bhat*x_i, residual e_i, and e_i^2, with a Sum row and a
## footnote showing Bhat and s_e^2; NULL otherwise
## col_labels --- named list (y, x, yhat) of RAW latex for those
## working-table columns, so this function can be reused
## (e.g. by cluster_ratio(), upswr_ratio()) with
## context-appropriate labels
## extra_col --- optional column spec, or list of column specs, adding
## RAW-latex-labelled columns of per-unit `values` right before
## the y column of the working table (e.g. the within-cluster
## means $\bar y_i$ used by cluster_ratio() to build $\hat
## t_i = \bar y_i M_i$); a spec is a list (label, values,
## formula, total): `formula` (optional) is shown as a footnote
## on the y column documenting that relationship, and
## `total = TRUE` puts the column's sum in the Sum row
## tail_col --- optional column spec, or list of specs, in the same format,
## appended AFTER the $e_i^2$ column (e.g. the per-cluster
## second-stage variance contributions $\hat v_i$); each spec's
## `formula` is footnoted on its own column
## var2 --- extra variance component to add to $\hat V(\hat B)\,\bar x^2$,
## i.e. $\hat V(\hat B) = [(1-n/N) s_e^2/n + \mathrm{var2}]/\bar x^2$;
## 0 (default) gives the usual one-stage ratio variance. Used by
## cluster_ratio() to add the within-cluster (second-stage) term
## var2_label --- RAW latex for that extra term, shown symbolically in the
## standard-error source note when var2 > 0
## xbar_label --- RAW latex naming $\bar x$ in the standard-error note
## (cluster_ratio() passes "\\bar M", the mean cluster size)
## returns list ($estimate = c(Est., S.E., ci.low, ci.upp), $table)
ratio_est <- function (ydata, xdata, xbarU = NULL, N = Inf,
estimate = c ("mean", "total", "model"), show.details = TRUE,
col_labels = list (y = "y_i", x = "x_i", yhat = "\\hat y_i"),
extra_col = NULL, B_label = "\\hat B",
tail_col = NULL, var2 = 0, var2_label = NULL,
xbar_label = "\\bar x")
{
estimate <- match.arg (estimate)
n <- length (xdata)
xbar <- mean (xdata)
ybar <- mean (ydata)
B_hat <- ybar / xbar
yhat <- B_hat * xdata
e <- ydata - yhat
var_e <- sum (e^2) / (n - 1)
sd_B_hat <- sqrt ((1 - n/N) * var_e / n + var2) / xbar
if (estimate == "model")
{
est <- B_hat
sd_est <- sd_B_hat
}
else
{
if (is.null (xbarU))
stop ("xbarU (population mean of x) is required when estimate is \"mean\" or \"total\"")
if (estimate == "total")
{
if (!is.finite (N))
stop ("N must be finite to compute estimate = \"total\"")
est <- B_hat * xbarU * N
sd_est <- sd_B_hat * xbarU * N
}
else
{
est <- B_hat * xbarU
sd_est <- sd_B_hat * xbarU
}
}
mem <- qt (0.975, df = n - 1) * sd_est
estimate_vec <- c (Est. = est, S.E. = sd_est, ci.low = est - mem, ci.upp = est + mem)
if (!show.details)
return (list (estimate = estimate_vec, table = NULL))
## a column spec is list (label, values, formula, total); accept either one
## spec or a list of them, so callers can add several columns at once
norm_specs <- function (x)
if (is.null (x)) list () else if (!is.null (x$label)) list (x) else x
extras <- norm_specs (extra_col)
tails <- norm_specs (tail_col)
nm <- function (specs, prefix) if (length (specs)) paste0 (prefix, seq_along (specs)) else character (0)
e_nms <- nm (extras, "extra")
t_nms <- nm (tails, "tail")
add_specs <- function (cols, specs, nms)
{
for (k in seq_along (specs))
cols[[nms[k]]] <- c (specs[[k]]$values,
if (isTRUE (specs[[k]]$total)) sum (specs[[k]]$values) else NA)
cols
}
working_cols <- add_specs (list (i = c (seq_len (n), "Sum")), extras, e_nms)
working_cols$y <- c (ydata, sum (ydata))
working_cols$x <- c (xdata, sum (xdata))
working_cols$yhat <- c (yhat, sum (yhat))
working_cols$e <- c (e, sum (e))
working_cols$e2 <- c (e^2, sum (e^2))
working_cols <- add_specs (working_cols, tails, t_nms)
working <- as.data.frame (working_cols, check.names = FALSE)
working <- truncate_working (working, id_col = "i")
fmt_cols <- c (e_nms, "y", "x", "yhat", "e", "e2", t_nms)
label_args <- list (i = md ("$i$"))
for (k in seq_along (extras))
label_args[[e_nms[k]]] <- md (sprintf ("$%s$", extras[[k]]$label))
label_args$y <- md (sprintf ("$%s$", col_labels$y))
label_args$x <- md (sprintf ("$%s$", col_labels$x))
label_args$yhat <- md (sprintf ("$%s$", col_labels$yhat))
label_args$e <- md ("$e_i$")
label_args$e2 <- md ("$e_i^2$")
for (k in seq_along (tails))
label_args[[t_nms[k]]] <- md (sprintf ("$%s$", tails[[k]]$label))
table <- gt (working) |>
tab_header (title = "Ratio Estimation: Working Table", subtitle = md (sprintf ("$n = %d$", n)))
table <- do.call (cols_label, c (list (table), label_args))
table <- table |>
cols_width (i ~ px (40), everything () ~ px (100)) |>
fmt_auto (data = working, columns = fmt_cols, decimals = 2) |>
sub_missing (missing_text = "...") |>
tab_style (
style = cell_text (weight = "bold"),
locations = cells_body (rows = i == "Sum")
) |>
tab_footnote (
footnote = md (sprintf ("$\\hat B = \\sum_i %s / \\sum_i %s = %s$, $%s = \\hat B\\, %s$, $s_e^2 = %s$.",
col_labels$y, col_labels$x, fmt_num (B_hat), col_labels$yhat, col_labels$x, fmt_num (var_e))),
locations = cells_column_labels (columns = yhat)
)
## an extra column's formula documents how it feeds the y column; a tail
## column's formula belongs on the tail column itself
for (k in seq_along (extras))
if (!is.null (extras[[k]]$formula))
table <- table |>
tab_footnote (
footnote = md (sprintf ("$%s$.", extras[[k]]$formula)),
locations = cells_column_labels (columns = y)
)
for (k in seq_along (tails))
if (!is.null (tails[[k]]$formula))
table <- table |>
tab_footnote (
footnote = md (sprintf ("$%s$.", tails[[k]]$formula)),
locations = cells_column_labels (columns = tidyselect::all_of (t_nms[k]))
)
table <- table |>
tab_source_note (
source_note = md (if (estimate == "model")
sprintf ("$%s = \\dfrac{\\sum_i %s}{\\sum_i %s} = \\dfrac{%s}{%s} = %s$",
B_label, col_labels$y, col_labels$x, fmt_num (sum (ydata)), fmt_num (sum (xdata)), fmt_num (B_hat))
else if (estimate == "total")
sprintf ("$\\hat t = \\hat B\\,\\bar x_U\\,N = %s \\times %s \\times %s = %s$",
fmt_num (B_hat), fmt_num (xbarU), fmt_num (N, 0), fmt_num (est))
else
sprintf ("$\\bar y_r = \\hat B\\,\\bar x_U = %s \\times %s = %s$",
fmt_num (B_hat), fmt_num (xbarU), fmt_num (est)))
) |>
tab_source_note (
source_note = md (if (estimate == "model")
(if (var2 > 0)
## two parts: the ratio variance over the sampled units, plus
## whatever second-stage variance the caller supplied
sprintf ("$\\mathrm{SE}(%s) = \\sqrt{\\left(1-\\dfrac{n}{N}\\right)\\dfrac{s_e^2}{n\\,%s^2} + %s} = \\sqrt{\\left(1-\\dfrac{%d}{%d}\\right)\\dfrac{%s}{%d \\times %s^2} + %s} = %s$",
B_label, xbar_label,
if (is.null (var2_label)) "\\mathrm{var}_2" else var2_label,
n, N, fmt_num (var_e), n, fmt_num (xbar),
fmt_num (var2 / xbar^2), fmt_num (sd_B_hat))
else if (is.finite (N))
sprintf ("$\\mathrm{SE}(%s) = \\dfrac{1}{\\bar x}\\sqrt{\\left(1-\\dfrac{n}{N}\\right)\\dfrac{s_e^2}{n}} = \\dfrac{1}{%s}\\sqrt{\\left(1-\\dfrac{%d}{%d}\\right)\\dfrac{%s}{%d}} = %s$",
B_label, fmt_num (xbar), n, N, fmt_num (var_e), n, fmt_num (sd_B_hat))
else
sprintf ("$\\mathrm{SE}(%s) = \\dfrac{1}{\\bar x}\\sqrt{\\dfrac{s_e^2}{n}} = \\dfrac{1}{%s}\\sqrt{\\dfrac{%s}{%d}} = %s$",
B_label, fmt_num (xbar), fmt_num (var_e), n, fmt_num (sd_B_hat)))
else if (estimate == "total")
sprintf ("$\\mathrm{SE}(\\hat t) = \\mathrm{SE}(\\hat B)\\,\\bar x_U\\,N = %s \\times %s \\times %s = %s$",
fmt_num (sd_B_hat), fmt_num (xbarU), fmt_num (N, 0), fmt_num (sd_est))
else
sprintf ("$\\mathrm{SE}(\\bar y_r) = \\mathrm{SE}(\\hat B)\\,\\bar x_U = %s \\times %s = %s$",
fmt_num (sd_B_hat), fmt_num (xbarU), fmt_num (sd_est)))
)
list (estimate = estimate_vec, table = table)
}
## ydata --- observations of the variable of interest
## xdata --- observations of the auxiliary variable
## xbarU --- population mean of the auxiliary variable; required when
## estimate is "mean" or "total"; ignored when estimate = "model"
## N --- population size
## estimate --- "mean" (default) returns the regression estimate of the
## population MEAN of y, Bhat0 + Bhat1*xbarU; "total" returns
## the regression estimate of the population TOTAL of y
## (requires finite N); "model" returns the raw fitted slope
## Bhat1 itself (with its own SE from the lm fit)
## show.details --- if TRUE (default), also return a gt table of the
## row-level working values y_i, x_i, fitted yhat_i, e_i,
## e_i^2, with a Sum row and a footnote showing Bhat0,
## Bhat1, and s_e^2; NULL otherwise
## returns list ($estimate = c(Est., S.E., ci.low, ci.upp), $table)
reg_est <- function (ydata, xdata, xbarU = NULL, N = Inf,
estimate = c ("mean", "total", "model"), show.details = TRUE)
{
estimate <- match.arg (estimate)
n <- length (ydata)
lmfit <- lm (ydata ~ xdata)
Bhat <- lmfit$coefficients
yhat <- lmfit$fitted.values
e <- lmfit$residuals
SSe <- sum (e^2) / (n - 2)
if (estimate == "model")
{
est <- unname (Bhat[2])
sd_est <- summary (lmfit)$coefficients[2, "Std. Error"]
}
else
{
if (is.null (xbarU))
stop ("xbarU (population mean of x) is required when estimate is \"mean\" or \"total\"")
yhat_reg <- unname (Bhat[1] + Bhat[2] * xbarU)
se_yhat_reg <- sqrt ((1-n/N) * SSe / n)
if (estimate == "total")
{
if (!is.finite (N))
stop ("N must be finite to compute estimate = \"total\"")
est <- yhat_reg * N
sd_est <- se_yhat_reg * N
}
else
{
est <- yhat_reg
sd_est <- se_yhat_reg
}
}
mem <- qt (0.975, df = n - 2) * sd_est
estimate_vec <- c (Est. = est, S.E. = sd_est, ci.low = est - mem, ci.upp = est + mem)
if (!show.details)
return (list (estimate = estimate_vec, table = NULL))
Sxx <- sum ((xdata - mean (xdata))^2)
Sxy <- sum ((xdata - mean (xdata)) * (ydata - mean (ydata)))
working <- data.frame (
i = c (seq_len (n), "Sum"),
y = c (ydata, sum (ydata)),
x = c (xdata, sum (xdata)),
yhat = c (yhat, sum (yhat)),
e = c (e, sum (e)),
e2 = c (e^2, sum (e^2))
)
working <- truncate_working (working, id_col = "i")
table <- gt (working) |>
tab_header (title = "Regression Estimation: Working Table", subtitle = md (sprintf ("$n = %d$", n))) |>
cols_label (
i = md ("$i$"),
y = md ("$y_i$"),
x = md ("$x_i$"),
yhat = md ("$\\hat y_i$"),
e = md ("$e_i$"),
e2 = md ("$e_i^2$")
) |>
cols_width (i ~ px (40), everything () ~ px (100)) |>
fmt_auto (data = working, columns = c ("y", "x", "yhat", "e", "e2"), decimals = 2) |>
sub_missing (missing_text = "...") |>
tab_style (
style = cell_text (weight = "bold"),
locations = cells_body (rows = i == "Sum")
) |>
tab_footnote (
footnote = md (sprintf ("$\\hat B_0 = %s$, $\\hat B_1 = %s$, $\\hat y_i = \\hat B_0 + \\hat B_1 x_i$, $s_e^2 = %s$.",
fmt_num (Bhat[1]), fmt_num (Bhat[2]), fmt_num (SSe))),
locations = cells_column_labels (columns = yhat)
) |>
tab_source_note (
source_note = md (if (estimate == "model")
sprintf ("$\\hat B_1 = \\dfrac{\\sum_i(x_i-\\bar x)(y_i-\\bar y)}{\\sum_i(x_i-\\bar x)^2} = \\dfrac{%s}{%s} = %s$",
fmt_num (Sxy), fmt_num (Sxx), fmt_num (unname (Bhat[2])))
else if (estimate == "total")
sprintf ("$\\hat t_{reg} = N\\bar y_{reg} = %s \\times %s = %s$",
fmt_num (N, 0), fmt_num (yhat_reg), fmt_num (est))
else
sprintf ("$\\bar y_{reg} = \\hat B_0+\\hat B_1\\bar x_U = %s + %s \\times %s = %s$",
fmt_num (Bhat[1]), fmt_num (Bhat[2]), fmt_num (xbarU), fmt_num (yhat_reg)))
) |>
tab_source_note (
source_note = md (if (estimate == "model")
sprintf ("$\\mathrm{SE}(\\hat B_1) = \\sqrt{\\dfrac{s_e^2}{\\sum_i(x_i-\\bar x)^2}} = \\sqrt{\\dfrac{%s}{%s}} = %s$",
fmt_num (SSe), fmt_num (Sxx), fmt_num (sd_est))
else if (estimate == "total")
sprintf ("$\\mathrm{SE}(\\hat t_{reg}) = N\\times\\mathrm{SE}(\\bar y_{reg}) = %s \\times %s = %s$",
fmt_num (N, 0), fmt_num (se_yhat_reg), fmt_num (sd_est))
else if (is.finite (N))
sprintf ("$\\mathrm{SE}(\\bar y_{reg}) = \\sqrt{\\left(1-\\dfrac{n}{N}\\right)\\dfrac{s_e^2}{n}} = \\sqrt{\\left(1-\\dfrac{%d}{%d}\\right)\\dfrac{%s}{%d}} = %s$",
n, N, fmt_num (SSe), n, fmt_num (se_yhat_reg))
else
sprintf ("$\\mathrm{SE}(\\bar y_{reg}) = \\sqrt{\\dfrac{s_e^2}{n}} = \\sqrt{\\dfrac{%s}{%d}} = %s$",
fmt_num (SSe), n, fmt_num (se_yhat_reg)))
)
list (estimate = estimate_vec, table = table)
}
## sdata --- a vector of original survey data in a domain
## N --- population size
## n --- total sample size (not the sample size in the domain)
## to find total, multiply domain size N_d to the estimate returned by this function
domain_est <- function (sdata, n, N = Inf)
{
n_d <- length (sdata)
ybar <- mean (sdata)
se.ybar <- sqrt((1 - n / N)) * sd (sdata) / sqrt(n_d)
mem <- qt (0.975, df = n_d - 1) * se.ybar
c (Est. = ybar, S.E. = se.ybar, ci.low = ybar - mem, ci.upp = ybar + mem)
}
# ---------------------------------------------------------
# Stratified Estimator
# ---------------------------------------------------------
## Given per-stratum summaries, compute the stratified (or post-stratified)
## mean or total.
## estimate --- "mean" (default) returns the stratified MEAN estimate;
## "total" returns the stratified TOTAL estimate
## show.details --- if TRUE (default), also return a gt table showing the
## stratum-by-stratum working, with a Total row
## for poststratification, use nh = n * Nh/N
## returns list ($estimate = c(Est., S.E., ci.low, ci.upp), $table)
str_est <- function (ybarh, sh, nh, Nh, estimate = c ("mean", "total"), show.details = TRUE)
{
estimate <- match.arg (estimate)
N <- sum (Nh)
Pi_h <- Nh / N
v_h <- (1 - nh / Nh) * Pi_h^2 * sh^2 / nh
ybar <- sum (ybarh * Pi_h)
seybar <- sqrt (sum (v_h))
if (estimate == "total")
{
est <- N * ybar
se.est <- N * seybar
}
else
{
est <- ybar
se.est <- seybar
}
mem <- 1.96 * se.est
estimate_vec <- c (Est. = est, S.E. = se.est, ci.low = est - mem, ci.upp = est + mem)
if (!show.details)
return (list (estimate = estimate_vec, table = NULL))
stratum <- names (Nh)
if (is.null (stratum)) stratum <- seq_along (Nh)
working <- data.frame (
Stratum = c (stratum, "Total"),
Nh = c (Nh, N),
nh = c (nh, sum (nh)),
Pi_h = c (Pi_h, sum (Pi_h)),
ybarh = c (ybarh, NA),
sh2 = c (sh^2, NA),
weighted_mean = c (Pi_h * ybarh, ybar),
v_h = c (v_h, seybar^2)
)
table <- gt (working) |>
tab_header (title = "Stratified Mean Estimation: Working Table") |>
cols_label (
Stratum = md ("Stratum"),
Nh = md ("$N_h$"),
nh = md ("$n_h$"),
Pi_h = md ("$\\pi_h$"),
ybarh = md ("$\\bar y_h$"),
sh2 = md ("$s_h^2$"),
weighted_mean = md ("$\\pi_h\\bar y_h$"),
v_h = md ("$v_h$")
) |>
cols_width (c (Nh, nh) ~ px (60), everything () ~ px (110)) |>
fmt_number (columns = c (Nh, nh), decimals = 0) |>
fmt_number (columns = Pi_h, decimals = 2) |>
fmt_auto (data = working, columns = c ("ybarh", "sh2", "weighted_mean", "v_h"), decimals = 2) |>
sub_missing (missing_text = "") |>
tab_style (
style = cell_text (weight = "bold"),
locations = cells_body (rows = Stratum == "Total")
) |>
tab_footnote (
footnote = md ("$\\pi_h = N_h/N$ is the stratum weight."),
locations = cells_column_labels (columns = Pi_h)
) |>
tab_footnote (
footnote = md ("$v_h = (1-n_h/N_h)\\,\\pi_h^2\\, s_h^2/n_h$ is stratum $h$'s contribution to $V(\\bar y_{str}) = \\sum_h v_h$."),
locations = cells_column_labels (columns = v_h)
) |>
tab_source_note (
source_note = md (if (estimate == "total")
sprintf ("$\\hat t_{str} = N\\bar y_{str} = %s \\times %s = %s$",
fmt_num (N, 0), fmt_num (ybar), fmt_num (est))
else
sprintf ("$\\bar y_{str} = \\sum_h \\pi_h\\bar y_h = %s$", fmt_num (ybar)))
) |>
tab_source_note (
source_note = md (if (estimate == "total")
sprintf ("$\\mathrm{SE}(\\hat t_{str}) = N \\times \\mathrm{SE}(\\bar y_{str}) = %s \\times %s = %s$",
fmt_num (N, 0), fmt_num (seybar), fmt_num (se.est))
else
sprintf ("$\\mathrm{SE}(\\bar y_{str}) = \\sqrt{\\sum_h v_h} = \\sqrt{%s} = %s$",
fmt_num (sum (v_h)), fmt_num (seybar)))
)
list (estimate = estimate_vec, table = table)
}
## this function finds statistical estimates given a dataset with sampling weight
# stratdata --- data.frame containing stratified sample
# y --- name of variable for which we want to estiamte population mean
# stratum --- name of variable that will be used as stratum variable
# weight --- name of variable indicating sampling weight
# note: from weights we can find Nh (see the code for formula)
# returns the same list ($estimate, $table) as str_est(), which
# this function calls after deriving ybarh/sh/nh/Nh from the raw data
str_est_data <- function (stratdata, y, stratum, weight, estimate = c ("mean", "total"), show.details = TRUE)
{
## compute stratum-wise data
sh <- tapply (stratdata[, y], stratdata[,stratum], sd)
ybarh <- tapply (stratdata[, y], stratdata[,stratum], mean)
## find population stratum size using sampling weight included in the data set
Nh <- tapply (stratdata[, weight], stratdata[,stratum], sum)
nh <- tapply (1:nrow(stratdata), stratdata[,stratum], length)
str_est (ybarh, sh, nh, Nh, estimate = estimate, show.details = show.details)
}
# ---------------------------------------------------------
# Cluster (ratio-to-size) Estimator
# ---------------------------------------------------------
## Cluster ratio estimation (one-stage if Mi = mi, i.e. every element in the
## sampled cluster is measured; two-stage/ratio-to-size if only mi < Mi
## elements are subsampled, so t_hat_i = Mi * ybari is an ESTIMATED total).
## data --- original data frame for holding data
## cname --- variable recording cluster (psu) identity
## csize --- variable recording cluster (psu) population size (not sample size)
## yvar --- variable of interest
## N --- total number of clusters (psus) in the population
## estimate --- "mean" (default) or "model" both return the raw ratio
## Bhat = sum(t_hat_i)/sum(Mi) --- for this ratio-to-size
## cluster estimator, that raw B IS the population MEAN per
## element, so the two coincide; "total" returns the
## population TOTAL of y, Bhat*Mtotal_U (requires Mtotal_U,
## the population total of the cluster-size variable)
## Mtotal_U --- population total of the cluster-size variable; required
## when estimate = "total"
## show.details --- if TRUE, also return a gt table of the cluster-level
## working values (m_i, ybar_i, s_i^2, t_hat_i, M_i,
## fitted, e_i, e_i^2, v_hat_i), via ratio_est()'s own
## working table
##
## Variance, in two readable parts (Lohr, Sampling: Design and Analysis):
## $\hat V(\bar y_r) = \hat V_{\mathrm{ratio}} + (n/N)\, \hat V_{\mathrm{str}}$.
## The first part is the variance of the RATIO estimate $\sum \hat t_i /
## \sum M_i$ formed from the n sampled clusters,
## $\hat V_{\mathrm{ratio}} = (1 - n/N)\, s_e^2 / (n \bar M^2)$, with
## $s_e^2 = \sum_i (\hat t_i - \bar y_r M_i)^2/(n-1)$,
## and is all a one-stage sample needs. The second part is the variance of a
## STRATIFIED sample in which the n sampled clusters are the strata, of sizes
## $M_i$ inside $M = \sum_{i \in S} M_i = n \bar M$ elements:
## $\hat V_{\mathrm{str}} = \sum_{i \in S} \hat v_i$,
## $\hat v_i = (M_i/M)^2 (1 - m_i/M_i)\, s_i^2 / m_i$,
## where $m_i$ is the number of elements measured in cluster $i$ and $s_i^2$
## their sample variance. Sub-sampling makes $\hat t_i = M_i \bar y_i$ an
## estimate rather than a total, and $s_e^2$ recovers only the fraction
## $1 - n/N$ of that noise; the $(n/N)\hat V_{\mathrm{str}}$ term adds back
## the rest. It vanishes for a one-stage sample ($m_i = M_i$) and for an
## unknown or infinite $N$, and at $n = N$ it is the whole variance --- a
## census of the clusters IS a stratified sample, and the formula says so.
## returns list ($estimate = c(Est., S.E., ci.low, ci.upp), $table)
cluster_ratio <- function (data, cname, csize, yvar, N = Inf,
estimate = c ("mean", "total", "model"),
show.details = TRUE, Mtotal_U = NULL)
{
estimate <- match.arg (estimate)
clust <- data[, cname]
ydata <- data[, yvar]
ybari <- tapply (ydata, clust, mean)
mi <- tapply (ydata, clust, length) # elements measured in cluster i
si2 <- tapply (ydata, clust, function (v) if (length (v) > 1) var (v) else 0)
Mi <- tapply (data[, csize], clust, function (x) x[1])
## same as the cluster total if all Mi elements were measured (mi = Mi)
t_hat_cls <- ybari * Mi
n <- length (Mi)
## Second-stage contributions: treat the n sampled clusters as strata of
## sizes M_i within M = sum_{i in S} M_i = n Mbar elements, so sum(vi) is
## exactly the stratified variance V_str and the term to add is (n/N) V_str.
## ratio_est() divides by xbar^2 = Mbar^2, so var2 carries that factor back.
M <- sum (Mi)
vi <- (Mi / M)^2 * (1 - mi / Mi) * si2 / mi
subsampled <- any (mi < Mi)
two_stage <- subsampled && is.finite (N)
var2 <- if (two_stage) (n / N) * sum (vi) * mean (Mi)^2 else 0
ybar_col <- list (label = "\\bar y_i", values = ybari, formula = "\\hat t_i = \\bar y_i \\, M_i")
m_col <- list (label = "m_i", values = mi, total = TRUE)
s2_col <- list (label = "s_i^2", values = si2)
v_col <- list (label = "\\hat v_i", values = vi, total = TRUE,
formula = paste0 ("\\hat v_i = \\left(\\dfrac{M_i}{M}\\right)^2",
"\\left(1-\\dfrac{m_i}{M_i}\\right)\\dfrac{s_i^2}{m_i},",
"\\quad M = \\sum_i M_i = n\\bar M,",
"\\quad \\sum_i \\hat v_i = \\hat V_{\\mathrm{str}}"))
## m_i is shown whenever the clusters were sub-sampled; s_i^2 and v_hat_i
## only when they actually enter the standard error
extras <- if (two_stage) list (m_col, ybar_col, s2_col)
else if (subsampled) list (m_col, ybar_col)
else list (ybar_col)
tails <- if (two_stage) list (v_col) else NULL
v2_lab <- "\\dfrac{n}{N}\\sum_i \\hat v_i"
if (estimate == "total")
{
if (is.null (Mtotal_U))
stop ("Mtotal_U (population total of the cluster-size variable) is required when estimate = \"total\"")
ratio_est (t_hat_cls, Mi, xbarU = Mtotal_U, N = N, estimate = "mean", show.details = show.details,
col_labels = list (y = "\\hat t_i", x = "M_i", yhat = "\\hat B M_i"),
extra_col = extras, tail_col = tails, var2 = var2, var2_label = v2_lab,
xbar_label = "\\bar M")
}
else
{
ratio_est (t_hat_cls, Mi, N = N, estimate = "model", show.details = show.details,
col_labels = list (y = "\\hat t_i", x = "M_i", yhat = "\\hat B M_i"),
extra_col = extras, tail_col = tails, var2 = var2, var2_label = v2_lab,
xbar_label = "\\bar M", B_label = "\\bar y_r")
}
}
# ---------------------------------------------------------
# UPS Estimators
# ---------------------------------------------------------
## Total estimation for datasets collected with UPS with Replacement
## (psi-estimator / Hansen-Hurwitz).
## estimate --- "total" (default) is this function's whole purpose, computed
## as the SAMPLE MEAN of the expanded pseudo-values t_i/psi_i
## (a with-replacement HT-type property); "mean" instead
## divides that total by N (requires finite N)
upswr_total <- function (total, psi, N = Inf, estimate = c ("total", "mean"), show.details = TRUE)
{
estimate <- match.arg (estimate)
## N = Inf: draws are with replacement, so no finite population correction
res <- srs_est (total/psi, N = Inf, estimate = "mean", show.details = show.details,
col_labels = list (y = "t_i/\\psi_i", yhat = "\\hat t_{\\mathrm{HH}}",
ybar = "\\hat t_{\\mathrm{HH}}", s = "s_u"),
title = "UPSWR Hansen-Hurwitz Estimation: Working Table")
if (estimate == "mean")
{
if (!is.finite (N)) stop ("N must be finite to compute estimate = \"mean\" (mean total per psu)")
res$estimate <- res$estimate / N
}
res
}
## Ratio estimation for datasets collected with UPS with Replacement.
## estimate --- "mean" (default) or "model" both return the raw ratio
## Bhat = sum(t_i/psi_i)/sum(M_i/psi_i) --- that raw B IS the
## population MEAN per element; "total" returns the population
## TOTAL of y, Bhat*Mtotal_U (requires Mtotal_U, the population
## total of the size variable)
upswr_ratio <- function (total, M, psi, N = Inf, estimate = c ("mean", "total", "model"),
show.details = TRUE, Mtotal_U = NULL)
{
estimate <- match.arg (estimate)
if (estimate == "total")
{
if (is.null (Mtotal_U))
stop ("Mtotal_U (population total of the size variable) is required when estimate = \"total\"")
ratio_est (total/psi, M/psi, xbarU = Mtotal_U, N = N, estimate = "mean", show.details = show.details,
col_labels = list (y = "t_i/\\psi_i", x = "M_i/\\psi_i", yhat = "\\hat B\\,M_i/\\psi_i"))
}
else
{
ratio_est (total/psi, M/psi, N = N, estimate = "model", show.details = show.details,
col_labels = list (y = "t_i/\\psi_i", x = "M_i/\\psi_i", yhat = "\\hat B\\,M_i/\\psi_i"))
}
}
## Cluster UPS with Replacement: aggregates to the cluster level (ybari, Mi,
## psi per sampled cluster), then delegates to upswr_ratio(), so the
## working table and estimate/total logic are shared with the base function.
## data --- data frame holding the sample
## sid --- variable recording cluster (psu) identity
## csize --- variable recording cluster (psu) population size
## cpsi --- variable recording the cluster's selection probability psi_i
## yvar --- variable of interest (element level)
## estimate, Mtotal_U --- see upswr_ratio()
cluster_upswr_ratio <- function (data, sid, csize, cpsi, yvar, N = Inf,
estimate = c ("mean", "total", "model"),
show.details = TRUE, Mtotal_U = NULL)
{
estimate <- match.arg (estimate)
clust <- data[, sid]
ydata <- data[, yvar]
ybari <- tapply (ydata, clust, mean)
Mi <- tapply (data[,csize], clust, function (x) x[1])
psi <- tapply (data[,cpsi], clust, function (x) x[1])
t_hat_cls <- ybari * Mi
upswr_ratio (t_hat_cls, Mi, psi, N = N, estimate = estimate,
show.details = show.details, Mtotal_U = Mtotal_U)
}
## Cluster UPS without Replacement: same idea as cluster_upswr_ratio(),
## but dividing by the cluster's inclusion probability pik instead of psi.
## data --- data frame holding the sample
## cname --- variable recording cluster identity
## csize --- variable recording cluster population size
## cpik --- variable recording the cluster's inclusion probability pi_k
## yvar --- variable of interest (element level)
cluster_upswo_ratio <- function (data, cname, csize, cpik, yvar, N = Inf,
estimate = c ("mean", "total", "model"),
show.details = TRUE, Mtotal_U = NULL)
{
estimate <- match.arg (estimate)
clust <- data[, cname]
ydata <- data[, yvar]
ybari <- tapply (ydata, clust, mean)
Mi <- tapply (data[,csize], clust, function (x) x[1])
pik <- tapply (data[,cpik], clust, function (x) x[1])
t_hat_cls <- ybari * Mi
if (estimate == "total")
{
if (is.null (Mtotal_U))
stop ("Mtotal_U (population total of the size variable) is required when estimate = \"total\"")
ratio_est (t_hat_cls/pik, Mi/pik, xbarU = Mtotal_U, N = N, estimate = "mean", show.details = show.details,
col_labels = list (y = "\\hat t_i/\\pi_i", x = "M_i/\\pi_i", yhat = "\\hat B\\,M_i/\\pi_i"))
}
else
{
ratio_est (t_hat_cls/pik, Mi/pik, N = N, estimate = "model", show.details = show.details,
col_labels = list (y = "\\hat t_i/\\pi_i", x = "M_i/\\pi_i", yhat = "\\hat B\\,M_i/\\pi_i"))
}
}We first manually define a small dataset describing hours spent studying, where the probability of sampling a specific student group (\(\psi_i\)) was proportional to the number of students in that group (\(M_i\)).
# 1. Define the raw data
psi <- c(24, 100, 100, 76, 44) / 647
Mi <- c(24, 100, 100, 76, 44)
hrs <- c(75, 203, 203, 191, 168)
# Display the data cleanly
data.frame(Mi, psi, hrs) |>
gt() |>
cols_label(
Mi = md("$M_i$"),
psi = md("$\\psi_i$"),
hrs = md("**Total Hours ($y_i$)**")
) |>
fmt_number(columns = psi, decimals = 3)| \(M_i\) | \(\psi_i\) | Total Hours (\(y_i\)) |
|---|---|---|
| 24 | 0.037 | 75 |
| 100 | 0.155 | 203 |
| 100 | 0.155 | 203 |
| 76 | 0.117 | 191 |
| 44 | 0.068 | 168 |
plot(hrs ~ Mi, main = "Study Hours vs Number of Students",
xlab = "Number of Students (Mi)", ylab = "Study Hours (hrs)")
abline(lm(hrs ~ Mi), col = "blue")
Estimate the total and mean study hours:
# Estimate the total hours spent by all students
res_hrs_total <- upswr_total(total = hrs, psi = psi, show.details = TRUE)
res_hrs_total$table| UPSWR Hansen-Hurwitz Estimation: Working Table | ||||
| \(n = 5\) | ||||
| \(i\) | \(t_i/\psi_i\) | \(\hat t_{\mathrm{HH}}\) | \(e_i\) | \((t_i/\psi_i-\hat t_{\mathrm{HH}})^2\)1 |
|---|---|---|---|---|
| 1 | 2,021.88 | 1,749.01 | 272.86 | 7.45 × 104 |
| 2 | 1,313.41 | 1,749.01 | −435.60 | 1.90 × 105 |
| 3 | 1,313.41 | 1,749.01 | −435.60 | 1.90 × 105 |
| 4 | 1,626.01 | 1,749.01 | −123.00 | 1.51 × 104 |
| 5 | 2,470.36 | 1,749.01 | 721.35 | 5.20 × 105 |
| Sum | 8,745.07 | 1,749.01 | 0.00 | 9.89 × 105 |
| 1 \(s_u^2 = \sum_i (t_i/\psi_i-\hat t_{\mathrm{HH}})^2/(n-1) = 2.47 \times 10^{5}\). | ||||
| \(\hat t_{\mathrm{HH}} = \dfrac{\sum_i t_i/\psi_i}{n} = \dfrac{8,745.07}{5} = 1,749.01\) | ||||
| \(\mathrm{SE}(\hat t_{\mathrm{HH}}) = \dfrac{s_u}{\sqrt n} = \dfrac{497.35}{\sqrt{5}} = 222.42\) | ||||
res_hrs_total$estimate |> format_est_gt("hh_total")| \(\hat{t}_{\mathrm{HH}}\) | \(\mathrm{SE}(\hat{t}_{\mathrm{HH}})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 1,749.01 | 222.42 | 1,131.47 | 2,366.56 |
# Estimate the mean hours spent by each student (SE and CI scale down too);
# under PPS, t_HH / M is exactly the HH ratio estimator
(res_hrs_total$estimate / 647) |> format_est_gt("hh_ratio_mean")| \(\hat{\bar{y}}_{\mathrm{HH},r}\) | \(\mathrm{SE}(\hat{\bar{y}}_{\mathrm{HH},r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 2.70 | 0.34 | 1.75 | 3.66 |
# Estimate the mean hours spent by each student
res_hrs_ratio <- upswr_ratio(total = hrs, M = Mi, psi = psi, estimate = "mean")
res_hrs_ratio$table| Ratio Estimation: Working Table | |||||
| \(n = 5\) | |||||
| \(i\) | \(t_i/\psi_i\) | \(M_i/\psi_i\) | \(\hat B\,M_i/\psi_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 2,021.88 | 647.00 | 1,749.01 | 272.86 | 7.45 × 104 |
| 2 | 1,313.41 | 647.00 | 1,749.01 | −435.60 | 1.90 × 105 |
| 3 | 1,313.41 | 647.00 | 1,749.01 | −435.60 | 1.90 × 105 |
| 4 | 1,626.01 | 647.00 | 1,749.01 | −123.00 | 1.51 × 104 |
| 5 | 2,470.36 | 647.00 | 1,749.01 | 721.35 | 5.20 × 105 |
| Sum | 8,745.07 | 3,235.00 | 8,745.07 | 0.00 | 9.89 × 105 |
| 1 \(\hat B = \sum_i t_i/\psi_i / \sum_i M_i/\psi_i = 2.70\), \(\hat B\,M_i/\psi_i = \hat B\, M_i/\psi_i\), \(s_e^2 = 2.47 \times 10^{5}\). | |||||
| \(\hat B = \dfrac{\sum_i t_i/\psi_i}{\sum_i M_i/\psi_i} = \dfrac{8,745.07}{3,235.00} = 2.70\) | |||||
| \(\mathrm{SE}(\hat B) = \dfrac{1}{\bar x}\sqrt{\dfrac{s_e^2}{n}} = \dfrac{1}{647.00}\sqrt{\dfrac{2.47 \times 10^{5}}{5}} = 0.34\) | |||||
res_hrs_ratio$estimate |> format_est_gt("hh_ratio_mean")| \(\hat{\bar{y}}_{\mathrm{HH},r}\) | \(\mathrm{SE}(\hat{\bar{y}}_{\mathrm{HH},r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 2.70 | 0.34 | 1.75 | 3.66 |
statepop.csv collected with one-stage UPSWRThe file statepop.csv contains data from an unequal-probability sample of 100 U.S. counties. Counties were chosen with probabilities proportional to their populations. The total U.S. population is \(M = 255,077,536\). Because sampling was done with replacement, large counties like Los Angeles occur multiple times in the sample.
statepop <- read.csv("data/statepop.csv")
# Constants
M <- 255077536
popn <- statepop$popn
psi <- popn / M
phys <- statepop$phys
# Paged display of the raw data: printing a data frame under `df-print: paged`
# draws one page at a time instead of a long scrolling box
statepop |>
select(state, county, popn, phys) |>
mutate(psi = popn / M)plot(popn, phys, log = "xy",
main = "Physicians vs Population (Log Scale)",
xlab = "Population", ylab = "Physicians")
Estimating physicians three ways. One sample answers three different questions, using the three estimators of Section 6.1: the population total \(t\), the mean per unit \(\bar{t}_U\) (physicians per county), and the mean per element \(\bar{y}_U\) (physicians per capita).
Total physicians in the US, \(\hat{t}_{\mathrm{HH}}\):
res_phys_total <- upswr_total(phys, psi, show.details = FALSE)
res_phys_total$estimate |> format_est_gt("hh_total")| \(\hat{t}_{\mathrm{HH}}\) | \(\mathrm{SE}(\hat{t}_{\mathrm{HH}})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 570,304.30 | 41,401.23 | 488,155.27 | 652,453.32 |
Output Note: Each draw contributes \(\mathrm{phys}_i/\psi_i\), an unbiased estimate of the national total on its own; their average puts about 570,304 physicians in the US, with a 95% interval of [488,155, 652,453].
Physicians per county, \(\hat{\bar{t}}_{\mathrm{HH}} = \hat{t}_{\mathrm{HH}}/N\):
N <- 3141 # counties in the US
res_phys_county <- upswr_total(phys, psi, N = N, estimate = "mean", show.details = FALSE)
res_phys_county$estimate |> format_est_gt("hh_mean")| \(\hat{\bar{t}}_{\mathrm{HH}}\) | \(\mathrm{SE}(\hat{\bar{t}}_{\mathrm{HH}})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 181.57 | 13.18 | 155.41 | 207.72 |
Output Note: This is the same estimate divided by the known \(N\), so its standard error shrinks by the same factor and the two statements carry identical information. It averages over counties, each weighted equally, so the many small rural counties pull it down to about 182 physicians per county.
Physicians per capita, \(\hat{\bar{y}}_{\mathrm{HH},r}\):
res_phys_ratio <- upswr_ratio(phys, popn, psi, estimate = "mean")
res_phys_ratio$table| Ratio Estimation: Working Table | |||||
| \(n = 100\) | |||||
| \(i\) | \(t_i/\psi_i\) | \(M_i/\psi_i\) | \(\hat B\,M_i/\psi_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 7.46 × 104 | 2.55 × 108 | 5.70 × 105 | −4.96 × 105 | 2.46 × 1011 |
| 2 | 4.99 × 105 | 2.55 × 108 | 5.70 × 105 | −7.16 × 104 | 5.13 × 109 |
| 3 | 4.99 × 105 | 2.55 × 108 | 5.70 × 105 | −7.16 × 104 | 5.13 × 109 |
| 4 | 1.29 × 105 | 2.55 × 108 | 5.70 × 105 | −4.41 × 105 | 1.95 × 1011 |
| 5 | 4.39 × 105 | 2.55 × 108 | 5.70 × 105 | −1.31 × 105 | 1.72 × 1010 |
| 6 | 2.22 × 105 | 2.55 × 108 | 5.70 × 105 | −3.48 × 105 | 1.21 × 1011 |
| ... | ... | ... | ... | ... | ... |
| Sum | 5.70 × 107 | 2.55 × 1010 | 5.70 × 107 | −2.10 × 10−9 | 1.70 × 1013 |
| 1 \(\hat B = \sum_i t_i/\psi_i / \sum_i M_i/\psi_i = 0.00\), \(\hat B\,M_i/\psi_i = \hat B\, M_i/\psi_i\), \(s_e^2 = 1.71 \times 10^{11}\). | |||||
| \(\hat B = \dfrac{\sum_i t_i/\psi_i}{\sum_i M_i/\psi_i} = \dfrac{5.70 \times 10^{7}}{2.55 \times 10^{10}} = 0.00\) | |||||
| \(\mathrm{SE}(\hat B) = \dfrac{1}{\bar x}\sqrt{\dfrac{s_e^2}{n}} = \dfrac{1}{2.55 \times 10^{8}}\sqrt{\dfrac{1.71 \times 10^{11}}{100}} = 1.62 \times 10^{-4}\) | |||||
res_phys_ratio$estimate |> format_est_gt("hh_ratio_mean")| \(\hat{\bar{y}}_{\mathrm{HH},r}\) | \(\mathrm{SE}(\hat{\bar{y}}_{\mathrm{HH},r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 0.00 | 0.00 | 0.00 | 0.00 |
Output Note: Here the size variable is the auxiliary one (\(x_i = \mathrm{popn}_i\)) and the design is PPS on it, so this is exactly the collapse of Section 2: every \(x_i/\psi_i\) equals \(M\), and the ratio estimator reduces to the plain average of the county rates \(\mathrm{phys}_i/\mathrm{popn}_i\). The working table shows those pairs. The estimate is about 0.00224 physicians per person, i.e. roughly 224 per 100,000 people — a population-weighted figure, unlike the per-county average above.
Let’s compare our UPSWR results against a naive Simple Random Sample of 100 counties.
# Read SRS sample, treating -99 as NA
counties <- read.csv("data/counties.csv", na.strings = "-99")
totpop <- 255077536
N <- 3141 # Total number of counties
# 1. Naive SRS estimate for total
res_counties_srs <- srs_est(counties$physician, N = N, estimate = "total", show.details = FALSE)
res_counties_srs$estimate |> format_est_gt("total")| \(\hat{t}\) | \(\mathrm{SE}(\hat{t})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 933,410.97 | 491,982.79 | −42,789.62 | 1,909,611.56 |
# 2. Ratio estimate for total using SRS data
res_counties_ratio <- ratio_est(counties$physician, counties$totpop, xbarU = totpop/N, N = N, estimate = "total")
res_counties_ratio$table| Ratio Estimation: Working Table | |||||
| \(n = 100\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 24.00 | 36,023.00 | 90.31 | −66.31 | 4,397.47 |
| 2 | 44.00 | 73,524.00 | 184.33 | −140.33 | 19,693.17 |
| 3 | 7.00 | 6,408.00 | 16.07 | −9.07 | 82.18 |
| 4 | 11.00 | 19,261.00 | 48.29 | −37.29 | 1,390.49 |
| 5 | 3.00 | 7,649.00 | 19.18 | −16.18 | 261.69 |
| 6 | 327.00 | 188,377.00 | 472.28 | −145.28 | 21,106.52 |
| ... | ... | ... | ... | ... | ... |
| Sum | 29,717.00 | 11,853,116.00 | 29,717.00 | 0.00 | 17,054,519.52 |
| 1 \(\hat B = \sum_i y_i / \sum_i x_i = 0.00\), \(\hat y_i = \hat B\, x_i\), \(s_e^2 = 1.72 \times 10^{5}\). | |||||
| \(\hat t = \hat B\,\bar x_U\,N = 0.00 \times 81,209.02 \times 3,141 = 6.40 \times 10^{5}\) | |||||
| \(\mathrm{SE}(\hat t) = \mathrm{SE}(\hat B)\,\bar x_U\,N = 3.45 \times 10^{-4} \times 81,209.02 \times 3,141 = 87,885.27\) | |||||
res_counties_ratio$estimate |> format_est_gt("total")| \(\hat{t}\) | \(\mathrm{SE}(\hat{t})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 639,506.03 | 87,885.27 | 465,122.60 | 813,889.46 |
# 3. Regression estimate for total using SRS data
res_counties_reg <- reg_est(counties$physician, counties$totpop, xbarU = totpop/N, N = N, estimate = "total")
res_counties_reg$table| Regression Estimation: Working Table | |||||
| \(n = 100\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 24.00 | 36,023.00 | 52.56 | −28.56 | 815.88 |
| 2 | 44.00 | 73,524.00 | 163.74 | −119.74 | 14,337.75 |
| 3 | 7.00 | 6,408.00 | −35.23 | 42.23 | 1,783.70 |
| 4 | 11.00 | 19,261.00 | 2.87 | 8.13 | 66.09 |
| 5 | 3.00 | 7,649.00 | −31.55 | 34.55 | 1,194.03 |
| 6 | 327.00 | 188,377.00 | 504.24 | −177.24 | 31,413.03 |
| ... | ... | ... | ... | ... | ... |
| Sum | 29,717.00 | 11,853,116.00 | 29,717.00 | 0.00 | 11,349,761.43 |
| 1 \(\hat B_0 = -54.23\), \(\hat B_1 = 0.00\), \(\hat y_i = \hat B_0 + \hat B_1 x_i\), \(s_e^2 = 1.16 \times 10^{5}\). | |||||
| \(\hat t_{reg} = N\bar y_{reg} = 3,141 \times 186.52 = 5.86 \times 10^{5}\) | |||||
| \(\mathrm{SE}(\hat t_{reg}) = N\times\mathrm{SE}(\bar y_{reg}) = 3,141 \times 33.49 = 1.05 \times 10^{5}\) | |||||
res_counties_reg$estimate |> format_est_gt("reg_total")| \(\hat{t}_{reg}\) | \(\mathrm{SE}(\hat{t}_{reg})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 585,870.60 | 105,177.42 | 377,149.44 | 794,591.76 |
Output Note: The true value is 532,638. Notice how poorly the naive SRS estimate performs (predicting over 800,000) because it ignores county sizes. Both the SRS Ratio and SRS Regression methods perform much better by incorporating auxiliary population data.
studenthrs.csv collected with two-stage UPSWRHere, we sample classes (Primary Sampling Units) with replacement based on class size, and then subsample students within them to ask about their study hours.
studenthrs <- read.table("data/studentshrs.txt", header = TRUE)
# Calculate sampling probability (cpsi) based on known total students (647)
studenthrs <- studenthrs |>
mutate(cpsi = class_size / 647)
# Run the 2-stage UPSWR function
res_hrs <- cluster_upswr_ratio(
data = studenthrs,
sid = "sid",
csize = "class_size",
cpsi = "cpsi",
yvar = "hrs",
estimate = "mean"
)
# Display the cluster-level working table
res_hrs$table| Ratio Estimation: Working Table | |||||
| \(n = 5\) | |||||
| \(i\) | \(t_i/\psi_i\) | \(M_i/\psi_i\) | \(\hat B\,M_i/\psi_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 2,393.90 | 647.00 | 1,617.50 | 776.40 | 6.03 × 105 |
| 2 | 1,811.60 | 647.00 | 1,617.50 | 194.10 | 3.77 × 104 |
| 3 | 1,552.80 | 647.00 | 1,617.50 | −64.70 | 4.19 × 103 |
| 4 | 1,035.20 | 647.00 | 1,617.50 | −582.30 | 3.39 × 105 |
| 5 | 1,294.00 | 647.00 | 1,617.50 | −323.50 | 1.05 × 105 |
| Sum | 8,087.50 | 3,235.00 | 8,087.50 | 0.00 | 1.09 × 106 |
| 1 \(\hat B = \sum_i t_i/\psi_i / \sum_i M_i/\psi_i = 2.50\), \(\hat B\,M_i/\psi_i = \hat B\, M_i/\psi_i\), \(s_e^2 = 2.72 \times 10^{5}\). | |||||
| \(\hat B = \dfrac{\sum_i t_i/\psi_i}{\sum_i M_i/\psi_i} = \dfrac{8,087.50}{3,235.00} = 2.50\) | |||||
| \(\mathrm{SE}(\hat B) = \dfrac{1}{\bar x}\sqrt{\dfrac{s_e^2}{n}} = \dfrac{1}{647.00}\sqrt{\dfrac{2.72 \times 10^{5}}{5}} = 0.36\) | |||||
Estimate the overall mean study hours:
res_hrs$estimate |> format_est_gt("hh_ratio_mean")| \(\hat{\bar{y}}_{\mathrm{HH},r}\) | \(\mathrm{SE}(\hat{\bar{y}}_{\mathrm{HH},r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 2.50 | 0.36 | 1.50 | 3.50 |
agpop dataLet’s simulate 2000 random samples to see how UPSWR empirically compares to Simple Random Sampling.
# Read population data
agpop <- acres_in_thousands(read.csv("data/agpop.csv", na.strings = "-99"))
agpop <- agpop |> filter(!is.na(acres92), !is.na(acres87))
N <- nrow(agpop)
total_acres92 <- sum(agpop$acres92)
total_acres87 <- sum(agpop$acres87)
acres87 <- agpop$acres87
# Storage vectors for simulation
est_srs300 <- rep(0, 2000)
est_srs300_reg_acres87 <- rep(0, 2000)
est_ups300_acres87 <- rep(0, 2000)
for (i in 1:2000) {
# SRS Simulation
srs300 <- sample(1:N, size = 300)
acres92_srs300 <- agpop$acres92[srs300]
acres87_srs300 <- agpop$acres87[srs300]
est_srs300[i] <- srs_est(acres92_srs300, N = N, estimate = "total", show.details = FALSE)$estimate[1]
est_srs300_reg_acres87[i] <- reg_est(acres92_srs300, acres87_srs300,
xbarU = total_acres87/N, N = N, estimate = "total", show.details = FALSE)$estimate[1]
# UPSWR Simulation (Probability proportional to acres87)
ups300 <- sample(1:N, size = 300, prob = acres87, replace = TRUE)
acres92_ups300 <- agpop$acres92[ups300]
psi <- (acres87 / sum(acres87))[ups300]
est_ups300_acres87[i] <- upswr_total(acres92_ups300, psi, show.details = FALSE)$estimate[1]
}par(mfrow = c(2,1))
# Compare Naive SRS, Regression SRS, and UPSWR
boxplot(data.frame(est_srs300, est_srs300_reg_acres87, est_ups300_acres87),
main = "Comparing SRS, SRS-Regression, and UPSWR",
names = c("SRS Naive", "SRS Regression", "UPSWR"),
col = c("lightblue", "lightgreen", "salmon"))
abline(h = total_acres92, col = "red", lwd = 2)
# Zoom in on the two high-precision estimators
boxplot(data.frame(est_srs300_reg_acres87, est_ups300_acres87),
main = "Zoomed: SRS-Regression vs. UPSWR",
names = c("SRS Regression", "UPSWR"),
col = c("lightgreen", "salmon"))
abline(h = total_acres92, col = "red", lwd = 2)
Output Note: The red line represents the true total acreage. The naive SRS estimator has massive variance. Both SRS Regression and UPSWR drastically reduce variance, centering tightly around the true mean.
The simulation for this section runs as a standalone app: SRS (average and ratio) with UPSWR (Hansen–Hurwitz). It opens in a new tab, together with its documentation and the full list of apps.
Both estimators are \(\bar{x}_U\) times an estimate of the slope \(B\). SRS uses a ratio of sample totals, with every unit equally likely. UPSWR uses a mean of ratios, over a sample tilted towards large \(x_i\). Their errors can be written in terms of the \(e_i\) of the sampled units:
UPSWR. Each draw contributes \(\bar{x}_U r_i\), and its error is exactly a scaled residual: \[ \bar{x}_U r_i - \bar{t}_U = \bar{x}_U\,(r_i - B) = \bar{x}_U\,\frac{t_i - Bx_i}{x_i} = \frac{\bar{x}_U}{x_i}\,e_i . \] A draw lands on unit \(i\) with probability \(\psi_i = x_i/t_x\), so these errors average to \(\sum_i \psi_i\,\bar{x}_U e_i/x_i = \sum_i e_i/N = 0\); this is why HH is unbiased. The \(n\) draws are independent, so the variance is \(1/n\) times the \(\psi\)-weighted mean square of the draw errors: \[ V(\hat{\bar{t}}_{\mathrm{HH}}) = \frac{1}{n}\sum_{i=1}^N \psi_i \left(\frac{\bar{x}_U}{x_i}\,e_i\right)^2 = \frac{1}{n}\sum_{i=1}^N \frac{x_i}{N\bar{x}_U}\cdot\frac{\bar{x}_U^2}{x_i^2}\,e_i^2 = \frac{1}{n}\cdot\frac{1}{N}\sum_{i=1}^N e_i^2\,\frac{\bar{x}_U}{x_i} . \] A unit is drawn in proportion to \(x_i\), but its error is scaled by \(1/x_i\). It squares to \(1/x_i^2\), so each unit’s net contribution is proportional to \(e_i^2/x_i\).
SRS + ratio. Since \(\sum_{i=1}^n t_i - B\sum_{i=1}^n x_i = \sum_{i=1}^n e_i\), \[ \hat{\bar{t}}_{r} - \bar{t}_U = \bar{x}_U\,(\hat{B} - B) = \frac{\bar{x}_U}{\bar{x}}\,\bar{e}, \qquad \bar{e} = \frac{1}{n}\sum_{i=1}^n e_i , \] so the error is (up to \(\bar{x}_U/\bar{x} \approx 1\)) the SRS mean of the residuals. Their population mean is \(0\), so \[ V(\hat{\bar{t}}_{r}) \approx \left(1-\frac{n}{N}\right)\frac{1}{n}\cdot\frac{1}{N-1}\sum_{i=1}^N e_i^2 . \]
Ignoring the finite population correction (and taking \(N - 1 \approx N\)), the two variances differ only by the weight \(\bar{x}_U/x_i\) on each squared residual: \[ V(\hat{\bar{t}}_{r}) \approx \frac{1}{n}\cdot\frac{1}{N}\sum_{i=1}^N e_i^2, \qquad V(\hat{\bar{t}}_{\mathrm{HH}}) = \frac{1}{n}\cdot\frac{1}{N}\sum_{i=1}^N e_i^2\,\frac{\bar{x}_U}{x_i}. \] Under UPSWR, units larger than average count less and smaller ones count more. UPSWR therefore wins when the large residuals belong to the units with \(x_i > \bar{x}_U\), i.e. when the spread of \(t_i\) about the line grows with \(x_i\), as it does for county acreages. When the spread is constant in \(x_i\), the SRS ratio estimator does better. The app’s “Variance reduction vs. SRS average” legend reports both variances, simulated and theoretical, relative to the plain SRS average \(\hat{\bar{t}}\).
When sampling Without Replacement with unequal probabilities, we rely on the inclusion probabilities (\(\pi_k\)) rather than single-draw probabilities (\(\psi_i\)).
classpps <- read.csv("data/classpps.csv")
res_upswo <- cluster_upswo_ratio(
data = classpps,
cname = "class",
csize = "clssize",
cpik = "clssize", # Assuming inclusion prob is proportional to class size
yvar = "hours",
estimate = "mean"
)
# Display the cluster-level working table
res_upswo$table| Ratio Estimation: Working Table | |||||
| \(n = 5\) | |||||
| \(i\) | \(\hat t_i/\pi_i\) | \(M_i/\pi_i\) | \(\hat B\,M_i/\pi_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 3.50 | 1.00 | 3.45 | 0.05 | 0.00 |
| 2 | 5.00 | 1.00 | 3.45 | 1.55 | 2.40 |
| 3 | 3.62 | 1.00 | 3.45 | 0.17 | 0.03 |
| 4 | 3.12 | 1.00 | 3.45 | −0.33 | 0.11 |
| 5 | 2.00 | 1.00 | 3.45 | −1.45 | 2.10 |
| Sum | 17.25 | 5.00 | 17.25 | 0.00 | 4.64 |
| 1 \(\hat B = \sum_i \hat t_i/\pi_i / \sum_i M_i/\pi_i = 3.45\), \(\hat B\,M_i/\pi_i = \hat B\, M_i/\pi_i\), \(s_e^2 = 1.16\). | |||||
| \(\hat B = \dfrac{\sum_i \hat t_i/\pi_i}{\sum_i M_i/\pi_i} = \dfrac{17.25}{5.00} = 3.45\) | |||||
| \(\mathrm{SE}(\hat B) = \dfrac{1}{\bar x}\sqrt{\dfrac{s_e^2}{n}} = \dfrac{1}{1.00}\sqrt{\dfrac{1.16}{5}} = 0.48\) | |||||
Estimate the overall mean study hours (Horvitz-Thompson ratio style):
res_upswo$estimate |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 3.45 | 0.48 | 2.11 | 4.79 |