# Binary outcomes (2x2 data)
study,year,events_int,n_int,events_ctrl,n_ctrl,subgroup
Smith,2020,10,50,5,50,early
Jones,2021,15,60,8,55,late

# Continuous outcomes
study,year,n_int,mean_int,sd_int,n_ctrl,mean_ctrl,sd_ctrl
Smith,2020,50,25.3,5.2,50,30.1,6.1

# Proportions (single-arm)
study,year,events,total
Smith,2020,42,50
Jones,2021,38,55

Run a meta-analysis of mortality using OR with random effects.
Input: extractions/pooled_data.csv

Run a meta-analysis of ODI scores using MD.
Input: extractions/continuous_data.csv

Calculate the pooled proportion of surgical success.
Input: extractions/case_series_data.csv

library(meta)
library(metafor)

data <- read.csv("{{INPUT}}")

ma <- metabin(
  event.e = events_int,
  n.e = n_int,
  event.c = events_ctrl,
  n.c = n_ctrl,
  studlab = paste(study, year),
  data = data,
  sm = "{{MEASURE}}",  # OR, RR, RD
  method = "MH",        # MH, Inverse, GLMM
  random = TRUE,
  fixed = FALSE,
  prediction = TRUE,
  hakn = TRUE           # Knapp-Hartung adjustment
)

# Forest plot
png("{{OUTPUT}}/forest_plot.png", width = 1200, height = 800, res = 150)
forest(ma,
       sortvar = TE,
       prediction = TRUE,
       print.tau2 = TRUE,
       print.I2 = TRUE,
       leftcols = c("studlab", "event.e", "n.e", "event.c", "n.c"),
       rightcols = c("effect", "ci", "w.random"))
dev.off()

# Comprehensive heterogeneity assessment
cat("=== HETEROGENEITY ASSESSMENT ===\n\n")

# Statistics
cat(sprintf("I² = %.1f%% [%.1f%%, %.1f%%]\n",
            ma$I2*100, ma$lower.I2*100, ma$upper.I2*100))
cat(sprintf("τ² = %.4f (τ = %.4f)\n", ma$tau2, sqrt(ma$tau2)))
cat(sprintf("H² = %.2f\n", ma$H^2))
cat(sprintf("Q = %.2f, df = %d, p = %.4f\n", ma$Q, ma$df.Q, ma$pval.Q))
cat(sprintf("\nPrediction interval: [%.2f, %.2f]\n",
            exp(ma$lower.predict), exp(ma$upper.predict)))

# Baujat plot (outlier detection)
png("{{OUTPUT}}/baujat_plot.png", width = 800, height = 600, res = 150)
baujat(ma)
dev.off()

# GOSH plot (if < 15 studies)
if (ma$k <= 15) {
  gosh_result <- gosh(ma)
  png("{{OUTPUT}}/gosh_plot.png", width = 800, height = 600, res = 150)
  plot(gosh_result)
  dev.off()
}

# Only reliable if k >= 10
if (ma$k >= 10) {
  # Funnel plot
  png("{{OUTPUT}}/funnel_plot.png", width = 800, height = 600, res = 150)
  funnel(ma, studlab = TRUE, cex.studlab = 0.7)
  dev.off()
  
  # Contour-enhanced funnel plot
  png("{{OUTPUT}}/funnel_contour.png", width = 800, height = 600, res = 150)
  funnel(ma, contour = c(0.9, 0.95, 0.99),
         col.contour = c("darkgray", "gray", "lightgray"))
  legend("topright", c("p < 0.10", "p < 0.05", "p < 0.01"),
         fill = c("darkgray", "gray", "lightgray"), bty = "n")
  dev.off()
  
  # Statistical tests
  cat("\n=== PUBLICATION BIAS TESTS ===\n\n")
  
  # Egger's test (continuous outcomes)
  egger <- metabias(ma, method = "Egger")
  cat(sprintf("Egger's test: t = %.2f, p = %.4f\n",
              egger$statistic, egger$p.value))
  
  # Peters' test (binary outcomes - preferred)
  peters <- metabias(ma, method = "Peters")
  cat(sprintf("Peters' test: t = %.2f, p = %.4f\n",
              peters$statistic, peters$p.value))
  
  # Trim-and-fill
  tf <- trimfill(ma)
  cat(sprintf("\nTrim-and-fill: %d studies imputed\n", tf$k0))
  cat(sprintf("Adjusted estimate: %.2f [%.2f, %.2f]\n",
              exp(tf$TE.random), exp(tf$lower.random), exp(tf$upper.random)))
  
  # Funnel with imputed studies
  png("{{OUTPUT}}/funnel_trimfill.png", width = 800, height = 600, res = 150)
  funnel(tf, studlab = TRUE)
  dev.off()
  
} else {
  cat("\nNote: k < 10, publication bias tests unreliable.\n")
  cat("Visual inspection of funnel plot only.\n")
}

# Leave-one-out analysis
l1o <- metainf(ma, pooled = "random")
png("{{OUTPUT}}/leave_one_out.png", width = 1000, height = 600, res = 150)
forest(l1o)
dev.off()

# Influence diagnostics
inf <- influence(ma)
png("{{OUTPUT}}/influence.png", width = 1000, height = 800, res = 150)
plot(inf)
dev.off()

# Cumulative meta-analysis (by year)
cum <- metacum(ma, sortvar = data$year)
png("{{OUTPUT}}/cumulative.png", width = 1000, height = 600, res = 150)
forest(cum)
dev.off()

# Report influential studies
cat("\n=== INFLUENTIAL STUDIES ===\n\n")
influential <- which(inf$is.infl)
if (length(influential) > 0) {
  cat("The following studies are potentially influential:\n")
  for (i in influential) {
    cat(sprintf("  - %s\n", ma$studlab[i]))
  }
} else {
  cat("No studies identified as influential.\n")
}

# Update with subgroup
ma_sub <- update(ma, subgroup = data${{SUBGROUP}})

# Forest plot with subgroups
png("{{OUTPUT}}/forest_subgroup.png", width = 1400, height = 1000, res = 150)
forest(ma_sub,
       sortvar = TE,
       prediction = TRUE,
       subgroup = TRUE,
       print.subgroup.labels = TRUE,
       subgroup.hetstat = TRUE)
dev.off()

# Test for subgroup differences
cat("\n=== SUBGROUP ANALYSIS ===\n\n")
cat(sprintf("Test for subgroup differences:\n"))
cat(sprintf("  Q = %.2f, df = %d, p = %.4f\n",
            ma_sub$Q.b.random, ma_sub$df.Q.b, ma_sub$pval.Q.b.random))

# Within-subgroup heterogeneity
for (sg in unique(data${{SUBGROUP}})) {
  ma_sg <- subset(ma_sub, subgroup == sg)
  cat(sprintf("\n%s (k = %d):\n", sg, ma_sg$k))
  cat(sprintf("  Effect: %.2f [%.2f, %.2f]\n",
              exp(ma_sg$TE.random), exp(ma_sg$lower.random), exp(ma_sg$upper.random)))
  cat(sprintf("  I² = %.1f%%\n", ma_sg$I2 * 100))
}

metabin(..., incr = 0.5, method.incr = "only0")

library(estmeansd)
bc.mean.sd(q1.val = Q1, med.val = median, q3.val = Q3, n = n)

log_hr <- log(hr)
se <- (log(ci_upper) - log(ci_lower)) / 3.92