2  Simple Random Sampling

2.1 Estimators for Simple Random Sampling

A simple random sample (SRS) of \(n\) units is drawn without replacement from a population of \(N\) units, so that every subset of \(n\) units is equally likely. Write \(y_i\) for the value of unit \(i\), \(\bar{y}_U = \frac{1}{N}\sum_{i=1}^N y_i\) for the population mean, \(t = N\bar{y}_U\) for the population total, and \[ S^2 = \frac{1}{N-1}\sum_{i=1}^N (y_i - \bar{y}_U)^2 \] for the population variance.

1. Estimating the Population Mean

The sample mean is unbiased for \(\bar{y}_U\): \[ \bar{y} = \frac{1}{n}\sum_{i=1}^n y_i . \] Its true variance and its estimate, which replaces \(S^2\) by the sample variance \(s^2\), are \[ V(\bar{y}) = \left(1 - \frac{n}{N}\right)\frac{S^2}{n}, \qquad \hat{V}(\bar{y}) = \left(1 - \frac{n}{N}\right)\frac{s^2}{n}, \qquad s^2 = \frac{1}{n-1}\sum_{i=1}^n (y_i - \bar{y})^2 , \] and \(\mathrm{SE}(\bar{y}) = \sqrt{\hat{V}(\bar{y})}\). The factor \(1 - n/N\) is the finite population correction (fpc). It is close to \(1\) when the sampling fraction \(n/N\) is small.

2. Estimating the Population Total

\[ \hat{t} = N\bar{y}, \qquad \hat{V}(\hat{t}) = N^2\,\hat{V}(\bar{y}) = N^2\left(1 - \frac{n}{N}\right)\frac{s^2}{n}, \qquad \mathrm{SE}(\hat{t}) = N\,\mathrm{SE}(\bar{y}). \]

3. Estimating a Proportion

For a 0/1 indicator \(y_i\), the population mean is the proportion \(p\), and the sample mean is the sample proportion \(\hat{p}\). Since then \(s^2 = \frac{n}{n-1}\hat{p}(1-\hat{p})\), \[ \hat{p} = \frac{1}{n}\sum_{i=1}^n y_i, \qquad \hat{V}(\hat{p}) = \left(1 - \frac{n}{N}\right)\frac{\hat{p}(1 - \hat{p})}{n-1}. \]

4. Confidence Intervals

By the central limit theorem for sample means, an approximate \(95\%\) confidence interval is \[ \bar{y} \pm t_{0.975,\,n-1}\,\mathrm{SE}(\bar{y}), \] and likewise for \(\hat{t}\) and \(\hat{p}\). With a skewed population, \(n\) must be larger for the normal approximation, and hence the nominal coverage, to hold. The simulation and the app later in this chapter show this.

---

2.2 Analysis of agsrs.csv Data

Throughout the agriculture examples in this book, county acreages are expressed in thousands of acres — the raw figures run into the millions, and every mean, variance, and total below reads far more easily on that scale. The helper acres_in_thousands() performs the conversion as the data is read.

2.2.1 Step by step calculation without using a function

To evaluate a Simple Random Sample (SRS), we calculate the point estimate, the standard error (incorporating the Finite Population Correction factor), and the margin of error to construct a 95% confidence interval.

Key Formulas:

• Sample Mean: \(\bar{y} = \frac{1}{n} \sum_{i=1}^n y_i\)
• Standard Error: \(\mathrm{SE}(\bar{y}) = \sqrt{1 - \frac{n}{N}} \left( \frac{s}{\sqrt{n}} \right)\)
• Margin of Error: \(\mathrm{ME} = t_{\alpha/2, n-1} \times \mathrm{SE}(\bar{y})\)
• Population Total: \(\hat{t} = N\bar{y}\)

## read survey data
agsrs <- read.csv ("data/agsrs.csv")
## work in thousands of acres, the unit used for every agriculture example here
agsrs[c("acres92", "acres87", "acres82")] <- agsrs[c("acres92", "acres87", "acres82")] / 1000
head(agsrs)

## extract the variable of interest
sdata <- agsrs$acres92
N <- 3078

## do calculation
n <- length (sdata)
ybar <- mean (sdata)
se.ybar <- sqrt((1 - n / N)) * sd (sdata) / sqrt(n)  
mem <- qt (0.975, df = n - 1) * se.ybar

