Cluster sampling divides the population into \(N\) clusters, out of which a simple random sample of \(n\) clusters is selected. This design is highly practical when a sampling frame of individual elements is unavailable, or when individuals are geographically grouped.
In a one-stage design, every element \(M_i\) within the sampled cluster \(i\) is surveyed.
Total in Cluster \(i\): \(t_i = \sum_{j=1}^{M_i} y_{ij}\)
Ratio Estimator of the Population Mean: The ratio estimator leverages cluster sizes to estimate the overall element mean: \[ \overline{y}_r = \frac{\sum_{i=1}^n t_i}{\sum_{i=1}^n M_i}. \]
Estimated Variance of the Mean: \[ \hat{V}(\overline{y}_r) = \left(1 - \frac{n}{N}\right) \frac{1}{n\overline{M}^2} \frac{\sum_{i=1}^n (t_i - \overline{y}_r M_i)^2}{n-1}. \]
In a two-stage design, we first sample \(n\) clusters, then take a subsample of \(m_i\) elements from the \(M_i\) elements within each selected cluster. Because we don’t observe the true cluster total \(t_i\), we estimate it:
Estimated Total in Cluster \(i\): \(\hat{t}_i = M_i \overline{y}_i\), where \(\overline{y}_i = \frac{1}{m_i} \sum_{j=1}^{m_i} y_{ij}\)
Ratio Estimator of the Population Mean: \[ \overline{y}_r = \frac{\sum_{i=1}^n \hat{t}_i}{\sum_{i=1}^n M_i} \]
Estimated Variance of the Mean: it splits into two recognisable pieces — the variance of a ratio estimate computed over the \(n\) sampled clusters, plus the cluster sampling fraction \(n/N\) times the variance of a stratified sample in which those same clusters are the strata: \[ \hat{V}(\overline{y}_r) = \underbrace{\left(1 - \frac{n}{N}\right) \frac{1}{n\overline{M}^2} \frac{\sum_{i=1}^n (\hat{t}_i - \overline{y}_r M_i)^2}{n-1}}_{\textstyle \hat{V}_{\mathrm{ratio}}} \;+\; \frac{n}{N} \underbrace{\sum_{i=1}^n \left(\frac{M_i}{M}\right)^{2} \left(1 - \frac{m_i}{M_i}\right)\frac{s_i^2}{m_i}}_{\textstyle \hat{V}_{\mathrm{str}}}, \] \[ \qquad M = \sum_{i=1}^n M_i = n\overline{M} . \] Here \(s_i^2\) is the sample variance of the \(m_i\) elements measured in cluster \(i\), and \(M\) is the number of elements the sampled clusters contain. \(\hat{V}_{\mathrm{ratio}}\) is exactly the ratio-estimate variance for \(\sum \hat{t}_i / \sum M_i\), and \(\hat{V}_{\mathrm{str}}\) is exactly the stratified variance \(\sum_h W_h^2(1 - n_h/N_h)s_h^2/n_h\) with the sampled clusters as strata: stratum sizes \(N_h = M_i\), sample sizes \(n_h = m_i\), and weights \(W_h = M_i/M\). Each term \[\hat{v}_i = \left(\frac{M_i}{M}\right)^{2}\left(1 - \frac{m_i}{M_i}\right)\frac{s_i^2}{m_i}\] is one stratum’s contribution; the working tables below carry them in a \(\hat{v}_i\) column whose Sum is \(\hat{V}_{\mathrm{str}}\).
Why is the second piece there at all? If we could see the true cluster totals \(t_i\), the ratio piece alone would be the whole variance. Sub-sampling replaces \(t_i\) by the estimate \(\hat{t}_i\), and the between-cluster sum of squares absorbs only the fraction \(1 - n/N\) of the noise that creates; the stratified piece pays the rest. Three special cases follow straight from the formula:
The estimating functions used throughout this book (including cluster_ratio, used below for one-stage and two-stage cluster samples) 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"))
}
}Consider a population of 187 high school algebra classes in a city. An investigator takes an SRS of 12 of those classes (the clusters). Because this is a one-stage sample, every student in the 12 sampled classes takes the test to evaluate their function knowledge.
We read the raw data and display it using gt’s interactive pagination feature to keep the document clean.
algebra <- read.csv("data/algebra.csv")
algebra |>
gt() |>
tab_header(title = "Raw Data: Algebra Classes") |>
cols_label(Mi = html("M<sub>i</sub>")) |>
opt_interactive(use_pagination = TRUE, page_size_default = 10)Output Note: The table above displays the raw, unaggregated survey data. The interactive pagination allows you to flip through the dataset compactly.
Next, we run our updated function to compute both the cluster-level working table and the overall ratio estimate for the average student score.
alg_res <- cluster_ratio(data = algebra,
cname = "class",
csize = "Mi",
yvar = "score",
N = 187)
# Display the cluster-level working table
alg_res$table| Ratio Estimation: Working Table | ||||||
| \(n = 12\) | ||||||
| \(i\) | \(\bar y_i\) | \(\hat t_i\)1 | \(M_i\) | \(\hat B M_i\)2 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|---|
| 1 | 61.50 | 1,230.00 | 20.00 | 1,251.37 | −21.37 | 456.73 |
| 2 | 64.23 | 1,670.00 | 26.00 | 1,626.78 | 43.22 | 1,867.74 |
| 3 | 58.42 | 1,402.00 | 24.00 | 1,501.65 | −99.65 | 9,929.22 |
| 4 | 58.00 | 1,972.00 | 34.00 | 2,127.33 | −155.33 | 24,127.75 |
| 5 | 58.00 | 1,508.00 | 26.00 | 1,626.78 | −118.78 | 14,109.31 |
| 6 | 64.86 | 1,816.00 | 28.00 | 1,751.92 | 64.08 | 4,106.28 |
| ... | ... | ... | ... | ... | ... | ... |
| Sum | ... | 18,708.00 | 299.00 | 18,708.00 | 0.00 | 194,827.04 |
| 1 \(\hat t_i = \bar y_i \, M_i\). | ||||||
| 2 \(\hat B = \sum_i \hat t_i / \sum_i M_i = 62.57\), \(\hat B M_i = \hat B\, M_i\), \(s_e^2 = 17,711.55\). | ||||||
| \(\bar y_r = \dfrac{\sum_i \hat t_i}{\sum_i M_i} = \dfrac{18,708.00}{299.00} = 62.57\) | ||||||
| \(\mathrm{SE}(\bar y_r) = \dfrac{1}{\bar x}\sqrt{\left(1-\dfrac{n}{N}\right)\dfrac{s_e^2}{n}} = \dfrac{1}{24.92}\sqrt{\left(1-\dfrac{12}{187}\right)\dfrac{17,711.55}{12}} = 1.49\) | ||||||
Output Note: The working table above displays, for each of the 12 randomly selected classes, the class size (\(M_i\)), the total class score (\(t_i\)), the fitted value \(\hat B M_i\), and the residual and squared residual used to estimate \(s_e^2\).
Finally, we extract the ratio estimate from the function output list to assess the overall population mean.
# Display the overall ratio estimate
alg_res$estimate |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 62.57 | 1.49 | 59.29 | 65.85 |
Output Note: Using the ratio estimator, the estimated mean score for a student across all 187 algebra classes is 62.57, with a 95% confidence interval of [59.29, 65.85]. The standard error accounts for the clustering effect, which typically increases variance compared to a pure SRS of students.
In two-stage cluster sampling, we sample clusters, and then subsample elements within them. In this dataset, researcher Arnold investigated coot eggs. He first selected a random sample of clutches (the clusters), and then measured the volume of a random subsample of eggs within each selected clutch.
coots <- read.csv("data/coots.csv")
cootsWe apply the exact same ratio estimator formulation, though mathematically, the cluster total \(t_i\) is now an estimated total \(\hat{t}_i\).
coots_res <- cluster_ratio(data = coots,
cname = "clutch",
csize = "csize",
yvar = "volume")
coots_res$table| Ratio Estimation: Working Table | |||||||
| \(n = 184\) | |||||||
| \(i\) | \(m_i\) | \(\bar y_i\) | \(\hat t_i\)1 | \(M_i\) | \(\hat B M_i\)2 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|---|---|
| 1 | 2.00 | 3.86 | 50.24 | 13.00 | 32.38 | 17.86 | 318.92 |
| 2 | 2.00 | 4.19 | 54.52 | 13.00 | 32.38 | 22.15 | 490.48 |
| 3 | 2.00 | 0.92 | 5.50 | 6.00 | 14.94 | −9.45 | 89.23 |
| 4 | 2.00 | 3.00 | 32.98 | 11.00 | 27.40 | 5.59 | 31.20 |
| 5 | 2.00 | 2.50 | 24.96 | 10.00 | 24.91 | 0.05 | 0.00 |
| 6 | 2.00 | 3.98 | 51.80 | 13.00 | 32.38 | 19.42 | 377.05 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| Sum | 368.00 | ... | 4,375.95 | 1,757.00 | 4,375.95 | 0.00 | 11,439.58 |
| 1 \(\hat t_i = \bar y_i \, M_i\). | |||||||
| 2 \(\hat B = \sum_i \hat t_i / \sum_i M_i = 2.49\), \(\hat B M_i = \hat B\, M_i\), \(s_e^2 = 62.51\). | |||||||
| \(\bar y_r = \dfrac{\sum_i \hat t_i}{\sum_i M_i} = \dfrac{4,375.95}{1,757.00} = 2.49\) | |||||||
| \(\mathrm{SE}(\bar y_r) = \dfrac{1}{\bar x}\sqrt{\dfrac{s_e^2}{n}} = \dfrac{1}{9.55}\sqrt{\dfrac{62.51}{184}} = 0.06\) | |||||||
Output Note: Unlike the previous example, we only measured a subset of eggs per clutch, and the working table now carries an \(m_i\) column recording how many. The function automatically estimates the total volume of each clutch (\(\hat{t}_i\)) by multiplying the true clutch size (\(M_i\)) by the observed subsample mean (\(\overline{y}_i\)); the working table (truncated to the first 6 rows + Sum) shows these estimated totals alongside the ratio fit and residuals. No within-cluster variance column appears here: with \(N\) unknown (N = Inf), the factor \(1/(nN)\) in front of the second-stage term is zero, and \(s_e^2\) already carries all of the sub-sampling noise on its own.
# Display the overall ratio estimate
coots_res$estimate |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 2.49 | 0.06 | 2.37 | 2.61 |
Output Note: The estimated overall average egg volume is 2.491 units. Because the total number of clutches in the population (\(N\)) is not explicitly provided (set to Inf by default), the finite population correction factor is safely omitted from the standard error calculation.
Neither example above actually exercises the within-cluster term — algebra is one-stage, and coots has an unknown \(N\). To see it at work we need a population whose clusters we can count. The agpop.csv file provides one: its 3,059 counties (dropping those coded -99) sit inside exactly \(N = 50\) states, so the states are clusters of known number and known sizes \(M_i\).
We take a two-stage sample in the obvious way: an SRS of \(n = 10\) states, followed by an SRS of \(m_i = 30\%\) of the counties in each selected state.
agpop <- acres_in_thousands(read.csv("data/agpop.csv"))
agpop <- agpop[!is.na(agpop$acres92), ]
## cluster sizes M_i: the number of counties in each state
Mi_all <- table(agpop$state)
agpop$Mi <- as.integer(Mi_all[agpop$state])
set.seed(2026)
states <- sample(names(Mi_all), 10) # stage 1: n = 10 clusters
agclus <- do.call(rbind, lapply(states, function(st) {
rows <- which(agpop$state == st) # stage 2: m_i = 0.3 M_i
agpop[sample(rows, max(1, round(0.30 * length(rows)))), ]
}))
c(clusters = length(states), counties = nrow(agclus))clusters counties
10 184 Everything the estimator needs is a cluster-level summary: the known size \(M_i\) of each sampled state, the number of counties actually measured in it (\(m_i\)), their mean \(\overline{y}_i\) and variance \(s_i^2\), and the estimated state total \(\hat{t}_i = M_i \overline{y}_i\).
n <- length(states) # clusters drawn
N <- 50 # clusters in all
Mi <- tapply(agclus$Mi, agclus$state, function(x) x[1]) # cluster sizes
mi <- tapply(agclus$acres92, agclus$state, length) # counties measured
ybari <- tapply(agclus$acres92, agclus$state, mean) # subsample means
si2 <- tapply(agclus$acres92, agclus$state, var) # within-cluster vars
thati <- Mi * ybari # estimated totals
clus_summary <- data.frame(Mi, mi, ybari, si2, thati)
clus_summary |>
gt(rownames_to_stub = TRUE) |>
cols_label(Mi = md("$M_i$"), mi = md("$m_i$"),
ybari = md("$\\bar y_i$"), si2 = md("$s_i^2$"),
thati = md("$\\hat t_i$")) |>
fmt_auto(data = clus_summary, columns = c("ybari", "si2", "thati"))| \(M_i\) | \(m_i\) | \(\bar y_i\) | \(s_i^2\) | \(\hat t_i\) | |
|---|---|---|---|---|---|
| CA | 58 | 17 | 443.89 | 203,331.78 | 25,745.46 |
| NC | 100 | 30 | 90.67 | 3,996.75 | 9,066.94 |
| NE | 93 | 28 | 489.72 | 100,362.46 | 45,543.98 |
| NJ | 21 | 6 | 40.47 | 2,483.35 | 849.86 |
| NV | 15 | 4 | 213.67 | 41,391.58 | 3,204.99 |
| NY | 59 | 18 | 123.28 | 9,047.88 | 7,273.64 |
| OK | 77 | 23 | 395.78 | 21,374.08 | 30,475.10 |
| PA | 65 | 20 | 104.15 | 5,516.48 | 6,769.90 |
| UT | 29 | 9 | 396.64 | 130,980.12 | 11,502.66 |
| VA | 96 | 29 | 67.00 | 2,568.08 | 6,431.95 |
Output Note: The states drawn range from 15 to 100 counties, of which 4 to 30 were visited. \(s_i^2\) is large throughout: county acreages vary enormously even inside a single state, and that is the variability the second-stage term will account for.
The ratio estimate itself divides the estimated totals by the cluster sizes.
Mbar <- mean(Mi) # average size of the sampled clusters
ybarr <- sum(thati) / sum(Mi) # ybar_r = sum(t_i hat) / sum(M_i)
c(Mbar = Mbar, ybar_r = ybarr) Mbar ybar_r
61.3000 239.5832 Its variance has the two pieces from the formula above. The ratio piece measures how the estimated totals scatter about the fitted line \(\overline{y}_r M_i\):
## ratio piece: the variance of the ratio estimate over the 10 sampled states
ei <- thati - ybarr * Mi # ratio residuals
se2 <- sum(ei^2) / (n - 1) # s_e^2
V_ratio <- (1 - n / N) * se2 / (n * Mbar^2)
c(s_e2 = se2, V_ratio = V_ratio) s_e2 V_ratio
1.650481e+08 3.513821e+03 The stratified piece treats those same 10 states as strata, of sizes \(M_i\) inside the \(M = \sum_i M_i\) counties they contain, and is then scaled by the cluster sampling fraction \(n/N\):
## stratified piece: strata = the sampled states, weights W_i = M_i / M
M <- sum(Mi) # elements in the 10 states
vi <- (Mi / M)^2 * (1 - mi / Mi) * si2 / mi # each stratum's contribution
V_str <- sum(vi) # the stratified variance
c(M = M, n_Mbar = n * Mbar, V_str = V_str, scaled = (n / N) * V_str) M n_Mbar V_str scaled
613.00000 613.00000 180.37073 36.07415 Adding the two gives the standard error, and the interval uses \(n - 1\) degrees of freedom — the clusters, not the counties, are the units that were randomised.
SE <- sqrt(V_ratio + (n / N) * V_str)
mem <- qt(0.975, df = n - 1) * SE
c(Est. = ybarr, S.E. = SE, ci.low = ybarr - mem, ci.upp = ybarr + mem) Est. S.E. ci.low ci.upp
239.58317 59.58099 104.80160 374.36475 Output Note: The estimated mean 1992 acreage per county is 239.6 thousand acres, with a 95% confidence interval of [104.8, 374.4] thousand acres. The interval is wide: only 10 states were visited, so \(n - 1 = 9\) degrees of freedom is all the between-cluster scatter provides.
cluster_ratioThe function performs exactly these steps, provided \(N\) is supplied.
clus_res <- cluster_ratio(data = agclus,
cname = "state",
csize = "Mi",
yvar = "acres92",
N = 50)
clus_res$table| Ratio Estimation: Working Table | |||||||||
| \(n = 10\) | |||||||||
| \(i\) | \(m_i\) | \(\bar y_i\) | \(s_i^2\) | \(\hat t_i\)1 | \(M_i\) | \(\hat B M_i\)2 | \(e_i\) | \(e_i^2\) | \(\hat v_i\)3 |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 17.00 | 443.89 | 203,331.78 | 25,745.46 | 58.00 | 13,895.82 | 11,849.63 | 1.40 × 108 | 75.69 |
| 2 | 30.00 | 90.67 | 3,996.75 | 9,066.94 | 100.00 | 23,958.32 | −14,891.37 | 2.22 × 108 | 2.48 |
| 3 | 28.00 | 489.72 | 100,362.46 | 45,543.98 | 93.00 | 22,281.23 | 23,262.75 | 5.41 × 108 | 57.66 |
| 4 | 6.00 | 40.47 | 2,483.35 | 849.86 | 21.00 | 5,031.25 | −4,181.39 | 1.75 × 107 | 0.35 |
| 5 | 4.00 | 213.67 | 41,391.58 | 3,204.99 | 15.00 | 3,593.75 | −388.75 | 1.51 × 105 | 4.54 |
| 6 | 18.00 | 123.28 | 9,047.88 | 7,273.64 | 59.00 | 14,135.41 | −6,861.77 | 4.71 × 107 | 3.24 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| Sum | 184.00 | ... | ... | 146,864.48 | 613.00 | 146,864.48 | 0.00 | 1.49 × 109 | 180.37 |
| 1 \(\hat t_i = \bar y_i \, M_i\). | |||||||||
| 2 \(\hat B = \sum_i \hat t_i / \sum_i M_i = 239.58\), \(\hat B M_i = \hat B\, M_i\), \(s_e^2 = 1.65 \times 10^{8}\). | |||||||||
| 3 \(\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}}\). | |||||||||
| \(\bar y_r = \dfrac{\sum_i \hat t_i}{\sum_i M_i} = \dfrac{1.47 \times 10^{5}}{613.00} = 239.58\) | |||||||||
| \(\mathrm{SE}(\bar y_r) = \sqrt{\left(1-\dfrac{n}{N}\right)\dfrac{s_e^2}{n\,\bar M^2} + \dfrac{n}{N}\sum_i \hat v_i} = \sqrt{\left(1-\dfrac{10}{50}\right)\dfrac{1.65 \times 10^{8}}{10 \times 61.30^2} + 36.07} = 59.58\) | |||||||||
Output Note: Because this sample is two-stage and \(N\) is finite, the working table carries the three extra columns computed above: \(m_i\), the number of counties measured in each state; \(s_i^2\), the variance among them; and \(\hat v_i = (M_i/M)^2 (1 - m_i/M_i) s_i^2/m_i\), that state’s contribution to the stratified variance. The Sum of the \(\hat v_i\) column is \(\hat{V}_{\mathrm{str}}\) itself, and the source note under the table spells the standard error out as \(\hat{V}_{\mathrm{ratio}} + (n/N)\hat{V}_{\mathrm{str}}\).
clus_res$estimate |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 239.58 | 59.58 | 104.80 | 374.36 |
Output Note: The estimate, standard error, and interval reproduce the direct computation above.
The simulation for this section runs as a standalone app: Two-Stage Cluster Sampling (ratio estimator) with SRS. It opens in a new tab, together with its documentation and the full list of apps.