Relative Risk (RR) Test

# Relative Risk: Smoking and Lung Cancer Incidence (Cohort Study)
# Prospective design with incident cases over 20-year follow-up

library(epitools)     # For riskratio()
library(ggplot2)
library(survival)     # For survival curves
library(survminer)    # For ggsurvplot

# Set seed for reproducibility
set.seed(2025)

# === Simulate prospective cohort data ===
# 20-year follow-up, RR ≈ 15 (strong effect)
# Incidence: 15% in smokers, 1% in never-smokers

n_smokers <- 2500
n_never <- 2500

# Incident lung cancer during 20-year follow-up
incidence_smokers <- 0.15    # 15% cumulative incidence
incidence_never <- 0.01      # 1% cumulative incidence

# Generate outcome data
smokers_cancer <- rbinom(1, n_smokers, incidence_smokers)
smokers_no_cancer <- n_smokers - smokers_cancer

never_cancer <- rbinom(1, n_never, incidence_never)
never_no_cancer <- n_never - never_cancer

# For realistic simulation, use typical values
smokers_cancer <- round(n_smokers * incidence_smokers)  # ~375 cases
never_cancer <- round(n_never * incidence_never)        # ~25 cases
smokers_no_cancer <- n_smokers - smokers_cancer
never_no_cancer <- n_never - never_cancer

# Create 2x2 table
# Rows: Exposure (Smoker, Never-smoker)
# Columns: Outcome (Cancer, No Cancer)
table_cohort <- matrix(c(
  smokers_cancer, smokers_no_cancer,
  never_cancer, never_no_cancer
), nrow=2, byrow=TRUE,
dimnames=list(
  Exposure = c("Smoker", "Never-smoker"),
  Outcome = c("Cancer", "No Cancer")
))

print("=== 2x2 Contingency Table(Cohort Data) ===")
print(table_cohort)
print(addmargins(table_cohort))

# === STEP 1: Check Assumptions ===

cat("\n=== Assumption Checks ===\n")

# 1. Study design
cat("Design: Prospective cohort study ✓\n")
cat("Sampling: By exposure status(smokers/never-smokers) ✓\n")
cat("Follow-up: 20 years(incident cases only) ✓\n")

# 2. Adequate events
min_events <- min(smokers_cancer, never_cancer)
cat("\nMinimum events per exposure group:", min_events, 
    ifelse(min_events >= 5, "✓", "✗ Consider exact methods"), "\n")

# 3. Temporal sequence
cat("Temporal sequence: Exposure(baseline smoking status) precedes outcome(incident cancer) ✓\n")

# === STEP 2: Calculate Cumulative Incidence (Risk) ===

cat("\n=== Cumulative Incidence(Risk) ===\n")

risk_smokers <- smokers_cancer / n_smokers
risk_never <- never_cancer / n_never

cat("Risk in smokers:", round(risk_smokers, 4), "(", 
    round(risk_smokers * 100, 2), "%)\n")
cat("Risk in never-smokers:", round(risk_never, 4), "(", 
    round(risk_never * 100, 2), "%)\n")

# === STEP 3: Calculate Relative Risk with 95% CI ===

cat("\n=== Relative Risk Calculation ===\n")

# Method 1: Manual calculation
RR_manual <- risk_smokers / risk_never
cat("RR = Risk(smokers) / Risk(never-smokers)\n")
cat("RR =", round(risk_smokers, 4), "/", round(risk_never, 4), "=", round(RR_manual, 2), "\n")

# Calculate 95% CI using log method
log_RR <- log(RR_manual)
SE_log_RR <- sqrt((1 - risk_smokers) / smokers_cancer + (1 - risk_never) / never_cancer)
CI_lower <- exp(log_RR - 1.96 * SE_log_RR)
CI_upper <- exp(log_RR + 1.96 * SE_log_RR)

cat("95% CI:", round(CI_lower, 2), "-", round(CI_upper, 2), "\n")

# Method 2: Using epitools package
rr_result <- riskratio(table_cohort, method="wald")
print(rr_result)

# === STEP 4: Calculate Risk Difference (Attributable Risk) ===

cat("\n=== Risk Difference(Attributable Risk) ===\n")

RD <- risk_smokers - risk_never
cat("RD = Risk(smokers) - Risk(never-smokers)\n")
cat("RD =", round(RD, 4), "(", round(RD * 100, 2), "percentage points)\n")