## return estimate vector for pop mean
c (Est. = ybar, S.E. = se.ybar, ci.low = ybar - mem, ci.upp = ybar + mem)

     Est.      S.E.    ci.low    ci.upp 
297.89705  18.89843 260.70626 335.08784 

## return estimate vector for pop total
c (Est. = ybar, S.E. = se.ybar, ci.low = ybar - mem, ci.upp = ybar + mem) * N

      Est.       S.E.     ci.low     ci.upp 
 916927.11   58169.38  802453.86 1031400.36 

2.2.2 Write a function for repeated use

2.2.2.1 A function for doing data analysis for srs sample

Wrap the previous manual calculations into a custom R function to quickly estimate means, standard errors, and confidence intervals for any variable of interest. To avoid repeating this code in every chapter, all estimating functions (and their gt-table formatting helpers) live in a single shared file, samplingestimate.r, which we source below.

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"))
  }
}

2.2.2.2 Apply srs_est to agsrs.csv data

Import Data

Load the survey sample data representing US counties.

agsrs <- acres_in_thousands (read.csv ("data/agsrs.csv"))

Estimating the mean of acre92

Calculate the estimated average farm acreage per county using the custom function. The function’s $table shows the row-level working values (truncated to the first 6 rows + Sum), and $estimate holds the final point estimate, SE, and CI, formatted with format_est_gt().

res_mean <- srs_est(agsrs[,"acres92"], N = 3078, estimate = "mean")
res_mean$table

SRS Mean Estimation: Working Table
\(n = 300\)
\(i\)	\(y_i\)	\(\hat y_i\)	\(e_i\)	\((y_i-\hat y_i)^2\)1
1	175.21	297.90	−122.69	15,052.36
2	138.13	297.90	−159.76	25,523.91
3	56.10	297.90	−241.80	58,464.84
4	199.12	297.90	−98.78	9,757.50
5	89.23	297.90	−208.67	43,542.77
6	96.19	297.90	−201.70	40,684.12
...	...	...	...	...
Sum	89,369.11	297.90	0.00	35,496,086.45
1 \(\hat y_i = \bar y\) for every unit under the SRS model; \(s_y^2 = \sum_i (y_i-\hat y_i)^2/(n-1) = 1.19 \times 10^{5}\).
\(\bar y = \dfrac{\sum_i y_i}{n} = \dfrac{89,369.11}{300} = 297.90\)
\(\mathrm{SE}(\bar y) = \sqrt{1-\dfrac{n}{N}}\,\dfrac{s_y}{\sqrt n} = \sqrt{1-\dfrac{300}{3078}}\,\dfrac{344.55}{\sqrt{300}} = 18.90\)

res_mean$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

Estimating the total of acre92

Set estimate = "total" to have the function itself scale the mean and confidence interval by the population size (\(N = 3078\)) to estimate the total farm acreage across all counties.

res_total <- srs_est(agsrs[,"acres92"], N = 3078, estimate = "total", show.details = FALSE)
res_total$estimate |> format_est_gt("total")

\(\hat{t}\)	\(\mathrm{SE}(\hat{t})\)	95% CI Lower	95% CI Upper
916,927.11	58,169.38	802,453.86	1,031,400.36

Estimating the proportion of counties with fewer than 200K acres for farming in 1992

To estimate a proportion under SRS, create an indicator variable (\(y_i = 1\) if true, \(0\) if false). The mean of this indicator variable serves as the sample proportion (\(\hat{p} = \bar{y}_{indicator}\)).

acres92.is.fewer.200k <- as.numeric (agsrs[,"acres92"] < 200)   # 200 thousand acres
head(acres92.is.fewer.200k)

[1] 1 1 1 1 1 1

res_prop <- srs_est(acres92.is.fewer.200k, N = 3078, estimate = "mean", show.details = FALSE)
res_prop$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
0.51	0.03	0.46	0.56

Estimating the total number of counties with fewer than 200K acres for farming in 1992

Set estimate = "total" to estimate the absolute count of counties meeting the criteria.

res_prop_total <- srs_est(acres92.is.fewer.200k, N = 3078, estimate = "total", show.details = FALSE)
res_prop_total$estimate |> format_est_gt("total")

\(\hat{t}\)	\(\mathrm{SE}(\hat{t})\)	95% CI Lower	95% CI Upper
1,569.78	84.54	1,403.42	1,736.14

