# =============================================================================
# HSCI 230: Evaluating Epidemiological Research  -  Lesson 12: Integrated Appraisal of Epidemiological Research
# Answer key for the in-lesson R activities
# Data file(s): none  (the three published ORs and CIs are typed directly into R)
# Packages: metafor  (cross-check only; install once with install.packages("metafor"))
# Reproduces every code block in the lesson, then answers each activity question.
# =============================================================================

# ==== Section 3: Red Flags, Quality Indicators, and Applied Synthesis / Activity: Pooling published odds ratios with inverse-variance meta-analysis ====

# Three hypothetical studies of exposure E and outcome Y
study  <- c("Study A", "Study B", "Study C")
or     <- c(1.40, 1.10, 1.55)
lci    <- c(1.05, 0.85, 1.10)
uci    <- c(1.85, 1.43, 2.18)

# SE of log-OR derived from the published 95% CI
log_or <- log(or)
log_se <- (log(uci) - log(lci)) / (2 * 1.96)
w      <- 1 / log_se^2                         # inverse-variance weights

# Pooled log-OR and 95% CI
pool_log <- sum(w * log_or) / sum(w)
pool_se  <- sqrt(1 / sum(w))
pool_or  <- exp(pool_log)
pool_ci  <- exp(pool_log + c(-1, 1) * 1.96 * pool_se)
round(c(pooled_OR = pool_or, lower = pool_ci[1], upper = pool_ci[2]), 2)

# Cross-check against metafor::rma() (fixed-effect, inverse-variance)
library(metafor)
fit_fe <- rma(yi = log_or, sei = log_se, method = "FE", slab = study)
fit_fe
round(exp(c(coef(fit_fe), fit_fe$ci.lb, fit_fe$ci.ub)), 2)   # same 1.30 (1.10 to 1.53)
forest(fit_fe, transf = exp, refline = 1, xlab = "Odds ratio")

# ---- Activity questions (rActivity-230-12-1) ----
wt <- data.frame(study, OR = or, CI_ratio = round(uci / lci, 2), log_se = round(log_se, 3),
                 weight = round(w, 1), share = round(100 * w / sum(w), 1))
cat("Q1: inverse-variance weights and their share of the pooled estimate:\n"); print(wt)
cat("   ", study[which.max(w)], "carries the most weight: its CI is the narrowest",
    "(upper/lower ratio", min(round(uci / lci, 2)), "), and weight = 1/SE^2.\n")

cat("Q2: pooled 95% CI", round(pool_ci[1], 2), "to", round(pool_ci[2], 2),
    "(ratio", round(pool_ci[2] / pool_ci[1], 2), ") vs individual CI ratios",
    paste(round(uci / lci, 2), collapse = ", "), ". Pooled variance 1/sum(w) =",
    round(1 / sum(w), 4), "is smaller than any single study variance (",
    paste(round(log_se^2, 4), collapse = ", "), "). Fair only if the studies estimate one",
    "common effect (fixed-effect assumption). Heterogeneity check: Q =", round(fit_fe$QE, 2),
    "on 2 df, p =", round(fit_fe$QEp, 3), ", I^2 =", round(fit_fe$I2, 1), "%.\n")

fit_re <- rma(yi = log_or, sei = log_se, method = "REML", slab = study)
cat("    Random-effects (REML) pooled OR =", round(exp(coef(fit_re)), 2), "(",
    round(exp(fit_re$ci.lb), 2), "to", round(exp(fit_re$ci.ub), 2), "), tau^2 =",
    round(fit_re$tau2, 4), "\n")

# Q3: pooling propagates precision, not validity. Sensitivity analysis: drop Study A
# (cross-sectional, weak confounder control) and pool the remaining two.
keep <- study != "Study A"
pool_log_bc <- sum(w[keep] * log_or[keep]) / sum(w[keep])
pool_se_bc  <- sqrt(1 / sum(w[keep]))
cat("Q3: without Study A the pooled OR is", round(exp(pool_log_bc), 2), "(",
    round(exp(pool_log_bc - 1.96 * pool_se_bc), 2), "to",
    round(exp(pool_log_bc + 1.96 * pool_se_bc), 2), "). A precise pooled CI says nothing about",
    "bias in the inputs: report pooled estimates by design and risk-of-bias category,",
    "add sensitivity analyses excluding high-risk studies, and grade certainty (ROBINS-I /",
    "GRADE) rather than relying on the CI alone.\n")