# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 8: Sampling, Selection, and External Validity
# Answer key for the in-lesson R activities
# Data file(s): none  (the 10,000-person population is simulated with set.seed(230))
# Packages: none beyond base R
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 1: Selection Bias Mechanisms / Simulate the healthy-worker effect with a convenience sample ====

set.seed(230)

# 10,000-person "population" with true mean BMI = 27
N        <- 10000
bmi      <- rnorm(N, mean = 27, sd = 5)
gym_goer <- rbinom(N, 1, prob = plogis(2 - 0.15*bmi))   # lower BMI -> more likely

mean(bmi)                              # truth: ~27

# (1) Simple random sample of 200
srs <- sample(bmi, size = 200)
mean(srs)

# (2) Convenience sample: 200 gym-goers
conv <- sample(bmi[gym_goer == 1], size = 200)
mean(conv)

# Stretch: 1000 replicates of each strategy
srs_means  <- replicate(1000, mean(sample(bmi, 200)))
conv_means <- replicate(1000, mean(sample(bmi[gym_goer == 1], 200)))

par(mfrow = c(1, 2))
hist(srs_means,  main = "Simple Random Sample", xlim = c(20, 30), xlab = "Mean BMI")
abline(v = 27, col = "red", lwd = 2)
hist(conv_means, main = "Convenience (gym)",    xlim = c(20, 30), xlab = "Mean BMI")
abline(v = 27, col = "red", lwd = 2)
par(mfrow = c(1, 1))

# ---- Activity questions (rActivity-230-8-1) ----
cat("Q1: true mean BMI =", round(mean(bmi), 2), "; mean(srs) =", round(mean(srs), 2),
    "(", round(mean(srs) - mean(bmi), 2), "from truth ); mean(conv) =", round(mean(conv), 2),
    "(", round(mean(conv) - mean(bmi), 2), "from truth, i.e.", round(mean(bmi) - mean(conv), 1),
    "BMI units too low ).\n")
cat("    Over 1000 replicates: SRS means centre at", round(mean(srs_means), 2),
    "(SD", round(sd(srs_means), 2), "); convenience means centre at", round(mean(conv_means), 2),
    "(SD", round(sd(conv_means), 2), "), about", round(mean(bmi) - mean(conv_means), 1),
    "units below the truth.\n")

p_sel <- plogis(2 - 0.15 * c(20, 27, 35))
cat("Q2: P(gym-goer) at BMI 20 =", round(p_sel[1], 3), ", at BMI 27 =", round(p_sel[2], 3),
    ", at BMI 35 =", round(p_sel[3], 3), ": lean people are about",
    round(p_sel[1] / p_sel[3]), "times as likely to enter the convenience frame, so the",
    "sampled subset (", sum(gym_goer), "gym-goers, mean BMI", round(mean(bmi[gym_goer == 1]), 2),
    ") is lighter than the population by construction.\n")

conv_means_400 <- replicate(1000, mean(sample(bmi[gym_goer == 1], 400)))
cat("Q3: doubling the convenience sample to 400 gives replicate means centred at",
    round(mean(conv_means_400), 2), "(SD", round(sd(conv_means_400), 2),
    "): narrower, same wrong centre. Selection bias is structural, like the",
    "healthy-worker SMR below 1.0; only a different sampling frame or an adjustment",
    "for the selection mechanism fixes it.\n")

# ==== Section 4: Final Assessment / Activity: Convenience sampling and the healthy-worker effect ====

set.seed(230)

# 10,000-person "population" with true mean BMI = 27
N        <- 10000
bmi      <- rnorm(N, mean = 27, sd = 5)
gym_goer <- rbinom(N, 1, prob = plogis(2 - 0.15*bmi))   # lower BMI -> more likely

mean(bmi)                              # truth: ~27

# (1) Simple random sample of 200
srs <- sample(bmi, size = 200)
mean(srs)

# (2) Convenience sample: 200 gym-goers
conv <- sample(bmi[gym_goer == 1], size = 200)
mean(conv)

## -----------------------------------------------------------------------------
## Stretch: replicate 1000 times and look at the sampling distribution
## -----------------------------------------------------------------------------
srs_means  <- replicate(1000, mean(sample(bmi, 200)))
conv_means <- replicate(1000, mean(sample(bmi[gym_goer == 1], 200)))

par(mfrow = c(1, 2))
hist(srs_means,  main = "Simple Random Sample", xlim = c(20, 30), xlab = "Mean BMI")
abline(v = 27, col = "red", lwd = 2)
hist(conv_means, main = "Convenience (gym)",     xlim = c(20, 30), xlab = "Mean BMI")
abline(v = 27, col = "red", lwd = 2)
par(mfrow = c(1, 1))