# Supplementary File 2
# Complete statistical-analysis and figure-generation code
# Manuscript: Slow Flow/No-Reflow After Primary Percutaneous Coronary Intervention in ST-Segment Elevation Myocardial Infarction: A Meta-Analysis
#
# Run this script from a directory containing Supplementary_File_1_Raw_Data.xlsx.
# The script reproduces the data preparation, primary and sensitivity meta-analyses, risk-of-bias and GRADE summaries, tables, figures, and validation checks.


# ============================================================
# lib/00_setup.R
# ============================================================
options(stringsAsFactors = FALSE, scipen = 999)

suppressPackageStartupMessages({
  library(readxl)
  library(openxlsx)
  library(dplyr)
  library(tidyr)
  library(purrr)
  library(ggplot2)
  library(patchwork)
  library(metafor)
  library(svglite)
})

set.seed(74063)

project_root <- normalizePath(file.path(getwd()), winslash = "/", mustWork = TRUE)
output_root <- normalizePath(file.path(project_root, "Reanalysis_Output"), winslash = "/", mustWork = FALSE)
input_workbook <- file.path(project_root, "Supplementary_File_1_Raw_Data.xlsx")

dirs <- c("DATA", "TABLES", "FIGURES", "VALIDATION")
walk(file.path(output_root, dirs), dir.create, recursive = TRUE, showWarnings = FALSE)

palette_contract <- c(
  neutral_dark = "#272727", neutral_mid = "#767676", neutral_light = "#D8D8D8",
  signal_blue = "#2F6B9A", signal_teal = "#2A9D8F", accent_orange = "#D9822B",
  accent_red = "#C44E52", pale_blue = "#DCEAF4", pale_orange = "#F6E4CF"
)

theme_meta <- function(base_size = 8, base_family = "Arial") {
  theme_classic(base_size = base_size, base_family = base_family) +
    theme(
      axis.line = element_line(linewidth = 0.35, colour = "black"),
      axis.ticks = element_line(linewidth = 0.35, colour = "black"),
      panel.grid = element_blank(),
      strip.background = element_blank(),
      strip.text = element_text(face = "bold", size = base_size),
      legend.title = element_text(size = base_size),
      legend.text = element_text(size = base_size - 0.5),
      plot.title = element_text(face = "bold", size = base_size + 1),
      plot.subtitle = element_text(size = base_size),
      plot.caption = element_text(size = base_size - 0.5, hjust = 0),
      plot.margin = margin(5, 7, 5, 7)
    )
}
theme_set(theme_meta())

ascii_clean <- function(x) {
  replacements <- c(
    "<U+00AD>" = "", "<U+00B1>" = "+/-", "<U+00B2>" = "2",
    "<U+00D7>" = "x", "<U+00ED>" = "i", "<U+00F6>" = "o",
    "<U+00F8>" = "o", "<U+00FC>" = "u", "<U+011F>" = "g",
    "<U+0131>" = "i", "<U+015E>" = "S", "<U+015F>" = "s",
    "<U+2002>" = " ", "<U+2013>" = "-", "<U+2014>" = "-",
    "<U+201C>" = "\"", "<U+201D>" = "\"", "<U+2264>" = "<=",
    "<U+2265>" = ">="
  )
  out <- as.character(x)
  for (needle in names(replacements)) out <- gsub(needle, replacements[[needle]], out, fixed = TRUE)
  unicode_codes <- c(
    "00AD" = "", "00B1" = "+/-", "00B2" = "2", "00D7" = "x",
    "00ED" = "i", "00F6" = "o", "00F8" = "o", "00FC" = "u",
    "011F" = "g", "0131" = "i", "015E" = "S", "015F" = "s",
    "2002" = " ", "2013" = "-", "2014" = "-", "201C" = "\"",
    "201D" = "\"", "2264" = "<=", "2265" = ">="
  )
  for (hex in names(unicode_codes)) {
    out <- gsub(intToUtf8(strtoi(hex, base = 16L)), unicode_codes[[hex]], out, fixed = TRUE)
  }
  out
}

write_csv_utf8 <- function(x, path) {
  x <- as.data.frame(x)
  x[] <- lapply(x, function(col) if (is.character(col)) ascii_clean(col) else col)
  write.csv(x, path, row.names = FALSE, na = "", fileEncoding = "UTF-8")
}

or_effect <- function(or, lo, hi) {
  tibble(yi = log(as.numeric(or)), sei = (log(as.numeric(hi)) - log(as.numeric(lo))) / (2 * 1.96)) %>%
    mutate(vi = sei^2)
}

pi_null_value <- function(measure) ifelse(measure == "OR", 1, 0)

pi_relation <- function(measure, lower, upper) {
  null <- pi_null_value(measure)
  includes <- is.finite(lower) && is.finite(upper) && lower <= null && upper >= null
  paste(if (includes) "includes" else "excludes", if (measure == "OR") "1" else "0")
}

pi_sentence <- function(measure, lower, upper) {
  paste0("Prediction interval ", pi_relation(measure, lower, upper), ".")
}

fit_meta <- function(dat, factor, analysis, domain, measure = "OR", interpretation = "") {
  dat <- dat %>% filter(is.finite(yi), is.finite(vi), vi > 0)
  if (nrow(dat) < 3) stop("Meta-analysis requires k >= 3: ", factor, " / ", analysis)
  model <- rma.uni(yi = dat$yi, vi = dat$vi, method = "REML", test = "knha")
  pred <- predict(model)
  transform <- if (measure == "OR") exp else identity
  pi_lower_value <- transform(pred$pi.lb)
  pi_upper_value <- transform(pred$pi.ub)
  pi_includes_null <- pi_lower_value <= pi_null_value(measure) && pi_upper_value >= pi_null_value(measure)
  interpretation <- paste(trimws(interpretation), pi_sentence(measure, pi_lower_value, pi_upper_value))
  tibble(
    domain = domain,
    factor = factor,
    analysis = analysis,
    measure = measure,
    k = model$k,
    estimate = transform(as.numeric(model$b)),
    ci_lower = transform(model$ci.lb),
    ci_upper = transform(model$ci.ub),
    tau2 = model$tau2,
    Q = model$QE,
    Q_p = model$QEp,
    I2 = max(0, (model$QE - (model$k - 1)) / model$QE) * 100,
    pi_lower = pi_lower_value,
    pi_upper = pi_upper_value,
    pi_includes_null = pi_includes_null,
    study_ids = paste(dat$Study_No, collapse = ","),
    interpretation = interpretation,
    model_object = list(model),
    input_data = list(dat)
  )
}

save_pub_r <- function(plot, filename, width_mm = 183, height_mm = 120, dpi = 600) {
  w <- width_mm / 25.4
  h <- height_mm / 25.4
  svglite::svglite(paste0(filename, ".svg"), width = w, height = h)
  print(plot)
  dev.off()
  grDevices::cairo_pdf(paste0(filename, ".pdf"), width = w, height = h, family = "Arial")
  print(plot)
  dev.off()
  grDevices::tiff(paste0(filename, ".tiff"), width = w, height = h, units = "in", res = dpi,
                  compression = "lzw", type = "cairo", family = "Arial")
  print(plot)
  dev.off()
  grDevices::png(paste0(filename, ".png"), width = w, height = h, units = "in", res = 300,
                 type = "cairo", family = "Arial")
  print(plot)
  dev.off()
}

format_effect <- function(measure, est, lo, hi) {
  digits <- 2
  sprintf(paste0("%.", digits, "f (%.", digits, "f to %.", digits, "f)"), est, lo, hi)
}

style_workbook <- function(wb, sheets) {
  header <- createStyle(fontName = "Arial", fontSize = 10, fontColour = "#FFFFFF", fgFill = "#274C77",
                        textDecoration = "bold", halign = "center", valign = "center", wrapText = TRUE,
                        border = "Bottom", borderColour = "#FFFFFF")
  body <- createStyle(fontName = "Arial", fontSize = 9, valign = "top")
  note <- createStyle(fontName = "Arial", fontSize = 9, fontColour = "#5A5A5A", textDecoration = "italic")
  for (s in sheets) {
    nr <- nrow(readWorkbook(wb, sheet = s, colNames = FALSE))
    nc <- ncol(readWorkbook(wb, sheet = s, colNames = FALSE))
    if (nr > 0 && nc > 0) {
      addStyle(wb, s, header, rows = 1, cols = seq_len(nc), gridExpand = TRUE, stack = TRUE)
      if (nr > 1) addStyle(wb, s, body, rows = 2:nr, cols = seq_len(nc), gridExpand = TRUE, stack = TRUE)
      setRowHeights(wb, s, rows = 1, heights = 34)
      freezePane(wb, s, firstActiveRow = 2)
      addFilter(wb, s, rows = 1, cols = seq_len(nc))
      setColWidths(wb, s, cols = seq_len(nc), widths = "auto")
      showGridLines(wb, s, showGridLines = FALSE)
    }
  }
  invisible(wb)
}


# ============================================================
# lib/01_import_and_validate.R
# ============================================================
sheet_names <- excel_sheets(input_workbook)
required_sheets <- c(
  "Studies", "Eligibility", "Age_Effects", "Effect_Sources", "Binary_Effects",
  "Killip_Effects", "Thrombus_Effects", "Risk_of_Bias", "GRADE", "Descriptive_Evidence"
)
missing_sheets <- setdiff(required_sheets, sheet_names)
if (length(missing_sheets)) stop("Missing required sheets: ", paste(missing_sheets, collapse = ", "))

