Code
cherry <- read.csv("data/cherry.csv", header = T)
head(cherry) |> gt()| diameter | height | volume |
|---|---|---|
| 8.3 | 70 | 10.3 |
| 8.6 | 65 | 10.3 |
| 8.8 | 63 | 10.2 |
| 10.5 | 72 | 16.4 |
| 10.7 | 81 | 18.8 |
| 10.8 | 83 | 19.7 |
4 Ratio and Regression Estimation
4.1 Ratio, Regression, and Domain Estimation
When auxiliary information (a variable \(x\)) is highly correlated with our target variable of interest (\(y\)), we can use Ratio or Regression estimation to improve precision without increasing the sample size. We can also use auxiliary population structures to adjust our estimates after the fact via Post-stratification, or estimate means for specific sub-populations (Domains).
Throughout, \((x_i, y_i)\), \(i = 1, \ldots, n\), is an SRS of size \(n\) from a population of \(N\) units, and \(\overline{x}_U\) (equivalently \(t_x = N\overline{x}_U\)) is known.
1. Ratio Estimation
If the relationship between \(y\) and \(x\) passes through the origin, we estimate the ratio \(B = \overline{y}_U/\overline{x}_U\) and use it to adjust the mean: \[ \hat{B} = \frac{\overline{y}}{\overline{x}} = \frac{\sum_{i=1}^n y_i}{\sum_{i=1}^n x_i}, \qquad \overline{y}_r = \hat{B}\,\overline{x}_U, \qquad \hat{t}_r = \hat{B}\,t_x = N\overline{y}_r . \] Its error is driven by the residuals \(e_i = y_i - \hat{B}x_i\) from the line through the origin: \[ \hat{V}(\hat{B}) = \frac{1}{\overline{x}^2}\left(1 - \frac{n}{N}\right)\frac{s_e^2}{n}, \qquad \hat{V}(\overline{y}_r) = \overline{x}_U^2\,\hat{V}(\hat{B}) = \left(\frac{\overline{x}_U}{\overline{x}}\right)^2 \left(1 - \frac{n}{N}\right) \frac{s_e^2}{n}, \qquad s_e^2 = \frac{1}{n-1}\sum_{i=1}^n \left(y_i - \hat{B}x_i\right)^2 , \] with \(\hat{V}(\hat{t}_r) = N^2\hat{V}(\overline{y}_r)\) and \(95\%\) intervals using \(t_{0.975,\,n-1}\). The ratio estimator is slightly biased, but the bias is negligible when \(n\) is large.
\(e_i\) is the vertical distance from the point \((x_i, y_i)\) to the fitted line through the origin, \(y = \hat{B}x\). No centring appears in \(s_e^2\) because \(\hat{B} = \overline{y}/\overline{x}\) forces \(\sum_i e_i = \sum_i y_i - \hat{B}\sum_i x_i = 0\) exactly, so the residuals already average to zero; the divisor is \(n-1\) for the one parameter \(\hat{B}\) estimated from the sample. Compare the variance with the SRS one, \(\left(1-\frac{n}{N}\right)s_y^2/n\): ratio estimation replaces the scatter of \(y\) about \(\overline{y}\) by the scatter about the line. That is the whole source of the gain — and the whole risk, since an \(x\) that predicts \(y\) poorly gives \(s_e^2 > s_y^2\) and a less precise estimate than the plain sample mean.
2. Regression Estimation
If the relationship does not pass through the origin, a linear regression is more appropriate. With the OLS coefficients \[ \hat{B}_1 = \frac{\sum_{i=1}^n (x_i - \overline{x})(y_i - \overline{y})}{\sum_{i=1}^n (x_i - \overline{x})^2}, \qquad \hat{B}_0 = \overline{y} - \hat{B}_1\overline{x}, \] the regression estimator of the mean is the fitted line evaluated at \(\overline{x}_U\): \[ \overline{y}_{reg} = \hat{B}_0 + \hat{B}_1 \overline{x}_U = \overline{y} + \hat{B}_1\left(\overline{x}_U - \overline{x}\right), \qquad \hat{t}_{reg} = N\overline{y}_{reg} . \] Its estimated variance uses the residual mean square of the regression: \[ \hat{V}(\overline{y}_{reg}) = \left(1 - \frac{n}{N}\right) \frac{s_e^2}{n}, \qquad s_e^2 = \frac{1}{n-2}\sum_{i=1}^n \left(y_i - \hat{B}_0 - \hat{B}_1 x_i\right)^2 , \] with \(\hat{V}(\hat{t}_{reg}) = N^2\hat{V}(\overline{y}_{reg})\) and \(95\%\) intervals using \(t_{0.975,\,n-2}\). The divisor is now \(n-2\), since both \(\hat{B}_0\) and \(\hat{B}_1\) are estimated. For the same data \(\sum_i e_i^2\) can only be smaller here than under the ratio fit — the ratio line is the special case \(\hat{B}_0 = 0\) — so the regression estimator is the safer of the two whenever the relationship misses the origin, at the cost of one degree of freedom. Note also that the factor \((\overline{x}_U/\overline{x})^2\) is absent: the fitted line is evaluated at \(\overline{x}_U\) directly, so no correction for \(\overline{x} \neq \overline{x}_U\) is needed.
3. Domain Estimation
To estimate the mean of a subgroup (domain \(d\)) when the domain sample size \(n_d\) is a random variable, use the \(n_d\) sampled units that fall in the domain: \[ \overline{y}_d = \frac{1}{n_d} \sum_{i \in S_d} y_i, \qquad \hat{V}(\overline{y}_d) = \left(1 - \frac{n}{N}\right) \frac{s_d^2}{n_d}, \qquad s_d^2 = \frac{1}{n_d - 1}\sum_{i \in S_d}\left(y_i - \overline{y}_d\right)^2 , \] where \(S_d\) is the set of sampled units in domain \(d\), \(n\) is the overall sample size, and \(N\) is the population size. \(95\%\) intervals use \(t_{0.975,\,n_d-1}\). The fpc uses the overall sampling fraction \(n/N\), because the domain’s population size \(N_d\) is usually unknown. \(\overline{y}_d\) is in fact a ratio estimator, \(\sum_{i=1}^n u_i/\sum_{i=1}^n x_i\) with \(x_i = 1\) if unit \(i\) is in the domain and \(u_i = x_i y_i\).
4. Post-stratification
If we know the true population sizes \(N_h\) of \(H\) groups (post-strata), but used a Simple Random Sample, we can post-stratify to reduce variance. We weight the group sample means by the known \(\pi_h = N_h/N\), as in stratified sampling: \[ \overline{y}_{post} = \sum_{h=1}^H \pi_h\,\overline{y}_h . \] The observed group sizes \(n_h\) are random, so for the variance we substitute their expected values under SRS, the proportional sample sizes \[ n_h^{post} = n\,\frac{N_h}{N} = n\pi_h , \] into the stratified variance formula. Every group then has sampling fraction \(n/N\), and the variance becomes \[ \hat{V}(\overline{y}_{post}) = \sum_{h=1}^H \pi_h^2\left(1 - \frac{n_h^{post}}{N_h}\right)\frac{s_h^2}{n_h^{post}} = \left(1 - \frac{n}{N}\right)\frac{1}{n}\sum_{h=1}^H \pi_h s_h^2 , \] the variance of proportional stratified sampling. With \(\hat{t}_{post} = N\overline{y}_{post}\), \(\hat{V}(\hat{t}_{post}) = N^2\hat{V}(\overline{y}_{post})\). Here \(n_h\), \(\overline{y}_h\) and \(s_h^2\) are the count, mean and variance of the units the sample happened to place in group \(h\), so every group needs \(n_h \geq 2\) for \(s_h^2\) to exist. Substituting \(n_h^{post}\) treats the \(n_h\) as fixed at \(n\pi_h\); carrying their randomness adds a second-order term \(\frac{1}{n^2}\sum_h (1-\pi_h)s_h^2\), negligible unless some \(n_h\) are small.
Post-stratification Viewed as a Special Regression Estimation
Post-stratification is not a separate technique: it is the regression estimator of Section 2 with the group indicators as the auxiliary variables. Take the \(H\) dummies \(x_{ih} = \mathbf{1}\{i \in \text{group } h\}\) and fit without an intercept. Then
- the fitted coefficients are the group means, \(\hat{B}_h = \overline{y}_h\);
- the population means of those dummies are the known weights, \(\overline{x}_{U,h} = N_h/N = \pi_h\);
so the regression estimator is the post-stratified one, \[ \overline{y}_{reg} = \sum_h \hat{B}_h\,\overline{x}_{U,h} = \sum_h \pi_h \overline{y}_h = \overline{y}_{post} . \]
Its residuals are each unit’s deviation from its own group mean, \(e_i = y_i - \overline{y}_{h(i)}\), and the general variance \(\left(1-\frac{n}{N}\right)s_e^2/n\) applies with \(p = H\) fitted parameters: \[ s_e^2 = \frac{1}{n-H}\sum_h \sum_{i \in h} \left(y_i - \overline{y}_h\right)^2 = \frac{1}{n-H}\sum_h (n_h - 1)\,s_h^2 . \]
This and the post-stratification variance are the same weighted average of the same \(s_h^2\); only the weights differ, \[ \underbrace{\frac{n_h-1}{n-H}}_{\text{sample proportions}} \qquad \text{versus} \qquad \underbrace{\pi_h = \frac{N_h}{N}}_{\text{population proportions}} , \] so the substitution \(n_h \to n\pi_h\) above is simply: pool the residuals with the weights we know, not with the ones the sample happened to deliver. It is the same reasoning that puts \(\overline{x}_U\), rather than \(\overline{x}\), in the ratio estimator.
Seen this way, the estimators above are one estimator with different auxiliary variables, and \(s_e^2 = \frac{1}{n-p}\sum_i e_i^2\) throughout, where \(p\) is the number of fitted parameters:
| Estimator | Auxiliary variables | Residual \(e_i\) | \(p\) |
|---|---|---|---|
| Ratio | one \(x\), line forced through the origin | \(y_i - \hat{B}x_i\) | 1 |
| Regression | intercept and \(x\) | \(y_i - \hat{B}_0 - \hat{B}_1 x_i\) | 2 |
| Post-stratification | \(H\) group indicators | \(y_i - \overline{y}_h\) | \(H\) |
Domain estimation stands slightly apart: there the divisor is the random \(n_d\) rather than the full \(n\), because only the units falling in the domain carry information about it.
4.2 Functions and packages for Analyzing Data
We begin by loading the necessary libraries, including gt for formatting mathematical tables, and defining our custom functions for Ratio, Regression, Domain, and Post-stratified estimation. We also define a versatile gt helper to render our estimation outputs cleanly.
4.2.1 General functions for ratio and regression estimation
The estimating functions used throughout this book (including reg_est, ratio_est, domain_est, srs_est, and the format_est_gt table-formatting helper used below) live in a single shared file, samplingestimate.r, which we source below.
For post-stratification, we reuse the shared str_est() stratified estimator, but we plug in the expected sample sizes \(n_h = n \times \frac{N_h}{N}\) instead of the actual observed random \(n_h\).
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"))
}
}4.3 Ratio Estimation for cherry.csv dataset
We want to estimate the volume of wood in cherry trees using their diameter as an auxiliary variable.
4.3.1 Importing data
First, let’s load the data and verify that a strong linear relationship passing near the origin exists between volume and diameter.
Code
plot(cherry$volume ~ cherry$diameter, main = "Volume vs Diameter")
Code
anova(lm(cherry$volume ~ 0 + cherry$diameter))Output Note: The plot and the linear model output (which forces a 0 intercept) confirm a very strong relationship, making this a prime candidate for Ratio estimation.
4.3.2 SRS estimate
To establish a baseline, we first estimate the mean and total volume using standard Simple Random Sampling (ignoring the auxiliary variable).
Code
N <- 2967
## estimating the mean of volume
res_srs_vol <- srs_est(cherry$volume, N = N, estimate = "mean", show.details = FALSE)
srs_mean_vol_raw <- res_srs_vol$estimate
srs_mean_vol_raw |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 30.17 | 2.94 | 24.17 | 36.17 |
Code
## estimating the total of volume
res_srs_total_vol <- srs_est(cherry$volume, N = N, estimate = "total", show.details = FALSE)
srs_total_vol_raw <- res_srs_total_vol$estimate
srs_total_vol_raw |> format_est_gt("total")| \(\hat{t}\) | \(\mathrm{SE}(\hat{t})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 89,517.26 | 8,713.67 | 71,721.58 | 107,312.94 |
Output Note: Under SRS, the standard error for the mean volume is 2.94.
4.3.3 Ratio Estimation Using the Function
We calculate the ratio \(\hat{B} = \overline{y}/\overline{x}\); the function’s $table shows the row-level working values \(x_i\), \(y_i\), fitted \(\hat{y}_i = \hat{B}x_i\), residuals \(e_i = y_i - \hat{y}_i\), and \(e_i^2\) (with a Sum row and a footnote giving \(\hat B\) and \(s_e^2\)), so we no longer need to build this table by hand.
Code
ydata <- cherry$volume
xdata <- cherry$diameter
N <- 2967
res_ratio_model <- ratio_est(ydata, xdata, N = N, estimate = "model")
res_ratio_model$table| Ratio Estimation: Working Table | |||||
| \(n = 31\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 10.30 | 8.30 | 18.90 | −8.60 | 73.99 |
| 2 | 10.30 | 8.60 | 19.59 | −9.29 | 86.21 |
| 3 | 10.20 | 8.80 | 20.04 | −9.84 | 96.84 |
| 4 | 16.40 | 10.50 | 23.91 | −7.51 | 56.43 |
| 5 | 18.80 | 10.70 | 24.37 | −5.57 | 31.00 |
| 6 | 19.70 | 10.80 | 24.60 | −4.90 | 23.96 |
| ... | ... | ... | ... | ... | ... |
| Sum | 935.30 | 410.70 | 935.30 | 0.00 | 2,821.59 |
| 1 \(\hat B = \sum_i y_i / \sum_i x_i = 2.28\), \(\hat y_i = \hat B\, x_i\), \(s_e^2 = 94.05\). | |||||
| \(\hat B = \dfrac{\sum_i y_i}{\sum_i x_i} = \dfrac{935.30}{410.70} = 2.28\) | |||||
| \(\mathrm{SE}(\hat B) = \dfrac{1}{\bar x}\sqrt{\left(1-\dfrac{n}{N}\right)\dfrac{s_e^2}{n}} = \dfrac{1}{13.25}\sqrt{\left(1-\dfrac{31}{2967}\right)\dfrac{94.05}{31}} = 0.13\) | |||||
Code
plot(cherry$volume ~ cherry$diameter, main="Ratio Estimation Fit")
abline(a = 0, b = res_ratio_model$estimate["Est."], col="blue", lwd=2)
Output Note: The blue line represents our estimated ratio \(\hat{B} \approx 2.28\). The working table above shows the individual residuals used to estimate \(s_e^2\).
Code
B_v2d <- res_ratio_model$estimate
B_v2d |> format_est_gt("ratio")| \(\hat{B}\) | \(\mathrm{SE}(\hat{B})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 2.28 | 0.13 | 2.01 | 2.54 |
Output Note: The ratio \(\hat{B}\) is estimated at 2.277 with a standard error of 0.1308.
4.3.3.1 Estimating the mean volume of wood
If we know the true total of the population diameters is 41,835, we can derive the population mean diameter \(\overline{x}_{\text{U}}\) and pass it as xbarU with estimate = "mean".
Code
mean_diameters <- 41835 / N
res_ratio_mean <- ratio_est(ydata, xdata, xbarU = mean_diameters, N = N, estimate = "mean")
res_ratio_mean$table| Ratio Estimation: Working Table | |||||
| \(n = 31\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 10.30 | 8.30 | 18.90 | −8.60 | 73.99 |
| 2 | 10.30 | 8.60 | 19.59 | −9.29 | 86.21 |
| 3 | 10.20 | 8.80 | 20.04 | −9.84 | 96.84 |
| 4 | 16.40 | 10.50 | 23.91 | −7.51 | 56.43 |
| 5 | 18.80 | 10.70 | 24.37 | −5.57 | 31.00 |
| 6 | 19.70 | 10.80 | 24.60 | −4.90 | 23.96 |
| ... | ... | ... | ... | ... | ... |
| Sum | 935.30 | 410.70 | 935.30 | 0.00 | 2,821.59 |
| 1 \(\hat B = \sum_i y_i / \sum_i x_i = 2.28\), \(\hat y_i = \hat B\, x_i\), \(s_e^2 = 94.05\). | |||||
| \(\bar y_r = \hat B\,\bar x_U = 2.28 \times 14.10 = 32.11\) | |||||
| \(\mathrm{SE}(\bar y_r) = \mathrm{SE}(\hat B)\,\bar x_U = 0.13 \times 14.10 = 1.84\) | |||||
Code
ratio_mean_volume_raw <- res_ratio_mean$estimate
ratio_mean_volume_raw |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 32.11 | 1.84 | 28.34 | 35.88 |
4.3.3.2 Estimating the total volume of wood
Using estimate = "total" with the known population total diameter gives the estimated population total volume directly.
Code
t_diameters <- 41835
res_ratio_total <- ratio_est(ydata, xdata, xbarU = t_diameters / N, N = N, estimate = "total")
res_ratio_total$table| Ratio Estimation: Working Table | |||||
| \(n = 31\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 10.30 | 8.30 | 18.90 | −8.60 | 73.99 |
| 2 | 10.30 | 8.60 | 19.59 | −9.29 | 86.21 |
| 3 | 10.20 | 8.80 | 20.04 | −9.84 | 96.84 |
| 4 | 16.40 | 10.50 | 23.91 | −7.51 | 56.43 |
| 5 | 18.80 | 10.70 | 24.37 | −5.57 | 31.00 |
| 6 | 19.70 | 10.80 | 24.60 | −4.90 | 23.96 |
| ... | ... | ... | ... | ... | ... |
| Sum | 935.30 | 410.70 | 935.30 | 0.00 | 2,821.59 |
| 1 \(\hat B = \sum_i y_i / \sum_i x_i = 2.28\), \(\hat y_i = \hat B\, x_i\), \(s_e^2 = 94.05\). | |||||
| \(\hat t = \hat B\,\bar x_U\,N = 2.28 \times 14.10 \times 2,967 = 95,272.16\) | |||||
| \(\mathrm{SE}(\hat t) = \mathrm{SE}(\hat B)\,\bar x_U\,N = 0.13 \times 14.10 \times 2,967 = 5,471.43\) | |||||
Code
ratio_total_volume <- res_ratio_total$estimate
ratio_total_volume |> format_est_gt("total")| \(\hat{t}\) | \(\mathrm{SE}(\hat{t})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 95,272.16 | 5,471.43 | 84,098.00 | 106,446.32 |
4.3.3.3 Percentage of Variance Reduction
How much did ratio estimation help compared to ignoring the auxiliary variable?
Code
variance_reduction <- 1 - (ratio_mean_volume_raw["S.E."] / srs_mean_vol_raw["S.E."])^2
unname(variance_reduction)[1] 0.6057237Output Note: Using the diameter as an auxiliary variable reduced the variance of our estimate by 60.6% compared to simple random sampling.
4.4 Simulation study with agpop.csv
Let’s test this practically by drawing 2000 random samples and seeing how SRS and Ratio estimators behave under repeated sampling.
4.4.1 Information of population data
Code
agpop <- acres_in_thousands(read.csv("data/agpop.csv"))
## acres87 is the auxiliary variable, so drop counties missing either one
agpop <- agpop[!is.na(agpop$acres92) & !is.na(agpop$acres87), ]
# sample size
n <- 300
# population size
N <- nrow(agpop)
# true values that we want to estimate
tyU <- sum(agpop[, "acres92"])
# suppose known for ratio estimate
txU <- sum(agpop[, "acres87"])
B <- tyU / txUCode
plot(agpop[, "acres87"], agpop[, "acres92"], main="Acres 1992 vs Acres 1987")
abline(a = 0, b = B, col="red")
Code
# linear model output
lm_gt <- format_lm_gt(lm(agpop$acres92 ~ agpop$acres87))
lm_gt$coef| Coefficient Estimates | ||||
| Estimate | S.E. | t value | p value | |
|---|---|---|---|---|
| (Intercept) | −2.4509 | 0.8517 | −2.8775 | 0.0040 |
| agpop$acres87 | 0.9882 | 0.0016 | 618.4281 | 0.0000 |
Code
lm_gt$fit| Model Fit Summary | ||||
| \(R^2\) | Adj. \(R^2\) | \(\hat\sigma\) | F-statistic | p-value |
|---|---|---|---|---|
| 0.9921 | 0.9921 | 37.8382 | 382,453.3759 | 0.0000 |
Code
# expected reduction of variance of ratio estimate to srs
exp_var_red <- 1 - var(agpop$acres92 - B * agpop$acres87) / var(agpop$acres92)
exp_var_red[1] 0.9920478Output Note: Past acreage (1987) strongly predicts current acreage (1992). The theoretically expected variance reduction is roughly 99.2%.
4.4.2 Simulation Studies
Code
nsim <- 2000
sim_rat <- data.frame(simple = rep(0, nsim), ratio = rep(0, nsim), B = rep(0, nsim))
for (i in 1:nsim) {
srs <- sample(N, n)
y_srs <- agpop[srs, "acres92"]
x_srs <- agpop[srs, "acres87"]
# SRS estimate (Total)
sim_rat$simple[i] <- mean(y_srs) * N
# ratio estimate (Total)
sim_rat$B[i] <- mean(y_srs) / mean(x_srs)
sim_rat$ratio[i] <- sim_rat$B[i] * txU
}
boxplot(sim_rat[, 1:2], main="Estimator Variance: SRS vs Ratio")
abline(h = tyU, col = "red", lwd=2)
We now calculate the actual variance across our 2000 simulations.
Code
sim_var <- sapply(sim_rat[, 1:2], var)
ratio_sim_var <- sim_var / sim_var[1]
percentage_reduction <- 1 - ratio_sim_var
data.frame(sim_var, ratio_sim_var, percentage_reduction) |>
gt(rownames_to_stub = TRUE) |>
cols_label(
sim_var = "Variance",
ratio_sim_var = "Relative Variance",
percentage_reduction = "Variance Reduction"
) |>
fmt_number(columns = c(sim_var), decimals = 0) |>
fmt_percent(columns = c(percentage_reduction), decimals = 2) |>
fmt_number(columns = c(ratio_sim_var), decimals = 4)| Variance | Relative Variance | Variance Reduction | |
|---|---|---|---|
| simple | 4,887,704,641 | 1.0000 | 0.00% |
| ratio | 38,193,285 | 0.0078 | 99.22% |
Output Note: The boxplot dramatically illustrates the precision gain. The variance of the ratio estimator is only 0.8% of the SRS variance, confirming our theoretical expectations.
4.5 Estimating Domain Means and Post-stratification with agsrs.csv example
4.5.1 Estimating Domain Means
When estimating the mean for a specific subgroup (like a region), the sample size in that subgroup (\(n_d\)) is a random variable, which increases our uncertainty.
Code
agsrs <- acres_in_thousands(read.csv("data/agsrs.csv"))
n <- nrow(agsrs)
## estimating domain means
region_agsrs_domain <- data.frame(rbind(
NC = domain_est(agsrs$acres92[agsrs$region=="NC"], n=n, N = 3078),
NE = domain_est(agsrs$acres92[agsrs$region=="NE"], n=n, N = 3078),
S = domain_est(agsrs$acres92[agsrs$region=="S"], n=n, N = 3078),
W = domain_est(agsrs$acres92[agsrs$region=="W"], n=n, N = 3078)
))
region_agsrs_domain |> format_est_gt("domain_mean")| \(\bar{y}_{d}\) | \(\mathrm{SE}(\bar{y}_{d})\) | 95% CI Lower | 95% CI Upper | |
|---|---|---|---|---|
| NC | 323.42 | 30.01 | 263.66 | 383.17 |
| NE | 106.55 | 28.96 | 45.45 | 167.64 |
| S | 246.19 | 23.50 | 199.77 | 292.61 |
| W | 518.98 | 70.30 | 377.21 | 660.75 |
If we incorrectly used the standard SRS formula (pretending \(n_d\) was fixed), our standard errors would be artificially low:
Code
## naively applying SRS estimation to each domain (incorrect analysis)
region_agsrs_srs <- data.frame(rbind(
NC = srs_est(agsrs$acres92[agsrs$region=="NC"], N = 1054, estimate = "mean", show.details = FALSE)$estimate,
NE = srs_est(agsrs$acres92[agsrs$region=="NE"], N = 220, estimate = "mean", show.details = FALSE)$estimate,
S = srs_est(agsrs$acres92[agsrs$region=="S"], N = 1382, estimate = "mean", show.details = FALSE)$estimate,
W = srs_est(agsrs$acres92[agsrs$region=="W"], N = 422, estimate = "mean", show.details = FALSE)$estimate
))
region_agsrs_srs |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper | |
|---|---|---|---|---|
| NC | 323.42 | 30.40 | 262.89 | 383.94 |
| NE | 106.55 | 29.21 | 44.93 | 168.17 |
| S | 246.19 | 23.26 | 200.24 | 292.13 |
| W | 518.98 | 70.03 | 377.74 | 660.21 |
Output Note: Comparing the two tables, you can see the correct domain_mean Standard Errors are larger, appropriately penalizing us for the uncertainty in how many counties from each region randomly fell into our sample.
4.5.2 Post-stratification Analysis
If we know the true regional population sizes (\(N_h\)), we can post-stratify our SRS data.
Code
nh <- tapply(agsrs[, "acres92"], agsrs[,"region"], length)
n <- sum(nh)
sh <- tapply(agsrs[, "acres92"], agsrs[,"region"], sd)
ybarh <- tapply(agsrs[, "acres92"], agsrs[,"region"], mean)
# create a vector with external information
Nh <- c(NC = 1054, NE = 220, S = 1382, W = 422)Create \(n_h^{\text{post-strat}} = n \times \frac{N_h}{N}\), instead of using the randomly observed \(n_h\).
Code
nh_post <- Nh / sum(Nh) * n
## regional summary for stratified sampling estimation
data.frame(
ybarh = ybarh,
sh = sh,
nh = nh,
pi_obs = nh/n,
Nh = Nh,
pi_pop = Nh/sum(Nh),
nh_post = nh_post
) |>
gt(rownames_to_stub = TRUE) |>
cols_label(
ybarh = md("$\\bar{y}_h$"),
sh = md("$s_h$"),
nh = md("$n_h^{\\text{obs}}$"),
pi_obs = md("$\\pi^{\\text{obs}}_h$"),
Nh = md("$N_h$"),
pi_pop = md("$\\pi^{\\text{pop}}_h$"),
nh_post = md("$n_h^{\\text{post}}$")
) |>
fmt_number(columns = c(ybarh, sh, pi_obs, pi_pop, nh_post), decimals = 3)| \(\bar{y}_h\) | \(s_h\) | \(n_h^{\text{obs}}\) | \(\pi^{\text{obs}}_h\) | \(N_h\) | \(\pi^{\text{pop}}_h\) | \(n_h^{\text{post}}\) | |
|---|---|---|---|---|---|---|---|
| NC | 323.416 | 278.970 | 78 | 0.260 | 1054 | 0.342 | 102.729 |
| NE | 106.549 | 129.320 | 18 | 0.060 | 220 | 0.071 | 21.442 |
| S | 246.186 | 312.941 | 160 | 0.533 | 1382 | 0.449 | 134.698 |
| W | 518.978 | 490.842 | 44 | 0.147 | 422 | 0.137 | 41.131 |
Output Note: The table above shows the difference between the sample proportion (\(\pi^{obs}_h\)) and the true population proportion (\(\pi^{pop}_h\)). Post-stratification forces the weights to reflect the true population structure.
Code
## strata mean estimate
res_post_strat <- str_est(ybarh, sh, nh_post, Nh, estimate = "mean")
res_post_strat$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 | 323.42 | 77,824.28 | 110.75 | 80.17 |
| NE | 220 | 21 | 0.07 | 106.55 | 16,723.71 | 7.62 | 3.60 |
| S | 1,382 | 135 | 0.45 | 246.19 | 97,932.26 | 110.54 | 132.28 |
| W | 422 | 41 | 0.14 | 518.98 | 240,925.91 | 71.15 | 99.37 |
| Total | 3,078 | 300 | 1.00 | 300.05 | 315.43 | ||
| 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 = 300.05\) | |||||||
| \(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{315.43} = 17.76\) | |||||||
Code
post_strat_est <- res_post_strat$estimate
post_strat_est |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 300.05 | 17.76 | 265.24 | 334.86 |
Code
## compare with SRS estimate
res_srs_est <- srs_est(agsrs$acres92, N = 3087, estimate = "mean", show.details = FALSE)
srs_est_raw <- res_srs_est$estimate
srs_est_raw |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 297.90 | 18.90 | 260.70 | 335.09 |
Code
## true mean check
mean(agpop$acres92)[1] 309.9004Output Note: The post-stratified estimate gives a Standard Error of 17.8, which is more precise than the unadjusted SRS standard error of 18.9.
4.6 Estimating Domain Means and Post-stratification with teacher example
4.6.1 Importing Data and Calculating Contingency Table
Code
teacher <- read.csv("data/college_teacher.csv")
table(teacher[, 2:1]) teacher
gender 0 1
1 156 84
2 120 40Code
addmargins(table(teacher[, 2:1]), 2) teacher
gender 0 1 Sum
1 156 84 240
2 120 40 160Code
# Proportions of students who want to become a teacher in two domains (female and male)
prop.table(table(teacher[, 2:1]), margin = "gender") teacher
gender 0 1
1 0.65 0.35
2 0.75 0.254.6.2 Domain Summary
Code
n <- nrow(teacher)
ybarh <- tapply(teacher$teacher, INDEX = teacher$gender, FUN = mean)
sh <- tapply(teacher$teacher, INDEX = teacher$gender, FUN = sd)
nh <- tapply(teacher$teacher, INDEX = teacher$gender, FUN = length)
prop_h_obs <- nh / sum(nh)
data.frame(
ybarh = ybarh,
sh = sh,
nh = nh,
pi_obs = prop_h_obs
) |>
gt(rownames_to_stub = TRUE) |>
cols_label(
ybarh = md("$\\bar{y}_h$"),
sh = md("$s_h$"),
nh = md("$n^{\\text{obs}}_h$"),
pi_obs = md("$\\pi^{\\text{obs}}_h$")
) |>
fmt_number(columns = everything(), decimals = 3)| \(\bar{y}_h\) | \(s_h\) | \(n^{\text{obs}}_h\) | \(\pi^{\text{obs}}_h\) | |
|---|---|---|---|---|
| 1 | 0.350 | 0.478 | 240.000 | 0.600 |
| 2 | 0.250 | 0.434 | 160.000 | 0.400 |
4.6.3 Estimating Domain Means
Code
Nh <- c(3000, 1000)
n <- nrow(teacher)
domain_teacher_est <- data.frame(rbind(
female = domain_est(subset(teacher, gender==1)$teacher, n=n, N=4000),
male = domain_est(subset(teacher, gender==2)$teacher, n=n, N=4000)
))
domain_teacher_est |> format_est_gt("domain_mean")| \(\bar{y}_{d}\) | \(\mathrm{SE}(\bar{y}_{d})\) | 95% CI Lower | 95% CI Upper | |
|---|---|---|---|---|
| female | 0.35 | 0.03 | 0.29 | 0.41 |
| male | 0.25 | 0.03 | 0.19 | 0.31 |
4.6.4 Post-stratification Analysis
Code
Nh <- c(3000, 1000)
prop_h_pop <- Nh / sum(Nh)
prop_h_pop[1] 0.75 0.25Create \(n_h^{\text{post-strat}} = n \times \frac{N_h}{N}\)
Code
nh_post <- Nh / sum(Nh) * n
## show the differences
data.frame(
phat = ybarh,
sh = sh,
nh = nh,
pi_obs = nh/n,
Nh = Nh,
pi_pop = Nh/sum(Nh),
nh_post = nh_post
) |>
gt(rownames_to_stub = TRUE) |>
cols_label(
phat = md("$\\hat{p}_h$"),
sh = md("$s_h$"),
nh = md("$n_h$"),
pi_obs = md("$\\pi^{\\text{obs}}_h$"),
Nh = md("$N_h$"),
pi_pop = md("$\\pi^{\\text{pop}}_h$"),
nh_post = md("$n_h^{\\text{post-strat}}$")
) |>
fmt_number(columns = c(phat, sh, pi_obs, pi_pop, nh_post), decimals = 3)| \(\hat{p}_h\) | \(s_h\) | \(n_h\) | \(\pi^{\text{obs}}_h\) | \(N_h\) | \(\pi^{\text{pop}}_h\) | \(n_h^{\text{post-strat}}\) | |
|---|---|---|---|---|---|---|---|
| 1 | 0.350 | 0.478 | 240 | 0.600 | 3000 | 0.750 | 300.000 |
| 2 | 0.250 | 0.434 | 160 | 0.400 | 1000 | 0.250 | 100.000 |
Calculate the post-stratification estimation of the mean:
Code
## poststratification estimation of mean
res_teacher_post_strat <- str_est(ybarh, sh, nh_post, Nh, estimate = "mean")
res_teacher_post_strat$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 | 3,000 | 300 | 0.75 | 0.35 | 0.23 | 0.26 | 3.86 × 10−4 |
| 2 | 1,000 | 100 | 0.25 | 0.25 | 0.19 | 0.06 | 1.06 × 10−4 |
| Total | 4,000 | 400 | 1.00 | 0.32 | 4.92 × 10−4 | ||
| 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 = 0.32\) | |||||||
| \(\mathrm{SE}(\bar y_{str}) = \sqrt{\sum_h v_h} = \sqrt{4.92 \times 10^{-4}} = 0.02\) | |||||||
Code
teacher_post_strat <- res_teacher_post_strat$estimate
teacher_post_strat |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 0.32 | 0.02 | 0.28 | 0.37 |
This post-stratification estimate explicitly corrects for the sampling imbalance using the true strata sizes:
Code
# 0.35 * 0.75 + 0.25 * 0.25
(ybarh[1] * prop_h_pop[1]) + (ybarh[2] * prop_h_pop[2]) 1
0.325 4.6.5 Comparing to the Analysis with SRS
If we didn’t post-stratify, the simple mean is heavily weighted toward males because they were over-sampled by chance (40% in sample vs 25% in population).
Code
res_teacher_srs <- srs_est(teacher$teacher, N=sum(Nh), estimate = "mean", show.details = FALSE)
res_teacher_srs$tableNULLCode
teacher_srs <- res_teacher_srs$estimate
teacher_srs |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 0.31 | 0.02 | 0.27 | 0.35 |
This equals the naive sample proportion:
Code
# 0.35 * 0.6 + 0.25 * 0.4
mean(teacher$teacher == 1)[1] 0.314.7 Regression Estimation for the cherry.csv dataset
When the relationship between \(y\) and \(x\) doesn’t pass exactly through the origin, we use Regression Estimation.
4.7.1 Importing data
Code
cherry <- read.csv("data/cherry.csv", header = T)
ydata <- cherry$volume
xdata <- cherry$diameter
t_diameters <- 41835
xbarU <- t_diameters / 2967
N <- 29674.7.2 Regression Estimation Using the Function
4.7.2.1 Fitting a linear regression model
Code
n <- length(ydata)
lmfit <- lm(ydata ~ xdata)
lmfit_gt <- format_lm_gt(lmfit)
lmfit_gt$coef| Coefficient Estimates | ||||
| Estimate | S.E. | t value | p value | |
|---|---|---|---|---|
| (Intercept) | −36.9435 | 3.3651 | −10.9783 | 0.0000 |
| xdata | 5.0659 | 0.2474 | 20.4783 | 0.0000 |
Code
lmfit_gt$fit| Model Fit Summary | ||||
| \(R^2\) | Adj. \(R^2\) | \(\hat\sigma\) | F-statistic | p-value |
|---|---|---|---|---|
| 0.9353 | 0.9331 | 4.2520 | 419.3603 | 0.0000 |
Code
plot(xdata, ydata, main="Linear Regression: Volume vs Diameter")
abline(lmfit, col="green", lwd=2)
The function’s $table shows the row-level working values \(x_i\), \(y_i\), fitted \(\hat{y}_i = \hat{B}_0 + \hat{B}_1 x_i\), residuals \(e_i = y_i - \hat{y}_i\), and \(e_i^2\) (with a Sum row and a footnote giving \(\hat B_0\), \(\hat B_1\), and \(s_e^2\)):
Code
res_reg_model <- reg_est(ydata, xdata, N = N, estimate = "model")
res_reg_model$table| Regression Estimation: Working Table | |||||
| \(n = 31\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 10.30 | 8.30 | 5.10 | 5.20 | 27.01 |
| 2 | 10.30 | 8.60 | 6.62 | 3.68 | 13.52 |
| 3 | 10.20 | 8.80 | 7.64 | 2.56 | 6.57 |
| 4 | 16.40 | 10.50 | 16.25 | 0.15 | 0.02 |
| 5 | 18.80 | 10.70 | 17.26 | 1.54 | 2.37 |
| 6 | 19.70 | 10.80 | 17.77 | 1.93 | 3.73 |
| ... | ... | ... | ... | ... | ... |
| Sum | 935.30 | 410.70 | 935.30 | 0.00 | 524.30 |
| 1 \(\hat B_0 = -36.94\), \(\hat B_1 = 5.07\), \(\hat y_i = \hat B_0 + \hat B_1 x_i\), \(s_e^2 = 18.08\). | |||||
| \(\hat B_1 = \dfrac{\sum_i(x_i-\bar x)(y_i-\bar y)}{\sum_i(x_i-\bar x)^2} = \dfrac{1,496.64}{295.44} = 5.07\) | |||||
| \(\mathrm{SE}(\hat B_1) = \sqrt{\dfrac{s_e^2}{\sum_i(x_i-\bar x)^2}} = \sqrt{\dfrac{18.08}{295.44}} = 0.25\) | |||||
4.7.2.2 Estimate the mean
We compute the regression estimate \(\overline{y}_{reg}\) by passing xbarU with estimate = "mean".
Code
res_reg_mean <- reg_est(ydata, xdata, xbarU = xbarU, N = N, estimate = "mean")
res_reg_mean$table| Regression Estimation: Working Table | |||||
| \(n = 31\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 10.30 | 8.30 | 5.10 | 5.20 | 27.01 |
| 2 | 10.30 | 8.60 | 6.62 | 3.68 | 13.52 |
| 3 | 10.20 | 8.80 | 7.64 | 2.56 | 6.57 |
| 4 | 16.40 | 10.50 | 16.25 | 0.15 | 0.02 |
| 5 | 18.80 | 10.70 | 17.26 | 1.54 | 2.37 |
| 6 | 19.70 | 10.80 | 17.77 | 1.93 | 3.73 |
| ... | ... | ... | ... | ... | ... |
| Sum | 935.30 | 410.70 | 935.30 | 0.00 | 524.30 |
| 1 \(\hat B_0 = -36.94\), \(\hat B_1 = 5.07\), \(\hat y_i = \hat B_0 + \hat B_1 x_i\), \(s_e^2 = 18.08\). | |||||
| \(\bar y_{reg} = \hat B_0+\hat B_1\bar x_U = -36.94 + 5.07 \times 14.10 = 34.49\) | |||||
| \(\mathrm{SE}(\bar y_{reg}) = \sqrt{\left(1-\dfrac{n}{N}\right)\dfrac{s_e^2}{n}} = \sqrt{\left(1-\dfrac{31}{2967}\right)\dfrac{18.08}{31}} = 0.76\) | |||||
Code
output_reg <- res_reg_mean$estimate
output_reg |> format_est_gt("reg_mean")| \(\bar{y}_{reg}\) | \(\mathrm{SE}(\bar{y}_{reg})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 34.49 | 0.76 | 32.93 | 36.04 |
Output Note: The Regression estimate for the mean volume is 34.486 with a standard error of 0.7597.
4.7.2.3 Estimate the total
Code
reg_total_volume_raw <- reg_est(ydata, xdata, xbarU = xbarU, N = N, estimate = "total")$estimate
reg_total_volume_raw |> format_est_gt("reg_total")| \(\hat{t}_{reg}\) | \(\mathrm{SE}(\hat{t}_{reg})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 102,318.86 | 2,253.97 | 97,708.98 | 106,928.74 |
4.7.2.4 Percentage of Variance Reduction
Code
reg_mean_volume_raw <- output_reg
var_reduction_reg <- 1 - (reg_mean_volume_raw["S.E."] / srs_mean_vol_raw["S.E."])^2
unname(var_reduction_reg)[1] 0.9330895Output Note: The variance reduction for the regression estimator is 93.3%, which is similar to the ratio estimator since the true intercept is very close to zero.
4.8 Regression estimation for photo counts of dead trees
To estimate the number of dead trees in an area, we divide the area into 100 square plots and count dead trees on photographs. We select an SRS of 25 plots for exact field counts. We know the true mean photo count is 11.3.
Code
photocounts <- c(10, 12, 7, 13, 13, 6, 17, 16, 15, 10, 14, 12, 10, 5, 12, 10,
10, 9, 6, 11, 7, 9, 11, 10, 10)
fieldcounts <- c(15, 14, 9, 14, 8, 5, 18, 15, 13, 15, 11, 15, 12, 8, 13,
9, 11, 12, 9, 12, 13, 11, 10, 9, 8)
lmfit <- lm(fieldcounts ~ photocounts)
lmfit_gt <- format_lm_gt(lmfit)
lmfit_gt$coef| Coefficient Estimates | ||||
| Estimate | S.E. | t value | p value | |
|---|---|---|---|---|
| (Intercept) | 5.0593 | 1.7635 | 2.8689 | 0.0087 |
| photocounts | 0.6133 | 0.1601 | 3.8316 | 0.0009 |
Code
lmfit_gt$fit| Model Fit Summary | ||||
| \(R^2\) | Adj. \(R^2\) | \(\hat\sigma\) | F-statistic | p-value |
|---|---|---|---|---|
| 0.3896 | 0.3631 | 2.4062 | 14.6815 | 0.0009 |
Code
plot(photocounts, fieldcounts, main="Field Counts vs Photo Counts")
abline(lmfit, col="purple", lwd=2)
Let’s compare the three different estimation methods to see which provides the lowest Standard Error for the population mean.
Code
# 1. Regression estimate for the mean
res_dt_reg <- reg_est(ydata = fieldcounts, xdata = photocounts, xbarU = 11.3, N = 100, estimate = "mean")
res_dt_reg$table| Regression Estimation: Working Table | |||||
| \(n = 25\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 15.00 | 10.00 | 11.19 | 3.81 | 14.50 |
| 2 | 14.00 | 12.00 | 12.42 | 1.58 | 2.50 |
| 3 | 9.00 | 7.00 | 9.35 | −0.35 | 0.12 |
| 4 | 14.00 | 13.00 | 13.03 | 0.97 | 0.94 |
| 5 | 8.00 | 13.00 | 13.03 | −5.03 | 25.32 |
| 6 | 5.00 | 6.00 | 8.74 | −3.74 | 13.98 |
| ... | ... | ... | ... | ... | ... |
| Sum | 289.00 | 265.00 | 289.00 | 0.00 | 133.16 |
| 1 \(\hat B_0 = 5.06\), \(\hat B_1 = 0.61\), \(\hat y_i = \hat B_0 + \hat B_1 x_i\), \(s_e^2 = 5.79\). | |||||
| \(\bar y_{reg} = \hat B_0+\hat B_1\bar x_U = 5.06 + 0.61 \times 11.30 = 11.99\) | |||||
| \(\mathrm{SE}(\bar y_{reg}) = \sqrt{\left(1-\dfrac{n}{N}\right)\dfrac{s_e^2}{n}} = \sqrt{\left(1-\dfrac{25}{100}\right)\dfrac{5.79}{25}} = 0.42\) | |||||
Code
dt_reg_raw <- res_dt_reg$estimate
dt_reg_raw |> format_est_gt("reg_mean")| \(\bar{y}_{reg}\) | \(\mathrm{SE}(\bar{y}_{reg})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 11.99 | 0.42 | 11.13 | 12.85 |
Code
# (And estimating the total)
reg_est(ydata = fieldcounts, xdata = photocounts, xbarU = 11.3, N = 100, estimate = "total")$estimate |>
format_est_gt("reg_total")| \(\hat{t}_{reg}\) | \(\mathrm{SE}(\hat{t}_{reg})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 1,198.93 | 41.68 | 1,112.72 | 1,285.14 |
Code
# 2. Simple random sampling estimate for the mean (Ignoring photo counts)
res_dt_srs <- srs_est(sdata = fieldcounts, N = 100, estimate = "mean", show.details = FALSE)
dt_srs_raw <- res_dt_srs$estimate
dt_srs_raw |> format_est_gt("mean")| \(\bar{y}\) | \(\mathrm{SE}(\bar{y})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 11.56 | 0.52 | 10.48 | 12.64 |
Code
# 3. Ratio estimate for the mean
res_dt_ratio <- ratio_est(ydata = fieldcounts, xdata = photocounts, xbarU = 11.3, N = 100, estimate = "mean")
res_dt_ratio$table| Ratio Estimation: Working Table | |||||
| \(n = 25\) | |||||
| \(i\) | \(y_i\) | \(x_i\) | \(\hat y_i\)1 | \(e_i\) | \(e_i^2\) |
|---|---|---|---|---|---|
| 1 | 15.00 | 10.00 | 10.91 | 4.09 | 16.76 |
| 2 | 14.00 | 12.00 | 13.09 | 0.91 | 0.83 |
| 3 | 9.00 | 7.00 | 7.63 | 1.37 | 1.87 |
| 4 | 14.00 | 13.00 | 14.18 | −0.18 | 0.03 |
| 5 | 8.00 | 13.00 | 14.18 | −6.18 | 38.16 |
| 6 | 5.00 | 6.00 | 6.54 | −1.54 | 2.38 |
| ... | ... | ... | ... | ... | ... |
| Sum | 289.00 | 265.00 | 289.00 | 0.00 | 184.64 |
| 1 \(\hat B = \sum_i y_i / \sum_i x_i = 1.09\), \(\hat y_i = \hat B\, x_i\), \(s_e^2 = 7.69\). | |||||
| \(\bar y_r = \hat B\,\bar x_U = 1.09 \times 11.30 = 12.32\) | |||||
| \(\mathrm{SE}(\bar y_r) = \mathrm{SE}(\hat B)\,\bar x_U = 0.05 \times 11.30 = 0.51\) | |||||
Code
dt_ratio_raw <- res_dt_ratio$estimate
dt_ratio_raw |> format_est_gt("ratio_mean")| \(\bar{y}_{r}\) | \(\mathrm{SE}(\bar{y}_{r})\) | 95% CI Lower | 95% CI Upper |
|---|---|---|---|
| 12.32 | 0.51 | 11.27 | 13.38 |
Output Note:
SRS Mean: SE = 0.522
Ratio Mean: SE = 0.512
Regression Mean: SE = 0.417
Because the intercept of the true line is non-zero, the Regression estimator provides the most precise estimate in this scenario.
4.9 Shinylive App Illustrating Ratio Estimation and Regression Estimation
The simulation for this section runs as a standalone app: Ratio Estimation and Regression Estimation. It opens in a new tab, together with its documentation and the full list of apps.
4.10 Interactive Demonstration: Post-stratifying a Simple Random Sample
The simulation for this section runs as a standalone app: Post-stratification of a Simple Random Sample. It opens in a new tab, together with its documentation and the full list of apps.