2.2.3 Comparing with true value

Load the complete census (population) dataset to calculate the true parameters. This allows us to assess the accuracy of our previous SRS estimates against the exact population values.

agpop <- acres_in_thousands (read.csv ("data/agpop.csv"))
#true mean (thousands of acres)
mean (agpop[, "acres92"], na.rm = T)

[1] 308.5824

# true total (thousands of acres)
sum (agpop[, "acres92"], na.rm = T)

[1] 943953.6

# true proportion of counties with less than 200K acres for farming
mean (agpop[, "acres92"] < 200, na.rm = T)

[1] 0.5145472

# true number of counties with less than 200K acres for farming
sum (agpop[, "acres92"] < 200, na.rm = T)

[1] 1574

2.3 A Simulation Demonstration of SRS Inference

Load the full population data and filter out missing values to serve as a complete, known finite population for resampling experiments.

# read population data, with acreages in thousands of acres
agpop <- acres_in_thousands (read.csv ("data/agpop.csv"))
# remove those counties with na
agpop <- subset( agpop, !is.na (acres92))

2.3.1 True Values

Define the target sample size and compute the true population mean (\(\mu\)) and theoretical standard error (\(\mathrm{SE}_{true}\)) using the known population standard deviation (\(S\)).

Theoretical Standard Error:

\[\mathrm{SE}_{true}(\bar{y}) = \sqrt{1 - \frac{n}{N}} \left(\frac{S}{\sqrt{n}}\right)\]

# sample size
n <- 300
# population size
N <- nrow (agpop); N

[1] 3059

# true value of population mean
ybarU <- mean (agpop[,"acres92"]); ybarU

[1] 308.5824

# true value of deviation of sample mean
true.se.ybar <- sqrt (1- n/N) * sd (agpop[,"acres92"]) / sqrt (n); true.se.ybar

[1] 23.32029

2.3.2 One SRS sampling

Draw a single random sample of size \(n=300\) from the population and run our estimator function to see a single instance of survey inference.

##
# srs sampling
srs <- sample (1:N,n)
head(agpop [srs, ])

# get data of variable "acres92"
sdata <- agpop [srs, "acres92"]
# analysis
res_one <- srs_est(sdata, N, estimate = "mean", show.details = FALSE)
res_one$estimate |> format_est_gt("mean")

\(\bar{y}\)	\(\mathrm{SE}(\bar{y})\)	95% CI Lower	95% CI Upper
313.20	22.15	269.62	356.78

2.3.3 Repeating SRS sampling 5000 times

To demonstrate the Central Limit Theorem and unbiasedness of the estimator, draw 5,000 independent samples, storing the estimates from each iteration. We plot the resulting sampling distribution against the true population mean.

(Note: We use show.details = FALSE and pull out just $estimate here, so the loop skips building 5,000 working tables and stays fast.)

nres <- 5000 # number of repeated sampling
simulation.results <- matrix (0, nres, 4) # matrix recording repeated results
colnames(simulation.results) <- c( "Est.",   "S.E.",   "ci.low", "ci.upp")

for (i in 1:nres)
{
    srs <- sample (N, n)
   sdata <- agpop [srs, "acres92"]
   simulation.results [i,] <- srs_est (sdata, N, estimate = "mean", show.details = FALSE)$estimate
}

head(simulation.results)

         Est.     S.E.   ci.low   ci.upp
[1,] 290.7471 21.78580 247.8741 333.6200
[2,] 300.6398 21.54840 258.2341 343.0456
[3,] 286.4379 17.31731 252.3586 320.5172
[4,] 333.7544 22.35838 289.7547 377.7541
[5,] 304.5797 18.34546 268.4771 340.6822
[6,] 328.4969 25.07866 279.1439 377.8499

# look at the distribution of sample mean
par (mfrow= c(2,2))
hist (agpop$acres92,main = "Population Distribution of acre92")
hist (simulation.results[,1], main = "Sampling Distribution of Sample Mean for acre92")
abline (v = ybarU, col = "red")
qqnorm (simulation.results[,1], main="QQ plot of Sample Mean"); qqline(simulation.results[,1])
boxplot (simulation.results[,1], main = "Boxplot of Sample Mean")
abline (h = ybarU, col = "red")

mean (simulation.results[,1])

[1] 309.1513

ybarU