read_sheet <- function(s) read_excel(input_workbook, sheet = s, .name_repair = "unique")

source_data <- list(
  EligibilityData = read_sheet("Eligibility"),
  BinaryEffects = read_sheet("Binary_Effects"),
  DescriptiveEvidence = read_sheet("Descriptive_Evidence"),
  GradeData = read_sheet("GRADE")
)
support <- list(
  Studies = read_sheet("Studies"),
  AgeEffects = read_sheet("Age_Effects"),
  EffectSources = read_sheet("Effect_Sources"),
  KillipEffects = read_sheet("Killip_Effects"),
  ThrombusEffects = read_sheet("Thrombus_Effects"),
  RiskOfBias = read_sheet("Risk_of_Bias")
)

elig <- source_data$EligibilityData
if (nrow(elig) != 44) stop("Eligibility data must retain all 44 assessed records.")
included <- elig %>% filter(grepl("^Included", Eligibility))
if (nrow(included) != 43) stop("Final review must contain exactly 43 studies.")
if (!any(elig$Study_No == 33 & elig$Eligibility != "Included")) stop("The mixed-population record exclusion was not found.")
if (!any(elig$Study_No == 41 & grepl("eligibility sensitivity", elig$Sensitivity_status, ignore.case = TRUE))) stop("The delayed-PCI sensitivity status was not found.")

validation_summary <- tibble(
  check = c("Assessed records", "Final included studies", "Mixed-population record excluded", "Delayed-PCI sensitivity flag",
            "Qu 2024 source present", "Zhang and He DOI present"),
  observed = c(
    nrow(elig), nrow(included),
    any(elig$Study_No == 33 & elig$Eligibility != "Included"),
    any(elig$Study_No == 41 & grepl("eligibility sensitivity", elig$Sensitivity_status, ignore.case = TRUE)),
    any(support$Studies$Study_No == 16 & grepl("10.62347/HQVS4428", support$Studies$DOI, fixed = TRUE)),
    any(support$Studies$Study_No == 44 & grepl("10.14744/AnatolJCardiol.2020.29267", support$Studies$DOI, fixed = TRUE))
  ),
  expected = c(44, 43, TRUE, TRUE, TRUE, TRUE)
) %>% mutate(status = ifelse(as.character(observed) == as.character(expected), "PASS", "ERROR"))
if (any(validation_summary$status != "PASS")) stop("Input validation failed.")


# ============================================================
# lib/02_build_effect_sizes.R
# ============================================================
study <- support$Studies %>%
  transmute(
    Study_No = as.integer(Study_No), Study, DOI, Country, Design,
    NR_n = as.numeric(NR_n), Normal_n = as.numeric(Normal_n), Total_N = as.numeric(Total_N),
    outcome_definition = Outcome_Definition, outcome_class = Outcome_Class,
    eligibility_note = Primary_PCI_Eligibility,
    boundary_note = Special_Population_or_Eligibility_Flag
  ) %>%
  left_join(elig %>% transmute(Study_No = as.integer(Study_No), final_eligibility = Eligibility,
                               primary_status = Primary_analysis_status, sensitivity_status = Sensitivity_status,
                               eligibility_reason = Reason), by = "Study_No")

final_study <- study %>% filter(grepl("^Included", final_eligibility))
write_csv_utf8(final_study, file.path(output_root, "DATA", "final_study_characteristics.csv"))

binary_es <- function(dat, factor_name = NULL, source_sheet = "") {
  if (!is.null(factor_name)) dat$Factor <- factor_name
  dat <- dat %>%
    mutate(across(c(Study_No, a, b, c, d), as.numeric)) %>%
    filter(is.finite(a), is.finite(b), is.finite(c), is.finite(d))
  es <- escalc(measure = "OR", ai = a, bi = b, ci = c, di = d, data = dat, add = 0.5, to = "only0")
  as_tibble(es) %>%
    mutate(sei = sqrt(vi), source_type = "raw-derived", source_sheet = source_sheet,
           effect = exp(yi), ci_lower = exp(yi - 1.96 * sei), ci_upper = exp(yi + 1.96 * sei),
           contrast_unit = "Exposed vs unexposed")
}

trace12 <- support$EffectSources %>%
  transmute(
    Factor, Study_No = as.integer(Study_No), Study, outcome_class = Outcome_Class,
    a = as.numeric(NR_events), nr_total = as.numeric(NR_total), c = as.numeric(Normal_events), normal_total = as.numeric(Normal_total),
    primary_use = Primary_analysis_use, adjusted_use = Adjusted_sensitivity_use, source_note = Notes,
    b = nr_total - a, d = normal_total - c
  )

hypertension <- binary_es(trace12 %>% filter(Factor == "Hypertension"), source_sheet = "Effect_Sources")
diabetes <- binary_es(trace12 %>% filter(Factor == "Diabetes mellitus"), source_sheet = "Effect_Sources")

final_bin <- source_data$BinaryEffects %>%
  transmute(Factor, Study_No = as.integer(Study_No), a = as.numeric(NR_event), c = as.numeric(Normal_event),
            trace_status = Extraction_basis, source_note = Source_excerpt) %>%
  left_join(study %>% select(Study_No, Study, NR_n, Normal_n, outcome_class, final_eligibility), by = "Study_No") %>%
  mutate(b = NR_n - a, d = Normal_n - c) %>%
  filter(grepl("^Included", final_eligibility), !grepl("exclude", trace_status, ignore.case = TRUE))

from_final_bin <- function(name, exclude_ids = integer()) {
  binary_es(final_bin %>% filter(Factor == name, !Study_No %in% exclude_ids), source_sheet = "Binary_Effects")
}

male <- from_final_bin("Male sex")
smoking <- from_final_bin("Smoking")
dyslipidemia <- from_final_bin("Dyslipidemia/hyperlipidemia", 24)
previous_mi <- from_final_bin("Previous myocardial infarction")
advanced_killip <- from_final_bin(unique(final_bin$Factor[grepl("Advanced Killip", final_bin$Factor)]))
multivessel <- from_final_bin("Multivessel disease", 4)
thrombus_aspiration <- from_final_bin("Thrombus aspiration", 43)
direct_stenting <- from_final_bin("Direct stenting")

killip <- support$KillipEffects %>%
  transmute(Factor = "Killip >=II / acute HF", Study_No = as.integer(Study_No),
            a = as.numeric(`NR exposed`), nr_total = as.numeric(`NR total`),
            c = as.numeric(`Normal exposed`), normal_total = as.numeric(`Normal total`),
            source_note = Notes, b = nr_total - a, d = normal_total - c) %>%
  left_join(study %>% select(Study_No, Study, outcome_class), by = "Study_No") %>%
  binary_es(source_sheet = "Killip_Effects")

high_thrombus <- support$ThrombusEffects %>%
  transmute(Factor = "High thrombus burden >=4", Study_No = as.integer(Study_No),
            a = as.numeric(`NR high-thrombus`), nr_total = as.numeric(`NR total`),
            c = as.numeric(`Normal high-thrombus`), normal_total = as.numeric(`Normal total`),
            source_note = `Adjusted-use note`, b = nr_total - a, d = normal_total - c) %>%
  left_join(study %>% select(Study_No, Study, outcome_class), by = "Study_No") %>%
  binary_es(source_sheet = "Thrombus_Effects")

age <- support$AgeEffects %>%
  transmute(Study_No = as.integer(Study_No), Study, Factor = "Age", yi = as.numeric(MD_years),
            ci_lower = as.numeric(CI_Lower), ci_upper = as.numeric(CI_Upper), source_note = Source) %>%
  mutate(sei = (ci_upper - ci_lower) / (2 * 1.96), vi = sei^2, source_type = "raw-derived",
         source_sheet = "Age_Effects", effect = yi, contrast_unit = "Mean difference, years") %>%
  left_join(study %>% select(Study_No, outcome_class, final_eligibility), by = "Study_No") %>%
  filter(grepl("^Included", final_eligibility))

stent_inputs <- tibble(
  Study_No = c(4, 6, 22, 28), Factor = "Stent length",
  m1 = c(25.00, 24.0, 33.98, 26.5), sd1 = c(10.30, 10.0, 18.12, 6.5),
  m2 = c(20.61, 21.0, 29.27, 24.8), sd2 = c(7.73, 8.3, 15.10, 8.3)
) %>% left_join(study %>% select(Study_No, Study, NR_n, Normal_n, outcome_class), by = "Study_No")
stent_length <- escalc(measure = "MD", m1i = m1, sd1i = sd1, n1i = NR_n, m2i = m2, sd2i = sd2, n2i = Normal_n,
                       data = stent_inputs) %>% as_tibble() %>%
  mutate(sei = sqrt(vi), source_type = "raw-derived", source_sheet = "Continuous source data",
         effect = yi, ci_lower = yi - 1.96 * sei, ci_upper = yi + 1.96 * sei,
         contrast_unit = "Mean difference, mm", source_note = "Source-confirmed mean, SD and group size")

