3  Stratified Random Sampling

3.1 Estimation for Stratified Sampling

Stratified random sampling divides a population of size \(N\) into \(H\) mutually exclusive strata. A simple random sample of size \(n_h\) is then drawn from each stratum of size \(N_h\). This design often yields higher precision than Simple Random Sampling (SRS).

Write \(N_h\), \(n_h\) for the population and sample sizes of stratum \(h\) (\(\sum_h N_h = N\), \(\sum_h n_h = n\)), \(\pi_h = N_h/N\) for the stratum weight, \(\bar{y}_h\) and \(s_h^2\) for the sample mean and variance in stratum \(h\), and \(S_h^2\) for the population variance in stratum \(h\).

1. Stratified Population Mean Estimator

The estimated population mean is a weighted average of the sample means from each stratum: \[ \overline{y}_{str} = \sum_{h=1}^H \pi_h \overline{y}_h, \qquad \pi_h = \frac{N_h}{N}. \] Because the strata are sampled independently, its variance is the weighted sum of the within-stratum SRS variances. Only the within-stratum variability counts; the differences between strata do not: \[ V(\overline{y}_{str}) = \sum_{h=1}^H \pi_h^2 \left(1 - \frac{n_h}{N_h}\right) \frac{S_h^2}{n_h}. \]

2. Estimated Variance, SE, and Confidence Interval

The estimated variance replaces each \(S_h^2\) by \(s_h^2\), keeping the finite population correction in each stratum: \[ \hat{V}(\overline{y}_{str}) = \sum_{h=1}^H v_h, \qquad v_h = \pi_h^2 \left(1 - \frac{n_h}{N_h}\right) \frac{s_h^2}{n_h}, \qquad \mathrm{SE}(\overline{y}_{str}) = \sqrt{\hat{V}(\overline{y}_{str})}, \] where \(v_h\) is stratum \(h\)’s contribution (the \(v_h\) column of the working tables below). An approximate \(95\%\) confidence interval is \(\overline{y}_{str} \pm 1.96\,\mathrm{SE}(\overline{y}_{str})\).

3. Estimating the Population Total

\[ \hat{t}_{str} = N\overline{y}_{str} = \sum_{h=1}^H N_h \overline{y}_h, \qquad \hat{V}(\hat{t}_{str}) = N^2\,\hat{V}(\overline{y}_{str}) = \sum_{h=1}^H N_h^2 \left(1 - \frac{n_h}{N_h}\right) \frac{s_h^2}{n_h}. \]

4. Proportional Allocation

Sample sizes are allocated proportional to the size of the stratum: \[ n_h = n \left( \frac{N_h}{N} \right) = n\pi_h . \] Every stratum then has the same sampling fraction \(n/N\), and the variance simplifies to \[ V_{prop}(\overline{y}_{str}) = \left(1 - \frac{n}{N}\right)\frac{1}{n}\sum_{h=1}^H \pi_h S_h^2 , \] which is never larger than the SRS variance \((1 - n/N)S^2/n\), up to terms of order \(1/N_h\).

5. Neyman (Optimal) Allocation

If surveying costs are equal across strata, Neyman allocation minimizes the variance of the estimator by sampling more heavily from larger and more highly variable strata: \[ n_h = n \frac{N_h S_h}{\sum_{h=1}^H N_h S_h}, \] which gives the smallest possible variance for a total sample size \(n\): \[ V_{opt}(\overline{y}_{str}) = \frac{1}{n}\left(\sum_{h=1}^H \pi_h S_h\right)^2 - \frac{1}{N}\sum_{h=1}^H \pi_h S_h^2 . \] In practice the \(S_h\) are unknown and are replaced by estimates from a previous survey or a pilot sample.

3.2 Functions for Analyzing Data

The estimating functions used throughout this book (including str_est, str_est_data, and srs_est, used below to compute the stratified mean, standard error, and confidence intervals) live in a single shared file, samplingestimate.r, which we source below. str_est_data computes these directly from a dataset, while str_est computes them from pre-calculated summary statistics.

