Shinylive App Comparing SRS (average and ratio) with UPSWR (Hansen–Hurwitz)
SRS (average and ratio) with UPSWR (Hansen–Hurwitz)
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| label: upswr-app
#| standalone: true
#| viewerHeight: 860
library(shiny)
## ---------------------------------------------------------------- population
data_url <- "https://raw.githubusercontent.com/longhaiSK/sampling/main/data/agpop.csv"
raw <- tryCatch({
download.file(data_url, "agpop.csv")
read.csv("agpop.csv")
}, error = function(e) NULL)
if (is.null(raw)) { # fallback so the app runs without network access
set.seed(1)
Nf <- 3000
x1 <- rgamma(Nf, shape = 1.2, scale = 1.8e5)
raw <- data.frame(
acres92 = pmax(0, 0.93 * x1 + rnorm(Nf, 0, 0.5 * x1^0.85)),
acres87 = x1,
farms92 = pmax(1, round(300 + x1 / 6000 + rnorm(Nf, 0, 250)))
)
}
raw <- raw[is.finite(raw$acres92) & raw$acres92 != -99, ]
yv <- raw$acres92
N <- length(yv)
## every numeric column other than y is a candidate auxiliary variable, both
## for the UPS sizes M and for the ratio estimator's x;
## -99 marks a missing entry and is replaced by the mean of the rest
num <- vapply(raw, is.numeric, logical(1))
X <- raw[, num & names(raw) != "acres92", drop = FALSE]
X <- as.data.frame(lapply(X, function(v) {
v[!is.finite(v) | v == -99] <- NA
v[is.na(v)] <- mean(v, na.rm = TRUE)
v
}))
keep <- vapply(X, function(v) all(is.finite(v)) && all(v >= 0) && var(v) > 0,
logical(1))
X <- X[, keep, drop = FALSE]
xvars <- names(X)
xdef <- if ("acres87" %in% xvars) "acres87" else xvars[1]
mdef <- xdef
ybarU <- mean(yv)
Vy <- var(yv) # S^2 of y, for the SRS average
ZCRIT <- 1.96 # nominal 95% intervals
set.seed(11) # display-only vertical jitter, fixed once
ypl <- yv + rnorm(N, 0, 0.004 * diff(range(yv)))
rgy <- range(ypl)
## ------------------------------------------------------------ estimator specs
EST <- c("srs", "ratio", "hh")
LAB <- c(srs = "SRS average", ratio = "SRS ratio", hh = "UPSWR HH")
COL <- c(srs = "#e07b39", ratio = "#2ca02c", hh = "#9467bd")
PCH <- c(srs = 17, ratio = 16, hh = 15)
col_samp <- "#d62728" # missed intervals and the population mean
col_pts <- "blue" # sampled observations
col_true <- "black"
col_pop <- "#9aa0a655"
col_band <- "#eef3f8"
## point size grows with the per-draw selection probability; the scale is
## shared by both designs, so SRS points (all 1/N) show what "equal" looks like
psize <- function(p, pmax) 0.25 + 2.4 * sqrt(p / pmax)
## ------------------------------------------------------- broken-axis mapping
## "top": linear from the minimum up to a cut, then the sparse upper tail
## squeezed into a thin strip (SRS view: every unit is equally likely).
## "low": the many small values squeezed into a thin strip at the bottom and
## the large values spread out (UPSWR view: large units dominate).
BRK_THIN <- 0.14
BRK_GAP <- 0.025
BRK_WIDE <- 1 - BRK_THIN - BRK_GAP
## SRS view: the population mean falls at mid-height
cut_at_centre <- function(rmin, rmax, m) {
cut <- rmin + (m - rmin) * BRK_WIDE / 0.5
if (!is.finite(cut) || cut >= rmax) rmin + 0.95 * (rmax - rmin) else cut
}
## UPSWR view: squeeze the values that PPS reaches with only `share` of its
## draws, i.e. the lower `share` of the psi-weighted distribution
cut_psi_share <- function(v, psi, share = 0.2) {
o <- order(v)
v[o][which(cumsum(psi[o]) >= share)[1]]
}
mk_break <- function(rmin, rmax, cut, mode = c("top", "low")) {
mode <- match.arg(mode)
eps <- 1e-9 * (rmax - rmin)
cut <- min(max(cut, rmin + eps), rmax - eps)
lo_h <- if (mode == "top") BRK_WIDE else BRK_THIN
hi_h <- 1 - lo_h - BRK_GAP
tf <- function(v) {
v <- pmin(pmax(v, rmin), rmax)
ifelse(v <= cut, (v - rmin) / (cut - rmin) * lo_h,
lo_h + BRK_GAP + (v - cut) / (rmax - cut) * hi_h)
}
list(tf = tf, rmin = rmin, rmax = rmax, cut = cut, mode = mode,
gap = c(lo_h, lo_h + BRK_GAP))
}
mk_plain <- function(rmin, rmax) {
list(tf = function(v) (pmin(pmax(v, rmin), rmax) - rmin) / (rmax - rmin),
rmin = rmin, rmax = rmax, cut = NA, gap = NULL)
}
## counts stay as counts, acreages are shown in thousands
unit_div <- function(rmax) if (rmax >= 1e5) 1000 else 1
unit_lab <- function(div) if (div == 1000) " (thousands)" else ""
axis_at <- function(bk, side, div, rg, nt) {
v <- pretty(rg, nt)
v <- v[v >= rg[1] & v <= rg[2]]
if (!length(v)) return(invisible())
d <- if (diff(rg) / div < 40) 1 else 0
axis(side, at = bk$tf(v), labels = formatC(v / div, format = "f", digits = d,
big.mark = ","), las = 1, cex.axis = 0.75)
}
draw_axis <- function(bk, side, div) {
if (is.null(bk$gap)) {
axis_at(bk, side, div, c(bk$rmin, bk$rmax), 5)
} else {
nt <- if (bk$mode == "top") c(5, 2) else c(2, 5)
axis_at(bk, side, div, c(bk$rmin, bk$cut), nt[1])
axis_at(bk, side, div, c(bk$cut + 1e-9, bk$rmax), nt[2])
}
}
mask_gap <- function(bk, side) {
if (is.null(bk$gap)) return(invisible())
usr <- par("usr")
g <- bk$gap
m <- mean(g)
if (side == 2) {
rect(usr[1], g[1], usr[2], g[2], col = "white", border = NA)
d <- 0.010 * (usr[2] - usr[1])
segments(usr[1] - d, c(m - 0.012, m + 0.004), usr[1] + d,
c(m - 0.004, m + 0.012), xpd = NA, lwd = 1.2)
} else {
rect(g[1], usr[3], g[2], usr[4], col = "white", border = NA)
d <- 0.010 * (usr[4] - usr[3])
segments(c(m - 0.012, m + 0.004), usr[3] - d, c(m - 0.004, m + 0.012),
usr[3] + d, xpd = NA, lwd = 1.2)
}
}
## interval drawn as a capped vertical bar with the estimate marked on it;
## the bar turns red when the interval misses the population mean
ci_bar <- function(x, est, ci, col, pch, covered, tf = identity, cap = 0.012) {
bcol <- if (covered) col else col_samp
segments(x, tf(ci[1]), x, tf(ci[2]), col = bcol, lwd = 2)
segments(x - cap, tf(ci), x + cap, tf(ci), col = bcol, lwd = 2)
points(x, tf(est), pch = pch, col = col, cex = 1.7)
}
## symmetric density polygon, drawn vertically about x = at
violin <- function(x, at, col, w = 0.32) {
if (length(x) < 10 || diff(range(x)) <= 0) return(invisible())
d <- density(x)
hw <- w * d$y / max(d$y)
polygon(c(at + hw, rev(at - hw)), c(d$x, rev(d$x)),
col = adjustcolor(col, 0.35), border = col)
m <- median(x)
segments(at - w * 0.6, m, at + w * 0.6, m, col = col, lwd = 2)
}
pct <- function(v) if (is.na(v)) "--" else sprintf("%.1f%%", 100 * v)
## ---- inline SVG button icons (no font or icon-library dependency) ----------
svg_icon <- function(paths)
HTML(sprintf(paste0('<svg width="16" height="16" viewBox="0 0 16 16" ',
'fill="currentColor" style="vertical-align:-2px">%s</svg>'), paths))
icon_play <- svg_icon('<path d="M3 2v12l10-6z"/>')
icon_pause <- svg_icon('<rect x="3" y="2" width="4" height="12"/><rect x="9" y="2" width="4" height="12"/>')
icon_skip <- svg_icon('<path d="M1 2v12l7-6zM8 2v12l6-6z"/><rect x="14" y="2" width="2" height="12"/>')
## ---------------------------------------------------------------------- ui
ui <- fluidPage(
titlePanel("Shinylive App for UPSWR"),
tags$style(HTML(
".well {padding: 10px 12px; margin-bottom: 10px;}
.form-group {margin-bottom: 6px;}
.checkbox, .radio {margin-top: 2px; margin-bottom: 2px;}
.irs {margin-bottom: 0;}
.btn {margin-bottom: 6px;}"
)),
wellPanel(
fluidRow(
column(
3,
selectInput("mvar", HTML("Aux. var. for UPS (M<sub>i</sub>)"),
choices = xvars, selected = mdef),
selectInput("xvar", HTML("Aux. var. for estimation (x<sub>i</sub>)"),
choices = xvars, selected = xdef),
checkboxInput("sync", HTML("Same variable (M<sub>i</sub> = x<sub>i</sub>)"), TRUE),
checkboxInput("brk", "Break the axes", TRUE)
),
column(
3,
checkboxGroupInput("show", "Estimators",
choices = c("SRS + average" = "srs",
"SRS + ratio" = "ratio",
"UPSWR + HH ratio" = "hh"),
selected = EST),
sliderInput("n", "Sample size (n)", 10, 500, 100, step = 10)
),
column(
3,
radioButtons("design", "Show sample from",
choices = c("SRS" = "srs", "UPSWR" = "ups"),
selected = "ups", inline = TRUE),
sliderInput("speed", "Samples per second", min = 1, max = 10, value = 5, step = 1),
numericInput("nrep", "Stop after this many samples",
value = 500, min = 10, max = 20000, step = 50)
),
column(
3,
actionButton("draw", "Draw new samples", class = "btn-primary", width = "100%"),
br(), br(),
fluidRow(
column(4, actionButton("play", icon_play, class = "btn-success", width = "100%",
title = "Play")),
column(4, actionButton("pause", icon_pause, class = "btn-warning", width = "100%",
title = "Pause")),
column(4, actionButton("skip", icon_skip, class = "btn-info", width = "100%",
title = "Skip to the end"))
),
actionButton("reset", "Reset", width = "100%"),
div(style = "font-size: 90%;", textOutput("counter"))
)
)
),
plotOutput("mainplot", height = "520px")
)
## ------------------------------------------------------------------- server
server <- function(input, output, session) {
rv <- reactiveValues(s_srs = NULL, s_ups = NULL, n = 100, ests = NULL)
running <- reactiveVal(FALSE)
nrep <- reactive({
v <- input$nrep
if (is.null(v) || is.na(v) || v < 1) 1L else min(as.integer(v), 20000L)
})
## everything that depends on the two auxiliary variables: M sets the UPS
## selection probabilities, x feeds the SRS ratio estimator
aux <- reactive({
mv <- X[[input$mvar]]
xv <- X[[input$xvar]]
psi <- mv / sum(mv) # per-draw selection probabilities
B <- sum(yv) / sum(xv)
pos <- psi > 0
list(m = mv, x = xv, psi = psi, B = B, xbar = mean(xv),
## per-draw variances, so the SE of each estimator with n draws is
## sqrt(fpc * V / n) for the SRS ones and sqrt(V / n) for HH
V = c(srs = Vy, ratio = var(yv - B * xv),
hh = sum(psi[pos] * (yv[pos] / psi[pos] - sum(yv))^2) / N^2),
## horizontal axis of each view: x for SRS, M for UPSWR
ax = list(
top = list(v = xv, rg = range(xv), slope = B, name = input$xvar,
cut = cut_at_centre(min(xv), max(xv), mean(xv)),
cuty = cut_at_centre(rgy[1], rgy[2], ybarU)),
low = list(v = mv, rg = range(mv), slope = sum(yv) / sum(mv),
name = input$mvar, cut = cut_psi_share(mv, psi),
cuty = cut_psi_share(ypl, psi))))
})
se_theory <- function(a, n)
sqrt(c(srs = 1 - n / N, ratio = 1 - n / N, hh = 1) * a$V / n)
## one replication: an SRS (without replacement) feeding the average and the
## ratio estimator, and an independent PPS-with-replacement sample feeding HH
one_draw <- function(a, n) {
s <- sample.int(N, n)
u <- sample.int(N, n, replace = TRUE, prob = a$psi)
ys <- yv[s]
xs <- a$x[s]
fpc <- 1 - n / N
Bhat <- if (sum(xs) > 0) sum(ys) / sum(xs) else 0
z <- yv[u] / a$psi[u] / N # each draw's estimate of the mean
est <- c(srs = mean(ys),
ratio = Bhat * a$xbar,
hh = mean(z))
se <- c(srs = sqrt(fpc * var(ys) / n),
ratio = (a$xbar / mean(xs)) * sqrt(fpc * var(ys - Bhat * xs) / n),
hh = sd(z) / sqrt(n))
e <- c(setNames(est, paste0("est_", EST)),
setNames(est - ZCRIT * se, paste0("lo_", EST)),
setNames(est + ZCRIT * se, paste0("hi_", EST)),
setNames(as.numeric(abs(est - ybarU) <= ZCRIT * se), paste0("cov_", EST)),
Bhat = Bhat, Bhh = mean(yv[u] / a$m[u]))
list(s = s, u = u, e = e)
}
## draw k replications; only the last one's samples are kept for the left
## panel. "Skip" calls this with the full remaining count.
draw_many <- function(k) {
a <- isolate(aux())
n <- isolate(input$n)
E <- vector("list", k)
for (i in seq_len(k)) {
d <- one_draw(a, n)
E[[i]] <- d$e
}
rv$s_srs <- d$s
rv$s_ups <- d$u
rv$n <- n
rv$ests <- rbind(rv$ests, do.call(rbind, E))
}
draw_one <- function() draw_many(1)
## a fresh history always starts from one sample, so the plot is never empty
reset_history <- function() {
running(FALSE)
rv$ests <- NULL
draw_one()
}
observeEvent(list(input$n, input$mvar, input$xvar), reset_history(), ignoreInit = FALSE)
## keep M and x the same while "Same variable" is ticked: a change to either
## box is copied to the other, and ticking the box copies M to x
sync_to <- function(id, v) {
if (isTRUE(input$sync) && !identical(input[[id]], v))
updateSelectInput(session, id, selected = v)
}
observeEvent(input$mvar, sync_to("xvar", input$mvar))
observeEvent(input$xvar, sync_to("mvar", input$xvar))
observeEvent(input$sync, sync_to("xvar", input$mvar))
observeEvent(input$reset, reset_history())
## one sample per click; Play repeats until Pause or until the cap is reached
observeEvent(input$draw, { running(FALSE); draw_one() })
observeEvent(input$play, running(TRUE))
observeEvent(input$pause, running(FALSE))
observeEvent(input$skip, {
running(FALSE)
k <- nrep() - NROW(rv$ests)
if (k > 0) draw_many(k)
})
observe({
if (!isTRUE(running())) return()
if (NROW(isolate(rv$ests)) >= isolate(nrep())) { running(FALSE); return() }
invalidateLater(1000 / isolate(input$speed), session)
isolate(draw_one())
})
output$counter <- renderText(
sprintf("%d of %d samples", NROW(rv$ests), nrep())
)
output$mainplot <- renderPlot({
req(rv$s_srs)
sel <- EST[EST %in% input$show]
validate(need(length(sel) > 0, "Select at least one estimator."))
a <- aux()
n <- rv$n
M <- rv$ests
nrep <- NROW(M)
cur <- M[nrep, ]
ups <- identical(input$design, "ups")
s <- if (ups) rv$s_ups else rv$s_srs
xpos <- setNames(0.02 + 0.04 * (seq_along(sel) - 1), sel) # CI bar positions
## violin window, wide enough for the least precise estimator shown
half <- 4 * max(se_theory(a, n)[sel])
lo <- ybarU - half
hi <- ybarU + half
md <- if (ups) "low" else "top"
h <- a$ax[[md]]
if (isTRUE(input$brk)) {
bky <- mk_break(rgy[1], rgy[2], h$cuty, md)
bkx <- mk_break(h$rg[1], h$rg[2], h$cut, md)
} else {
bky <- mk_plain(rgy[1], rgy[2])
bkx <- mk_plain(h$rg[1], h$rg[2])
}
tfy <- bky$tf
tfx <- bkx$tf
divy <- unit_div(rgy[2])
divx <- unit_div(h$rg[2])
xg <- seq(h$rg[1], h$rg[2], length.out = 1500)
draw_line <- function(yl, ...) { # stop the line at the frame
yl[yl < rgy[1] | yl > rgy[2]] <- NA
lines(tfx(xg), tfy(yl), ...)
}
## point sizes: per-draw selection probability under the design shown
pmax_ <- max(a$psi)
cexp <- if (ups) psize(a$psi, pmax_) else rep(psize(1 / N, pmax_), N)
layout(matrix(1:2, nrow = 1), widths = c(1, 1))
## ---- left: population sized by selection probability, with the sample ----
par(mar = c(4.2, 4.8, 3.0, 1.0))
plot(NA, xlim = c(0, 1), ylim = c(0, 1), axes = FALSE,
xlab = paste0(if (ups) "M: " else "x: ", h$name, unit_lab(divx)),
ylab = paste0("t: acres92", unit_lab(divy)),
main = sprintf("Population and %s sample, n = %d",
if (ups) "UPSWR" else "SRS", n))
rect(0, tfy(lo), 1, tfy(hi), col = col_band, border = NA)
points(tfx(h$v), tfy(ypl), pch = 16, cex = cexp, col = col_pop)
draw_line(h$slope * xg, col = col_true, lwd = 1.8)
if (ups) {
if ("hh" %in% sel)
draw_line(cur[["Bhh"]] * xg, col = COL["hh"], lwd = 2, lty = 2)
} else {
if ("ratio" %in% sel)
draw_line(cur[["Bhat"]] * xg, col = COL["ratio"], lwd = 2, lty = 2)
if ("srs" %in% sel)
abline(h = tfy(cur[["est_srs"]]), col = COL["srs"], lwd = 2, lty = 2)
}
## a unit drawn more than once is drawn once, ringed by its draw count
tab <- table(s)
su <- as.integer(names(tab))
points(tfx(h$v[su]), tfy(ypl[su]), pch = 16, cex = cexp[su], col = col_pts)
rep_u <- su[tab > 1]
if (length(rep_u))
text(tfx(h$v[rep_u]), tfy(ypl[rep_u]), labels = tab[tab > 1],
pos = 4, cex = 0.8, col = col_pts, font = 2)
abline(h = tfy(ybarU), col = col_samp, lty = 3)
for (k in sel) ci_bar(xpos[[k]], cur[[paste0("est_", k)]],
c(cur[[paste0("lo_", k)]], cur[[paste0("hi_", k)]]),
COL[[k]], PCH[[k]], cur[[paste0("cov_", k)]] == 1, tfy)
mask_gap(bky, 2); mask_gap(bkx, 1)
box(); draw_axis(bky, 2, divy); draw_axis(bkx, 1, divx)
leg <- data.frame(l = c(if (ups) "size ~ psi = M / M0" else "size ~ 1/N (equal)",
"sampled t"),
c = c("grey55", col_pts), lt = NA, p = 16, w = NA)
leg <- rbind(leg, data.frame(l = if (ups) "true t = (t/M0) M" else "true t = (t/tx) x",
c = col_true, lt = 1, p = NA, w = 1.8))
if (ups && "hh" %in% sel) leg <- rbind(leg, data.frame(
l = "HH slope mean(t/M)", c = COL[["hh"]], lt = 2, p = NA, w = 2))
if (!ups && "ratio" %in% sel) leg <- rbind(leg, data.frame(
l = "ratio fit", c = COL[["ratio"]], lt = 2, p = NA, w = 2))
if (!ups && "srs" %in% sel) leg <- rbind(leg, data.frame(
l = "sample average", c = COL[["srs"]], lt = 2, p = NA, w = 2))
leg <- rbind(leg, data.frame(l = "population mean", c = col_samp, lt = 3, p = NA, w = 1))
legend("topleft", bty = "n", cex = 0.9, inset = c(0.15, 0.01),
legend = leg$l, col = leg$c, lty = leg$lt, pch = leg$p, lwd = leg$w)
## ---- right: sampling distributions as vertical violins, with current CIs ----
par(mar = c(4.2, 4.8, 3.0, 1.0))
k_n <- length(sel)
at <- setNames(seq_len(k_n), sel)
plot(NA, xlim = c(0.4, k_n + 0.6), ylim = c(lo, hi), axes = FALSE,
xlab = "", ylab = paste0("estimate of mean t (acres92)", unit_lab(divy)),
main = "Sampling distributions of the estimates")
usr <- par("usr")
rect(usr[1], usr[3], usr[2], usr[4], col = col_band, border = NA)
abline(h = ybarU, col = col_samp, lty = 3)
for (k in sel) {
violin(M[, paste0("est_", k)], at[[k]], COL[[k]])
ci <- c(cur[[paste0("lo_", k)]], cur[[paste0("hi_", k)]])
bc <- if (cur[[paste0("cov_", k)]] == 1) COL[[k]] else col_samp
segments(at[[k]], ci[1], at[[k]], ci[2], col = bc, lwd = 2.5)
segments(at[[k]] - 0.08, ci, at[[k]] + 0.08, ci, col = bc, lwd = 2.5)
points(at[[k]], cur[[paste0("est_", k)]], pch = PCH[[k]], col = COL[[k]], cex = 1.7)
}
box()
vy <- pretty(c(lo, hi), 5)
vy <- vy[vy >= lo & vy <= hi]
axis(2, at = vy, labels = formatC(vy / divy, format = "f",
digits = if ((hi - lo) / divy < 40) 1 else 0,
big.mark = ","), las = 1, cex.axis = 0.75)
axis(1, at = at, labels = LAB[sel], tick = FALSE, cex.axis = 0.9)
cov_txt <- function(k) if (nrep >= 10) pct(mean(M[, paste0("cov_", k)])) else "--"
legend("topleft", bty = "n", cex = 0.9, inset = c(0.02, 0.01),
title = "95% CI coverage", fill = COL[sel], border = NA,
legend = paste(LAB[sel], vapply(sel, cov_txt, "")))
## reductions are relative to the SRS average, which is always simulated
cmp <- setdiff(sel, "srs")
if (length(cmp)) {
se_th <- se_theory(a, n)
red <- lapply(cmp, function(k) {
mc <- if (nrep > 5) 1 - var(M[, paste0("est_", k)]) / var(M[, "est_srs"]) else NA
th <- 1 - se_th[[k]]^2 / se_th[["srs"]]^2
list(txt = sprintf("%s %s (theory %s)", LAB[[k]], pct(mc), pct(th)),
neg = if (is.na(mc)) th < 0 else mc < 0)
})
legend("bottomleft", bty = "n", cex = 0.9, inset = c(0.02, 0.01),
title = "Variance reduction vs. SRS average", title.col = "black",
fill = COL[cmp], border = NA,
legend = vapply(red, `[[`, "", "txt"),
text.col = ifelse(vapply(red, `[[`, TRUE, "neg"), col_samp, "black"))
}
})
}
shinyApp(ui, server)About this app
The app compares two ways of using a size variable \(x_i\) to estimate the mean per unit \(\bar{t}_U = t/N\) (here \(t_i\) is acres92 for county \(i\)): at the estimation stage (SRS + ratio) or at the design stage (UPSWR with \(\psi_i = x_i/t_x\)). Choosing the same variable for \(M_i\) and \(x_i\) in the app gives exactly this comparison. The plain SRS average is included as a baseline. Throughout, write \[
B = \frac{t}{t_x}, \qquad r_i = \frac{t_i}{x_i}, \qquad e_i = t_i - B x_i ,
\] so \(e_i\) is unit \(i\)’s residual from the line through the origin with the population slope. Note that \(\sum_{i=1}^N e_i = t - B t_x = 0\).
SRS + ratio. Use \(x_i\) only at the estimation stage: \[ \hat{\bar{t}}_{r} = \hat{B}\,\bar{x}_U, \quad \hat{B} = \frac{\sum_{i=1}^n t_i}{\sum_{i=1}^n x_i}, \qquad \widehat{\mathrm{SE}}(\hat{\bar{t}}_{r}) = \frac{\bar{x}_U}{\bar{x}}\sqrt{\left(1-\frac{n}{N}\right)\frac{s_e^2}{n}}, \qquad s_e^2 = \frac{1}{n-1}\sum_{i=1}^n \left(t_i - \hat{B}x_i\right)^2 . \]
UPSWR + ratio. Use \(x_i\) at the design stage, with \(\psi_i = x_i/t_x\). Every draw then has \(x_i/\psi_i = t_x\). So the HH ratio estimator, scaled by \(\bar{x}_U\), is exactly the HH estimator of the mean, \(\hat{t}_{\mathrm{HH}}/N\), and it reduces to \(\bar{x}_U\) times the average of the sampled ratios \(r_i\): \[ \hat{\bar{t}}_{\mathrm{HH}} = \bar{x}_U\,\hat{\bar{y}}_{\mathrm{HH},r} = \bar{x}_U\,\frac{\sum_{i=1}^n t_i/\psi_i}{\sum_{i=1}^n x_i/\psi_i} = \frac{\hat{t}_{\mathrm{HH}}}{N} = \bar{x}_U\,\frac{1}{n}\sum_{i=1}^n r_i , \] \[ \widehat{\mathrm{SE}}(\hat{\bar{t}}_{\mathrm{HH}}) = \bar{x}_U\,\frac{s_r}{\sqrt{n}}, \qquad s_r^2 = \frac{1}{n-1}\sum_{i=1}^n \left(r_i - \frac{1}{n}\sum_{j=1}^n r_j\right)^2 . \]
This app accompanies Shinylive App Comparing SRS and UPSWR in Elements of Sampling Survey.