# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 5: Cohort Studies
# Answer key for the in-lesson R activities
# Data file(s): none  (the cohort summary counts are typed directly into R)
# Packages: none beyond base R
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 1: Introduction & Cohort Study Design / Risk ratio and incidence rate from cohort data ====

# 1000 exposed and 1000 unexposed individuals followed for up to 5 years.
#   exposed: 80 events in 4500 person-years
# unexposed: 30 events in 4900 person-years

events <- c(exposed = 80,   unexposed = 30)
n      <- c(exposed = 1000, unexposed = 1000)
py     <- c(exposed = 4500, unexposed = 4900)

risk <- events / n
rate <- events / py * 1000          # per 1000 person-years

RR   <- risk["exposed"] / risk["unexposed"]   # risk ratio
IRR  <- rate["exposed"] / rate["unexposed"]   # incidence rate ratio

round(data.frame(risk, rate, RR = RR, IRR = IRR), 3)

# ---- Activity questions (rActivity-230-5-1) ----
cat("Q1: 5-year risk =", risk["exposed"], "(exposed) vs", risk["unexposed"], "(unexposed);",
    "RR =", round(RR, 3), ": exposed people had", round(RR, 2), "times the cumulative risk;",
    "risk difference =", risk["exposed"] - risk["unexposed"], "(", 1000 * (risk["exposed"] - risk["unexposed"]),
    "per 1000 over 5 years).\n")

cat("Q2: IRR =", round(IRR, 3), "> RR =", round(RR, 3), "because the exposed group accrued",
    py["exposed"], "person-years vs", py["unexposed"], ": mean follow-up",
    round(py["exposed"] / n["exposed"], 2), "vs", round(py["unexposed"] / n["unexposed"], 2),
    "years per person. The exposed group lost more person-time (earlier events and/or",
    "earlier censoring), so dividing by person-years raises its rate relative to its risk.\n")

# Q3: closed cohort, everyone followed exactly 5 years, no losses.
#     Person-time = 5 years for non-cases plus the time to event for cases;
#     with events spread evenly over follow-up, cases contribute ~2.5 years each.
py_closed <- (n - events) * 5 + events * 2.5
rate_closed <- events / py_closed * 1000
IRR_closed <- rate_closed["exposed"] / rate_closed["unexposed"]
cat("Q3: closed cohort with no losses: person-years would be", py_closed["exposed"], "and",
    py_closed["unexposed"], "; IRR =", round(IRR_closed, 3), "vs RR =", round(RR, 3),
    ". The two measures converge (a small gap remains because the group with more events",
    "loses a little more person-time). Report the risk ratio in that setting; use rates",
    "(IRR/HR) once follow-up differs between groups.\n")

# ==== Section 5: Final Review & Assessment / Activity, Risk, rate, RR/IRR, and a 95% CI for the rate ratio ====

# 1000 exposed and 1000 unexposed individuals followed for up to 5 years.
#   exposed:  80 events in 4500 person-years
# unexposed:  30 events in 4900 person-years

events <- c(exposed = 80,   unexposed = 30)
n      <- c(exposed = 1000, unexposed = 1000)
py     <- c(exposed = 4500, unexposed = 4900)

risk <- events / n
rate <- events / py * 1000          # per 1000 person-years

RR  <- risk["exposed"] / risk["unexposed"]   # risk ratio
IRR <- rate["exposed"] / rate["unexposed"]   # incidence rate ratio

round(data.frame(risk, rate, RR = RR, IRR = IRR), 3)

## -----------------------------------------------------------------------------
## Stretch: 95% CI for the rate ratio (Wald approximation on the log scale)
## -----------------------------------------------------------------------------
log_irr   <- log(IRR)
se_logirr <- sqrt(1/events["exposed"] + 1/events["unexposed"])
ci_irr    <- exp(log_irr + c(-1, 1) * 1.96 * se_logirr)
round(c(IRR = IRR, lower = ci_irr[1], upper = ci_irr[2]), 3)
cat("Stretch: IRR =", round(IRR, 3), "with 95% CI", round(ci_irr[1], 3), "to",
    round(ci_irr[2], 3), "; the CI excludes 1.\n")