Odds Ratio (OR) Test

# Odds Ratio: Smoking and Lung Cancer (Case-Control Study)
# Simulates classic epidemiological design establishing smoking-cancer link

library(epitools)   # For oddsratio() with CI
library(vcd)        # For fourfold plot
library(ggplot2)
library(dplyr)

# Set seed for reproducibility
set.seed(2025)

# === Simulate realistic case-control data ===
# OR ≈ 15 (strong association based on Doll & Hill findings)
# Exposure prevalence: 80% in cases, 30% in controls

n_cases <- 200
n_controls <- 200

# Cases: 80% exposed (160 smokers, 40 non-smokers)
cases_exposed <- 160
cases_unexposed <- 40

# Controls: 30% exposed (60 smokers, 140 non-smokers)
controls_exposed <- 60
controls_unexposed <- 140

# Create 2x2 contingency table
# Rows: Exposure (Smoker, Non-smoker)
# Columns: Outcome (Case, Control)
table_data <- matrix(c(cases_exposed, controls_exposed,
                       cases_unexposed, controls_unexposed),
                     nrow = 2, byrow = TRUE,
                     dimnames = list(
                       Exposure = c("Smoker", "Non-smoker"),
                       Outcome = c("Case", "Control")
                     ))

print("=== 2x2 Contingency Table ===")
print(table_data)
print(addmargins(table_data))  # Add row/column totals

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

# 1. Independence: Verify study design (case-control, no matching)
cat("\n=== Assumption Checks ===\n")
cat("Independence: Case-control design with independent sampling ✓\n")

# 2. Adequate cell counts (all ≥ 5)
min_cell <- min(table_data)
cat("Minimum cell count:", min_cell, ifelse(min_cell >= 5, "✓", "✗ Use exact methods"), "\n")

# 3. No zero cells
zero_cells <- sum(table_data == 0)
cat("Zero cells:", zero_cells, ifelse(zero_cells == 0, "✓", "✗ Add continuity correction"), "\n")

# === STEP 2: Calculate Odds Ratio with 95% CI ===

cat("\n=== Odds Ratio Calculation ===\n")

# Method 1: Manual calculation
a <- table_data[1,1]  # Cases exposed
b <- table_data[1,2]  # Controls exposed
c <- table_data[2,1]  # Cases unexposed
d <- table_data[2,2]  # Controls unexposed

OR_manual <- (a * d) / (b * c)
cat("Manual OR = (a×d)/(b×c) =", OR_manual, "\n")

# Calculate 95% CI using Woolf's method (log scale)
log_OR <- log(OR_manual)
SE_log_OR <- sqrt(1/a + 1/b + 1/c + 1/d)
CI_lower <- exp(log_OR - 1.96 * SE_log_OR)
CI_upper <- exp(log_OR + 1.96 * SE_log_OR)

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

# Method 2: Using epitools package
or_result <- oddsratio(table_data, method = "wald")
print(or_result)

# Method 3: Fisher's exact test (provides exact p-value)
fisher_result <- fisher.test(table_data)
cat("\nFisher's Exact Test:")
cat("\nOR =", fisher_result$estimate)
cat("\n95% CI:", fisher_result$conf.int[1], "-", fisher_result$conf.int[2])
cat("\np-value =", fisher_result$p.value, "\n")

# === STEP 3: Hypothesis Test ===
cat("\n=== Hypothesis Test ===\n")
cat("H₀: OR = 1 (no association)\n")
cat("Hₐ: OR ≠ 1 (association exists)\n")

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

# === STEP 4: Visualizations ===

# 4.1 Fourfold plot (association display)
png("fourfold_plot.png", width=800, height=600)
fourfoldplot(table_data, color=c("#E69F00", "#56B4E9"),
             conf.level=0, margin=1,
             main="Smoking and Lung Cancer\n(Case-Control Study)")
dev.off()

# 4.2 Mosaic plot
png("mosaic_plot.png", width=800, height=600)
mosaicplot(table_data, color=TRUE, shade=TRUE,
           main="Association: Smoking → Lung Cancer",
           xlab="Smoking Status", ylab="Disease Status")
dev.off()