# SE for risk difference
SE_RD <- sqrt(risk_smokers * (1 - risk_smokers) / n_smokers + 
              risk_never * (1 - risk_never) / n_never)
RD_CI_lower <- RD - 1.96 * SE_RD
RD_CI_upper <- RD + 1.96 * SE_RD

cat("95% CI for RD:", round(RD_CI_lower, 4), "-", round(RD_CI_upper, 4), "\n")
cat("\nInterpretation:", round(RD * 100, 2), 
    "% excess lung cancer risk in smokers attributable to smoking\n")

# === STEP 5: Hypothesis Test ===

cat("\n=== Hypothesis Test ===\n")
cat("H₀: RR = 1 (no association)\n")
cat("Hₐ: RR ≠ 1 (association exists)\n\n")

if (CI_lower > 1) {
  cat("Result: Reject H₀. RR significantly > 1 (p < .05)\n")
  cat("Smoking is associated with INCREASED risk of lung cancer\n")
} else if (CI_upper < 1) {
  cat("Result: Reject H₀. RR significantly < 1 (p < .05)\n")
  cat("Exposure associated with DECREASED risk\n")
} else {
  cat("Result: Fail to reject H₀. 95% CI includes 1.0\n")
  cat("No significant association detected\n")
}

# Chi-square test
chi2_result <- chisq.test(table_cohort, correct=FALSE)
cat("\nChi-square test: χ² =", round(chi2_result$statistic, 2), 
    ", p =", format.pval(chi2_result$p.value, digits=3), "\n")

# === STEP 6: Visualizations ===

# 6.1 Forest plot (RR with CI)
forest_data <- data.frame(
  Comparison = "Smokers vs Never-Smokers",
  RR = RR_manual,
  CI_lower = CI_lower,
  CI_upper = CI_upper
)

ggplot(forest_data, aes(y=Comparison, x=RR)) +
  geom_point(size=5, color="darkred") +
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), height=0.2, linewidth=1.5) +
  geom_vline(xintercept=1, linetype="dashed", color="blue", linewidth=1.2) +
  scale_x_log10(breaks=c(0.5, 1, 2, 5, 10, 20, 30)) +
  labs(title="Relative Risk: Smoking and Lung Cancer Incidence",
       subtitle="20-year Prospective Cohort Study(n=5,000)",
       x="Relative Risk(log scale) with 95% CI",
       y="") +
  theme_classic(base_size=14) +
  theme(plot.title = element_text(hjust=0.5, face="bold"),
        plot.subtitle = element_text(hjust=0.5)) +
  annotate("text", x=RR_manual, y=1.3, 
           label=paste0("RR = ", round(RR_manual, 2), 
                       "\n95% CI: [", round(CI_lower, 2), ", ", round(CI_upper, 2), "]"),
           size=5, fontface="bold")

ggsave("forest_plot_rr.png", width=12, height=6)

# 6.2 Bar plot: Cumulative incidence by exposure
incidence_data <- data.frame(
  Group = c("Smokers", "Never-Smokers"),
  Incidence = c(risk_smokers * 100, risk_never * 100),
  SE = c(sqrt(risk_smokers * (1 - risk_smokers) / n_smokers) * 100,
         sqrt(risk_never * (1 - risk_never) / n_never) * 100)
)

ggplot(incidence_data, aes(x=Group, y=Incidence, fill=Group)) +
  geom_bar(stat="identity", width=0.6, alpha=0.8) +
  geom_errorbar(aes(ymin=Incidence - 1.96*SE, ymax=Incidence + 1.96*SE),
                width=0.2, linewidth=1) +
  geom_text(aes(label=paste0(round(Incidence, 2), "%")), 
            vjust=-2, size=6, fontface="bold") +
  scale_fill_manual(values=c("Smokers"="#E69F00", "Never-Smokers"="#56B4E9")) +
  labs(title="20-Year Cumulative Incidence of Lung Cancer",
       subtitle="Prospective Cohort Study",
       y="Cumulative Incidence(%) ± 95% CI",
       x="Smoking Status") +
  theme_classic(base_size=14) +
  theme(legend.position="none",
        plot.title = element_text(hjust=0.5, face="bold"),
        plot.subtitle = element_text(hjust=0.5)) +
  ylim(0, 20)

ggsave("incidence_by_exposure.png", width=10, height=8)

