# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 2: Systematic Reviews and Meta-Analysis
# Answer key for the in-lesson R activities
# Data file(s): none  (the eight trial 2x2 counts are typed directly into R)
# Packages: metafor  (install once with install.packages("metafor"))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 3: Forest Plots & Heterogeneity / Forest plot & meta-analysis with the metafor package ====

# install.packages("metafor")
library(metafor)

# Eight hypothetical RCTs of a smoking-cessation drug (events / total per arm)
dat <- data.frame(study  = paste("Trial", 1:8),
  ai = c(48, 54, 31, 84, 21, 64, 36, 15),  # events on drug
  n1i = c(120,160,100,220, 85,240,130,90),
  ci = c(30, 32, 21, 48, 18, 60, 28, 14),
  n2i = c(120,160,100,220, 85,240,130,90)
)

# Compute log-OR per study with continuity correction
es <- escalc(measure = "OR",
             ai = ai, n1i = n1i,
             ci = ci, n2i = n2i, data = dat)
es

# Random-effects pooled OR (DerSimonian-Laird estimator)
fit <- rma(yi, vi, data = es, method = "DL")
fit

# Forest plot
forest(fit, slab = es$study, transf = exp, refline = 1,
       xlab = "Odds ratio (smoking cessation)")

# ---- Activity questions (rActivity-230-2-1) ----
pooled <- exp(c(est = unname(coef(fit)), lower = fit$ci.lb, upper = fit$ci.ub))
cat("Q1: pooled log-OR =", round(coef(fit), 2), "-> OR =", round(pooled["est"], 2),
    "(95% CI", round(pooled["lower"], 2), "to", round(pooled["upper"], 2), ").",
    "Odds of quitting are about", round(pooled["est"], 1),
    "times higher on the drug; cessation is the desired outcome, so OR > 1 is beneficial.\n")

cat("Q2: I^2 =", round(fit$I2, 1), "%, tau^2 =", round(fit$tau2, 4),
    "; Q =", round(fit$QE, 2), "on", fit$k - 1, "df, p =", round(fit$QEp, 3),
    "-> low heterogeneity by the 25/50/75 benchmarks, and the Q-test agrees (not significant).\n")

# Q3: fixed-effect model for comparison (the diamond narrows)
fit_fe <- rma(yi, vi, data = es, method = "FE")
fit_fe
forest(fit_fe, slab = es$study, transf = exp, refline = 1,
       xlab = "Odds ratio (smoking cessation), fixed-effect model")
pooled_fe <- exp(c(est = unname(coef(fit_fe)), lower = fit_fe$ci.lb, upper = fit_fe$ci.ub))
cat("Q3: fixed-effect OR =", round(pooled_fe["est"], 2), "(95% CI",
    round(pooled_fe["lower"], 2), "to", round(pooled_fe["upper"], 2),
    "); random-effects CI width on the log scale =", round(fit$ci.ub - fit$ci.lb, 3),
    "vs fixed-effect", round(fit_fe$ci.ub - fit_fe$ci.lb, 3),
    ". FE pools only within-study variance; RE adds tau^2 to every weight, so its CI is wider.",
    "The narrower FE CI misleads whenever the trials do not share one true effect.\n")