# 4.3 Forest plot (OR with CI)
forest_data <- data.frame(
  Study = "Smoking vs Lung Cancer",
  OR = OR_manual,
  CI_lower = CI_lower,
  CI_upper = CI_upper
)

library(ggplot2)
ggplot(forest_data, aes(y=Study, x=OR)) +
  geom_point(size=4, color="darkblue") +
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), height=0.2, linewidth=1) +
  geom_vline(xintercept=1, linetype="dashed", color="red", linewidth=1) +
  scale_x_log10(breaks=c(0.5, 1, 2, 5, 10, 20, 30)) +
  labs(title="Odds Ratio: Smoking and Lung Cancer",
       subtitle="Case-Control Study(n=400)",
       x="Odds Ratio(log scale) with 95% CI",
       y="") +
  theme_classic(base_size=14) +
  annotate("text", x=OR_manual, y=1.3, 
           label=paste0("OR = ", round(OR_manual, 2), 
                       "\n95% CI: ", round(CI_lower, 2), "-", round(CI_upper, 2)),
           size=5, fontface="bold")

ggsave("forest_plot_or.png", width=10, height=6)

# 4.4 Bar plot comparing exposure prevalence
exposure_prev <- data.frame(
  Group = c("Cases\n(Lung Cancer)", "Controls\n(No Cancer)"),
  Prevalence = c(cases_exposed/n_cases * 100, 
                 controls_exposed/n_controls * 100)
)

ggplot(exposure_prev, aes(x=Group, y=Prevalence, fill=Group)) +
  geom_bar(stat="identity", width=0.6, alpha=0.8) +
  geom_text(aes(label=paste0(round(Prevalence, 1), "%")), 
            vjust=-0.5, size=6, fontface="bold") +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) +
  labs(title="Smoking Prevalence: Cases vs Controls",
       y="Smoking Prevalence(%)",
       x="") +
  theme_classic(base_size=14) +
  theme(legend.position="none") +
  ylim(0, 100)

ggsave("exposure_prevalence.png", width=8, height=6)

# === STEP 5: Sensitivity Analysis (Continuity Correction) ===
cat("\n=== Sensitivity Analysis ===\n")

# Add 0.5 to all cells (Haldane-Anscombe correction)
table_corrected <- table_data + 0.5
OR_corrected <- (table_corrected[1,1] * table_corrected[2,2]) / 
                (table_corrected[1,2] * table_corrected[2,1])

cat("OR with continuity correction(+0.5):", round(OR_corrected, 2), "\n")
cat("Change from uncorrected OR:", round(OR_corrected - OR_manual, 2), "\n")

# === APA-Style Reporting ===
cat("\n=== APA-Style Report ===\n")
cat(sprintf(
  "A case-control study(n=400) examined the association between smoking and lung cancer.\n
Among 200 lung cancer cases, 160 (80%%) were smokers compared to 60 (30%%) of 200 controls.\n
The odds of lung cancer were %.1f times higher in smokers than non-smokers\n(OR = %.2f, 95%% CI [%.2f, %.2f], p < .001, Fisher's exact test).\n\nThis strong association(OR > 10) is consistent with the causal relationship between\nsmoking and lung cancer established in landmark epidemiological studies(Doll & Hill, 1950).\n\nThe 95%% confidence interval excludes 1.0, indicating statistical significance at α = .05.\n\nInterpretation: Smokers have %.0f times the odds of lung cancer compared to non-smokers\nin this case-control sample. This OR magnitude aligns with meta-analytic estimates\n(pooled OR ≈ 8-20 for current smokers).",
  OR_manual, OR_manual, CI_lower, CI_upper, OR_manual
))

# === STEP 6: Calculate Additional Measures ===
cat("\n\n=== Additional Effect Measures ===\n")

# Attributable fraction among exposed (AFe)
AFe <- (OR_manual - 1) / OR_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)
p_exposed <- (cases_exposed + controls_exposed) / (n_cases + n_controls)
PAF <- p_exposed * (OR_manual - 1) / (1 + p_exposed * (OR_manual - 1)) * 100
cat("\nPopulation attributable fraction:", round(PAF, 1), "%\n")
cat("Interpretation:", round(PAF, 1), "% of all lung cancer could be prevented by eliminating smoking\n")