[1] 308.5824

sd (simulation.results [,1])

[1] 23.5814

true.se.ybar

[1] 23.32029

Figure 2.1: Population distribution and sampling distribution of the SRS sample mean over 5000 repeated samples (n = 300)

2.3.4 Empirical Coverage Rate of CIs

Verify the reliability of the 95% confidence intervals by checking whether the true population mean falls within the calculated bounds for each of the 5,000 samples. The empirical coverage rate should closely approximate 0.95. The plot highlights the first 100 intervals, flagging those that fail to capture the true mean in a distinct color.

library(dplyr)
library(plotrix)

# Convert to data frame (if it was a matrix) and use mutate to evaluate coverage
simulation.results <- as.data.frame(simulation.results) |>
  mutate(
    `Covered?` = as.numeric(ci.low < ybarU & ybarU < ci.upp)
  )

head(simulation.results)

par(mfrow=c(1,1))
plotCI(
  x = 1:100,
  y = simulation.results$Est.[1:100], 
  li = simulation.results$ci.low[1:100],
  ui = simulation.results$ci.upp[1:100],
  col = 2 - simulation.results$`Covered?`[1:100]
)
abline(h = ybarU, col = "blue")

# Empirical coverage rate
mean(simulation.results$`Covered?`)

[1] 0.9312

Figure 2.2: First 100 simulated 95% confidence intervals for the SRS sample mean, with coverage failures highlighted

2.4 Shinylive App for Illustrating SRS Theory

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

2.5 Prediction-Based Theory for Sampling (Optional)

Sampling is a prediction problem. The sample \(S\) plays the role of a training set, the non-sampled units \(R=U\setminus S\) play the role of test cases, and a population total or mean is a sum of predictions for the test cases plus the observed training values. This section states the framework in general, with a model for every unit and a generic estimator of the model parameters, and then specializes it to simple random sampling, where the mean of every unit is the same and is estimated by \(\bar y_S\).

2.5.1 Setup

Let \(U=\{1,\dots,N\}\) index a finite population with study variable \(Y_i\), let \(S\subset U\) with \(|S|=n\) denote the observed sample, and let \(R=U\setminus S\) denote the \(N-n\) non-sampled units. The quantities of interest are \[ T=\sum_{i\in U}Y_i=\sum_{i\in S}Y_i+\sum_{i\in R}Y_i,\qquad \bar Y_U=T/N . \] Only the second sum is unknown, so the task is to predict \(Y_i\) for each \(i\in R\).

The prediction approach treats \(Y_1,\dots,Y_N\) as random variables and \(S\) as fixed. All expectations are taken with respect to the model \(\xi\) rather than the randomization distribution. The sampling design enters only through the requirement that selection be uninformative, so that \(\xi\) governs the units in \(R\) as well as those in \(S\). The general working model is \[ \xi:\quad E(Y_i)=\mu_i(\theta),\quad \operatorname{Var}(Y_i)=\sigma_i^2,\quad \operatorname{Cov}(Y_i,Y_j)=0\ (i\neq j), \tag{2.1}\] where \(\theta\) is an unknown parameter and \(\mu_i(\theta)\) is the mean function, for example \(\mu_i(\theta)=x_i^{\top}\beta\) in regression. The model errors are \[ e_i=Y_i-\mu_i(\theta),\qquad E(e_i)=0,\quad E(e_i^2)=\sigma_i^2,\quad E(e_ie_j)=0\ (i\neq j). \tag{2.2}\]

A prediction rule estimates \(\theta\) from the sample by \(\hat\theta=\hat\theta(Y_S)\) and predicts \[ \hat Y_i= \begin{cases} Y_i, & i\in S,\\ \mu_i(\hat\theta), & i\in R, \end{cases} \qquad \hat T=\sum_{i\in U}\hat Y_i=\sum_{i\in S}Y_i+\sum_{i\in R}\mu_i(\hat\theta). \tag{2.3}\] Observed values are kept as they are, and the fitted mean function is used only for the unobserved units. Because \(\hat\theta\) is a function of \(Y_S\) alone and the errors are uncorrelated, \(\hat\theta\) is uncorrelated with \(e_i\) for every \(i\in R\). This is the only property of the estimator used below.

2.5.2 Two sources of error