baseline_timi_inputs <- tibble(
  Study_No = c(4, 5, 20, 21, 22, 29, 36, 41), Factor = "Impaired baseline TIMI flow",
  a = c(225, 307, 164, 405, 125, 55, 68, 28), nr_total = c(226, 324, 182, 491, 139, 98, 80, 28),
  c = c(480, 2283, 574, 847, 240, 68, 300, 131), normal_total = c(537, 2881, 776, 1167, 321, 119, 374, 201),
  source_note = c("Baseline TIMI <3", "Baseline TIMI <=1", "Initial MBG <2", "Preoperative TIMI <=1",
                  "Initial TIMI <3", "Initial TIMI 0/1", "Initial TIMI 0/1", "Pre-TIMI <=2")
) %>%
  mutate(b = nr_total - a, d = normal_total - c) %>%
  left_join(study %>% select(Study_No, Study, outcome_class), by = "Study_No") %>%
  binary_es(source_sheet = "Study-level source extraction; threshold sensitivity analyses reported")

adjusted_manual <- tribble(
  ~Factor, ~Study_No, ~or, ~lo, ~hi, ~contrast_unit, ~source_note,
  "Age", 5, 1.020, 1.000, 1.040, "per year", "Table 3 / Results",
  "Age", 6, 1.034, 1.014, 1.053, "per year", "Table 3 Model 1",
  "Age", 12, 1.036, 0.997, 1.077, "per year", "Detailed Table II",
  "Age", 34, 1.115, 1.032, 1.205, "per year", "Table 4",
  "Hypertension", 30, 1.432, 0.936, 2.191, "yes vs no", "Table 2",
  "Hypertension", 36, 1.820, 1.040, 3.180, "yes vs no", "Table 5",
  "Hypertension", 42, 1.783, 1.332, 2.388, "yes vs no", "Table 3",
  "Diabetes mellitus", 23, 1.847, 1.209, 2.822, "yes vs no", "Table 3",
  "Diabetes mellitus", 30, 2.173, 1.463, 3.227, "yes vs no", "Detailed Table 2",
  "Diabetes mellitus", 36, 1.730, 0.950, 3.170, "yes vs no", "Table 5",
  "Diabetes mellitus", 42, 2.233, 1.643, 3.035, "yes vs no", "Table 3",
  "Killip >=II / acute HF", 5, 1.990, 1.300, 3.040, ">=II vs lower", "Table 3",
  "Killip >=II / acute HF", 8, 3.590, 1.090, 11.830, "AHF vs no AHF", "Table 5",
  "Killip >=II / acute HF", 25, 1.893, 1.004, 3.570, ">I vs I", "Multivariable model",
  "Killip >=II / acute HF", 30, 1.381, 0.903, 2.112, "II-IV vs I", "Table 2",
  "High thrombus burden >=4", 8, 2.850, 1.100, 7.380, "heavy vs lower", "Table 5",
  "High thrombus burden >=4", 11, 2.330, 1.600, 3.390, ">=4 vs lower", "Table 3",
  "High thrombus burden >=4", 21, 2.033, 1.370, 3.018, ">=4 vs lower", "Development model",
  "High thrombus burden >=4", 32, 3.424, 1.857, 6.312, "high vs lower", "Multivariable model",
  "eGFR", 4, 0.991, 0.984, 0.999, "per unit", "Table 3",
  "eGFR", 10, 0.970, 0.940, 0.990, "per unit", "Table 2",
  "eGFR", 22, 0.960, 0.940, 0.980, "per unit", "Detailed Table 5 used as the prespecified source",
  "SYNTAX score", 1, 1.114, 1.050, 1.183, "per point", "Table 4 / Results",
  "SYNTAX score", 12, 1.015, 0.974, 1.087, "per point", "Detailed Table II",
  "SYNTAX score", 17, 1.065, 1.025, 1.107, "per point", "Multivariable model",
  "CRP / hs-CRP", 1, 1.693, 1.126, 2.547, "per mg/dL", "Table 4",
  "CRP / hs-CRP", 12, 1.115, 0.890, 1.396, "per mg/dL", "Detailed Table II",
  "CRP / hs-CRP", 35, 1.080, 1.000, 1.170, "per mg/dL", "Table 3",
  "Large reference/target vessel diameter", 5, 2.830, 1.870, 4.270, "study cutoff", "Table 3",
  "Large reference/target vessel diameter", 20, 4.686, 2.475, 8.872, "study cutoff", "Table 3",
  "Large reference/target vessel diameter", 21, 1.511, 1.087, 2.100, "study cutoff", "Table 3"
) %>%
  mutate(Study_No = as.integer(Study_No)) %>%
  bind_cols(or_effect(.$or, .$lo, .$hi)) %>%
  left_join(study %>% select(Study_No, Study, outcome_class), by = "Study_No") %>%
  mutate(source_type = "adjusted", source_sheet = "Published multivariable models",
         effect = or, ci_lower = lo, ci_upper = hi)

analysis_sets <- list(
  age = age, male = male, hypertension = hypertension, diabetes = diabetes, smoking = smoking,
  dyslipidemia = dyslipidemia, previous_mi = previous_mi, killip = killip,
  advanced_killip = advanced_killip, high_thrombus = high_thrombus, multivessel = multivessel,
  baseline_timi = baseline_timi_inputs, stent_length = stent_length,
  thrombus_aspiration = thrombus_aspiration, direct_stenting = direct_stenting
)

effect_source <- bind_rows(
  map_dfr(analysis_sets, ~ .x %>% select(any_of(c("Study_No", "Study", "Factor", "source_type", "contrast_unit", "effect", "ci_lower", "ci_upper", "outcome_class", "source_sheet", "source_note")))),
  adjusted_manual %>% select(Study_No, Study, Factor, source_type, contrast_unit, effect, ci_lower, ci_upper,
                             outcome_class, source_sheet, source_note)
) %>%
  mutate(
    included_in_primary = ifelse(source_type == "raw-derived", "Yes", "No"),
    included_in_adjusted_sensitivity = ifelse(source_type == "adjusted", "Yes", "No")
  )

write_csv_utf8(effect_source, file.path(output_root, "DATA", "final_analysis_dataset.csv"))
write_csv_utf8(adjusted_manual, file.path(output_root, "DATA", "adjusted_effects_final.csv"))
write_csv_utf8(source_data$DescriptiveEvidence, file.path(output_root, "DATA", "descriptive_effects.csv"))
write_csv_utf8(support$RiskOfBias, file.path(output_root, "DATA", "rob_grade_inputs.csv"))


# ============================================================
# 01_prepare_data.R
# ============================================================

ascii_study_labels <- c(
  `1`="Safak 2024", `2`="Sen 2024", `3`="Sondergaard 2023", `4`="Suleymanoglu 2020",
  `5`="Tasar 2019", `6`="Bayramoglu 2023", `7`="Chang 2022", `8`="Elkammash 2023",
  `9`="Huang 2021", `10`="Kai 2021", `11`="Kumar 2022", `12`="Koprulu 2024",
  `13`="Li 2022", `14`="Liu 2022", `15`="Mohamed 2021", `16`="Qu 2024",
  `17`="Quisi 2018", `18`="Rashed 2022", `19`="Ruiz-Avalos 2022", `20`="Aggarwal 2021",
  `21`="Yang 2020", `22`="Bamarinejad 2024", `23`="Sun 2023", `24`="Li 2024",
  `25`="Khalfallah 2022a", `26`="Wang 2024", `27`="Wu 2024", `28`="Yang 2022",
  `29`="Yildirim 2021", `30`="Zhang 2022", `31`="Zhang Qinyao 2022", `32`="Khalfallah 2022b",
  `33`="Ma 2024", `34`="Refaat 2021", `35`="Zhao 2019", `36`="Zhang 2020",
  `37`="Liu 2025", `38`="Gao 2018", `39`="Hu 2022", `40`="Wang 2018",
  `41`="Gong 2021", `42`="Zhang 2018", `43`="Simsek 2021", `44`="Zhang 2021"
)
study_label <- function(id) unname(ascii_study_labels[as.character(id)])

parse_mean_sd <- function(x) {
  x <- gsub("<U+00B1>", "|", x, fixed = TRUE)
  x <- gsub(intToUtf8(177), "|", x, fixed = TRUE)
  x <- gsub(intToUtf8(167), "|", x, fixed = TRUE)
  vals <- strsplit(x, "|", fixed = TRUE)
  tibble(
    mean = vapply(vals, function(v) suppressWarnings(as.numeric(trimws(v[1]))), numeric(1)),
    sd = vapply(vals, function(v) if (length(v) > 1) suppressWarnings(as.numeric(trimws(v[2]))) else NA_real_, numeric(1))
  )
}

# Rebuild age effects from the final workbook's group-level mean, SD and N.
age_source <- support$AgeEffects %>%
  select(Study_No, original_study = Study, NR_n, Normal_n) %>%
  left_join(support$Studies %>% transmute(Study_No = as.integer(Study_No), Age_NR, Age_Normal,
                                               outcome_class = Outcome_Class), by = "Study_No")
