# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 1: Foundations of Epidemiology
# Answer key for the in-lesson R activities
# Data file(s): none  (every dataset in this lesson is typed directly into R)
# Packages: none beyond base R
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 1: History of Epidemiology / Meet R: re-create Graunt's life table in modern code ====

# Graunt's 1662 life table: of 100 people born, how many survive to each age?
# These are Graunt's original (rough) figures, written into two vectors.
age       <- c(0,   6,  16, 26, 36, 46, 56, 66, 76, 80)
survivors <- c(100, 64, 40, 25, 16, 10,  6,  3,  1,  0)

# Bundle them into a data.frame, R's table-like object.
graunt <- data.frame(age = age, survivors = survivors)

# Probability of surviving from birth to each age:
graunt$survival_prob <- graunt$survivors / 100

# Print and plot the life table:
print(graunt)
plot(graunt$age, graunt$survival_prob,
     type = "b", pch = 19,
     xlab = "Age (years)",
     ylab = "Probability of surviving from birth",
     main = "Graunt's 1662 Life Table for London")

# Stretch challenge: conditional survival, the probability of surviving each
# age interval given that you reached the start of it (the building block of
# a survival curve).
graunt$cond_survival <- c(NA, graunt$survivors[-1] / graunt$survivors[-nrow(graunt)])
print(graunt)

# ---- Activity questions (rActivity-230-1-1) ----
cat("Q1: survival_prob at age 16 =", graunt$survival_prob[graunt$age == 16],
    "; at age 26 =", graunt$survival_prob[graunt$age == 26], "\n")
cat("    ", graunt$survivors[graunt$age == 16] - graunt$survivors[graunt$age == 26],
    "of every 100 births died between 16 and 26, but",
    100 - graunt$survivors[graunt$age == 16], "of 100 were already dead by 16:",
    "mortality was front-loaded into childhood.\n")

drops <- data.frame(from = age[-length(age)], to = age[-1],
                    drop = -diff(graunt$survival_prob))
cat("Q2: absolute drop in survival probability per interval:\n"); print(drops)
cat("    Steepest drop:", drops$from[which.max(drops$drop)], "to",
    drops$to[which.max(drops$drop)], "(", max(drops$drop),
    "); the first years of life were the riskiest.\n")

cs <- graunt$cond_survival[-1]
cat("Q3: conditional survival by interval:\n")
print(round(setNames(cs, paste0(age[-length(age)], "-", age[-1])), 3))
cat("    Lowest overall: 76-80 (", cs[length(cs)], "); among the full-length",
    "intervals 66-76 is lowest at", round(cs[8], 3), "(one chance in three).\n")

# ==== Section 1: History of Epidemiology / Try it: visualise the decline in global child mortality ====

# Approximate global under-5 mortality rate (% of children dying before age 5)
# Source: Gapminder / UN IGME estimates, rounded for illustration.
year    <- c(1800, 1900, 1950, 1980, 2000, 2020)
u5_pct  <- c(43.0, 36.0, 22.5, 12.0,  7.6,  3.7)

# A quick visual: black line + red points.
plot(year, u5_pct, type = "l", lwd = 2,
     ylim = c(0, 50),
     xlab = "Year",
     ylab = "Under-5 mortality (%)",
     main = "Global Child Mortality, 1800-2020")
points(year, u5_pct, pch = 19, col = "firebrick")

# Average annual reduction (rough): how many percentage points per decade?
total_drop_pct <- u5_pct[1] - u5_pct[length(u5_pct)]
years_span     <- max(year) - min(year)
cat("Drop of", total_drop_pct, "percentage points over", years_span, "years.\n")

# ---- Activity questions (rActivity-230-1-2) ----
cat("Q1: X =", total_drop_pct, "percentage points, Y =", years_span, "years;",
    "average drop per decade =", round(total_drop_pct / (years_span / 10), 2), "pp\n")

