# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 11: Confounding and Statistical Inference
# Answer key for the in-lesson R activities
# Data file(s): none  (the 5,000-person dataset 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: Confounding / Stratified analysis: confounding revealed and removed ====

set.seed(230)
n <- 5000
ses <- rbinom(n, 1, 0.5)                                  # 1 = high SES
hrt <- rbinom(n, 1, prob = ifelse(ses == 1, 0.6, 0.2)) # high SES uses HRT more
# CVD risk: lower in high SES, NOT affected by HRT (true null)
cvd <- rbinom(n, 1, prob = ifelse(ses == 1, 0.05, 0.15))

# Crude (unadjusted) OR -- looks "protective"
exp(coef(glm(cvd ~ hrt, family = binomial))["hrt"])

# Within each SES stratum -- the true effect: ~1.0
tapply(seq_len(n), ses, function(i) {
  exp(coef(glm(cvd[i] ~ hrt[i], family = binomial))[2])
})

# Adjusted OR, controlling for SES
exp(coef(glm(cvd ~ hrt + ses, family = binomial))["hrt"])

# ---- Activity questions (rActivity-230-11-1) ----
or_ci <- function(fit, term = 2) {
  est <- coef(fit)[term]; se <- sqrt(diag(vcov(fit)))[term]
  round(exp(c(OR = unname(est), lower = unname(est - 1.96 * se), upper = unname(est + 1.96 * se))), 2)
}
crude    <- or_ci(glm(cvd ~ hrt, family = binomial))
strata   <- tapply(seq_len(n), ses, function(i) or_ci(glm(cvd[i] ~ hrt[i], family = binomial)))
adjusted <- or_ci(glm(cvd ~ hrt + ses, family = binomial))

cat("Q1: crude OR =", crude["OR"], "(95% CI", crude["lower"], "to", crude["upper"], ").",
    "Conditions of confounding, each built into the simulation:\n")
cat("    (1) SES -> exposure: P(HRT | high SES) =", round(mean(hrt[ses == 1]), 2),
    "vs", round(mean(hrt[ses == 0]), 2), "in low SES.\n")
cat("    (2) SES -> outcome, independent of HRT: P(CVD | high SES) =", round(mean(cvd[ses == 1]), 3),
    "vs", round(mean(cvd[ses == 0]), 3), "in low SES (HRT never enters the cvd line).\n")
cat("    (3) Not on the causal pathway: SES is generated before HRT and is not caused by it.\n")

cat("Q2: stratum-specific ORs: low SES", strata[["0"]]["OR"], "(", strata[["0"]]["lower"], "to",
    strata[["0"]]["upper"], "), high SES", strata[["1"]]["OR"], "(", strata[["1"]]["lower"], "to",
    strata[["1"]]["upper"], "); SES-adjusted OR", adjusted["OR"], "(", adjusted["lower"], "to",
    adjusted["upper"], "). All three are close to 1 and to each other, and every CI includes 1:",
    "SES is a confounder, not an effect modifier (modification would show clearly different",
    "ORs across strata).\n")
tab <- table(hrt, cvd, ses)
mh <- mantelhaen.test(tab)
cat("    Mantel-Haenszel common OR =", round(unname(mh$estimate), 2), "(95% CI",
    paste(round(mh$conf.int, 2), collapse = " to "), "), p =", round(mh$p.value, 2), "\n")

cat("Q3: the crude OR of", crude["OR"], "arises with HRT doing nothing, purely from an",
    "SES confounder; observational data can only adjust for confounders that were measured,",
    "so a trial that randomises HRT (making exposure independent of all baseline",
    "characteristics, measured or not) was needed to settle the question.\n")

# ==== Section 3: Final Assessment / Activity: Crude vs stratified vs adjusted: an HRT/SES confounding demo ====

set.seed(230)
n   <- 5000
ses <- rbinom(n, 1, 0.5)                                   # 1 = high SES
hrt <- rbinom(n, 1, prob = ifelse(ses == 1, 0.6, 0.2))     # high SES uses HRT more

# CVD risk: lower in high SES, NOT affected by HRT (true null)
cvd <- rbinom(n, 1, prob = ifelse(ses == 1, 0.05, 0.15))

# Crude (unadjusted) OR -- looks "protective"
exp(coef(glm(cvd ~ hrt, family = binomial))["hrt"])

# Within each SES stratum -- the true effect: ~1.0
tapply(seq_len(n), ses, function(i) {
  exp(coef(glm(cvd[i] ~ hrt[i], family = binomial))[2])
})

# Adjusted OR, controlling for SES
exp(coef(glm(cvd ~ hrt + ses, family = binomial))["hrt"])