Shinylive App Demonstrating Post-stratification of a Simple Random Sample
Post-stratification of a Simple Random Sample
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 840
## 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; drop the counties missing any of
## the variables offered as post-stratifiers, so one population serves them all
for (v in c("acres92", "acres87", "acres82", "farms92"))
dat[[v]][dat[[v]] == -99] <- NA
dat <- dat[complete.cases(dat[, c("acres92", "acres87", "acres82", "farms92")]), ]
pop <- dat$acres92
N <- length(pop)
ybarU <- mean(pop)
S <- sd(pop)
## The post-stratifier is chosen from a list: the real regions, quantile groups
## of a numeric auxiliary, or random labels (a variable with no signal at all).
GROUP_CHOICES <- c("region (NC, NE, S, W)" = "region",
"acres82Quantiles" = "acres82",
"acres87Quantiles" = "acres87",
"farms92Quantiles" = "farms92",
"random labels (no signal)" = "random")
make_groups <- function(key) {
if (key == "region") {
g <- factor(dat$region)
} else if (key == "random") {
set.seed(7)
g <- factor(sample(rep(paste0("R", 1:4), length.out = N)))
} else {
## The same uneven quantile cut as the stratified-sampling app: the top
## 5% of counties, where the acreages run away, get a group of their own.
## Splitting at 25/75/95 rather than into equal quarters roughly doubles
## the between-group share of the variance, so the gain is unmistakable.
x <- dat[[key]]
g <- cut(x, breaks = quantile(x, probs = c(0, 0.25, 0.75, 0.95, 1)),
include.lowest = TRUE)
levels(g) <- paste0("Q", seq_len(nlevels(g)))
}
H <- nlevels(g)
idx <- split(seq_len(N), g)
Nh <- sapply(idx, length)
W <- Nh / N # pi_h, the known population weights
Sh2 <- sapply(idx, function(i) var(pop[i]))
ybarh <- sapply(idx, function(i) mean(pop[i]))
set.seed(42)
x_pos <- as.integer(g) + runif(N, -0.30, 0.30)
list(lab = levels(g), H = H, g = g, idx = idx, Nh = Nh, W = W,
Sh2 = Sh2, ybarh = ybarh, x_pos = x_pos,
R2 = 1 - sum(W * Sh2) / S^2) # between-group share of the variance
}
## ---- the "biased SRS" option -----------------------------------------------
## Selection weights: every county carries weight 1 except those in the top 5%
## of acres82, which carry weight B. Drawing n units WITHOUT replacement with
## these weights keeps the sample size exactly n while over-representing that
## group, so the plain average drifts upward. Crucially every county inside the
## group keeps the same weight, so the sample within a group is still a simple
## random sample of it -- which is exactly what lets post-stratification undo
## the damage. The mechanism is always tied to acres82, whatever variable is
## chosen for post-stratification.
bias_grp <- cut(dat$acres82,
breaks = quantile(dat$acres82, probs = c(0, 0.25, 0.75, 0.95, 1)),
include.lowest = TRUE)
bias_top <- which(as.integer(bias_grp) == nlevels(bias_grp))
draw_sample <- function(n, B) {
if (B <= 1) return(sample(N, n))
w <- rep(1, N); w[bias_top] <- B
sample(N, n, prob = w)
}
EST <- c("Average", "Post-strat (pps)", "Regression")
MCOL <- c("#e07b39", "#2ca02c", "#9467bd") # orange average, as in the ratio app
MAR <- c(5.2, 4.5, 3.0, 1.0)
## All three estimators from ONE simple random sample.
## Average y-bar, the intercept-only fit (p = 1)
## Post-strat sum pi_h ybar_h, variance pooled with the POPULATION weights
## Regression the same estimate as the GREG on group indicators (p = H),
## variance pooled with the SAMPLE weights (the ANOVA MSE)
## Groups with fewer than 2 sampled units carry no s_h^2, so they are dropped
## and the weights renormalised -- the usual collapsing rule.
estimate_all <- function(s, ii, a) {
y <- pop[ii]
n <- length(y)
gh <- s$g[ii]
nh <- as.vector(table(gh))
ybh <- as.vector(tapply(y, gh, mean))
sh2 <- as.vector(tapply(y, gh, var))
keep <- which(!is.na(nh) & nh >= 2)
Hk <- length(keep)
Wk <- s$W[keep] / sum(s$W[keep])
fpc <- 1 - n / N
m_avg <- mean(y)
v_avg <- fpc * var(y) / n
m_ps <- sum(Wk * ybh[keep])
v_ps <- fpc * sum(Wk * sh2[keep]) / n
## pooled residual mean square of the dummy-variable regression
s2e <- sum((nh[keep] - 1) * sh2[keep]) / (sum(nh[keep]) - Hk)
m_reg <- m_ps
v_reg <- fpc * s2e / n
m <- c(m_avg, m_ps, m_reg)
se <- sqrt(c(v_avg, v_ps, v_reg))
df <- c(n - 1, n - Hk, n - Hk)
d <- qt(a, df = df) * se
list(m = m, se = se, ii = ii, nh = nh, ybh = ybh, keep = keep,
dropped = s$H - Hk,
est = cbind(m, m - d, m + d))
}
## Anticipated variances from the full population, for the "theory" figures.
theory <- function(s, n) {
fpc <- 1 - n / N
c(avg = fpc * S^2 / n,
ps = fpc * sum(s$W * s$Sh2) / n, # post-strat and regression share this
reg = fpc * sum(s$W * s$Sh2) / n)
}
## ---- broken y-axis for the left panel --------------------------------------
## Linear over [0, cut], then the sparse upper tail squeezed into a thin lane,
## so the population mean sits at mid-height.
lane_frac <- 0.10
cut <- 2 * ybarU / (1 + lane_frac)
lane_w <- lane_frac * cut
ymax_top <- cut + lane_w
py <- pop
tail_id <- which(pop > cut)
if (length(tail_id) > 0) {
r <- rank(pop[tail_id], ties.method = "first") / (length(tail_id) + 1)
py[tail_id] <- cut + lane_w * (0.10 + 0.80 * r)
}
violin <- function(x, at, col, w = 0.30) {
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)
segments(at - w * 0.6, median(x), at + w * 0.6, median(x), col = col, lwd = 2)
}
## ---- fixed window for the sampling distributions ---------------------------
## Deliberately NOT reactive. The axis has to stay put as n and the
## post-stratifier change, so that a smaller standard error shows up as a
## narrower violin rather than being rescaled away. Same convention as the
## stratified-sampling app (five standard errors of an n = 100 SRS), so the
## two can be read against one another.
se100 <- sqrt(1 - 100 / N) * S / sqrt(100)
est_range <- ybarU + c(-1, 1) * 5 * se100
pct <- function(v) if (is.na(v)) "--" else sprintf("%.1f%%", 100 * v)
fmt <- function(x) format(round(x), big.mark = ",")
## ---- left panel: the population by post-stratum, with the SRS on it ---------
## j selects which estimator's FITTED VALUES are drawn: the overall mean for
## the average (an intercept-only fit), the group means for the other two.
panel_pop <- function(s, cur, n, lab, B = 1) {
H <- s$H
par(mar = MAR, yaxs = "i")
plot(NA, xlim = c(0.3, H + 0.7), ylim = c(0, ymax_top),
xlab = "", ylab = "acres92 (thousands)", xaxt = "n", yaxt = "n", bty = "n",
main = sprintf("Population by %s, with one %s and the fitted values",
lab, if (B > 1) "biased sample" else "SRS"))
for (h in seq_len(H))
rect(h - 0.38, 0, h + 0.38, ymax_top, col = adjustcolor("grey92", 0.6), border = NA)
points(s$x_pos, py, pch = 16, cex = 0.42, col = adjustcolor("grey55", 0.5))
if (!is.null(cur)) {
points(s$x_pos[cur$ii], py[cur$ii], pch = 16, cex = 0.55,
col = adjustcolor("red", 0.85))
## Both fits at once: the average is one flat line, while the other two
## share a fitted value per group and therefore share the same estimate.
segments(0.3, cur$m[1], H + 0.7, cur$m[1], col = MCOL[1], lwd = 3)
ok <- which(!is.na(cur$ybh))
segments(ok - 0.38, cur$ybh[ok], ok + 0.38, cur$ybh[ok], col = MCOL[2], lwd = 3)
abline(h = cur$m[2], col = MCOL[2], lwd = 1.5, lty = 4)
legend("topleft", bty = "o", bg = "white", box.col = NA, cex = 0.72,
inset = c(0.005, 0.01),
legend = c("Population mean",
sprintf("%s: one fitted value", EST[1]),
"Post-strat / Regression: group means",
"their (shared) estimate"),
col = c("black", MCOL[1], MCOL[2], MCOL[2]),
lty = c(2, 1, 1, 4), lwd = c(2, 3, 3, 1.5))
}
cx <- if (H > 8) 0.55 else 0.75
for (h in seq_len(H))
mtext(s$lab[h], side = 1, at = h, line = 0.2, cex = cx, col = "grey30")
if (!is.null(cur))
mtext(sprintf("N=%d\nn=%d (exp %.0f)", s$Nh, cur$nh, n * s$W),
side = 1, at = seq_len(H), line = 2.3, cex = 0.62, col = "grey30")
mtext(paste("post-strata: the observed n_h are random; post-stratification",
"reweights them to the known pi_h"),
side = 1, line = 4.1, 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.6, cut, col = "grey60", lty = 3)
text(H + 0.6, 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 panel: the three sampling distributions -------------------------
panel_dist <- function(s, rv, est_range, th, n, conf, nrep, B = 1) {
par(mar = MAR, yaxs = "i")
plot(NA, xlim = c(0.4, 3.6), ylim = est_range, xaxt = "n", yaxt = "n", bty = "n",
xlab = "", ylab = "Estimated 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 = EST, 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 sample' or Play to begin", col = "grey50")
return(invisible())
}
k <- NROW(rv$M)
for (jj in 1:3) {
violin(rv$M[, jj], at = jj, col = MCOL[jj])
ci_col <- if (rv$HIT[k, jj]) "forestgreen" else "red"
xb <- jj + 0.40
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)
## the mean estimated SE: where Post-strat and Regression actually differ
mtext(sprintf("mean SE %.1fK", mean(rv$SE[, jj]) / 1000), side = 1, at = jj,
line = 1.1, cex = 0.7, col = MCOL[jj])
}
legend("topleft", bty = "n", cex = 0.85, inset = c(0.02, 0.01),
title = sprintf("%.0f%% CI coverage", 100 * conf), fill = MCOL, border = NA,
legend = paste(EST, vapply(1:3, function(jj) pct(mean(rv$HIT[, jj])), "")))
if (k > 1) {
v0 <- var(rv$M[, 1])
emp <- vapply(2:3, function(jj) pct(1 - var(rv$M[, jj]) / v0), "")
## the theoretical figures assume an SRS, so drop them once selection is biased
thr <- if (B > 1) c("", "") else
c(sprintf(" (theory %s)", pct(1 - th["ps"] / th["avg"])), " (same theory)")
legend("bottomleft", bty = "n", cex = 0.85, inset = c(0.02, 0.01),
title = "Variance reduction vs. the average", fill = MCOL[2:3], border = NA,
legend = c(sprintf("%s %s%s", EST[2], emp[1], thr[1]),
sprintf("%s %s%s", EST[3], emp[2], thr[2])))
if (B > 1)
legend("bottomright", bty = "n", cex = 0.8, inset = c(0.02, 0.01), text.col = "red",
legend = c("biased selection:", "the average is not estimating mu"))
}
sub <- sprintf("n = %d | between-group R2 = %s | %d of %d draws", n, pct(s$R2), k, nrep)
if (B > 1) sub <- sprintf("n = %d | acres82 top 5%% over-sampled x%.1f | R2 = %s | %d of %d draws",
n, B, pct(s$R2), k, nrep)
outside <- sum(rv$M < est_range[1] | rv$M > est_range[2])
if (outside > 0) sub <- paste0(sub, sprintf(" | %d off axis", outside))
mtext(sub, side = 3, line = -0.2, cex = 0.72, col = "grey30")
}
## ---- 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 Simulation for Post-stratification of an SRS"),
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("groupvar", "Post-stratification variable",
choices = GROUP_CHOICES, selected = "acres82")
),
column(
3,
selectInput("n", "Sample size n",
choices = c(40, 60, 100, 200, 300, 500, 1000), selected = 300),
selectInput("bias", "Biased SRS: over-sample acres82 top 5% by", choices = c("1 (none)" = 1, "1.5" = 1.5, "2" = 2,
"3" = 3, "4" = 4, "5" = 5), selected = 1)
),
column(
3,
sliderInput("conf", "Confidence level", min = 0.80, max = 0.99,
value = 0.95, step = 0.01),
sliderInput("speed", "Samples per second", min = 1, max = 10, value = 3, step = 1)
),
column(
3,
numericInput("nrep", "Stop after this many samples",
value = 500, min = 10, max = 20000, step = 50),
actionButton("draw", "Draw new sample", 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%")
)
)
),
plotOutput("mainplot", height = "540px")
)
server <- function(input, output, session) {
rv <- reactiveValues(
cur = NULL, # the last sample and its three estimates
est = NULL, # matrix 3 x 3: estimate, lower, upper (last draw)
M = NULL, # accumulated estimates, one column per estimator
SE = NULL, # accumulated standard errors
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)
})
nn <- reactive(as.numeric(input$n)) # selectInput returns character
bb <- reactive(as.numeric(input$bias))
si <- reactive(make_groups(input$groupvar))
glab <- reactive(names(GROUP_CHOICES)[GROUP_CHOICES == input$groupvar])
th <- reactive(theory(si(), nn()))
draw_many <- function(k) {
n <- isolate(nn())
a <- 1 - (1 - isolate(input$conf)) / 2
s <- isolate(si())
B <- isolate(bb())
Mnew <- matrix(0, k, 3)
SEnew <- matrix(0, k, 3)
HITnew <- matrix(FALSE, k, 3)
for (i in seq_len(k)) {
r <- estimate_all(s, draw_sample(n, B), a)
Mnew[i, ] <- r$est[, 1]
SEnew[i, ] <- r$se
HITnew[i, ] <- r$est[, 2] <= ybarU & ybarU <= r$est[, 3]
if (i == k) { rv$cur <- r; rv$est <- r$est }
}
rv$M <- rbind(rv$M, Mnew)
rv$SE <- rbind(rv$SE, SEnew)
rv$HIT <- rbind(rv$HIT, HITnew)
}
draw_one <- function() draw_many(1)
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$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$cur <- NULL; rv$est <- NULL; rv$M <- NULL; rv$SE <- NULL; rv$HIT <- NULL
}
observeEvent(input$reset, reset_history())
observeEvent(input$n, reset_history())
observeEvent(input$conf, reset_history())
observeEvent(input$bias, reset_history())
observeEvent(input$groupvar, reset_history())
output$mainplot <- renderPlot({
layout(matrix(1:2, nrow = 1), widths = c(1.15, 1))
panel_pop(si(), rv$cur, nn(), glab(), bb())
panel_dist(si(), rv, est_range, th(), nn(), input$conf, nrep(), bb())
})
}
shinyApp(ui, server)About this app
This app draws one simple random sample of \(n\) counties from the full agpop.csv population and estimates the mean of acres92 from it three ways — the three rows of the table in the chapter’s formula section, all fitted to the same sample:
- Average — the plain \(\overline{y}\). In regression terms this is the intercept-only fit: one fitted value for every county, \(p = 1\).
- Post-strat (pps) — \(\overline{y}_{post} = \sum_h \pi_h \overline{y}_h\), with the variance pooled using the known population weights, \(\hat{V} = \left(1-\frac{n}{N}\right)\frac{1}{n}\sum_h \pi_h s_h^2\).
- Regression — the same estimator written as a regression on the group indicators, with the variance pooled using the sample weights, that is the ANOVA residual mean square \(s_e^2 = \frac{1}{n-H}\sum_h (n_h-1)s_h^2\).
The last two are the same estimator, so their point estimates coincide exactly and their violins are identical; they differ only in how \(s_e^2\) is pooled, which is visible in the “mean SE” printed under each column and in the coverage rates. That is the identity of the previous subsection, drawn.
A drop-down selects the post-stratification variable: the four real regions, the quantile groups of acres82, acres87 or farms92, or random labels carrying no signal at all. The quantile groups are cut at the 25th, 75th and 95th percentiles, as in the stratified-sampling app, so the runaway top 5% of counties forms a group of its own — that uneven split lifts the between-group \(R^2\) of acres82 from 46% to 74%, and the variance reduction follows it. The remaining controls set \(n\), the confidence level, and the simulation speed.
The left panel shows the population split into post-strata, the current sample in red, and the fitted values of the chosen estimator — one flat line for the average, one line per group for the other two. Under each group it prints \(N_h\), the observed \(n_h\), and the \(n\pi_h\) that a stratified design would have fixed in advance: post-stratification exists precisely because those two differ.
The right panel accumulates the three estimators over repeated samples, with the current confidence interval beside each violin (green when it covers \(\mu\)), the running coverage, and the variance reduction relative to the plain average.
Try acres82Quantiles against random labels: the reduction tracks the between-group \(R^2\) almost exactly, and a variable that explains nothing buys nothing.
The biased-SRS option
The Biased SRS slider turns the simple random sample into a deliberately unrepresentative one. Every county keeps selection weight 1 except those in the top 5% of acres82, which get weight \(B\); the \(n\) units are then drawn without replacement with those weights, so the sample size stays exactly \(n\) while that group is over-represented. This is the standard picture of a survey that reaches the largest units too easily. Two details make it work:
- the over-sampling is applied to whole groups, and every county within the group keeps the same weight, so the units observed in a group are still a simple random sample of it — which is precisely the condition under which post-stratification repairs the damage;
- the mechanism is always tied to
acres82, whichever variable you choose to post-stratify on.
At \(B = 3\) and \(n = 300\) the plain average is biased upward by roughly 113,000 acres and its 95% intervals cover \(\mu\) about 4% of the time, while post-stratifying on acres82Quantiles removes the bias almost exactly and restores coverage to 96%. The average is no longer estimating \(\mu\) at all, so the app says so in red and drops the “theory” figures, which assume an SRS.
Now change the post-stratification variable with the bias left on. acres87Quantiles also repairs it (that variable is nearly a copy of acres82); region repairs about two-thirds of it; farms92Quantiles and random labels repair none of it and leave coverage at 4%. Post-stratification only removes the bias it can see — it corrects a distorted sample along the variable you post-stratify on, and is blind to any distortion orthogonal to it.
This app accompanies Interactive Demonstration: Post-stratifying a Simple Random Sample in Elements of Sampling Survey.