Individual prediction. In machine learning a predictor is judged by its mean squared error on a test case. For \(i\in R\) add and subtract \(\mu_i(\theta)\): \[ \hat Y_i-Y_i=\underbrace{\mu_i(\hat\theta)-\mu_i(\theta)}_{\delta_i}-\underbrace{\big(Y_i-\mu_i(\theta)\big)}_{e_i}. \tag{2.4}\] The first part, \(\delta_i\), is the error in the fitted mean caused by the uncertainty in \(\hat\theta\), and the second part, \(e_i\), is the gap between the unobserved \(Y_i\) and its true mean. The cross term \(E(\delta_ie_i)\) vanishes, so \[ E\big(\hat Y_i-Y_i\big)^2=\underbrace{\sigma_i^2}_{\text{irreducible}}+\underbrace{E\big[\mu_i(\hat\theta)-\mu_i(\theta)\big]^2}_{\text{estimation error}} . \tag{2.5}\] The irreducible term is present even when \(\theta\) is known, because no predictor removes the scatter of \(Y_i\) about its mean. The estimation term is the mean squared error of \(\mu_i(\hat\theta)\) as an estimator of \(\mu_i(\theta)\) and decreases as the training sample grows. This is the usual decomposition of test error into noise and estimation variance, plus squared bias if the estimator is biased.

The total. The observed units contribute no error, so \[ \hat T-T=\sum_{i\in R}(\hat Y_i-Y_i)=\sum_{i\in R}\delta_i-\sum_{i\in R}e_i . \tag{2.6}\] The two sources of error appear again, now summed over the test cases. Squaring, the cross term between \(\sum_R\delta_i\) and \(\sum_R e_i\) vanishes by the same argument as before, and the \(e_i\) are uncorrelated, so \[ \operatorname{MSPE}(\hat T)=E\big[(\hat T-T)^2\big] =\underbrace{\sum_{i\in R}\sigma_i^2}_{\text{scatter of the unobserved }Y_i} +\underbrace{E\Big[\Big(\sum_{i\in R}\delta_i\Big)^{2}\Big]}_{\text{uncertainty in }\hat\theta} . \tag{2.7}\] Dividing by \(N^2\) gives \(\operatorname{MSPE}(\hat{\bar y}_U)\).

Comparing Equation 2.5 and Equation 2.7 shows what changes when the target moves from a single case to the total. The irreducible parts simply add, \(\sum_{i\in R}\sigma_i^2\). The estimation parts do not add, because the errors \(\delta_i\) all come from the same \(\hat\theta\) and are strongly correlated. Expanding the square gives \[ E\Big[\Big(\sum_{i\in R}\delta_i\Big)^{2}\Big]=\sum_{i\in R}E(\delta_i^2)+\sum_{i\neq j\in R}E(\delta_i\delta_j), \] and the double sum is typically of larger order than the first. A total is a more aggregated target than an individual prediction, and its error is dominated by the shared error in \(\hat\theta\), which does not average out across units.

2.5.3 Specialization to simple random sampling

For simple random sampling there are no auxiliary variables, so every unit has the same mean, \[ \mu_i(\theta)=\mu,\quad \sigma_i^2=\sigma^2,\qquad e_i=Y_i-\mu,\quad \bar e_A=|A|^{-1}\sum_{i\in A}e_i, \] and \(\theta=\mu\) is estimated by the sample mean \(\bar y_S=n^{-1}\sum_{i\in S}Y_i\). The predictor Equation 2.3 is \[ \hat Y_i= \begin{cases} Y_i, & i\in S,\\ \bar y_S, & i\in R, \end{cases} \qquad \hat T=\sum_{i\in S}Y_i+(N-n)\bar y_S=N\bar y_S , \tag{2.8}\] so the expansion estimator is a prediction rule that fills in every unobserved unit with the sample mean. Here the estimation error is common to all units, \(\delta_i=\bar y_S-\mu=\bar e_S\), and Equation 2.6 becomes \[ \hat T-T=(N-n)\,\bar e_S-\sum_{i\in R}e_i . \tag{2.9}\] Both terms have mean zero, so \(\hat T\) is model-unbiased.

Individual units. From Equation 2.5, with \(E(\bar e_S^{\,2})=\sigma^2/n\), \[ E\big(\hat Y_i-Y_i\big)^2=\sigma^2+\frac{\sigma^2}{n}=\sigma^2\Big(1+\frac1n\Big),\qquad i\in R . \tag{2.10}\]