Show the shared estimating functions (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"))
  }
}

3.3 Example: Analysis of “agstrat.csv” data

We will estimate the average 1992 farm acreage (acres92) across US counties, using geographic region (region) as the stratification variable.

3.3.1 Importing agstrat.csv data

agstrat <- acres_in_thousands(read.csv("data/agstrat.csv"))
head(agstrat) |> gt()

Table 3.1: Sample survey data for agstrat.csv, including stratum and weight

county	state	acres92	acres87	acres82	farms92	farms87	farms82	largef92	largef87	largef82	smallf92	smallf87	smallf82	region	rn	weight
PIERCE COUNTY	NE	297.326	332.862	319.619	725	857	865	54	54	42	58	67	48	NC	805	10.23301
JENNINGS COUNTY	IN	124.694	131.481	139.111	658	671	751	14	13	14	42	36	38	NC	241	10.23301
WAYNE COUNTY	OH	246.938	263.457	268.434	1582	1734	1866	20	19	16	175	186	184	NC	913	10.23301
VAN BUREN COUNTY	MI	206.781	190.251	197.055	1164	1278	1464	23	17	9	56	66	55	NC	478	10.23301
OZAUKEE COUNTY	WI	78.772	85.201	89.331	448	483	527	6	5	5	56	49	48	NC	1028	10.23301
CLEARWATER COUNTY	MN	210.897	229.537	213.105	583	699	693	34	32	23	8	19	13	NC	496	10.23301

Output Note: The table above shows the sample survey data, including the stratum (region) and weights for each county.

3.3.2 Spreadsheet calculation

3.3.2.1 Summarizing acre92 in each stratum

We first manually compute the sample size (\(n_h\)), standard deviation (\(s_h\)), and mean (\(\overline{y}_h\)) for each region, matching them with the known total counties in each region (\(N_h\)).

nh <- tapply (agstrat[, "acres92"], agstrat[,"region"], length)
sh <- tapply (agstrat[, "acres92"], agstrat[,"region"], sd)
ybarh <- tapply (agstrat[, "acres92"], agstrat[,"region"], mean)
# create a vector with external information (true population size per stratum)
Nh <- c(NC = 1054, NE = 220, S= 1382, W = 422)

data.frame(Nh, nh, ybarh, sh) |> 
  gt(rownames_to_stub = TRUE) |>
  cols_label(
    Nh = md("$N_h$"),
    nh = md("$n_h$"),
    ybarh = md("$\\bar{y}_h$"),
    sh = md("$s_h$")
  ) |>
  fmt_number(columns = c(ybarh, sh), decimals = 2)

Table 3.2: Sample summary statistics for acres92 by region stratum

	\(N_h\)	\(n_h\)	\(\bar{y}_h\)	\(s_h\)
NC	1054	103	300.50	172.10
NE	220	21	97.63	87.45
S	1382	135	211.32	231.49
W	422	41	662.30	629.43

Output Note: This displays the summary statistics per stratum. Note the significantly high variance (\(s_h\)) in the W region compared to the others.

3.3.2.2 Estimates

Rather than building the “spreadsheet-style” working table by hand, we call str_est(), whose $table shows the same intermediate stratum-level values (\(N_h\), \(n_h\), \(\pi_h\), \(\bar y_h\), \(s_h^2\), \(\pi_h\bar y_h\), \(v_h\)) needed for the final stratified mean and variance, with a Total row and footnotes for \(\pi_h\) and \(v_h\).

res_str_mean <- str_est(ybarh, sh, nh, Nh, estimate = "mean")
res_str_mean$table

Stratified Mean Estimation: Working Table
Stratum	\(N_h\)	\(n_h\)	\(\pi_h\)1	\(\bar y_h\)	\(s_h^2\)	\(\pi_h\bar y_h\)	\(v_h\)2
NC	1,054	103	0.34	300.50	29,618.18	102.90	30.42
NE	220	21	0.07	97.63	7,647.47	6.98	1.68
S	1,382	135	0.45	211.32	53,587.49	94.88	72.20
W	422	41	0.14	662.30	396,185.95	90.80	163.99
Total	3,078	300	1.00			295.56	268.30
1 \(\pi_h = N_h/N\) is the stratum weight.
2 \(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\).
\(\bar y_{str} = \sum_h \pi_h\bar y_h = 295.56\)
\(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{268.30} = 16.38\)