nr_age <- parse_mean_sd(age_source$Age_NR)
normal_age <- parse_mean_sd(age_source$Age_Normal)
age_inputs <- age_source %>%
  mutate(
    Study = study_label(Study_No),
    nr_mean = nr_age$mean, nr_sd = nr_age$sd,
    normal_mean = normal_age$mean, normal_sd = normal_age$sd
  ) %>%
  filter(is.finite(nr_mean), is.finite(nr_sd), is.finite(normal_mean), is.finite(normal_sd),
         is.finite(NR_n), is.finite(Normal_n), NR_n > 0, Normal_n > 0)
age_all <- escalc(measure = "MD", m1i = nr_mean, sd1i = nr_sd, n1i = NR_n,
              m2i = normal_mean, sd2i = normal_sd, n2i = Normal_n,
              data = age_inputs) %>% as_tibble() %>%
  mutate(
    Factor = "Age", sei = sqrt(vi), source_type = "raw-derived",
    source_sheet = "FINAL study-level age means, SDs and group sizes",
    source_note = "Mean difference from final group-level inputs",
    effect = yi, ci_lower = yi - 1.96 * sei, ci_upper = yi + 1.96 * sei,
    contrast_unit = "Mean difference, years"
  )
age <- age_all %>% filter(Study_No != 33)
analysis_sets$age <- age

age_input_summary <- age_all %>%
  transmute(Study_No, Study, outcome_class, NR_mean = nr_mean, NR_SD = nr_sd, NR_N = NR_n,
            Normal_mean = normal_mean, Normal_SD = normal_sd, Normal_N = Normal_n,
            MD = yi, SE = sei, Included_STEMI_primary = Study_No != 33,
            Included_TIMI_only = Study_No != 33 & outcome_class == "D1 TIMI-based",
            Included_excluding_41 = Study_No != 33 & Study_No != 41)
write_csv_utf8(age_input_summary, file.path(output_root, "DATA", "age_study_level_summary.csv"))

final_study <- final_study %>%
  mutate(Study = study_label(Study_No), DOI = ifelse(Study_No == 44, "10.14744/AnatolJCardiol.2020.29267", DOI))
write_csv_utf8(final_study, file.path(output_root, "DATA", "final_study_characteristics.csv"))

final_exclusion_log <- tibble(
  Study_No = c(33, 41), Study = study_label(c(33, 41)),
  final_status = c("Excluded from final review and all main analyses", "Included in primary set; excluded in eligibility sensitivity"),
  reason = c("Mixed STEMI/NSTEMI population with inseparable STEMI-specific effect",
             "PCI performed 2-7 days after symptom onset; boundary eligibility case"),
  analytical_handling = c("Excluded from the final review set", "Reported in a prespecified exclusion sensitivity")
)
write_csv_utf8(final_exclusion_log, file.path(output_root, "DATA", "final_exclusion_log.csv"))

effect_source <- effect_source %>% filter(Study_No != 33) %>% mutate(Study = study_label(Study_No))
write_csv_utf8(effect_source, file.path(output_root, "DATA", "final_analysis_dataset.csv"))


# ============================================================
# lib/03_primary_meta_analysis.R
# ============================================================
primary_specs <- tribble(
  ~key, ~factor, ~analysis, ~domain, ~measure, ~interpretation,
  "age", "Age", "STEMI-only crude primary", "Clinical", "MD", "Average age difference with substantial between-study heterogeneity.",
  "male", "Male sex", "Uniform crude primary", "Clinical", "OR", "No consistent association.",
  "hypertension", "Hypertension", "Uniform crude primary", "Clinical", "OR", "Crude association is not robust.",
  "diabetes", "Diabetes mellitus", "Uniform crude primary", "Clinical", "OR", "Positive pooled average with marked heterogeneity.",
  "smoking", "Smoking", "Uniform crude primary", "Clinical", "OR", "No consistent association; high heterogeneity.",
  "dyslipidemia", "Dyslipidemia/hyperlipidemia", "Strict crude exploratory", "Clinical", "OR", "Exploratory because definitions vary and heterogeneity is high.",
  "previous_mi", "Previous myocardial infarction", "Crude exploratory", "Clinical", "OR", "Sparse three-study evidence.",
  "killip", "Killip >=II / acute HF", "Strict crude primary", "Clinical", "OR", "Positive association with definition and setting uncertainty.",
  "advanced_killip", "Advanced Killip >=III/IV", "Strict crude", "Clinical", "OR", "Consistent observational association across five studies.",
  "high_thrombus", "High thrombus burden >=4", "Strict crude primary", "Angio/Procedural", "OR", "Positive average association with substantial between-study heterogeneity.",
  "multivessel", "Multivessel disease", "Strict crude exploratory", "Angio/Procedural", "OR", "Nonrobust exploratory association.",
  "baseline_timi", "Impaired baseline TIMI flow", "Broad crude primary", "Angio/Procedural", "OR", "Thresholds and perfusion constructs differ; see definition-restricted sensitivity analyses.",
  "stent_length", "Stent length", "Strict continuous MD", "Angio/Procedural", "MD", "Longer stent length on average; observational association.",
  "thrombus_aspiration", "Thrombus aspiration", "Crude observational association", "Angio/Procedural", "OR", "Confounding by indication; do not interpret causally.",
  "direct_stenting", "Direct stenting", "Crude observational association", "Angio/Procedural", "OR", "Anatomy-dependent treatment selection; do not interpret as benefit."
)

primary_models <- list()
primary_results_list <- vector("list", nrow(primary_specs))
for (i in seq_len(nrow(primary_specs))) {
  spec <- primary_specs[i, ]
  dat <- analysis_sets[[spec$key]]
  res <- fit_meta(dat, spec$factor, spec$analysis, spec$domain, spec$measure, spec$interpretation)
  primary_results_list[[i]] <- res
  primary_models[[spec$key]] <- res$model_object[[1]]
}
primary_results <- bind_rows(primary_results_list)
write_csv_utf8(primary_results %>% select(-model_object, -input_data),
               file.path(output_root, "DATA", "primary_meta_results.csv"))


# ============================================================
# 02_primary_meta.R
# ============================================================



# ============================================================
# lib/04_adjusted_sensitivity.R
# ============================================================
adjusted_specs <- tribble(
  ~factor, ~analysis, ~domain, ~interpretation,
  "Age", "Adjusted-only sensitivity", "Clinical", "Small per-year adjusted association.",
  "Hypertension", "Adjusted-only sensitivity", "Clinical", "Positive adjusted association based on three studies.",
  "Diabetes mellitus", "Adjusted-only sensitivity", "Clinical", "Consistent adjusted association after excluding one incompatible estimate.",
  "Killip >=II / acute HF", "Adjusted-only sensitivity", "Clinical", "Adjusted association remains uncertain across settings.",
  "High thrombus burden >=4", "Adjusted-only sensitivity", "Angio/Procedural", "Consistent adjusted association in compatible models.",
  "eGFR", "Adjusted-only", "Laboratory", "Heterogeneous and not robust under HKSJ.",
  "SYNTAX score", "Adjusted per-point", "Angio/Procedural", "Imprecise and heterogeneous per-point association.",
  "CRP / hs-CRP", "Adjusted exploratory", "Laboratory", "Exploratory; assay definitions vary and evidence is imprecise.",
  "Large reference/target vessel diameter", "Adjusted exploratory", "Angio/Procedural", "Threshold/reference categories differ; exploratory."
)

adjusted_results_list <- vector("list", nrow(adjusted_specs))
adjusted_models <- list()
for (i in seq_len(nrow(adjusted_specs))) {
  spec <- adjusted_specs[i, ]
  dat <- adjusted_manual %>% filter(Factor == spec$factor)
  res <- fit_meta(dat, spec$factor, spec$analysis, spec$domain, "OR", spec$interpretation)
  adjusted_results_list[[i]] <- res
  adjusted_models[[spec$factor]] <- res$model_object[[1]]
}
adjusted_results <- bind_rows(adjusted_results_list)
write_csv_utf8(adjusted_results %>% select(-model_object, -input_data),
               file.path(output_root, "DATA", "adjusted_sensitivity_results.csv"))


# ============================================================
# lib/05_definition_sensitivity.R
# ============================================================
timi_specs <- tribble(
  ~key, ~factor, ~domain, ~measure,
  "age", "Age", "Clinical", "MD",
  "hypertension", "Hypertension", "Clinical", "OR",
  "diabetes", "Diabetes mellitus", "Clinical", "OR"
)

timi_results_list <- vector("list", nrow(timi_specs))
definition_models <- list()
for (i in seq_len(nrow(timi_specs))) {
  spec <- timi_specs[i, ]
  dat <- analysis_sets[[spec$key]] %>% filter(outcome_class == "D1 TIMI-based")
  res <- fit_meta(dat, spec$factor, "TIMI-only (D1) sensitivity", spec$domain, spec$measure,
                  "Definition-restricted sensitivity; TIMI, perfusion-grade and cTFC phenotypes are not assumed identical.")
  timi_results_list[[i]] <- res
  definition_models[[spec$key]] <- res$model_object[[1]]
}
definition_results <- bind_rows(timi_results_list)
write_csv_utf8(definition_results %>% select(-model_object, -input_data),
               file.path(output_root, "DATA", "definition_sensitivity_results.csv"))