The total. Since \(\delta_i=\bar e_S\) for all \(i\in R\), the double sum in Equation 2.7 is as large as it can be: \(E[(\sum_R\delta_i)^2]=(N-n)^2E(\bar e_S^{\,2})=(N-n)^2\sigma^2/n\), whereas independent errors would give \((N-n)\sigma^2/n\). The scatter term is \((N-n)\sigma^2\), and therefore \[ \operatorname{MSPE}(\hat T)=\frac{(N-n)^2}{n}\sigma^2+(N-n)\sigma^2 =\frac{N(N-n)}{n}\,\sigma^2 , \tag{2.11}\] \[ \operatorname{MSPE}(\hat{\bar y}_U)=\frac{\sigma^2}{n}\left(1-\frac{n}{N}\right). \tag{2.12}\]

Expression Equation 2.12 is the classical simple-random-sampling variance with \(\sigma^2\) in place of the population variance \(S^2\), obtained without appeal to randomization. Four consequences follow.

1. The two terms of Equation 2.11 are the two sources of error. The term \((N-n)\sigma^2\) is the irreducible scatter of the unobserved units; the term \((N-n)^2\sigma^2/n\) is the cost of estimating \(\mu\) from \(n\) observations and alone decreases with \(n\).
2. The per-unit error \(\sigma^2(1+1/n)\) is barely above \(\sigma^2\) for moderate \(n\), so the estimation error looks negligible for individual prediction. For the total it is not: the estimation error is shared by all \(N-n\) unobserved units and accumulates to a factor \((N-n)^2\), whereas the scatter accumulates only to \(N-n\). The ratio of the two terms is \((N-n)/n\).
3. The finite-population correction is not an artifact of sampling without replacement. The \(n\) observed units contribute exactly zero to Equation 2.9, leaving \(N-n\) terms to predict; at \(n=N\) the MSPE is zero.
4. If \(\mu\) is known, \(\delta_i=0\) and the predictor \(\hat Y_i=\mu\) for \(i\in R\) has \(\operatorname{MSPE}(\hat T)=(N-n)\sigma^2\). This is the irreducible error, and the cost of not knowing \(\mu\) is the factor \(N/n\) applied to it. For the population mean the MSPE rises from \(\sigma^2(N-n)/N^2\) to approximately \(\sigma^2/n\), so precision is governed by the sample size rather than the population size.

Other sampling problems fit the same framework by changing \(\mu_i(\theta)\): a ratio or regression estimator uses \(\mu_i(\theta)=\beta x_i\) or \(\beta_0+\beta_1x_i\), and a stratified estimator uses a separate mean \(\mu_h\) in each stratum. Only \(\hat\theta\) and the shared estimation error \(\sum_R\delta_i\) change.

2.5.4 Illustration

Figure 2.3 displays one realization of a population under Equation 2.1 with \(\mu_i=\mu\). Since the labels are uninformative, the sample is placed on the left, shaded grey, and the units to be predicted on the right.

Figure 2.3: One realized population under the model. The sample occupies the shaded region on the left; the remaining units are to be predicted. Horizontal lines: the model mean \(\mu\) (orange), the estimand \(\bar Y_U\) (red), and the predictor \(\bar y_S\) (blue, dashed).

The scatter of the open points about the orange line is the irreducible error \(\sigma_i^2\) in Equation 2.5: it is the error of predicting an individual unit even with \(\mu\) known. The gap between the orange and blue lines is the estimation error \(\delta_i=\bar e_S\), with root mean square \(\sigma/\sqrt{n}\), which is the same for every unit in \(R\). The gap between the blue and red lines, scaled by \(N\), is the realized prediction error of the total; Equation 2.11 states that its root mean square is \(\sigma\sqrt{N(N-n)/n}\).

c(mu = mu, Ybar_U = ybar_U, ybar_S = ybar_S,
  error_total  = N * (ybar_S - ybar_U),
  rmspe_total  = sigma * sqrt(N * (N - n) / n))

         mu      Ybar_U      ybar_S error_total rmspe_total 
   50.00000    48.11714    59.15883  1104.16893   600.00000 

Estimating \(\sigma^2\) by the sample variance \(s_S^2\), which satisfies \(E(s_S^2)=\sigma^2\) under Equation 2.1, yields the usual variance estimator and recovers the design-based formulas.