Output Note: The table provides a detailed breakdown of how each individual stratum contributes to the final mean and variance.

Extracting $estimate gives our overall result.

res_str_mean$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
295.56	16.38	263.46	327.67

Output Note: The overall result indicates an estimated mean of 295.6 thousand acres, with a 95% confidence interval of [263.5, 327.7] thousand acres.

To find the population total (total acreage across all US counties), set estimate = "total".

## population total estimate
res_str_total <- str_est(ybarh, sh, nh, Nh, estimate = "total")
res_str_total$estimate |> format_est_gt("total")

\(\hat{t}\)	\(\mathrm{SE}(\hat{t})\)	95% CI Lower	95% CI Upper
909,736.04	50,417.25	810,918.23	1,008,553.84

Output Note: This scales the estimates up to the total area of 909,736 thousand acres.

3.3.3 Using the function “str_est_data”

Instead of pre-calculating the summary stats, this function extracts everything directly from the dataset using the pre-assigned sampling weights.

## In the function "str_est_data", we can find the stratum size with the variable "weights":
tapply(agstrat$weight, agstrat$region, sum)

  NC   NE    S    W 
1054  220 1382  422 

Output Note: This reconstructs the total population sizes (\(N_h\)) per region.

## if the dataset contains a variable "weight"
res_agstrat_acres92 <- str_est_data(agstrat, y="acres92", stratum="region", weight="weight", estimate = "mean")
res_agstrat_acres92$table

Stratified Mean Estimation: Working Table
Stratum	\(N_h\)	\(n_h\)	\(\pi_h\)1	\(\bar y_h\)	\(s_h^2\)	\(\pi_h\bar y_h\)	\(v_h\)2
NC	1,054	103	0.34	300.50	29,618.18	102.90	30.42
NE	220	21	0.07	97.63	7,647.47	6.98	1.68
S	1,382	135	0.45	211.32	53,587.49	94.88	72.20
W	422	41	0.14	662.30	396,185.95	90.80	163.99
Total	3,078	300	1.00			295.56	268.30
1 \(\pi_h = N_h/N\) is the stratum weight.
2 \(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\).
\(\bar y_{str} = \sum_h \pi_h\bar y_h = 295.56\)
\(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{268.30} = 16.38\)

res_agstrat_acres92$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
295.56	16.38	263.46	327.67

Output Note: Once again, this perfectly matches the manual calculation.

## estimating the mean of the number of small farms in 1992
res_agstrat_smallf92 <- str_est_data(agstrat, y="smallf92", stratum="region", weight="weight", estimate = "mean")
res_agstrat_smallf92$table

Stratified Mean Estimation: Working Table
Stratum	\(N_h\)	\(n_h\)	\(\pi_h\)1	\(\bar y_h\)	\(s_h^2\)	\(\pi_h\bar y_h\)	\(v_h\)2
NC	1,054	103	0.34	44.26	1,286.43	15.16	1.32
NE	220	21	0.07	47.24	2,364.79	3.38	0.52
S	1,382	135	0.45	47.39	6,205.45	21.28	8.36
W	422	41	0.14	124.39	100,640.94	17.05	41.66
Total	3,078	300	1.00			56.86	51.86
1 \(\pi_h = N_h/N\) is the stratum weight.
2 \(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\).
\(\bar y_{str} = \sum_h \pi_h\bar y_h = 56.86\)
\(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{51.86} = 7.20\)

res_agstrat_smallf92$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
56.86	7.20	42.75	70.98

Output Note: This demonstrates the function’s flexibility: we can easily estimate a different variable (an average of 56.9 small farms per county).

3.3.4 Comparing with SRS estimate

To see the benefit of stratifying, we compare our stratified estimates against a Simple Random Sample approach on the same population.

