# =============================================================================
# HSCI 410 Public Health Assessment and Analysis  -  Lesson 1: A Structured
# Approach to Data Analysis
# Answer key for the in-lesson R activities
# Data file(s): cohort.csv, bp01.csv  (download from the lesson page; save in
#               your R working directory. The script copies each file into the
#               data/raw/ folder that the lesson's project skeleton creates.)
# Packages: mediation, tidyverse (dplyr, readr, tidyr, ggplot2), here
#           (install once with install.packages(c("mediation", "tidyverse", "here")))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================


# ==== Section 1: Introduction & Data Collection / Fitting a mediation model in R (Baron & Kenny + the mediation package) ====

# install.packages(c("mediation", "dagitty"))
library(mediation)

# 1. Simulate data that match the DAG: education -> income -> health,
#    with a smaller direct path education -> health.
set.seed(410)
n         <- 800
education <- rnorm(n)
income    <- 0.6 * education + rnorm(n)              # path a
health    <- 0.3 * education + 0.5 * income + rnorm(n)  # direct + path b
dat       <- data.frame(education, income, health)

# 2. Baron & Kenny by hand --------------------------------------------------
#    Step 1: total effect c  (health on education)
coef(lm(health ~ education, data = dat))["education"]

#    Step 2: a (income on education)
fit_M  <- lm(income ~ education, data = dat)

#    Step 3: direct c' (education) and b (income), from health on both
fit_Y  <- lm(health ~ education + income, data = dat)
coef(fit_Y)                     # c' on education, b on income

# Indirect effect = a * b  (or equivalently c - c')
a <- coef(fit_M)["education"]
b <- coef(fit_Y)["income"]
a * b

# 3. Same answer, with bootstrap CIs, via the mediation package -------------
med <- mediate(fit_M, fit_Y,
                treat    = "education",
                mediator = "income",
                boot     = TRUE, sims = 1000)
summary(med)

# ---- Activity questions, Section 1 ------------------------------------------
c_total  <- unname(coef(lm(health ~ education, data = dat))["education"])
c_direct <- unname(coef(fit_Y)["education"])
cat("\nQ1: total effect c =", round(c_total, 3),
    "; direct effect c' =", round(c_direct, 3),
    "; difference c - c' =", round(c_total - c_direct, 3), "\n")
cat("    The total effect is the larger one. The gap between c and c' is the\n",
    "   indirect effect that travels through income (about half of the total).\n")

cat("\nQ2: a =", round(unname(a), 3), "; b =", round(unname(b), 3),
    "; a * b =", round(unname(a * b), 3), "\n")
cat("    Bootstrapped ACME =", round(med$d0, 3),
    " 95% CI (", round(med$d0.ci[1], 3), ",", round(med$d0.ci[2], 3), ")\n")
cat("    The product a*b equals the ACME to three decimals (the point estimate is\n",
    "   the same quantity; the bootstrap only adds the CI), and the CI excludes 0.\n")

cat("\nQ3: Prop. Mediated =", round(med$n0, 3),
    " (95% CI", round(med$n0.ci[1], 3), ",", round(med$n0.ci[2], 3), ")\n")
cat("    About", round(100 * med$n0), "% of education's total effect on health is\n",
    "   carried by income. If the ACME CI crossed zero we could not rule out that\n",
    "   income contributes nothing to the education-health link in these data.\n")


# ==== Section 1: Introduction & Data Collection / A reproducible project skeleton in RStudio ====

# NOTE (answer key): the lesson reads data/raw/cohort.csv. Download cohort.csv
# from the lesson page into your working directory; the two lines after the
# dir.create() calls copy it into the data/raw/ folder that the skeleton
# creates, so the pipeline below runs exactly as printed on the page.

# Create directories from R (or by hand). Run once at project start.
dir.create("data/raw",        recursive = TRUE)
dir.create("data/processed",  recursive = TRUE)
dir.create("R");  dir.create("output/figures", recursive = TRUE)

# NOTE (answer key): place the downloaded raw file where the skeleton expects it
file.copy("cohort.csv", "data/raw/cohort.csv", overwrite = TRUE)

# tidyverse: dplyr (manipulation), ggplot2 (graphics), readr (file IO),
# tidyr (reshape), stringr (text). Install once.
# install.packages("tidyverse")
library(tidyverse)
library(here)                                # robust file paths

# A canonical pipeline: read -> clean -> save -> analyse
raw <- read_csv(here("data/raw/cohort.csv"))
clean <- raw |>
  filter(!is.na(outcome)) |>
  mutate(age_grp = cut(age, c(0, 30, 50, 70, Inf)),
         smoker  = factor(smoker, levels = c("No", "Yes")))
write_csv(clean, here("data/processed/cohort_clean.csv"))

