# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 4: Case-Control Studies
# Answer key for the in-lesson R activities
# Data file(s): none  (the 2x2 table is typed directly into R)
# Packages: epitools  (stretch block only; install once with install.packages("epitools"))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 1: Introduction & The Study Base / Compute the odds ratio from a 2x2 case-control table ====

# Hypothetical study: 50 lung cancer cases + 150 controls; smoking status known.
#                  Smoker   Nonsmoker
# Cases (lung Ca)     45         5
# Controls            60        90

tab <- matrix(c(45, 60,
                5,  90),
              nrow = 2, byrow = FALSE,
              dimnames = list(Status = c("Case", "Control"),
                              Smoke  = c("Yes", "No")))
tab

# Odds ratio = (a*d) / (b*c)
a <- tab["Case",    "Yes"]; b <- tab["Case",    "No"]
c <- tab["Control", "Yes"]; d <- tab["Control", "No"]
OR <- (a * d) / (b * c)
OR

# 95% CI from the log-OR (Woolf method)
log_se <- sqrt(1/a + 1/b + 1/c + 1/d)
ci <- exp(log(OR) + c(-1, 1) * 1.96 * log_se)
round(c(OR = OR, lower = ci[1], upper = ci[2]), 2)

# ---- Activity questions (rActivity-230-4-1) ----
cat("Q1: OR =", round(OR, 2), ", 95% CI", round(ci[1], 2), "to", round(ci[2], 2),
    ": cases had", round(OR, 1), "times the odds of being smokers; the CI excludes 1,",
    "so the association is statistically significant at alpha = 0.05.\n")

recip <- c(`1/a` = 1/a, `1/b` = 1/b, `1/c` = 1/c, `1/d` = 1/d)
cat("Q2: reciprocals in the Woolf SE:\n"); print(round(recip, 3))
cat("    CI ratio upper/lower =", round(ci[2] / ci[1], 1), "; cell b =", b,
    "(non-smoking cases) contributes", round(100 * (1/b) / sum(recip)), "% of the variance",
    "and drives the imprecision.\n")

cat("Q3: rows (case/control) were fixed by sampling on outcome, so a/(a+c) and b/(b+d)",
    "are sampling fractions, not risks; risk ratio and risk difference are not identifiable.",
    "The OR is invariant to outcome-based sampling, so it is the valid measure here.\n")
cat("    Check: the 'risk' of being a case among smokers, a/(a+c) =", round(a / (a + c), 2),
    ", depends only on how many controls were recruited (150 here).\n")

# ==== Section 5: Final Review & Assessment / Activity, Odds ratios and 95% CIs from a 2x2 table ====

tab <- matrix(c(45, 60,
                5,  90),
              nrow = 2, byrow = FALSE,
              dimnames = list(Status = c("Case", "Control"),
                              Smoke  = c("Yes", "No")))
tab

# Odds ratio = (a*d) / (b*c)
a <- tab["Case",    "Yes"]; b <- tab["Case",    "No"]
c <- tab["Control", "Yes"]; d <- tab["Control", "No"]
OR <- (a * d) / (b * c)
OR

# 95% CI from the log-OR (Woolf method)
log_se <- sqrt(1/a + 1/b + 1/c + 1/d)
ci <- exp(log(OR) + c(-1, 1) * 1.96 * log_se)
round(c(OR = OR, lower = ci[1], upper = ci[2]), 2)

## -----------------------------------------------------------------------------
## Stretch: same analysis with epitools::oddsratio()
## -----------------------------------------------------------------------------
# install.packages("epitools")    # uncomment if not already installed
library(epitools)
# NOTE (answer key): the lesson leaves these two lines commented out; they are run
# here so you can compare. epitools treats rows as exposure levels and columns as
# outcome levels (first row / first column = reference). Our table is laid out the
# other way round (rows = case status), but the 2x2 odds ratio is symmetric, so
# the estimate and Wald CI printed for the second row equal the by-hand result:
# 13.5 (5.07 to 35.97).
oddsratio(tab, method = "wald")