agsrs <- acres_in_thousands(read.csv("data/agsrs.csv"))
res_agsrs_srs <- srs_est(agsrs[, "acres92"], N=3078, estimate = "mean", show.details = FALSE)
res_agsrs_srs$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
297.90	18.90	260.71	335.09

Output Note: The Simple Random Sample standard error is noticeably higher at 18.9 thousand acres.

## the ratio of estimated variance (Stratified SE / SRS SE)^2
ratio_var <- (res_agstrat_acres92$estimate["S.E."] / res_agsrs_srs$estimate["S.E."])^2
unname(ratio_var)

[1] 0.7512239

Output Note: The stratified variance is approximately 75.1% of the SRS variance.

## percentage of reduction of variance of str estimates from that of SRS estimate
unname(1 - ratio_var)

[1] 0.2487761

Output Note: Stratification by region reduced the variance by 24.9%. This is a substantial gain in precision for no extra sampling cost.

3.4 Allocation of stratum sample size

Deciding how to distribute the total sample \(n\) among the strata significantly impacts the precision of the estimator.

3.4.1 Analyzing seals.csv collected with stratified sampling

This dataset investigates ringed seal holes in different zones. The researchers used Proportional Allocation, taking 20% of the areas in each zone.

## load data
seals <- read.csv("data/seals.csv")

# survey data summary in each stratum
nh <- as.vector(table(seals[,"zone"]))
sh <- tapply(seals[, "holes"], seals[,"zone"], sd)
ybarh <- tapply(seals[, "holes"], seals[,"zone"], mean)
Nh <- c(68, 84, 48)

res_seals <- str_est(ybarh, sh, nh, Nh, estimate = "mean")
res_seals$table

Stratified Mean Estimation: Working Table
Stratum	\(N_h\)	\(n_h\)	\(\pi_h\)1	\(\bar y_h\)	\(s_h^2\)	\(\pi_h\bar y_h\)	\(v_h\)2
1	68	17	0.34	1.76	3.32	0.60	0.02
2	84	12	0.42	4.42	11.54	1.85	0.15
3	48	11	0.24	10.55	46.07	2.53	0.19
Total	200	40	1.00			4.99	0.35
1 \(\pi_h = N_h/N\) is the stratum weight.
2 \(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\).
\(\bar y_{str} = \sum_h \pi_h\bar y_h = 4.99\)
\(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{0.35} = 0.59\)

Output Note: Zone 3 has the highest standard deviation (\(s_h \approx 6.8\)), indicating a high degree of variability in seal density in that region.

res_seals$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
4.99	0.59	3.83	6.14

Output Note: This provides the estimated mean number of seal holes per square kilometer (5).

3.4.2 Neyman allocation of stratum sample size

Because Zone 3 was highly variable, we could have achieved a lower overall variance if we had sampled it more heavily. Neyman allocation formally calculates this optimal distribution.

## assuming the cost of surveying each zone is the same
Ch <- rep(1, 3)
Sh <- sh ## estimate population Sh with sample sh
NhShDCh <- Nh * Sh / sqrt(Ch)

## the optimal allocation proportion (Lh) for estimating the population total or mean:
Lh <- NhShDCh / sum(NhShDCh)

data.frame(Ch, Nh, Sh, NhShDCh, Lh) |> 
  gt(rownames_to_stub = TRUE) |> 
  cols_label(NhShDCh = md("$N_h S_h / \\sqrt{C_h}$"), Lh = md("$L_h$ (Prop)"), Ch=md("$C_h$"), Sh=md("$S_h$"), Nh=md("$N_h$")) |>
  fmt_number(columns = c(Sh, NhShDCh, Lh), decimals = 3)

Table 3.3: Neyman-optimal allocation proportions for the seals survey zones

	\(C_h\)	\(N_h\)	\(S_h\)	\(N_h S_h / \sqrt{C_h}\)	\(L_h\) (Prop)
1	1	68	1.821	123.831	0.168
2	1	84	3.397	285.327	0.388
3	1	48	6.788	325.809	0.443