# ============================================================
# lib/06_eligibility_sensitivity.R
# ============================================================
eligibility_specs <- tribble(
  ~key, ~factor, ~domain, ~measure,
  "age", "Age", "Clinical", "MD",
  "smoking", "Smoking", "Clinical", "OR",
  "killip", "Killip >=II / acute HF", "Clinical", "OR",
  "baseline_timi", "Impaired baseline TIMI flow", "Angio/Procedural", "OR"
)

eligibility_results_list <- vector("list", nrow(eligibility_specs))
eligibility_models <- list()
for (i in seq_len(nrow(eligibility_specs))) {
  spec <- eligibility_specs[i, ]
  dat <- analysis_sets[[spec$key]] %>% filter(Study_No != 41)
  res <- fit_meta(dat, spec$factor, "Eligibility sensitivity excluding #41", spec$domain, spec$measure,
                  "Excludes the delayed-PCI boundary study (2-7 days after symptom onset).")
  eligibility_results_list[[i]] <- res
  eligibility_models[[spec$key]] <- res$model_object[[1]]
}
eligibility_results <- bind_rows(eligibility_results_list)
write_csv_utf8(eligibility_results %>% select(-model_object, -input_data),
               file.path(output_root, "DATA", "eligibility_sensitivity_results.csv"))

sensitivity_results <- bind_rows(
  adjusted_results %>% mutate(sensitivity_type = "adjusted-only"),
  definition_results %>% mutate(sensitivity_type = "outcome-definition"),
  eligibility_results %>% mutate(sensitivity_type = "eligibility")
)
write_csv_utf8(sensitivity_results %>% select(-model_object, -input_data),
               file.path(output_root, "DATA", "sensitivity_results.csv"))


# ============================================================
# lib/07_leave_one_out_and_small_study.R
# ============================================================
loo_rows <- list()
for (i in seq_len(nrow(primary_results))) {
  row <- primary_results[i, ]
  if (row$k < 5 || row$I2 < 50) next
  dat <- row$input_data[[1]]
  measure <- row$measure
  for (j in seq_len(nrow(dat))) {
    keep <- dat[-j, , drop = FALSE]
    model <- rma.uni(yi = keep$yi, vi = keep$vi, method = "REML", test = "knha")
    tr <- if (measure == "OR") exp else identity
    loo_rows[[length(loo_rows) + 1]] <- tibble(
      factor = row$factor, omitted_study_id = dat$Study_No[j], omitted_study = dat$Study[j],
      k = model$k, measure = measure, estimate = tr(as.numeric(model$b)),
      ci_lower = tr(model$ci.lb), ci_upper = tr(model$ci.ub),
      I2 = max(0, (model$QE - (model$k - 1)) / model$QE) * 100,
      direction_changed = sign(as.numeric(model$b)) != sign(as.numeric(row$model_object[[1]]$b)),
      significance_changed = ((model$ci.lb <= 0 && model$ci.ub >= 0) !=
        (row$model_object[[1]]$ci.lb <= 0 && row$model_object[[1]]$ci.ub >= 0))
    )
  }
}
leave_one_out_results <- bind_rows(loo_rows)
write_csv_utf8(leave_one_out_results, file.path(output_root, "DATA", "leave_one_out_results.csv"))

small_rows <- list()
for (i in seq_len(nrow(primary_results))) {
  row <- primary_results[i, ]
  dat <- row$input_data[[1]]
  if (row$k >= 10) {
    eg <- regtest(row$model_object[[1]], model = "rma", predictor = "sei")
    small_rows[[length(small_rows) + 1]] <- tibble(
      factor = row$factor, k = row$k, test = "Egger regression", statistic = unname(eg$zval),
      p_value = eg$pval,
      interpretation = ifelse(eg$pval < 0.05, "Asymmetry signal detected; interpret cautiously.",
                              "No asymmetry signal detected; this does not establish absence of publication bias.")
    )
  } else {
    small_rows[[length(small_rows) + 1]] <- tibble(
      factor = row$factor, k = row$k, test = "Not performed", statistic = NA_real_, p_value = NA_real_,
      interpretation = "Fewer than 10 studies; funnel/Egger testing is not reliable."
    )
  }
}
small_study_results <- bind_rows(small_rows)
write_csv_utf8(small_study_results, file.path(output_root, "DATA", "small_study_effect_results.csv"))


# ============================================================
# 03_sensitivity_meta.R
# ============================================================

# Resolve baseline-flow threshold heterogeneity with explicit nested sets.
baseline_timi_only <- fit_meta(
  baseline_timi_inputs %>% filter(Study_No != 20),
  "Impaired baseline TIMI flow", "TIMI-flow-only sensitivity", "Angio/Procedural", "OR",
  "Excludes the MBG-based study; TIMI thresholds still range from <3 to <=1/<=2."
)
baseline_strict_01 <- fit_meta(
  baseline_timi_inputs %>% filter(Study_No %in% c(5, 21, 29, 36)),
  "Impaired baseline TIMI flow", "Strict baseline TIMI 0-1 sensitivity", "Angio/Procedural", "OR",
  "Uses only studies defining impaired baseline flow as TIMI 0-1 or <=1."
)
definition_results <- bind_rows(definition_results, baseline_timi_only, baseline_strict_01)
write_csv_utf8(
  definition_results %>% select(-model_object, -input_data),
  file.path(output_root, "DATA", "definition_sensitivity_results.csv")
)

all_meta_results <- bind_rows(primary_results, adjusted_results, definition_results, eligibility_results)



# ============================================================
# lib/08_risk_of_bias.R
# ============================================================
rob <- support$RiskOfBias %>%
  transmute(
    Study_No = as.integer(Study_No), Study,
    Selection, Confounding,
    `Outcome ascertainment`, `Missing data`, `Selective analysis/reporting`, Overall,
    `Key factual basis`
  ) %>% filter(Study_No != 33)

rob_long <- rob %>%
  pivot_longer(cols = c(Selection, Confounding, `Outcome ascertainment`, `Missing data`, `Selective analysis/reporting`, Overall),
               names_to = "domain", values_to = "judgment") %>%
  mutate(
    judgment = case_when(
      grepl("high", judgment, ignore.case = TRUE) ~ "High concern",
      grepl("some", judgment, ignore.case = TRUE) ~ "Some concerns",
      TRUE ~ "Low concern"
    ),
    judgment = factor(judgment, levels = c("Low concern", "Some concerns", "High concern"))
  )

rob_summary <- rob_long %>% count(domain, judgment) %>% group_by(domain) %>% mutate(percent = n / sum(n) * 100) %>% ungroup()

write_csv_utf8(rob_long, file.path(output_root, "DATA", "rob_domain_long.csv"))
write_csv_utf8(rob_summary, file.path(output_root, "DATA", "rob_summary.csv"))


# ============================================================
# lib/09_grade_inputs.R
# ============================================================
grade_final <- source_data$GradeData
grade_change_log <- tibble(
  association = character(), prior_rating = character(), revised_rating = character(), reason = character()
)
write_csv_utf8(grade_change_log, file.path(output_root, "VALIDATION", "GRADE_CHANGE_LOG.csv"))
writeLines(c(
  "# GRADE change log",
  "",
  "No certainty rating changed after consistency checks against the submitted GRADE input table."
), file.path(output_root, "VALIDATION", "GRADE_CHANGE_LOG.md"), useBytes = TRUE)


# ============================================================
# 04_rob_grade.R
# ============================================================

rob_long <- rob_long %>% filter(Study_No != 33) %>% mutate(Study = study_label(Study_No))
rob_summary <- rob_long %>%
  count(domain, judgment, name = "n") %>%
  group_by(domain) %>% mutate(percent = 100 * n / sum(n)) %>% ungroup()
write_csv_utf8(rob_long, file.path(output_root, "DATA", "rob_domain_long.csv"))
write_csv_utf8(rob_summary, file.path(output_root, "DATA", "rob_summary.csv"))

result_lookup <- bind_rows(primary_results, adjusted_results, definition_results, eligibility_results)
get_result <- function(factor_name, analysis_name) {
  result_lookup %>% filter(factor == factor_name, analysis == analysis_name) %>% slice(1)
}

grade_final <- source_data$GradeData
grade_final$Association <- ascii_clean(grade_final$Association)
age_main <- get_result("Age", "STEMI-only crude primary")
age_timi <- get_result("Age", "TIMI-only (D1) sensitivity")
age_adj <- get_result("Age", "Adjusted-only sensitivity")
diabetes_main <- get_result("Diabetes mellitus", "Uniform crude primary")
diabetes_adj <- get_result("Diabetes mellitus", "Adjusted-only sensitivity")
hypertension_main <- get_result("Hypertension", "Uniform crude primary")
hypertension_adj <- get_result("Hypertension", "Adjusted-only sensitivity")
killip_main <- get_result("Killip >=II / acute HF", "Strict crude primary")
killip_adj <- get_result("Killip >=II / acute HF", "Adjusted-only sensitivity")
advanced_killip <- get_result("Advanced Killip >=III/IV", "Strict crude")
ht_main <- get_result("High thrombus burden >=4", "Strict crude primary")
ht_adj <- get_result("High thrombus burden >=4", "Adjusted-only sensitivity")
bt_main <- get_result("Impaired baseline TIMI flow", "Broad crude primary")
egfr_adj <- get_result("eGFR", "Adjusted-only")
syntax_adj <- get_result("SYNTAX score", "Adjusted per-point")
stent_main <- get_result("Stent length", "Strict continuous MD")
aspiration_main <- get_result("Thrombus aspiration", "Crude observational association")
direct_main <- get_result("Direct stenting", "Crude observational association")

