SD01 — Statistical Thinking
Course: Statistics and Probability for Data Science (DS14270).
Source deck: Samer El Kababji, SD01 — Statistical Thinking. Equations and teaching simulations below extend the deck’s concepts; they are not additional observations from the deck. Each file runs independently using base R, with knitr/R Markdown needed only for rendering.
A sample is a subset of a population. The drug question in the deck targets people with severe heart disease; it is not restricted to the one anecdotal patient. For a computational sampling example, use a fictional league of players, where teams are groups and salary is measured in thousands of dollars. No real league salaries are used.
For a finite population of size \(N\), its mean is
\[ \mu=\frac{1}{N}\sum_{i=1}^N y_i. \]
Worked example — a tiny fictional population. If the entire population has four salaries \(y=(60,70,80,90)\) (in thousands of dollars), then
\[ \mu=\frac{60+70+80+90}{4}=75. \]
The sum adds every population value. The simulation below uses a larger fictional population.
set.seed(1801)
team_sizes <- rep(40L, 10)
team <- rep(seq_along(team_sizes), team_sizes)
salary <- 50 + 10 * team + rgamma(length(team), shape = 4, scale = 5)
population <- data.frame(id = seq_along(team), team, salary)
N <- nrow(population)
mu <- mean(population$salary)
c(population_size = N, population_mean = mu)population_size population_mean
400.0000 125.6629 This known population lets us compare samples with a truth that is usually unknown in practice.
A simple random sample of \(n\) distinct people gives each person inclusion probability \(\pi_i=n/N\). Its mean estimates \(\mu\):
\[ \bar y=\frac{1}{n}\sum_{i\in s}y_i,\qquad E(\bar y)=\mu. \]
Worked example — one sample. From the four-player illustration above, suppose the selected salaries are 60 and 80. The sum includes only these sampled players:
\[ \bar y=\frac{60+80}{2}=70. \]
One sample mean can differ from \(\mu=75\); \(E(\bar y)=\mu\) describes the average over all possible simple random samples.
n <- 80
srs <- population[sample.int(N, n, replace = FALSE), ]
c(sample_size = nrow(srs), sample_mean = mean(srs$salary), population_mean = mu) sample_size sample_mean population_mean
80.0000 120.0922 125.6629 Under this design, the estimated standard error includes a finite-population correction:
\[ \widehat{\operatorname{SE}}(\bar y) =\sqrt{\left(1-\frac{n}{N}\right)\frac{s^2}{n}},\qquad s^2=\frac{\sum_{i\in s}(y_i-\bar y)^2}{n-1}. \]
Worked example — the same sample. With \(N=4\), \(n=2\), and \(\bar y=70\), expand the squared deviations first:
\[ s^2=\frac{(60-70)^2+(80-70)^2}{2-1}=200, \qquad \widehat{\operatorname{SE}}(\bar y) =\sqrt{\left(1-\frac{2}{4}\right)\frac{200}{2}} \approx7.07. \]
The standard error is in thousands of dollars.
sqrt((1 - n / N) * var(srs$salary) / n)[1] 3.141283Treat team as a stratum and sample within every team. If stratum \(h\) has population size \(N_h\), its weight is \(W_h=N_h/N\). The estimator is
\[ \widehat\mu_{\mathrm{strat}}=\sum_{h=1}^H W_h\bar y_h. \]
Worked example — one stratified estimate. In a separate small example, suppose two teams contain 40 and 60 players, with sampled mean salaries of 70 and 90. Expand the sum into one weighted term per team:
\[ \widehat\mu_{\mathrm{strat}} =\frac{40}{100}\times70+\frac{60}{100}\times90 =28+54=82. \]
The weights come from population team sizes, and salaries are in thousands of dollars.
strata <- split(population, population$team)
stratified <- do.call(rbind, lapply(strata, function(d) {
d[sample.int(nrow(d), 8, replace = FALSE), ]
}))
stratum_means <- tapply(stratified$salary, stratified$team, mean)
weights <- table(population$team) / N
stratified_mean <- sum(as.numeric(weights) * as.numeric(stratum_means))
c(sample_size = nrow(stratified), weighted_mean = stratified_mean) sample_size weighted_mean
80.0000 124.7856 With equal team sizes and eight sampled per team, the ordinary sample mean also equals the weighted mean. With unequal team sizes and equal allocation, use the population weights rather than averaging strata equally.
Now treat team as a cluster. Select \(m\) teams from \(H\) teams and include all players on those teams. Each person’s inclusion probability is \(m/H\). For known population size \(N\), a design-unbiased estimator is
\[ \widehat\mu_{\mathrm{cluster}} =\frac{H}{mN}\sum_{h\in C}\sum_{i\in h}y_i. \]
Worked example — one selected cluster. In a separate tiny league with \(H=3\) teams of two players each, \(N=6\). Suppose \(m=1\) team is selected, and its two salaries are 60 and 80. The inner sum adds salaries within that team; the outer sum has only one selected team:
\[ \widehat\mu_{\mathrm{cluster}} =\frac{3}{1\times6}(60+80)=70. \]
H <- length(strata)
m <- 2
selected_teams <- sample(seq_len(H), m, replace = FALSE)
cluster_sample <- population[population$team %in% selected_teams, ]
cluster_mean <- H / (m * N) * sum(cluster_sample$salary)
c(selected_teams = selected_teams, sample_size = nrow(cluster_sample),
estimated_mean = cluster_mean)selected_teams1 selected_teams2 sample_size estimated_mean
3.0000 10.0000 80.0000 136.2494 Equal cluster sizes make this estimator equal to the sample mean here. Players on one team share a salary component, so 80 clustered observations need not contain as much information as 80 simple-random-sample observations. The SRS standard-error formula above does not apply to this clustered design.
Select \(m\) teams, then select \(k\) players from each selected team of size \(N_h\). The person’s inclusion probability and an inverse-probability estimator are
\[ \pi_{hi}=\frac{m}{H}\frac{k}{N_h},\qquad \widehat\mu_{\mathrm{multi}} =\frac{1}{N}\sum_{h\in C}\sum_{i\in s_h}\frac{y_{hi}}{\pi_{hi}}. \]
Worked example — one selected team. Use the tiny six-player league above, but now sample \(k=1\) player from the selected team of size \(N_h=2\). If that player’s salary is 80, then
\[ \pi_{hi}=\frac{1}{3}\times\frac{1}{2}=\frac{1}{6}, \qquad \widehat\mu_{\mathrm{multi}} =\frac{1}{6}\left(\frac{80}{1/6}\right)=80. \]
Here both sums have one term. Dividing the observed salary by its inclusion probability gives its estimated contribution to the population total, which is then divided by \(N\). The code below uses the larger simulated league.
k <- 10
multistage <- do.call(rbind, lapply(selected_teams, function(h) {
d <- strata[[as.character(h)]]
sampled <- d[sample.int(nrow(d), k, replace = FALSE), ]
sampled$inclusion_probability <- (m / H) * (k / nrow(d))
sampled
}))
multistage_mean <- sum(multistage$salary / multistage$inclusion_probability) / N
c(sample_size = nrow(multistage), estimated_mean = multistage_mean) sample_size estimated_mean
20.0000 134.8209 Taking players only from the first two teams gives zero inclusion probability to everyone else. Repeating that procedure cannot remove its systematic bias. For comparison, repeated SRS estimates fluctuate around the population mean:
\[ \operatorname{Bias}(\widehat\mu)=E(\widehat\mu)-\mu. \]
Worked example — average over all possible samples. Suppose a tiny population has mean \(\mu=10\), and a sampling procedure has two equally likely possible samples whose estimates are 11 and 13. Then
\[ E(\widehat\mu)=\tfrac12(11)+\tfrac12(13)=12,\qquad \operatorname{Bias}(\widehat\mu)=12-10=2. \]
The bias is the average error over the sampling procedure, rather than the error of just one observed sample.
convenience <- population[population$team %in% c(1, 2), ]
repeated_srs <- replicate(1000,
mean(population$salary[sample.int(N, 80, replace = FALSE)]))
c(population_mean = mu, convenience_mean = mean(convenience$salary),
average_srs_mean = mean(repeated_srs),
simulated_srs_bias = mean(repeated_srs) - mu) population_mean convenience_mean average_srs_mean simulated_srs_bias
125.66288137 85.17447250 125.65113548 -0.01174589 hist(repeated_srs, breaks = 25, xlab = "Mean from an SRS of 80 players",
main = "Sampling variation with a known population")
abline(v = mu, col = "steelblue", lwd = 2)
abline(v = mean(convenience$salary), col = "tomato", lwd = 2)
legend("topright", c("Population mean", "Convenience mean"),
col = c("steelblue", "tomato"), lwd = 2)
Random selection reduces selection bias under the sampling design, but nonresponse or incomplete population coverage can still bias the final data.
Try it: select only one cluster. Explain why its mean can differ greatly from the population mean even though the cluster was selected randomly.