# Sketch a DAG to anchor the analysis (see the earlier DAG course)
# library(dagitty)
# g <- dagitty("dag { smoker -> outcome ; age -> smoker ; age -> outcome }")

# ---- Activity questions, Section 2 ------------------------------------------
cat("\nQ1: folders created by the four dir.create() calls:\n")
print(list.dirs(c("data", "R", "output"), recursive = TRUE))
cat("    data/raw/ holds the irreplaceable original; data/processed/ holds files\n",
    "   that can be rebuilt from raw + code, so a cleaning mistake is never fatal.\n")

cat("\nQ2: raw has", nrow(raw), "rows and", ncol(raw), "columns; clean has",
    nrow(clean), "rows and", ncol(clean), "columns\n")
cat("    filter() dropped", nrow(raw) - nrow(clean), "rows with a missing outcome\n")
cat("    brand-new column(s):", setdiff(names(clean), names(raw)), "\n")
print(table(clean$age_grp))
cat("    re-encoded in place: smoker (character -> factor with levels",
    paste(levels(clean$smoker), collapse = " < "), ")\n")

cat("\nQ3: here() resolved the raw file to:", here("data/raw/cohort.csv"), "\n")
cat("    The path is built from the project root, so the same script runs on any\n",
    "   machine, user account or operating system once the project is unzipped.\n")


# ==== Section 2: Data Coding, Entry & File Management / A reproducible recoding pipeline (no overwrites, ever) ====

# NOTE (answer key): download bp01.csv from the lesson page into your working
# directory; this line copies it into data/raw/ (created above).
file.copy("bp01.csv", "data/raw/bp01.csv", overwrite = TRUE)

library(tidyverse)

# Read raw, never modify in place
bp_raw <- read_csv("data/raw/bp01.csv")

# Tidy: drop incomplete rows, build derived variables, lock factors
bp_clean <- bp_raw |>
  drop_na(systolic, diastolic, age) |>
  mutate(
    age_ct    = age - mean(age),                              # centred
    age_ctsq  = age_ct^2,                                       # quadratic term
    age_c3    = cut(age, c(0, 35, 55, Inf),
                    labels = c("young", "middle", "older")),
    htn       = factor(systolic >= 140 | diastolic >= 90,
                       levels = c(FALSE, TRUE),
                       labels = c("normotensive", "hypertensive"))
  )

# Persist as a new versioned file - and a small log line
write_csv(bp_clean, "data/processed/bp02.csv")
cat("bp02.csv", format(Sys.Date()), nrow(bp_clean), "obs",
    "\n", file = "data/file_log.txt", append = TRUE)

## At any time you can rebuild bp02 from bp01 by re-running this script.

# ---- Activity questions, Section 3 ------------------------------------------
cat("\nQ1: bp01 had", nrow(bp_raw), "obs and", ncol(bp_raw), "vars; bp02 has",
    nrow(bp_clean), "obs and", ncol(bp_clean), "vars (",
    nrow(bp_raw) - nrow(bp_clean), "incomplete rows dropped )\n")
cat("    The four new variables and their classes:\n")
str(bp_clean[, c("age_ct", "age_ctsq", "age_c3", "htn")])
cat("    age_ct: continuous (age minus the mean age of", round(mean(bp_clean$age), 1), ")\n",
    "   age_ctsq: derived continuous (age_ct squared, for a curved age effect)\n",
    "   age_c3: categorical factor with three age bands\n",
    "   htn: derived binary factor (140/90 rule)\n")
print(table(bp_clean$age_c3))
print(table(bp_clean$htn))

prev_140 <- mean(bp_clean$systolic >= 140 | bp_clean$diastolic >= 90)
prev_150 <- mean(bp_clean$systolic >= 150 | bp_clean$diastolic >= 90)
cat("\nQ2: prevalence of 'hypertensive' with systolic cutoff 140:",
    round(100 * prev_140, 1), "%\n")
cat("    prevalence with systolic cutoff 150:", round(100 * prev_150, 1), "%\n")
cat("    Raising the cutoff lowers the prevalence by",
    round(100 * (prev_140 - prev_150), 1), "percentage points; a categorical\n",
    "   recode is only as defensible as its threshold.\n")

cat("\nQ3: files now on disk:\n")
cat("    data/raw/bp01.csv exists:", file.exists("data/raw/bp01.csv"),
    "; data/processed/bp02.csv exists:", file.exists("data/processed/bp02.csv"), "\n")
cat("    file log contents:\n")
cat(readLines("data/file_log.txt"), sep = "\n")
cat("    Because bp01.csv is untouched, a recoding mistake in bp02 (for example\n",
    "   reversed htn labels) is repaired by fixing the script and re-running it;\n",
    "   an overwritten SPSS working file cannot be re-derived.\n")