Output Note: This displays the optimal allocation proportions (\(L_h\)). Even though Zone 2 is a larger area (\(N_h=84\)), Zone 3’s high variance means we should allocate the largest share (44%) of our sampling effort there.

3.5 Simulation to Study the Efficiency of Stratified Sampling with Different Allocation

We will use the full agpop.csv dataset (acting as our true population) to run 2000 repeated samples. We will compare SRS, Proportional (P2S), and Optimal (Neyman) allocations.

3.5.1 Using “Region” to Stratify

3.5.1.1 Read Population Data

agpop <- acres_in_thousands(read.csv("data/agpop.csv"))
agpop <- agpop[!is.na(agpop$acres92), ] ## remove those counties with na
N <- nrow(agpop)

3.5.1.2 Define Stratum Variable

agpop$stratum <- agpop$region
## reorder agpop for the ease of using strata of "sampling" (very important)
agpop <- agpop[order(agpop$stratum), ]

A stratification variable is most effective when it strongly explains the variance in the target variable. We check this via ANOVA and R-squared.

boxplot(acres92 ~ stratum, data = agpop)

Figure 3.1: 1992 farm acreage by region stratum

lm_region <- lm(agpop$acres92 ~ agpop$stratum)
anova(lm_region)

lm_region_gt <- format_lm_gt(lm_region)
lm_region_gt$fit

Model Fit Summary
\(R^2\)	Adj. \(R^2\)	\(\hat\sigma\)	F-statistic	p-value
0.1821	0.1813	384.8318	226.7299	0.0000

r_squared_region <- summary(lm_region)$r.squared

Output Note: The R-squared is 0.182. The region variable only explains about 18.2% of the variance in acreage, which is helpful but not extremely strong.

3.5.1.3 Stratified Sampling with P2S allocation (Single Run)

# doing one stratified sampling
Nh <- tapply(1:nrow(agpop), agpop$stratum, length) 
nh <- round(Nh/sum(Nh)*300)

## doing stratified sampling
strsample <- strata(agpop, "stratum", size = nh, method = "srswor")

# collecting data on sampled counties
agstrat <- agpop[strsample$ID_unit, ]
agstrat$weight <- 1/strsample$Prob

str_est_data(agstrat, "acres92", "stratum", "weight", estimate = "mean")$estimate |>
  format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
310.51	17.31	276.59	344.43

Output Note: This provides a single point estimate for demonstration.

3.5.1.4 Repeat stratified sampling with P2S allocation 2000 times

Nh <- tapply(1:nrow(agpop), agpop$stratum, length) 
nh <- round(Nh/sum(Nh)*300) ## make sure the order matches Nh

nres <- 2000
str_p2s_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    strsample <- strata(agpop, "stratum", size = nh, method = "srswor")
    agstrat <- agpop[strsample$ID_unit, ]
    agstrat$weight <- 1/strsample$Prob
    str_p2s_simulated[i,] <- str_est_data(agstrat, "acres92", "stratum", "weight", estimate = "mean", show.details = FALSE)$estimate
}

3.5.1.5 Repeat stratified sampling with optimal allocation 2000 times

Nh <- tapply(agpop$acres92, agpop$stratum, length) 
Sh <- tapply(agpop$acres92, agpop$stratum, sd)
nh_opt <- round((Nh*Sh)/sum(Nh*Sh) * 300)

data.frame(Nh, Sh, nh_opt) |> 
  gt(rownames_to_stub = TRUE) |> 
  fmt_number(columns = Sh, decimals = 2)

Table 3.4: Population sizes and Neyman-optimal sample sizes by region stratum

	Nh	Sh	nh_opt
NC	1052	271.19	87
NE	213	78.91	5
S	1376	244.13	102
W	418	836.61	106

Output Note: This table compares the population sizes to the optimal sample sizes. Notice that the West region (high variance) gets a disproportionately larger sample under Neyman allocation.

