Shinylive App Comparing Two-Stage Cluster Sampling (ratio estimator) with SRS
Two-Stage Cluster Sampling (ratio estimator) with SRS
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 820
## 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 the
## counties with no acres92, so that a single population underlies both
## clustering choices below.
dat$acres92[dat$acres92 == -99] <- NA
dat <- dat[!is.na(dat$acres92), ]
pop <- dat$acres92
Nel <- length(pop) # number of elements (counties) in the population
ybarU <- mean(pop)
S <- sd(pop)
## Build the clustering: "state" groups the counties as they really are,
## "random" permutes the state labels, so the clusters keep exactly the same
## sizes M_i but lose every trace of within-cluster homogeneity. Comparing the
## two isolates the effect of homogeneity inside clusters, and nothing else.
make_clusters <- function(method) {
lab <- factor(dat$state)
if (method == "random") {
set.seed(7)
lab <- sample(lab) # same sizes, scrambled membership
levels(lab) <- sprintf("C%02d", seq_len(nlevels(lab)))
}
H <- nlevels(lab) # N, the number of clusters
idx <- split(seq_len(Nel), lab)
Mi <- sapply(idx, length)
ti <- sapply(idx, function(i) sum(pop[i]))
Si2 <- sapply(idx, function(i) if (length(i) > 1) var(pop[i]) else 0)
set.seed(42)
x_pos <- as.integer(lab) + runif(Nel, -0.32, 0.32)
list(lev = levels(lab), H = H, idx = idx, Mi = Mi, ti = ti, Si2 = Si2,
M0 = sum(Mi), Mbar = mean(Mi), x_pos = x_pos,
R2 = 1 - sum(sapply(idx, function(i) sum((pop[i] - mean(pop[i]))^2))) /
sum((pop - ybarU)^2))
}
CLUST_CHOICES <- c("state (real clusters)" = "state",
"Random (scrambled labels)" = "random")
## subsample size in cluster i: m_i = p * M_i, at least 1 and at most M_i
subsize <- function(Mi, p) pmax(1, pmin(Mi, round(p * Mi)))
## Anticipated variances from the FULL population, for the "theory" figures.
## Two-stage ratio estimator (between-cluster term + within-cluster term):
## V = (1 - n/N) S_t^2 / (n Mbar^2)
## + sum_U M_i^2 (1 - m_i/M_i) S_i^2 / m_i / (n N Mbar^2),
## with S_t^2 the population variance of the ratio residuals t_i - ybarU M_i.
## SRS is compared at the sample size the cluster design delivers on average.
theory <- function(s, n, p) {
mi <- subsize(s$Mi, p)
St2 <- sum((s$ti - ybarU * s$Mi)^2) / (s$H - 1)
Vc <- (1 - n / s$H) * St2 / (n * s$Mbar^2) +
sum(s$Mi^2 * (1 - mi / s$Mi) * s$Si2 / mi) / (n * s$H * s$Mbar^2)
nel <- max(2, round(n * mean(mi)))
Vs <- (1 - nel / s$M0) * S^2 / nel
list(Vc = Vc, Vs = Vs, nel = nel)
}
## One two-stage cluster sample: an SRS of n clusters, then an SRS of
## m_i = p M_i counties inside each sampled cluster. The estimate is the ratio
## estimator with M_i as the auxiliary variable, ybar_r = sum(t_i hat)/sum(M_i),
## and its variance is estimated exactly as cluster_ratio() does: the ratio
## variance over the sampled clusters, plus n/N times the stratified variance
## V_str = sum_i (M_i/M)^2 (1 - m_i/M_i) s_i^2/m_i, M = sum_i M_i,
## that treats the sampled clusters as strata of sizes M_i.
clu_draw <- function(s, n, p, a) {
ci <- sort(sample(s$H, n))
Mis <- s$Mi[ci]
mi <- subsize(Mis, p)
ii <- unlist(lapply(seq_len(n), function(k) sample(s$idx[[ci[k]]], mi[k])))
yl <- split(pop[ii], rep(seq_len(n), mi))
ybari <- sapply(yl, mean)
si2 <- sapply(yl, function(v) if (length(v) > 1) var(v) else 0)
that <- Mis * ybari # t_i hat = M_i ybar_i
Mbar <- mean(Mis)
m <- sum(that) / sum(Mis)
Mtot <- sum(Mis)
v1 <- (1 - n / s$H) * sum((that - m * Mis)^2) / (n - 1) / (n * Mbar^2)
v2 <- (n / s$H) * sum((Mis / Mtot)^2 * (1 - mi / Mis) * si2 / mi)
list(ci = ci, mi = mi, Mi = Mis, ii = ii, ybari = ybari,
est = c(m, qt(a, df = n - 1) * sqrt(v1 + v2)))
}
## One simple random sample of counties, of the same size as the cluster
## sample just drawn, estimated by the plain sample average.
srs_draw <- function(nel, a) {
ii <- sample(Nel, nel)
y <- pop[ii]
m <- mean(y)
se <- sqrt(1 - nel / Nel) * sd(y) / sqrt(nel)
list(ii = ii, est = c(m, qt(a, df = nel - 1) * se))
}
## ---- 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
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)
}
METHOD <- c("SRS average", "Cluster ratio")
MCOL <- c("grey35", "#2ca02c")
MAR <- c(5.0, 4.5, 3.0, 1.0)
## symmetric density polygon, drawn vertically about x = at
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)
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)
fmt <- function(x) format(round(x), big.mark = ",")
## ---- left panel: the population, cluster by cluster -------------------------
## clu is the last cluster draw, srs the last SRS draw (either may be NULL);
## j = 1 highlights the SRS sample, j = 2 the two-stage cluster sample.
panel_pop <- function(s, j, clu, srs, p, lab) {
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 the %s sample", lab, METHOD[j]))
## one band per cluster; the clusters that made it into the two-stage
## sample are shaded, which is the whole point of the design
inS <- rep(FALSE, H)
if (j == 2 && !is.null(clu)) inS[clu$ci] <- TRUE
for (h in seq_len(H))
rect(h - 0.38, 0, h + 0.38, ymax_top,
col = if (inS[h]) "#fff0c2" else adjustcolor("grey92", 0.6), border = NA)
points(s$x_pos, py, pch = 16, cex = 0.42, col = adjustcolor("grey55", 0.5))
smp <- if (j == 1) srs else clu
if (!is.null(smp))
points(s$x_pos[smp$ii], py[smp$ii], pch = 16, cex = 0.55,
col = adjustcolor("red", 0.85))
## sample mean(s): one bar per sampled cluster for the two-stage design
## (the ybar_i that build t_i hat), a single bar for SRS
if (!is.null(smp)) {
if (j == 1) {
segments(0.3, srs$est[1], H + 0.7, srs$est[1], col = MCOL[1], lwd = 3)
} else {
segments(clu$ci - 0.38, clu$ybari, clu$ci + 0.38, clu$ybari,
col = MCOL[2], lwd = 3)
abline(h = clu$est[1], col = MCOL[2], lwd = 1.5, lty = 4)
}
leg <- if (j == 1)
list(txt = c("Population mean", "Sample mean"),
col = c("black", MCOL[1]), lty = c(2, 1), lwd = c(2, 3))
else
list(txt = c("Population mean", "Cluster means (ybar_i)",
"Ratio estimate (ybar_r)"),
col = c("black", MCOL[2], MCOL[2]), lty = c(2, 1, 4),
lwd = c(2, 3, 1.5))
legend("topleft", bty = "o", bg = "white", box.col = NA, cex = 0.75,
inset = c(0.005, 0.01), legend = leg$txt,
col = leg$col, lty = leg$lty, lwd = leg$lwd)
}
## cluster labels: sampled clusters in red, with m_i / M_i when there is room
cx <- if (H > 40) 0.45 else 0.6
for (h in seq_len(H)) {
mtext(s$lev[h], side = 1, at = h, line = 0.2, las = 2, cex = cx,
col = if (inS[h]) "red" else "grey45", font = if (inS[h]) 2 else 1)
}
if (j == 2 && !is.null(clu) && length(clu$ci) <= 20)
mtext(sprintf("%d/%d", clu$mi, clu$Mi), side = 1, at = clu$ci, line = 1.9,
las = 2, cex = 0.5, col = "red")
mtext(if (j == 2) sprintf("clusters (sampled ones shaded%s; p = %.0f%%)",
if (length(clu$ci) <= 20) ", labelled m_i/M_i" else "",
100 * p)
else "clusters",
side = 1, line = 3.6, cex = 0.72, 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.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 two sampling distributions ----------------------------
panel_dist <- function(s, rv, est_range, th, n, p, conf, nrep) {
par(mar = MAR, yaxs = "i")
plot(NA, xlim = c(0.4, 2.6), ylim = est_range, xaxt = "n", yaxt = "n", bty = "n",
xlab = "", ylab = "Estimated mean of acres92 (thousands)",
main = "Sampling distributions of the two estimators")
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
axis(1, at = 1:2, 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(1.5, ybarU, "Click 'Draw new samples' or Play to begin", col = "grey50")
return(invisible())
}
k <- NROW(rv$M)
for (jj in 1:2) {
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)
}
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(METHOD, vapply(1:2, function(jj) pct(mean(rv$HIT[, jj])), "")))
## design effect: variance of the cluster estimator relative to SRS at the
## same sample size, simulated so far and anticipated from the population
if (k > 1) {
deff <- var(rv$M[, 2]) / var(rv$M[, 1])
legend("bottomleft", bty = "n", cex = 0.85, inset = c(0.02, 0.01),
title = "Variance of cluster ratio vs. SRS", fill = MCOL[2], border = NA,
text.col = if (deff > 1) "red" else "black",
legend = sprintf("deff = %.2f (theory %.2f)", deff, th$Vc / th$Vs))
}
outside <- sum(rv$M < est_range[1] | rv$M > est_range[2])
sub <- sprintf("n = %d clusters, %s counties | R2 between clusters = %s | %d of %d draws%s",
n, fmt(rv$nel), pct(s$R2), k, nrep,
if (outside > 0) sprintf(" | %d off axis", outside) else "")
mtext(sub, side = 3, line = -0.2, cex = 0.7, 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 Two-Stage Cluster 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,
radioButtons("clustvar", "Cluster variable", choices = CLUST_CHOICES,
selected = "state"),
radioButtons("show", "Highlight sample from", choiceNames = METHOD,
choiceValues = as.character(1:2), selected = "2")
),
column(
3,
sliderInput("n", "Number of clusters drawn (n)", min = 2, max = 50,
value = 10, step = 1),
sliderInput("p", "Percentage sampled per cluster (p = m_i/M_i)",
min = 0.05, max = 1, value = 0.30, step = 0.05)
),
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 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%")
)
)
),
plotOutput("mainplot", height = "530px")
)
server <- function(input, output, session) {
rv <- reactiveValues(
clu = NULL, # last two-stage cluster draw
srs = NULL, # last simple random sample
est = NULL, # matrix 2 x 3: estimate, lower, upper (last draw)
nel = NA, # number of counties in the last sample
M = NULL, # accumulated estimates, 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_clusters(input$clustvar))
clab <- reactive(names(CLUST_CHOICES)[CLUST_CHOICES == input$clustvar])
th <- reactive(theory(si(), input$n, input$p))
## a window wide enough for whichever design is the more variable
est_range <- reactive({
t0 <- th()
ybarU + c(-1, 1) * 4 * sqrt(max(t0$Vc, t0$Vs))
})
## draw k replications; each one draws a two-stage cluster sample and an SRS
## of the same realised size. Only the last replication's sampled units are
## kept, for the left panel; "Skip" calls this with the full remaining count.
draw_many <- function(k) {
n <- isolate(input$n)
p <- isolate(input$p)
conf <- isolate(input$conf)
a <- 1 - (1 - conf) / 2
s <- isolate(si())
Mnew <- matrix(0, k, 2)
HITnew <- matrix(FALSE, k, 2)
for (i in seq_len(k)) {
cl <- clu_draw(s, n, p, a)
sr <- srs_draw(length(cl$ii), a)
est <- rbind(c(sr$est[1], sr$est[1] - sr$est[2], sr$est[1] + sr$est[2]),
c(cl$est[1], cl$est[1] - cl$est[2], cl$est[1] + cl$est[2]))
Mnew[i, ] <- est[, 1]
HITnew[i, ] <- est[, 2] <= ybarU & ybarU <= est[, 3]
if (i == k) {
rv$clu <- cl; rv$srs <- sr; rv$est <- est; rv$nel <- length(cl$ii)
}
}
rv$M <- rbind(rv$M, Mnew)
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$clu <- NULL; rv$srs <- NULL; rv$est <- NULL; rv$nel <- NA
rv$M <- NULL; rv$HIT <- NULL
}
observeEvent(input$reset, reset_history())
observeEvent(input$n, reset_history())
observeEvent(input$p, reset_history())
observeEvent(input$conf, reset_history())
observeEvent(input$clustvar, reset_history())
output$mainplot <- renderPlot({
layout(matrix(1:2, nrow = 1), widths = c(1.3, 1))
panel_pop(si(), as.integer(input$show), rv$clu, rv$srs, input$p, clab())
panel_dist(si(), rv, est_range(), th(), input$n, input$p, input$conf, nrep())
})
}
shinyApp(ui, server)About this app
This app repeats the chapter’s worked example many times over. At every step it draws a two-stage cluster sample — an SRS of \(n\) clusters, then an SRS of \(m_i = p \, M_i\) counties inside each — together with an ordinary SRS of counties of the same size, from the full agpop.csv population (with the -99 codes dropped). It estimates the mean of acres92 from each: the ratio estimator \(\overline{y}_r = \sum \hat{t}_i / \sum M_i\) for the cluster sample, and the plain sample average for the SRS.
A radio button chooses the clustering variable:
- state — the 50 real states, the natural clusters of
agpop.csv. - Random — the state labels randomly permuted, which keeps the cluster sizes \(M_i\) exactly as they are but scrambles which counties belong together.
Two sliders set the design: \(n\), the number of clusters drawn, and \(p = m_i/M_i\), the percentage of counties subsampled from each selected cluster (\(p = 100\%\) gives a one-stage sample). Play advances the simulation one draw at a time, Pause stops it, Skip jumps straight to the end, and changing any setting clears the accumulated history.
The left panel shows the population split into its clusters, with the highlighted sample in red. For the cluster design the selected clusters are shaded and labelled \(m_i/M_i\), a bar marks each subsample mean \(\overline{y}_i\), and the ratio estimate \(\overline{y}_r\) runs across the panel; for the SRS a single bar marks the sample mean.
The right panel accumulates the estimates of both designs, one column each, drawn as a violin once enough draws have accumulated. The bar beside each column is the current confidence interval, green when it covers \(\mu\) and red when it does not. The legends report the running CI coverage and the design effect — the variance of the ratio estimator relative to the SRS average, simulated and theoretical — and the subtitle gives the between-cluster \(R^2\) of acres92.
This app accompanies Interactive Demonstration: Clustering by State or at Random in Elements of Sampling Survey.