fmt_num <- function(x, digits) formatC(x, format = "f", digits = digits)
grade_pi_fragment <- function(res, digits = ifelse(res$measure == "MD", 2, 3)) {
  paste0("PI ", fmt_num(res$pi_lower, digits), " to ", fmt_num(res$pi_upper, digits),
         " (", pi_relation(res$measure, res$pi_lower, res$pi_upper), ")")
}
grade_result_text <- function(res, label, digits = ifelse(res$measure == "MD", 2, 3)) {
  paste0(label, " ", res$measure, " ", fmt_num(res$estimate, digits),
         " (95% CI ", fmt_num(res$ci_lower, digits), " to ", fmt_num(res$ci_upper, digits), "); ",
         grade_pi_fragment(res, digits))
}

grade_final <- grade_final %>%
  mutate(
    `Effect / finding` = case_when(
      Association == "Age" ~ paste(
        grade_result_text(age_main, "Primary"),
        grade_result_text(age_timi, "TIMI-only"),
        grade_result_text(age_adj, "Adjusted per-year"), sep = "; "),
      Association == "Diabetes" ~ paste(
        grade_result_text(diabetes_main, "Crude"),
        grade_result_text(diabetes_adj, "Adjusted"), sep = "; "),
      Association == "Hypertension" ~ paste(
        grade_result_text(hypertension_main, "Crude"),
        grade_result_text(hypertension_adj, "Adjusted"), sep = "; "),
      Association == "Killip >=II / acute HF" ~ paste(
        grade_result_text(killip_main, "Crude"),
        grade_result_text(killip_adj, "Adjusted"), sep = "; "),
      Association == "Advanced Killip >=III/IV" ~ grade_result_text(advanced_killip, "Crude"),
      grepl("High thrombus", Association) ~ paste(
        grade_result_text(ht_main, "Crude"),
        grade_result_text(ht_adj, "Adjusted"), sep = "; "),
      grepl("Baseline", Association) ~ paste0(
        grade_result_text(bt_main, "Broad crude"), "; thresholds heterogeneous"),
      Association == "eGFR" ~ grade_result_text(egfr_adj, "Adjusted"),
      Association == "SYNTAX score" ~ grade_result_text(syntax_adj, "Adjusted per-point"),
      Association == "Stent length" ~ grade_result_text(stent_main, "Crude"),
      Association == "Thrombus aspiration" ~ grade_result_text(aspiration_main, "Crude observational"),
      Association == "Direct stenting" ~ grade_result_text(direct_main, "Crude observational"),
      TRUE ~ `Effect / finding`
    ),
    Certainty = case_when(
      grepl("Diabetes|High thrombus", Association) ~ "Low",
      TRUE ~ Certainty
    )
  )
write_csv_utf8(grade_final, file.path(output_root, "DATA", "grade_summary_final.csv"))

grade_summary_row <- function(association, role, res, digits = ifelse(res$measure == "MD", 2, 3)) {
  tibble(
    Association = association, result_role = role, factor = res$factor, analysis = res$analysis,
    measure = res$measure, pi_lower = res$pi_lower, pi_upper = res$pi_upper,
    expected_relation = pi_relation(res$measure, res$pi_lower, res$pi_upper),
    expected_fragment = grade_pi_fragment(res, digits)
  )
}
grade_pi_summary <- bind_rows(
  grade_summary_row("Age", "primary", age_main),
  grade_summary_row("Age", "TIMI-only", age_timi),
  grade_summary_row("Age", "adjusted", age_adj),
  grade_summary_row("Diabetes", "crude", diabetes_main),
  grade_summary_row("Diabetes", "adjusted", diabetes_adj),
  grade_summary_row("Hypertension", "crude", hypertension_main),
  grade_summary_row("Hypertension", "adjusted", hypertension_adj),
  grade_summary_row("Killip >=II / acute HF", "crude", killip_main),
  grade_summary_row("Killip >=II / acute HF", "adjusted", killip_adj),
  grade_summary_row("Advanced Killip >=III/IV", "crude", advanced_killip),
  grade_summary_row("High thrombus burden", "crude", ht_main),
  grade_summary_row("High thrombus burden", "adjusted", ht_adj),
  grade_summary_row("Baseline impaired TIMI flow", "broad crude", bt_main),
  grade_summary_row("eGFR", "adjusted", egfr_adj),
  grade_summary_row("SYNTAX score", "adjusted", syntax_adj),
  grade_summary_row("Stent length", "crude", stent_main),
  grade_summary_row("Thrombus aspiration", "crude observational", aspiration_main),
  grade_summary_row("Direct stenting", "crude observational", direct_main)
) %>%
  left_join(grade_final %>% select(Association, grade_wording = `Effect / finding`), by = "Association") %>%
  rowwise() %>%
  mutate(text_matches_pi = grepl(expected_fragment, grade_wording, fixed = TRUE)) %>%
  ungroup()
write_csv_utf8(grade_pi_summary, file.path(output_root, "VALIDATION", "grade_pi_consistency_summary.csv"))


# ============================================================
# 05_figures.R
# ============================================================
fig_dir <- file.path(output_root, "FIGURES")

classify_analysis <- function(x) case_when(
  grepl("exploratory", x, ignore.case = TRUE) ~ "Exploratory",
  grepl("sensitivity|adjusted|TIMI|excluding", x, ignore.case = TRUE) ~ "Sensitivity",
  grepl("primary|Strict crude$|Broad crude", x, ignore.case = TRUE) ~ "Primary",
  TRUE ~ "Exploratory"
)
analysis_colours <- c(Primary = "#24557A", Sensitivity = "#2B8C7E", Exploratory = "#777777")
symbol_ge <- intToUtf8(8805)
symbol_le <- intToUtf8(8804)
symbol_sup2 <- intToUtf8(178)
symbol_endash <- intToUtf8(8211)
display_text <- function(x) {
  x <- gsub(">=", symbol_ge, x, fixed = TRUE)
  x <- gsub("<=", symbol_le, x, fixed = TRUE)
  x <- gsub("_", " ", x, fixed = TRUE)
  x <- gsub(symbol_ge, paste0(" ", symbol_ge, " "), x, fixed = TRUE)
  x <- gsub(symbol_le, paste0(" ", symbol_le, " "), x, fixed = TRUE)
  x <- gsub("  +", " ", x)
  x
}

forest_block <- function(dat, title, measure = unique(dat$measure), show_legend = FALSE) {
  stopifnot(length(measure) == 1)
  dat <- dat %>%
    mutate(
      analysis_class = factor(classify_analysis(analysis), levels = names(analysis_colours)),
      row_label = display_text(paste0(factor, "\n", analysis)),
      row_label = factor(row_label, levels = rev(unique(row_label))),
      effect_text = sprintf("%.2f (%.2f%s%.2f)", estimate, ci_lower, symbol_endash, ci_upper),
      stat_text = sprintf("k = %d   I%s = %.1f%%", k, symbol_sup2, I2)
    )
  left <- ggplot(dat, aes(y = row_label, x = 0, colour = analysis_class)) +
    geom_text(aes(label = as.character(row_label)), hjust = 0, lineheight = 0.93, size = 2.9) +
    scale_colour_manual(values = analysis_colours, drop = FALSE) +
    xlim(0, 1) + labs(title = title) + theme_void(base_family = "Arial", base_size = 8) +
    theme(plot.title = element_text(face = "bold", size = 8.5, margin = margin(b = 5)),
          legend.position = "none", plot.margin = margin(2, 2, 2, 2))
  mid <- ggplot(dat, aes(y = row_label, x = estimate, xmin = ci_lower, xmax = ci_upper,
                         colour = analysis_class)) +
    geom_vline(xintercept = ifelse(measure == "OR", 1, 0), colour = "#999999", linetype = 2, linewidth = 0.35) +
    geom_errorbar(orientation = "y", width = 0, linewidth = 0.7) +
    geom_point(shape = 21, fill = "white", stroke = 0.75, size = 2.3) +
    scale_colour_manual(values = analysis_colours, drop = FALSE) +
    labs(x = ifelse(measure == "OR", "Odds ratio (95% CI)", "Mean difference (95% CI)"), y = NULL) +
    theme_meta(8) + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
                          axis.line.y = element_blank(), legend.position = "none",
                          plot.margin = margin(2, 3, 2, 3))
  if (show_legend) {
    legend_dat <- tibble(
      analysis_class = factor(names(analysis_colours), levels = names(analysis_colours)),
      estimate = Inf,
      row_label = dat$row_label[1]
    )
    mid <- mid +
      geom_point(
        data = legend_dat,
        aes(x = estimate, y = row_label, colour = analysis_class),
        inherit.aes = FALSE, alpha = 0, show.legend = TRUE
      ) +
      guides(colour = guide_legend(
      title = NULL,
      override.aes = list(
        colour = unname(analysis_colours), shape = 21, fill = "white",
        alpha = 1, linetype = 1, linewidth = 0.8, size = 2.3
      )
    )) +
      theme(legend.position = "top", legend.direction = "horizontal", legend.key.width = unit(7, "mm"))
  }
  if (measure == "OR") mid <- mid + scale_x_log10()
  right <- ggplot(dat, aes(y = row_label, x = 0)) +
    geom_text(aes(label = effect_text), hjust = 0, size = 2.85, family = "Arial") +
    geom_text(aes(x = 1, label = stat_text), hjust = 0, size = 2.85, family = "Arial", colour = "#444444") +
    xlim(0, 2.1) + theme_void(base_family = "Arial", base_size = 8) +
    theme(plot.margin = margin(2, 2, 2, 2))
  (left | mid | right) + plot_layout(widths = c(2.6, 1.8, 2.2))
}