nres <- 2000
str_neyman_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    strsample <- strata(agpop, "stratum", size = nh_opt, method = "srswor")
    agstrat <- agpop[strsample$ID_unit, ]
    agstrat$weight <- 1/strsample$Prob
    str_neyman_simulated[i,] <- str_est_data(agstrat, "acres92", "stratum", "weight", estimate = "mean", show.details = FALSE)$estimate
}

3.5.1.6 Repeat simple random sampling 2000 times

nres <- 2000
srs_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    srs <- sample(sum(Nh), sum(nh))
    agsrs <- agpop[srs, ]
    srs_simulated[i,] <- srs_est(agsrs[, "acres92"], N = sum(Nh), estimate = "mean", show.details = FALSE)$estimate
}

3.5.1.7 Compare the efficiency of different methods

sim_results_str_region <- data.frame(
  "SRS" = srs_simulated[,1], 
  "Prop2size" = str_p2s_simulated[,1], 
  "Neyman" = str_neyman_simulated[,1]
)
boxplot(sim_results_str_region)
abline(h = mean(agpop$acres92), col = "red")

Figure 3.2: Sampling distribution of the mean under SRS, proportional, and Neyman allocation (region strata, 2000 simulations)

Output Note: The plot visually demonstrates that Neyman allocation has the tightest spread around the true mean (the red line).

sim_means <- sapply(sim_results_str_region, mean)
sim_var <- sapply(sim_results_str_region, var)
sim_var_relative <- sim_var / sim_var[1]
sim_var_reduction <- 1 - sim_var_relative

data.frame(
  "Mean" = sim_means, 
  "Variance" = sim_var, 
  "Relative Variance" = sim_var_relative, 
  "Percentage of Variance Reduction" = sim_var_reduction, 
  check.names = FALSE
) |> 
  gt(rownames_to_stub = TRUE) |> 
  fmt_number(columns = c(Mean, Variance), decimals = 1) |> 
  fmt_percent(columns = c(`Percentage of Variance Reduction`), decimals = 2) |> 
  fmt_number(columns = c(`Relative Variance`), decimals = 3)

Table 3.5: Simulated mean, variance, and variance reduction by allocation method (region strata)

	Mean	Variance	Relative Variance	Percentage of Variance Reduction
SRS	309.0	549.0	1.000	0.00%
Prop2size	308.7	452.6	0.824	17.56%
Neyman	308.8	310.9	0.566	43.37%

Output Note: The table summarizes the simulation. Proportional allocation reduces variance by 17.6%, while Neyman allocation reduces it by 43.4% compared to Simple Random Sampling.

3.5.2 Using “acres82” to Define Strata

What if we stratify by a variable highly correlated with our target, like the acreage from a decade prior (acres82)?

3.5.2.1 Read Population Data

agpop <- acres_in_thousands(read.csv("data/agpop.csv"))
## acres82 is the stratifier here, so drop the counties missing either variable
agpop <- agpop[!is.na(agpop$acres92) & !is.na(agpop$acres82), ] 
N <- nrow(agpop)

3.5.2.2 Define Stratum Variable with Quantiles of “acres82”

We convert the continuous acres82 variable into four categorical strata by cutting it at the 25th, 75th and 95th percentiles. The cut is deliberately uneven: the top 5% of counties, where the acreages run away, get a stratum of their own, which is what makes the stratification so effective below.

agpop$stratum <- cut(agpop$acres82,
                     breaks = quantile(agpop$acres82, probs = c(0,0.25,0.75,0.95,1)),
                     include.lowest = T)
levels(agpop$stratum) <- paste0("acres82", c("Q1", "Q2", "Q3", "Q4"))
agpop <- agpop[order(agpop$stratum), ]

Let’s check the explanatory power of this new stratification variable.

boxplot(acres92 ~ stratum, data = agpop)

Figure 3.3: 1992 farm acreage by acres82 quantile stratum

lm_acres82 <- lm(agpop$acres92 ~ agpop$stratum)
lm_acres82_gt <- format_lm_gt(lm_acres82)
lm_acres82_gt$fit

Model Fit Summary
\(R^2\)	Adj. \(R^2\)	\(\hat\sigma\)	F-statistic	p-value
0.7402	0.7400	216.9580	2,895.9129	0.0000

