# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 10: Design-Specific and Temporal Biases
# Answer key for the in-lesson R activities
# Data file(s): none  (the 2,000-patient cohort is simulated with set.seed(230))
# Packages: none beyond base R
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 2: Immortal Time & Lead-Time Bias / Immortal time bias in 30 lines of code ====

set.seed(230)
N <- 2000

# Each patient survives Exp(rate=1/3) years; about 60% eventually fill a Rx.
# If they die before filling, they go in the unexposed group (no benefit).
death_t <- rexp(N, rate = 1/3)
rx_t    <- rexp(N, rate = 1/2)
ever_rx <- rx_t < death_t

## (1) WRONG analysis: classify by ever_rx, ignore immortal time
mean_surv_wrong <- tapply(death_t, ever_rx, mean)
mean_surv_wrong   # exposed group looks much better

## (2) RIGHT analysis: split person-time at the prescription date
unexp_pt <- pmin(death_t, rx_t)              # every pt's unexposed person-time
exp_pt   <- pmax(death_t - rx_t, 0) * ever_rx # only after Rx, only if they got it

unexp_events <- sum(!ever_rx)
exp_events   <- sum(ever_rx)
c(unexp_rate = unexp_events / sum(unexp_pt),
  exp_rate   = exp_events   / sum(exp_pt))

# ---- Activity questions (rActivity-230-10-1) ----
rates <- c(unexp_rate = unexp_events / sum(unexp_pt), exp_rate = exp_events / sum(exp_pt))
immortal_pt <- sum(rx_t[ever_rx])       # pre-prescription time of the ever-treated
cat("Q1: wrong analysis: mean survival", round(mean_surv_wrong["FALSE"], 2), "years (never Rx) vs",
    round(mean_surv_wrong["TRUE"], 2), "years (ever Rx): an apparent",
    round(mean_surv_wrong["TRUE"] - mean_surv_wrong["FALSE"], 2), "-year advantage.",
    round(100 * mean(ever_rx)), "% ever filled a prescription. The", round(immortal_pt),
    "person-years the treated spent waiting for their prescription are immortal by",
    "construction (dying first would have put them in the untreated group), and the",
    "wrong analysis credits that time to the treated group.\n")

cat("Q2: right analysis: unexposed rate", round(rates["unexp_rate"], 3), "vs exposed rate",
    round(rates["exp_rate"], 3), "per person-year; both sit at the simulated 1/3 =",
    round(1/3, 3), ". Rate ratio =", round(rates["exp_rate"] / rates["unexp_rate"], 2),
    "(the naive comparison of mean survival made the treated look",
    round(mean_surv_wrong["TRUE"] / mean_surv_wrong["FALSE"], 1), "times better).",
    "Time-dependent classification removes the artefact.\n")

cat("Q3: check first how exposure and person-time are aligned at time zero: ever-use",
    "classification, exposure defined from prescriptions filled during follow-up while",
    "follow-up starts at cohort entry, no new-user design, no risk-set sampling.",
    "If patients had to survive to fill a prescription before being labelled exposed,",
    "the reported mortality reduction contains immortal time.\n")

# ==== Section 4: Final Assessment / Activity: Simulating Immortal Time Bias ====

set.seed(230)
N <- 2000

# Each patient survives Exp(rate=1/3) years; about 60% eventually fill a Rx.
# If they die before filling, they end up classified as unexposed.
death_t <- rexp(N, rate = 1/3)
rx_t    <- rexp(N, rate = 1/2)
ever_rx <- rx_t < death_t

## (1) WRONG analysis: classify by ever_rx, ignore immortal time
mean_surv_wrong <- tapply(death_t, ever_rx, mean)
mean_surv_wrong   # exposed group looks much better

## (2) RIGHT analysis: split person-time at the prescription date
unexp_pt <- pmin(death_t, rx_t)              # unexposed person-time
exp_pt   <- pmax(death_t - rx_t, 0) * ever_rx # exposed person-time

unexp_events <- sum(!ever_rx)
exp_events   <- sum( ever_rx)
c(unexp_rate = unexp_events / sum(unexp_pt),
  exp_rate   = exp_events   / sum(exp_pt))