# Figure 1: PRISMA 2020-style study-selection flow.
prisma_nodes <- tribble(
  ~x, ~y, ~w, ~h, ~label,
  0, 10.2, 2.65, 0.78, "Records identified\n(n = 1,278)",
  3.8, 10.2, 2.65, 0.78, "Duplicates removed\n(n = 853)",
  0, 8.5, 2.65, 0.78, "Records screened\n(n = 425)",
  3.8, 8.5, 2.65, 0.78, "Records excluded\n(n = 311)",
  0, 6.6, 2.65, 0.86, "Full-text reports assessed\nfor eligibility (n = 114)",
  3.8, 6.6, 4.35, 1.72, paste0("Full-text reports excluded (n = 71)\n",
    "Non-primary PCI (n = 53)\nAnimal studies (n = 17)\n",
    "Mixed STEMI/NSTEMI population without\nseparable STEMI-specific data (n = 1)"),
  0, 4.45, 2.65, 0.90, "Studies included in review\n(n = 43)"
)
prisma_edges <- tribble(
  ~x, ~y, ~xend, ~yend,
  0, 9.81, 0, 8.89,
  1.33, 10.2, 2.47, 10.2,
  1.33, 8.5, 2.47, 8.5,
  0, 8.11, 0, 7.03,
  1.33, 6.6, 1.62, 6.6,
  0, 6.17, 0, 4.90
)
fig1 <- ggplot() +
  geom_segment(data = prisma_edges, aes(x, y, xend = xend, yend = yend),
               arrow = arrow(length = unit(2.1, "mm"), type = "closed"), linewidth = 0.38, colour = "#4A4A4A") +
  geom_rect(data = prisma_nodes, aes(xmin = x-w/2, xmax = x+w/2, ymin = y-h/2, ymax = y+h/2),
            fill = "white", colour = "#4A4A4A", linewidth = 0.45) +
  geom_text(data = prisma_nodes, aes(x, y, label = label), family = "Arial", size = 2.8, lineheight = 1.05) +
  annotate("text", x = -2.05, y = 10.2, label = "Identification", angle = 90, family = "Arial", fontface = "bold", size = 2.8) +
  annotate("text", x = -2.05, y = 8.5, label = "Screening", angle = 90, family = "Arial", fontface = "bold", size = 2.8) +
  annotate("text", x = -2.05, y = 6.6, label = "Eligibility", angle = 90, family = "Arial", fontface = "bold", size = 2.8) +
  annotate("text", x = -2.05, y = 4.45, label = "Included", angle = 90, family = "Arial", fontface = "bold", size = 2.8) +
  coord_cartesian(xlim = c(-2.5, 6.3), ylim = c(3.75, 10.8), clip = "off") +
  labs(title = "Study selection") + theme_void(base_family = "Arial", base_size = 8) +
  theme(plot.title = element_text(face = "bold", size = 10, hjust = 0), plot.margin = margin(6, 6, 6, 6))
save_pub_r(fig1, file.path(fig_dir, "Figure_1_PRISMA"), width_mm = 178, height_mm = 145)

# Figure 2: clinical evidence hierarchy, with OR and MD kept on separate axes.
clinical_binary <- bind_rows(primary_results, adjusted_results) %>%
  filter(factor %in% c("Hypertension", "Diabetes mellitus", "Killip >=II / acute HF", "Advanced Killip >=III/IV"), measure == "OR") %>%
  filter(!(factor == "Hypertension" & analysis == "Adjusted-only sensitivity"))
age_md <- bind_rows(primary_results, definition_results, eligibility_results) %>%
  filter(factor == "Age", measure == "MD", analysis %in% c("STEMI-only crude primary", "TIMI-only (D1) sensitivity", "Eligibility sensitivity excluding #41")) %>%
  mutate(analysis = ifelse(analysis == "Eligibility sensitivity excluding #41",
                           "Delayed-PCI exclusion sensitivity", analysis))
lab_adjusted <- adjusted_results %>% filter(factor %in% c("eGFR", "CRP / hs-CRP"))
p2a <- forest_block(clinical_binary, "A  Clinical binary associations", "OR", show_legend = TRUE)
p2b <- forest_block(age_md, "B  Age", "MD")
p2c <- forest_block(lab_adjusted, "C  Laboratory adjusted evidence", "OR")
fig2 <- p2a / p2b / p2c + plot_layout(heights = c(1.45, 0.9, 0.75)) +
  plot_annotation(title = "Clinical and laboratory associations") &
  theme(plot.title = element_text(face = "bold", size = 10))
save_pub_r(fig2, file.path(fig_dir, "Figure_2_Clinical"), width_mm = 180, height_mm = 225)

# Figure 3: angiographic and procedural associations separated by effect scale and causal risk.
angio_or <- bind_rows(primary_results, adjusted_results, definition_results, eligibility_results) %>%
  filter(factor %in% c("High thrombus burden >=4", "Impaired baseline TIMI flow", "Multivessel disease"), measure == "OR") %>%
  filter(!(factor == "Impaired baseline TIMI flow" & analysis == "Eligibility sensitivity excluding #41"))
procedural_or <- primary_results %>% filter(factor %in% c("Thrombus aspiration", "Direct stenting"), measure == "OR")
stent_md <- primary_results %>% filter(factor == "Stent length", measure == "MD")
p3a <- forest_block(angio_or, "A  Angiographic associations", "OR", show_legend = TRUE)
p3b <- forest_block(stent_md, "B  Stent length", "MD")
p3c <- forest_block(procedural_or, "C  Procedural associations (observational)", "OR")
fig3 <- p3a / p3b / p3c + plot_layout(heights = c(1.5, 0.55, 0.7)) +
  plot_annotation(title = "Angiographic and procedural associations") &
  theme(plot.title = element_text(face = "bold", size = 10))
save_pub_r(fig3, file.path(fig_dir, "Figure_3_Angiographic_Procedural"), width_mm = 180, height_mm = 210)

# Figure 4: final included set only, author-year labels, robvis-style traffic light.
rob_colours <- c("Low concern" = "#4C9F70", "Some concerns" = "#E3B341", "High concern" = "#C44E52")
rob_domain_order <- c("Selection", "Confounding", "Outcome ascertainment", "Missing data", "Selective analysis/reporting", "Overall")
rob_plot_data <- rob_long %>%
  mutate(Study = factor(Study, levels = rev(study_label(sort(unique(Study_No))))),
         domain = factor(domain, levels = rob_domain_order))
fig4 <- ggplot(rob_plot_data, aes(x = domain, y = Study, fill = judgment)) +
  geom_point(shape = 21, size = 3.1, colour = "white", stroke = 0.35) +
  scale_fill_manual(values = rob_colours, drop = FALSE) +
  labs(title = "Risk of bias in the final included studies", x = NULL, y = NULL, fill = NULL) +
  theme_meta(8) +
  theme(axis.text.x = element_text(angle = 32, hjust = 1, size = 8), axis.text.y = element_text(size = 8),
        legend.position = "top", axis.line = element_blank(), axis.ticks = element_blank(),
        panel.grid.major.x = element_line(colour = "#E6E6E6", linewidth = 0.3))
save_pub_r(fig4, file.path(fig_dir, "Figure_4_RoB"), width_mm = 180, height_mm = 205)

# Supplementary age study-level forest.
age_pool <- primary_results %>% filter(factor == "Age", analysis == "STEMI-only crude primary")
age_plot <- age %>% filter(Study_No != 33) %>%
  transmute(Study, estimate = yi, ci_lower = yi - 1.96*sqrt(vi), ci_upper = yi + 1.96*sqrt(vi), type = "Study") %>%
  bind_rows(tibble(Study = "REML-HKSJ pooled", estimate = age_pool$estimate,
                   ci_lower = age_pool$ci_lower, ci_upper = age_pool$ci_upper, type = "Pooled")) %>%
  mutate(Study = factor(Study, levels = rev(Study)))
supp_age <- ggplot(age_plot, aes(y = Study, x = estimate, xmin = ci_lower, xmax = ci_upper, colour = type)) +
  geom_vline(xintercept = 0, colour = "#999999", linetype = 2, linewidth = 0.35) +
  geom_errorbar(orientation = "y", width = 0, linewidth = 0.55) + geom_point(size = 2) +
  scale_colour_manual(values = c(Study = "#777777", Pooled = "#24557A")) +
  labs(title = "Age: study-level and pooled mean differences", x = "Mean difference, years (95% CI)", y = NULL,
       caption = sprintf("Pooled MD %.2f (%.2f to %.2f); I%s %.1f%%; PI %.2f to %.2f.",
                         age_pool$estimate, age_pool$ci_lower, age_pool$ci_upper, symbol_sup2,
                         age_pool$I2, age_pool$pi_lower, age_pool$pi_upper)) +
  theme_meta(8) + theme(axis.text.y = element_text(size = 8), legend.position = "none")