r_squared_acres82 <- summary(lm_acres82)$r.squared

Output Note: The R-squared is 0.74. The past acreage explains 74% of the variance in 1992 acreage, making this a vastly superior stratification variable.

3.5.2.3 Stratified sampling with P2S allocation

Nh <- tapply(1:nrow(agpop), agpop$stratum, length) 
nh <- round(Nh/sum(Nh)*300)

nres <- 2000
str_p2s_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    strsample <- strata(agpop, "stratum", size = nh, method = "srswor")
    agstrat <- agpop[strsample$ID_unit, ]
    agstrat$weight <- 1/strsample$Prob
    str_p2s_simulated[i,] <- str_est_data(agstrat, "acres92", "stratum", "weight", estimate = "mean", show.details = FALSE)$estimate
}

3.5.2.4 Repeat stratified sampling with optimal allocation 2000 times

Nh <- tapply(1:nrow(agpop), agpop$stratum, length) 
Sh <- tapply(agpop$acres92, agpop$stratum, sd)
nh_opt <- round((Nh*Sh)/sum(Nh*Sh) * 300)

nres <- 2000
str_neyman_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    strsample <- strata(agpop, "stratum", size = nh_opt, method = "srswor")
    agstrat <- agpop[strsample$ID_unit, ]
    agstrat$weight <- 1/strsample$Prob
    str_neyman_simulated[i,] <- str_est_data(agstrat, "acres92", "stratum", "weight", estimate = "mean", show.details = FALSE)$estimate
}

3.5.2.5 Repeat simple random sampling 2000 times

nres <- 2000
srs_simulated <- matrix(0, nres, 4)
for (i in 1:nres) {
    srs <- sample(sum(Nh), sum(nh))
    agsrs <- agpop[srs, ]
    srs_simulated[i,] <- srs_est(agsrs[, "acres92"], N = sum(Nh), estimate = "mean", show.details = FALSE)$estimate
}

3.5.2.6 Compare the efficiency of different methods

sim_results_str_acres82 <- data.frame(
  "SRS" = srs_simulated[,1], 
  "Prop2size" = str_p2s_simulated[,1], 
  "Neyman" = str_neyman_simulated[,1]
)
boxplot(sim_results_str_acres82)
abline(h = mean(agpop$acres92), col = "red")

Figure 3.4: Sampling distribution of the mean under SRS, proportional, and Neyman allocation (acres82 strata, 2000 simulations)

Output Note: The boxplot visually highlights a massive reduction in variance compared to SRS, especially for the Neyman allocation.

sim_means <- sapply(sim_results_str_acres82, mean)
sim_var <- sapply(sim_results_str_acres82, var)
sim_var_relative <- sim_var / sim_var[1]
sim_var_reduction <- 1 - sim_var_relative

data.frame(
  "Mean" = sim_means, 
  "Variance" = sim_var, 
  "Relative Variance" = sim_var_relative, 
  "Percentage of Variance Reduction" = sim_var_reduction, 
  check.names = FALSE
) |> 
  gt(rownames_to_stub = TRUE) |> 
  fmt_number(columns = c(Mean, Variance), decimals = 1) |> 
  fmt_percent(columns = c(`Percentage of Variance Reduction`), decimals = 2) |> 
  fmt_number(columns = c(`Relative Variance`), decimals = 3)

Table 3.6: Simulated mean, variance, and variance reduction by allocation method (acres82 strata)

	Mean	Variance	Relative Variance	Percentage of Variance Reduction
SRS	309.0	546.4	1.000	0.00%
Prop2size	309.0	139.1	0.255	74.54%
Neyman	309.1	35.5	0.065	93.50%

Output Note: The table shows dramatic results: Proportional allocation reduces variance by 74.5%, and Neyman allocation reduces it by 93.5%. Stratifying by a highly correlated variable is immensely powerful.

3.6 Interactive Demonstration: Stratification by acres82 or region

The simulation for this section runs as a standalone app: Stratified Sampling. It opens in a new tab, together with its documentation and the full list of apps.