# 6.3 Simulate survival curves (for visualization)
# Generate individual-level data for survival analysis
time_smokers <- rexp(n_smokers, rate=-log(1 - incidence_smokers)/20)
time_never <- rexp(n_never, rate=-log(1 - incidence_never)/20)

# Cap at 20 years
time_smokers <- pmin(time_smokers, 20)
time_never <- pmin(time_never, 20)

# Event indicator (1 = cancer, 0 = censored)
event_smokers <- as.numeric(time_smokers < 20 & runif(n_smokers) < incidence_smokers)
event_never <- as.numeric(time_never < 20 & runif(n_never) < incidence_never)

# Combine into data frame
cohort_data <- data.frame(
  time = c(time_smokers, time_never),
  event = c(event_smokers, event_never),
  exposure = c(rep("Smoker", n_smokers), rep("Never-Smoker", n_never))
)

# Fit survival model
fit_surv <- survfit(Surv(time, event) ~ exposure, data=cohort_data)

# Plot survival curves (cancer-free survival)
ggsurvplot(fit_surv, data=cohort_data,
           risk.table=TRUE,
           pval=TRUE,
           conf.int=TRUE,
           title="Cancer-Free Survival by Smoking Status",
           xlab="Years of Follow-up",
           ylab="Probability of Remaining Cancer-Free",
           legend.title="Exposure",
           legend.labs=c("Never-Smoker", "Smoker"),
           palette=c("#56B4E9", "#E69F00"),
           ggtheme=theme_classic())

# === STEP 7: Calculate Population Impact ===

cat("\n=== Population Impact Measures ===\n")

# Attributable fraction among exposed (AFe)
AFe <- (RR_manual - 1) / RR_manual * 100
cat("Attributable fraction(exposed):", round(AFe, 1), "%\n")
cat("Interpretation:", round(AFe, 1), 
    "% of lung cancer in smokers is attributable to smoking\n")

# Population attributable fraction (PAF)
# Assume 50% prevalence of smoking in population
p_exposed_pop <- 0.50
PAF <- p_exposed_pop * (RR_manual - 1) / (p_exposed_pop * (RR_manual - 1) + 1) * 100
cat("\nPopulation attributable fraction(assuming 50% smoking prevalence):", 
    round(PAF, 1), "%\n")
cat("Interpretation:", round(PAF, 1), 
    "% of all lung cancer in population could be prevented by eliminating smoking\n")

# Number needed to harm (NNH)
NNH <- 1 / RD
cat("\nNumber needed to harm(NNH):", round(NNH, 1), "\n")
cat("Interpretation: For every", round(NNH, 0), 
    "smokers, 1 excess lung cancer case occurs over 20 years\n")

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat(sprintf(
  "A 20-year prospective cohort study(n=5,000) examined the association between\nsmoking and incident lung cancer.\n\nAmong 2,500 smokers, %d(%%. 1f%%) developed lung cancer compared to %d(%.1f%%)\nof 2,500 never-smokers.\n\nThe relative risk was RR = %.2f (95%% CI [%.2f, %.2f], p < .001), indicating\nsmokers had %.1f times the risk of lung cancer compared to never-smokers.\n\nThe risk difference was %.2f percentage points(95%% CI [%.2f, %.2f]), meaning\nan excess %.1f%% of smokers developed lung cancer attributable to smoking.\n\nThe attributable fraction among exposed was %.1f%%, indicating %.1f%% of lung\ncancer cases in smokers would not have occurred without smoking.\n\nAssuming 50%% smoking prevalence, the population attributable fraction was %.1f%%,\nsuggesting %.1f%% of all lung cancer could be prevented by eliminating smoking.\n\nThe number needed to harm was %.0f, meaning for every %.0f smokers over 20 years,\n1 excess lung cancer case occurred due to smoking.\n\nInterpretation: This very large relative risk(RR > 10) provides strong evidence\nfor a causal relationship between smoking and lung cancer, consistent with\nepidemiological consensus and biological plausibility(Doll & Hill, 1954).",
  smokers_cancer, risk_smokers * 100,
  never_cancer, risk_never * 100,
  RR_manual, CI_lower, CI_upper, RR_manual,
  RD * 100, RD_CI_lower * 100, RD_CI_upper * 100, RD * 100,
  AFe, AFe,
  PAF, PAF,
  NNH, NNH
))