# =============================================================================
# HSCI 410 Public Health Assessment and Analysis  -  Lesson 5: Survival Data
# Answer key for the in-lesson R activities
# Data file(s): phaa_followup.csv  (download from the lesson page; save in your R working directory)
# Packages: survival, survminer  (install once with install.packages(c("survival", "survminer")))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================
# phaa_followup.csv follows the 800 survey participants for up to 10 years:
#   fu_years      years observed (to the event, loss to follow-up, or 10 years)
#   cv_event      1 = cardiovascular event, 0 = censored
#   cv_event_time identical to fu_years (time to event or censoring)
# plus age, gender, smoker, bmi, systolic_bp and hypertension at baseline.
# If survminer will not install, every plot below has a base-R equivalent in
# the comments; the model output does not depend on survminer.

# ==== Section 3: Cox Proportional Hazards Model / Activity: Kaplan-Meier curves and a Cox model ====
library(survival);  library(survminer)
phaa <- read.csv("phaa_followup.csv", stringsAsFactors = FALSE,
                 na.strings = c("", "NA"))
phaa$smoker <- factor(phaa$smoker, levels = c("No","Yes"))

# How many events? Always check before fitting.
table(phaa$cv_event)          # 648 censored, 152 events
mean(phaa$cv_event)           # crude cumulative incidence = 0.19

# 1. Build the Surv object: time + event indicator
y <- Surv(time = phaa$fu_years, event = phaa$cv_event)
head(y, 10)                   # "3.2+" is censored at 3.2 years, "5.7" is an event

# 2. Kaplan-Meier: event-free probability over time, by smoking status
km <- survfit(y ~ smoker, data = phaa)
km_all <- survfit(y ~ 1, data = phaa)   # overall event-free curve (ignores smoking)
summary(km_all, times = c(1, 3, 5, 7, 10))   # event-free probability at landmark years
summary(km, times = c(1, 3, 5, 7, 10))       # the same, separately for non-smokers and smokers
ggsurvplot(km, data = phaa, conf.int = TRUE,
           pval = TRUE, risk.table = TRUE,
           xlab = "Years of follow-up")
# Base-R fallback if survminer is unavailable:
# plot(km, col = c("steelblue", "firebrick"), conf.int = TRUE,
#      xlab = "Years of follow-up", ylab = "Event-free probability")
# legend("bottomleft", c("Non-smoker", "Smoker"),
#        col = c("steelblue", "firebrick"), lty = 1)

# 3. Log-rank test: do the curves differ overall?
survdiff(y ~ smoker, data = phaa)

# 4. Cox model with multiple covariates -> hazard ratios
cox <- coxph(y ~ smoker + age + gender + bmi + hypertension,
             data = phaa)
summary(cox)

# 5. Test the proportional-hazards assumption
cox.zph(cox)
plot(cox.zph(cox))            # flat smooths support proportional hazards

# ---- Activity questions (box rActivity-410-8-1) -----------------------------
s_all <- summary(km_all, times = 5)
s_smk <- summary(km, times = 5)
cat("\nQ1: overall event-free probability at 5 years =", round(s_all$surv, 3),
    " (95% CI", round(s_all$lower, 3), "to", round(s_all$upper, 3), ")\n")
cat("    at 5 years: non-smokers", round(s_smk$surv[s_smk$strata == "smoker=No"], 3),
    " smokers", round(s_smk$surv[s_smk$strata == "smoker=Yes"], 3), "\n")
s10 <- summary(km, times = 10)
cat("    at 10 years: non-smokers", round(s10$surv[s10$strata == "smoker=No"], 3),
    " smokers", round(s10$surv[s10$strata == "smoker=Yes"], 3), "\n")
cat("    The smokers' curve drops faster; the gap is small in the first years and\n",
    "   widens to about 12 percentage points by year 10.\n")

lr <- survdiff(y ~ smoker, data = phaa)
cs <- summary(cox)
cat("\nQ2: log-rank chi-square =", round(lr$chisq, 2), "on 1 df, p =", signif(lr$pvalue, 2), "\n")
cat("    adjusted HR for smokerYes =", round(cs$conf.int["smokerYes", "exp(coef)"], 2),
    " 95% CI", round(cs$conf.int["smokerYes", "lower .95"], 2), "to",
    round(cs$conf.int["smokerYes", "upper .95"], 2),
    " (p =", signif(cs$coefficients["smokerYes", "Pr(>|z|)"], 2), ")\n")
cat("    At any moment during follow-up, a smoker's instantaneous risk of a\n",
    "   cardiovascular event is about 78% higher than a non-smoker's, holding age,\n",
    "   gender, BMI and hypertension constant. Age is the other clear predictor\n",
    "   (HR 1.05 per year); gender, BMI and hypertension are not significant here.\n")

zp <- cox.zph(cox)
cat("\nQ3: cox.zph global p =", round(zp$table["GLOBAL", "p"], 2), "\n")
print(round(zp$table, 3))
cat("    No predictor has p < .05 (the smallest is hypertension, p = 0.11), so the\n",
    "   proportional-hazards assumption is defensible for every term and a single\n",
    "   HR per predictor is an adequate summary. A small p-value would mean that\n",
    "   predictor's HR drifts over follow-up (for example strong early, weaker\n",
    "   later), and the fix is a time interaction, a stratified baseline hazard,\n",
    "   or a time-varying coefficient model.\n")

# ---- Stretch (not required): adjusted survival curves from the Cox model ----
new_dat <- data.frame(
  smoker       = factor(c("No", "Yes"), levels = c("No", "Yes")),
  age          = mean(phaa$age),
  gender       = "Woman",
  bmi          = mean(phaa$bmi),
  hypertension = "No"
)
plot(survfit(cox, newdata = new_dat),
     col = c("steelblue", "firebrick"),
     xlab = "Years", ylab = "Adjusted event-free probability",
     main = "Cox-adjusted curves: smoker vs non-smoker")
legend("bottomleft", c("Non-smoker", "Smoker"),
       col = c("steelblue", "firebrick"), lty = 1)