seg <- data.frame(from = year[-length(year)], to = year[-1],
                  drop_pp = -diff(u5_pct),
                  pp_per_decade = round(-diff(u5_pct) / diff(year) * 10, 2))
cat("Q2: fall between consecutive data points:\n"); print(seg)
cat("    Largest absolute fall:", seg$from[which.max(seg$drop_pp)], "-",
    seg$to[which.max(seg$drop_pp)], "(", max(seg$drop_pp), "pp );",
    "fastest rate of decline:", seg$from[which.max(seg$pp_per_decade)], "-",
    seg$to[which.max(seg$pp_per_decade)], "(", max(seg$pp_per_decade),
    "pp per decade ). The decline is not linear.\n")

u5_alt <- replace(u5_pct, length(u5_pct), 1.0)
cat("Q3: with 2020 = 1.0 the drop becomes", u5_alt[1] - u5_alt[length(u5_alt)],
    "pp (was", total_drop_pct, "): the headline gets more dramatic.\n")

# ==== Section 2: Ways of Knowing / Quantitative thinking in practice: a quick descriptive summary ====

# A toy dataset: 10 adults' systolic blood pressure (mmHg) and smoking status
sbp     <- c(118, 132, 145, 128, 155, 120, 142, 138, 125, 160)
smoker  <- c("no", "yes", "yes", "no", "yes", "no", "yes", "no", "no", "yes")

# Centre & spread, the standard quantitative summary.
mean(sbp)            # average
median(sbp)          # middle value
sd(sbp)              # standard deviation
range(sbp)           # minimum and maximum

# Group comparison: does mean SBP differ by smoking status?
tapply(sbp, smoker, mean)

# ---- Activity questions (rActivity-230-1-3) ----
cat("Q1: mean =", mean(sbp), ", median =", median(sbp), ", SD =", round(sd(sbp), 1),
    "mmHg; mean and median within ~1 mmHg -> roughly symmetric distribution.\n")
grp <- tapply(sbp, smoker, mean)
cat("Q2: smokers", grp["yes"], "vs non-smokers", grp["no"], "mmHg; gap =",
    grp["yes"] - grp["no"], "mmHg (large in clinical terms).\n")
cat("Q3: n =", table(smoker)["yes"], "per group, observational, no confounders measured:",
    "the gap is not causal evidence; the numbers cannot say why these people smoke,",
    "what else differs between the groups, or when the reading was taken.\n")

# ==== Section 4: Final Assessment / Activity, Foundations in code: history, ways of knowing, and reproducibility ====

# PART A -- Graunt's 1662 life table for London (same vectors as Section 1)
graunt <- data.frame(age = age, survivors = survivors)
graunt$survival_prob <- graunt$survivors / 100
print(graunt)

plot(graunt$age, graunt$survival_prob,
     type = "b", pch = 19,
     xlab = "Age (years)",
     ylab = "Probability of surviving from birth",
     main = "Graunt's 1662 Life Table for London")

# PART B -- ways of knowing: a quantitative summary of SBP by smoking
mean(sbp)     # average
median(sbp)   # middle value
sd(sbp)       # standard deviation
range(sbp)    # minimum and maximum

tapply(sbp, smoker, mean)            # mean SBP by smoking status

# PART C -- reproducibility: set.seed() + a bootstrap 95% CI for the mean
set.seed(230)                                # lock in random draws

x <- rnorm(100, mean = 10, sd = 2)         # Normal(10, 2)
boot_means <- replicate(1000, mean(sample(x, replace = TRUE)))

quantile(boot_means, c(0.025, 0.975))     # 95% bootstrap CI
sessionInfo()                                  # record package versions

cat("Part C check: sample mean =", round(mean(x), 3), "; bootstrap 95% CI =",
    paste(round(quantile(boot_means, c(0.025, 0.975)), 3), collapse = " to "),
    "(identical on every run because of set.seed(230)).\n")