Shinylive App for Illustrating Stratified Sampling
Stratified Sampling
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 800
## file: app.R
library(shiny)
data_url <- "https://raw.githubusercontent.com/longhaiSK/sampling/main/data/agpop.csv"
download.file(data_url, "agpop.csv")
dat <- read.csv("agpop.csv")
## -99 is the missing-value code in agpop.csv; recode to NA and drop
## counties missing the response (acres92) or any other numeric variable, so
## that the same population underlies every stratification choice below.
num_vars <- names(dat)[sapply(dat, is.numeric)]
for (v in num_vars) dat[[v]][dat[[v]] == -99] <- NA
dat <- dat[complete.cases(dat[, num_vars]), ]
## every numeric variable other than the response defines a stratification
## variable through its quantiles; region is the one categorical alternative
strat_vars <- setdiff(num_vars, "acres92")
pop <- dat$acres92
N <- length(pop)
ybarU <- mean(pop)
S <- sd(pop)
## build the stratum factor, population index list, and per-stratum summaries
## for any stratification variable (quantile groups of a numeric variable, or
## "region"); x_pos
## jitters each unit about its stratum's integer position, 1 (left) to H (right)
make_strata <- function(method) {
if (method %in% strat_vars) {
## unique() guards against tied quantiles (e.g. many zeros in largef*)
brk <- unique(quantile(dat[[method]], probs = c(0, 0.25, 0.75, 0.95, 1)))
stratum <- cut(dat[[method]], breaks = brk, include.lowest = TRUE)
levels(stratum) <- paste0("Q", seq_along(levels(stratum)))
} else {
stratum <- factor(dat$region)
}
H <- nlevels(stratum)
idx_h <- split(seq_len(N), stratum)
Nh <- sapply(idx_h, length)
Sh <- sapply(idx_h, function(i) sd(pop[i]))
W <- Nh / N
set.seed(42)
x_pos <- as.integer(stratum) + runif(N, -0.30, 0.30)
list(lev = levels(stratum), H = H, idx_h = idx_h,
Nh = Nh, Sh = Sh, W = W, x_pos = x_pos)
}
STRAT_CHOICES <- c(setNames(strat_vars, paste(strat_vars, "quantiles")),
"region" = "region")
## ---- broken y-axis, shared by the left panel --------------------------------
## Drawn region [0, 2 * ybarU]: visible scale [0, cut] plus a compressed lane
## for the right tail, so the population mean sits at the midpoint of the panel.
lane_frac <- 0.10
cut <- 2 * ybarU / (1 + lane_frac)
lane_w <- lane_frac * cut
ymax_top <- cut + lane_w
py <- pop # plotting positions on the y axis
tail_id <- which(pop > cut)
n_tail <- length(tail_id)
if (n_tail > 0) {
r <- rank(pop[tail_id], ties.method = "first") / (n_tail + 1)
py[tail_id] <- cut + lane_w * (0.10 + 0.80 * r)
}
## ---- fixed window for the sampling distributions ---------------------------
se100 <- sqrt(1 - 100 / N) * S / sqrt(100)
est_range <- ybarU + c(-1, 1) * 5 * se100
METHOD <- c("SRS", "Proportional", "Neyman")
MCOL <- c("grey35", "steelblue", "darkorange")
MAR <- c(6.5, 4.5, 3.0, 1.0) # identical in both panels, so the scales align
## allocation with at least 2 units per stratum, summing to n
alloc <- function(w, n, Nh) {
nh <- pmin(pmax(2, round(w / sum(w) * n)), Nh)
while (sum(nh) > n && any(nh > 2)) {
i <- which.max(nh - 2); nh[i] <- nh[i] - 1
}
while (sum(nh) < n && any(nh < Nh)) {
i <- which.max(Nh - nh); nh[i] <- nh[i] + 1
}
nh
}
## 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 <- fluidPage(
titlePanel("Shinylive App for Stratified Sampling"),
tags$style(HTML(
".well {padding: 10px 12px; margin-bottom: 10px;}
.form-group {margin-bottom: 6px;}
.radio {margin-top: 2px; margin-bottom: 2px;}
.irs {margin-bottom: 0;}
.btn {margin-bottom: 6px;}"
)),
wellPanel(
fluidRow(
column(
3,
selectInput("stratvar", "Stratification variable", choices = STRAT_CHOICES,
selected = "acres82"),
selectInput("show", "Highlight sample from",
choices = setNames(as.character(1:3), METHOD), selected = "3")
),
column(
3,
sliderInput("n", "Total sample size (n)", min = 40, max = 1000, value = 300, step = 10),
sliderInput("conf", "Confidence level", min = 0.80, max = 0.99, value = 0.95, step = 0.01)
),
column(
3,
sliderInput("speed", "Samples per second", min = 1, max = 10, value = 3, 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"))
),
br(),
actionButton("reset", "Reset", width = "100%")
)
)
),
plotOutput("mainplot", height = "520px")
)
server <- function(input, output, session) {
rv <- reactiveValues(
idx = NULL, # list of sampled indices, one element per design (last draw)
est = NULL, # matrix 3 x 3: mean, lower, upper (last draw)
mh_prop = NULL, # per-stratum sample means, proportional design (last draw)
mh_neyman = NULL, # per-stratum sample means, Neyman design (last draw)
M = NULL, # accumulated means, one column per design
HIT = NULL # accumulated coverage indicators
)
running <- reactiveVal(FALSE)
nrep <- reactive({
v <- input$nrep
if (is.null(v) || is.na(v) || v < 1) 1L else min(as.integer(v), 20000L)
})
si <- reactive(make_strata(input$stratvar))
strat_label <- reactive(names(STRAT_CHOICES)[STRAT_CHOICES == input$stratvar])
nh_prop <- reactive({ s <- si(); alloc(s$Nh, input$n, s$Nh) })
nh_neyman <- reactive({ s <- si(); alloc(s$Nh * s$Sh, input$n, s$Nh) })
# draw k replications; every replication draws its own SRS, proportional, and
# Neyman samples of total size n. Only the last replication's sampled indices
# are kept, for highlighting in the left panel; "Skip" calls this with the
# full remaining count so it jumps straight to the final violin plots.
draw_many <- function(k) {
n <- isolate(input$n)
conf <- isolate(input$conf)
a <- 1 - (1 - conf) / 2
s <- isolate(si())
H <- s$H
idx_h <- s$idx_h
Nh <- s$Nh
W <- s$W
npv <- isolate(nh_prop())
npy <- isolate(nh_neyman())
# stratified: ybar_str = sum W_h ybar_h, V = sum W_h^2 (1 - n_h/N_h) s_h^2 / n_h
# returns both the stratum-level means (mh) and the design estimate (est)
str_draw <- function(nh) {
ii <- unlist(lapply(seq_len(H), function(h) sample(idx_h[[h]], nh[h])))
yh <- split(pop[ii], rep(seq_len(H), nh))
mh <- sapply(yh, mean)
sh <- sapply(yh, sd)
m <- sum(W * mh)
v <- sum(W^2 * (1 - nh / Nh) * sh^2 / nh)
list(mh = mh, est = c(m, qt(a, df = sum(nh) - H) * sqrt(v)))
}
take <- function(nh) unlist(lapply(seq_len(H), function(h) sample(idx_h[[h]], nh[h])))
Mnew <- matrix(0, k, 3)
HITnew <- matrix(FALSE, k, 3)
for (i in seq_len(k)) {
i_srs <- sample(N, n)
y <- pop[i_srs]
m_srs <- mean(y)
se <- sqrt(1 - n / N) * sd(y) / sqrt(n)
d_srs <- qt(a, df = n - 1) * se
p <- str_draw(npv)
q <- str_draw(npy)
est <- rbind(c(m_srs, m_srs - d_srs, m_srs + d_srs),
c(p$est[1], p$est[1] - p$est[2], p$est[1] + p$est[2]),
c(q$est[1], q$est[1] - q$est[2], q$est[1] + q$est[2]))
Mnew[i, ] <- est[, 1]
HITnew[i, ] <- est[, 2] <= ybarU & ybarU <= est[, 3]
if (i == k) {
rv$idx <- list(i_srs, take(npv), take(npy))
rv$est <- est
rv$mh_prop <- p$mh
rv$mh_neyman <- q$mh
}
}
rv$M <- rbind(rv$M, Mnew)
rv$HIT <- rbind(rv$HIT, HITnew)
}
draw_one <- function() draw_many(1)
# 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))
# Skip: draw every remaining replication up to the cap in one go, landing
# directly on the final violin plots of the three estimators
observeEvent(input$skip, {
running(FALSE)
k <- nrep() - NROW(rv$M)
if (k > 0) draw_many(k)
})
observe({
if (!isTRUE(running())) return()
if (NROW(isolate(rv$M)) >= isolate(nrep())) { running(FALSE); return() }
invalidateLater(1000 / isolate(input$speed), session)
isolate(draw_one())
})
reset_history <- function() {
running(FALSE)
rv$idx <- NULL; rv$est <- NULL; rv$M <- NULL; rv$HIT <- NULL
rv$mh_prop <- NULL; rv$mh_neyman <- NULL
}
observeEvent(input$reset, reset_history())
observeEvent(input$n, reset_history())
observeEvent(input$conf, reset_history())
observeEvent(input$stratvar, reset_history())
fmt <- function(x) format(round(x), big.mark = ",")
output$mainplot <- renderPlot({
s <- si()
H <- s$H; Nh <- s$Nh; lev <- s$lev; x_pos <- s$x_pos
j <- as.integer(input$show)
nh <- switch(j, rep(NA_integer_, H), nh_prop(), nh_neyman())
layout(matrix(1:2, nrow = 1), widths = c(1, 1))
## ---- left: population by stratum, with the highlighted sample ----
par(mar = MAR, yaxs = "i")
plot(NA, xlim = c(0.3, H + 0.9), ylim = c(0, ymax_top),
xlab = "", ylab = "acres92 (thousands)", xaxt = "n", yaxt = "n", bty = "n",
main = sprintf("Population by %s stratum, with the %s sample", strat_label(), METHOD[j]))
for (h in seq_len(H)) {
rect(h - 0.36, 0, h + 0.36, ymax_top,
col = adjustcolor("grey92", 0.6), border = NA)
}
points(x_pos, py, pch = 16, cex = 0.5, col = adjustcolor("grey55", 0.5))
if (!is.null(rv$idx))
points(x_pos[rv$idx[[j]]], py[rv$idx[[j]]], pch = 16, cex = 0.6,
col = adjustcolor("red", 0.85))
# sample mean(s) of the highlighted design: one bar per stratum for the two
# stratified designs, or a single bar spanning all strata for SRS
if (!is.null(rv$idx)) {
if (j == 1) {
segments(0.3, rv$est[1, 1], H + 0.9, rv$est[1, 1], col = MCOL[1], lwd = 3)
} else {
mh <- if (j == 2) rv$mh_prop else rv$mh_neyman
segments(seq_len(H) - 0.36, mh, seq_len(H) + 0.36, mh, col = MCOL[j], lwd = 3)
}
legend("topleft", bty = "o", bg = "white", box.col = NA, cex = 0.8, inset = c(0.01, 0.01),
legend = c("Population mean", if (j == 1) "Sample mean" else "Stratum sample means"),
col = c("black", MCOL[j]), lty = c(2, 1), lwd = c(2, 3))
}
axis(1, at = seq_len(H), labels = lev, tick = FALSE, line = -0.5)
for (h in seq_len(H)) {
lab <- if (j == 1) sprintf("N = %d", Nh[h])
else sprintf("N = %d\nn = %d", Nh[h], nh[h])
mtext(lab, side = 1, at = h, line = 2.0, cex = 0.7, col = "grey30")
}
at <- pretty(c(0, cut), 6); at <- at[at <= cut]
axis(2, at = at, labels = format(at / 1000, big.mark = ",", trim = TRUE), las = 1)
segments(0.3, cut, H + 0.7, cut, col = "grey60", lty = 3)
text(H + 0.65, cut + lane_w / 2, sprintf("compressed: > %s", fmt(cut / 1000)),
adj = c(1, 0), col = "grey40", cex = 0.7)
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
## ---- right: sampling distributions of the three estimates, as violins ----
par(mar = MAR, yaxs = "i")
plot(NA, xlim = c(0.4, 3.6), ylim = est_range, xaxt = "n", yaxt = "n", bty = "n",
xlab = "", ylab = "Sample mean of acres92 (thousands)",
main = "Sampling distributions of the three estimators")
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
axis(1, at = 1:3, labels = METHOD, tick = FALSE, line = -0.5)
at_r <- pretty(est_range, 6); at_r <- at_r[at_r >= est_range[1] & at_r <= est_range[2]]
axis(2, at = at_r, labels = format(at_r / 1000, big.mark = ",", trim = TRUE), las = 1)
if (is.null(rv$M)) {
text(2, ybarU, "Click 'Draw new samples' or Play to begin", col = "grey50")
return(invisible())
}
k <- NROW(rv$M)
v_srs <- if (k > 1) var(rv$M[, 1]) else NA_real_
for (jj in 1:3) {
xc <- jj
violin(rv$M[, jj], at = xc, col = MCOL[jj])
# transient interval for the current draw
ci_col <- if (rv$HIT[k, jj]) "forestgreen" else "red"
xb <- xc + 0.42
segments(xb, rv$est[jj, 2], xb, rv$est[jj, 3], col = ci_col, lwd = 3)
segments(xb - 0.06, rv$est[jj, 2:3], xb + 0.06, rv$est[jj, 2:3], col = ci_col, lwd = 3)
points(xb, rv$est[jj, 1], pch = 19, cex = 1.1, col = ci_col)
}
legend("topleft", bty = "n", cex = 0.85, inset = c(0.02, 0.01),
title = "95% CI coverage", fill = MCOL, border = NA,
legend = paste(METHOD, vapply(1:3, function(jj) pct(mean(rv$HIT[, jj])), "")))
if (k > 1) {
# theoretical reduction, from the true population Nh/Sh and the current
# allocation, alongside the empirical reduction observed so far
V_srs_th <- (1 - input$n / N) * S^2 / input$n
nh_des <- list(nh_prop(), nh_neyman())
th_txt <- vapply(nh_des, function(nh) {
V_th <- sum(s$W^2 * (1 - nh / s$Nh) * s$Sh^2 / nh)
pct(1 - V_th / V_srs_th)
}, "")
emp_txt <- vapply(2:3, function(jj) pct(1 - var(rv$M[, jj]) / v_srs), "")
legend("bottomleft", bty = "n", cex = 0.85, inset = c(0.02, 0.01),
title = "Variance reduction vs. SRS", fill = MCOL[2:3], border = NA,
legend = sprintf("%s %s (theory %s)", METHOD[2:3], emp_txt, th_txt))
}
outside <- sum(rv$M < est_range[1] | rv$M > est_range[2])
sub <- sprintf("n = %d per design | samples drawn = %d of %d | nominal %.0f%%",
input$n, k, nrep(), 100 * input$conf)
if (outside > 0) sub <- paste0(sub, sprintf(" | %d mean(s) beyond the axis", outside))
mtext(sub, side = 3, line = -0.2, cex = 0.75, col = "grey30")
})
}
shinyApp(ui, server)
About this app
This app draws three samples of the same total size \(n\) at every step: a simple random sample, a stratified sample with proportional allocation, and a stratified sample with Neyman allocation. It uses the full agpop.csv population (with the -99 missing-value codes recoded to NA and counties with any missing value dropped), and a dropdown lets you choose the stratification variable: the quantile groups of any numeric variable other than the response acres92 (acres87, acres82, farms*, largef*, smallf*; cut at the 25th, 75th and 95th percentiles and labelled Q1 to Q4, with tied cut points merged), or the four region categories (NC, NE, S, W).
The two panels sit side by side. The left panel shows the population split into its strata along the x-axis, each drawn as a column of jittered points on a broken acres92 axis; the second dropdown chooses which of the three samples is highlighted in red, and each stratum is labelled with its \(N_h\) and the \(n_h\) that the chosen allocation assigns to it. A horizontal bar marks the sample mean of the highlighted design: one bar per stratum for the two stratified designs, or a single bar spanning the whole population for SRS, which has no stratum-level breakdown.
The right panel accumulates the sample means of all three designs on a common acres92 axis, one column per design (SRS, Proportional, Neyman), drawn as a violin once enough draws have accumulated. The bar beside each column is the current confidence interval for that design, green when it covers \(\mu\) and red when it does not, and replaced at every draw. The running CI coverage rate and the variance reduction relative to SRS are shown as legends, the latter alongside its theoretical counterpart (in brackets), computed from the true population \(N_h\) and \(S_h\) at the current allocation.
Play advances the simulation one draw at a time; Pause stops it; Skip jumps straight to the final violin plots by drawing every remaining replication up to the cap in a single step.
This app accompanies Interactive Demonstration: Stratification by acres82 or region in Elements of Sampling Survey.