save_pub_r(supp_age, file.path(fig_dir, "Supplement_Age_Sensitivity"), width_mm = 178, height_mm = 205)

supp_adj <- forest_block(adjusted_results %>% filter(factor != "Large reference/target vessel diameter"),
                         "Adjusted-only sensitivity analyses", "OR", show_legend = TRUE)
save_pub_r(supp_adj, file.path(fig_dir, "Supplement_Adjusted_Sensitivity"), width_mm = 180, height_mm = 150)

supp_rob <- ggplot(rob_summary, aes(x = percent, y = factor(domain, levels = rev(rob_domain_order)), fill = judgment)) +
  geom_col(width = 0.64) + scale_fill_manual(values = rob_colours, drop = FALSE) +
  scale_x_continuous(labels = function(x) paste0(x, "%"), limits = c(0,100), expand = c(0,0)) +
  labs(title = "Risk-of-bias summary", x = "Studies", y = NULL, fill = NULL) +
  theme_meta(8) + theme(legend.position = "top")
save_pub_r(supp_rob, file.path(fig_dir, "Supplement_RoB_Summary"), width_mm = 150, height_mm = 95)

figure_source_data <- bind_rows(
  clinical_binary %>% mutate(figure = "Figure 2A"), age_md %>% mutate(figure = "Figure 2B"),
  lab_adjusted %>% mutate(figure = "Figure 2C"), angio_or %>% mutate(figure = "Figure 3A"),
  stent_md %>% mutate(figure = "Figure 3B"), procedural_or %>% mutate(figure = "Figure 3C")
) %>% select(figure, domain, factor, analysis, measure, k, estimate, ci_lower, ci_upper, I2, tau2, pi_lower, pi_upper, study_ids)
write_csv_utf8(figure_source_data, file.path(output_root, "DATA", "figure_source_data.csv"))


# ============================================================
# 06_tables.R
# ============================================================
table_source_dir <- file.path(output_root, "VALIDATION", "table_sources")
dir.create(table_source_dir, recursive = TRUE, showWarnings = FALSE)

table1 <- final_study %>%
  select(Study_No, Study, Country, Design, NR_n, Normal_n, Total_N,
         DOI, outcome_definition, outcome_class, eligibility_note)

table2 <- support$Studies %>%
  transmute(
    Study_No = as.integer(Study_No), Study = study_label(Study_No),
    outcome_definition = Outcome_Definition,
    outcome_class = Outcome_Class,
    mechanical_obstruction_excluded = Mechanical_Obstruction_Excluded,
    primary_PCI_eligibility = Primary_PCI_Eligibility,
    adjusted_model_present = Adjusted_Model_Present,
    final_analytical_role = case_when(
      Study_No == 41 ~ "Included in primary analysis; excluded in eligibility sensitivity analysis",
      TRUE ~ "Included; contributes factor-specific quantitative or descriptive evidence"
    )
  ) %>% filter(Study_No %in% final_study$Study_No)

table3 <- bind_rows(primary_results, definition_results, eligibility_results) %>%
  mutate(
    effect_95CI = format_effect(measure, estimate, ci_lower, ci_upper),
    prediction_interval = format_effect(measure, estimate, pi_lower, pi_upper),
    PI_null_relation = mapply(pi_relation, measure, pi_lower, pi_upper)
  ) %>%
  select(domain, factor, analysis, measure, k, effect_95CI, prediction_interval,
         PI_null_relation, I2_percent = I2, tau2, Q, Q_p, study_ids, interpretation)

table4 <- adjusted_results %>%
  mutate(
    effect_95CI = format_effect(measure, estimate, ci_lower, ci_upper),
    prediction_interval = format_effect(measure, estimate, pi_lower, pi_upper),
    PI_null_relation = mapply(pi_relation, measure, pi_lower, pi_upper)
  ) %>%
  select(domain, factor, analysis, measure, k, effect_95CI, prediction_interval,
         PI_null_relation, I2_percent = I2, tau2, Q, Q_p, study_ids, interpretation)

final_effect_source_data <- effect_source %>%
  mutate(Study = study_label(Study_No)) %>%
  filter(Study_No != 33) %>%
  select(Study_No, Study, Factor, source_type, contrast_unit, effect, ci_lower, ci_upper,
         outcome_class, source_sheet, source_note, included_in_primary, included_in_adjusted_sensitivity)

baseline_input_summary <- baseline_timi_inputs %>%
  transmute(Study_No, Study = study_label(Study_No), threshold = source_note,
            outcome_class, exposed_NR = a, total_NR = nr_total,
            exposed_normal = c, total_normal = normal_total,
            OR = exp(yi), CI_lower = exp(yi - 1.96*sqrt(vi)), CI_upper = exp(yi + 1.96*sqrt(vi)),
            broad_primary = TRUE, timi_flow_only = Study_No != 20,
            strict_TIMI_0_1 = Study_No %in% c(5,21,29,36),
            eligibility_excluding_41 = Study_No != 41)

high_thrombus_adjusted_summary <- adjusted_manual %>%
  filter(Factor == "High thrombus burden >=4") %>%
  transmute(Study_No, Study = study_label(Study_No), exact_adjusted_OR = or,
            CI_lower = lo, CI_upper = hi, contrast_unit, source_note,
            final_decision = "Included in adjusted-only sensitivity",
            exclusion_note = "One incompatible-direction estimate was excluded; one additional study did not use the strict thrombus >=4 definition.")

prediction_intervals <- bind_rows(primary_results, adjusted_results, definition_results, eligibility_results) %>%
  filter(k >= 3) %>%
  mutate(PI_null_relation = mapply(pi_relation, measure, pi_lower, pi_upper)) %>%
  select(domain, factor, analysis, measure, k, estimate, ci_lower, ci_upper, I2, tau2,
         pi_lower, pi_upper, PI_null_relation)

write_csv_utf8(table1, file.path(table_source_dir, "Table_1.csv"))
write_csv_utf8(table2, file.path(table_source_dir, "Table_2.csv"))
write_csv_utf8(table3, file.path(table_source_dir, "Table_3.csv"))
write_csv_utf8(table4, file.path(table_source_dir, "Table_4.csv"))
write_csv_utf8(rob_long, file.path(table_source_dir, "Supplement_RoB.csv"))
write_csv_utf8(grade_final, file.path(table_source_dir, "Supplement_GRADE.csv"))
write_csv_utf8(final_effect_source_data, file.path(table_source_dir, "Effect_Source.csv"))
write_csv_utf8(source_data$DescriptiveEvidence, file.path(table_source_dir, "Nonpooled.csv"))
write_csv_utf8(eligibility_results %>% select(-model_object, -input_data), file.path(table_source_dir, "Eligibility.csv"))
write_csv_utf8(definition_results %>% select(-model_object, -input_data), file.path(table_source_dir, "Outcome_Definition.csv"))
write_csv_utf8(prediction_intervals, file.path(table_source_dir, "Prediction_Intervals.csv"))
write_csv_utf8(small_study_results, file.path(table_source_dir, "Small_Study_Effects.csv"))
write_csv_utf8(age_input_summary, file.path(table_source_dir, "Age_Input_Summary.csv"))
write_csv_utf8(baseline_input_summary, file.path(table_source_dir, "Baseline_TIMI_Summary.csv"))
write_csv_utf8(high_thrombus_adjusted_summary, file.path(table_source_dir, "Thrombus_Adjusted_Summary.csv"))
write_csv_utf8(grade_pi_summary, file.path(table_source_dir, "PI_Text_Consistency.csv"))

write_csv_utf8(final_effect_source_data, file.path(output_root, "DATA", "final_effect_source_data.csv"))


# ============================================================
# 07_validate.R
# ============================================================
pi_check <- bind_rows(primary_results, adjusted_results, definition_results, eligibility_results) %>%
  filter(k >= 3) %>%
  mutate(
    expected = mapply(pi_relation, measure, pi_lower, pi_upper),
    observed = case_when(
      grepl("includes 1", interpretation, fixed = TRUE) ~ "includes 1",
      grepl("excludes 1", interpretation, fixed = TRUE) ~ "excludes 1",
      grepl("includes 0", interpretation, fixed = TRUE) ~ "includes 0",
      grepl("excludes 0", interpretation, fixed = TRUE) ~ "excludes 0",
      TRUE ~ NA_character_
    )
  )

if (nrow(final_study) != 43) stop("Final study count is not 43.")
if (sum(final_study$Total_N) != 27405) stop("Total sample size is not 27,405.")
if (sum(final_study$NR_n) != 5631) stop("Slow-flow/no-reflow sample size is not 5,631.")
if (sum(final_study$Normal_n) != 21774) stop("Normal-flow sample size is not 21,774.")
if (any(is.na(pi_check$observed)) || any(pi_check$expected != pi_check$observed)) {
  stop("Prediction-interval interpretation is inconsistent with the calculated interval.")
}
message("Validation completed: study counts and prediction-interval wording are consistent.")