# Network meta-analysis: Retzius-sparing vs Ultra-preservation vs Conventional RARP
# PROSPERO CRD420251241363
#
# Publication-ready analysis version: consolidates linked publications into
# unique studies and uses the audited analysis-input CSV. The original source
# extraction file is retained unchanged alongside the correction log.
# The CSV must retain the original 89-column, two-header-row layout.
# Explicitly reported medians with interquartile ranges are converted in the
# validated input before analysis: Luo et al.'s S2 estimator is used for the
# mean when Q1 and Q3 are available, and Wan et al.'s finite-sample estimator
# is used for the SD. Median-range and ambiguous values are retained in the QC
# log and excluded from inverse-variance analyses.
#
# Recommended use in RStudio:
#   1. Save this script beside Supplementary_Data_2_Validated_analysis_input.csv.
#   2. Review qc_flags.csv and Data_corrections_log.csv after the run.
#   3. Run the complete script.

rm(list = ls())
options(stringsAsFactors = FALSE)

# ----------------------------- SETTINGS ------------------------------------

REFERENCE_TREATMENT <- "Conventional"
TREATMENT_ORDER <- c("Conventional", "Retzius", "Ultra")
RUN_ALL_CONDITIONAL_OUTCOMES <- TRUE
RUN_RCT_SENSITIVITY <- TRUE
RUN_RISK_OF_BIAS_SENSITIVITY <- TRUE
RUN_SMALL_STUDY_TESTS <- TRUE
MIN_STUDIES_FOR_EGGER <- 10
SUCRA_SIMULATIONS <- 20000
RANDOM_SEED <- 20260815

# Eligibility exclusions and duplicate-publication exclusions confirmed after
# the study-level audit. Gaona 2022 used nonrobotic laparoscopic radical
# prostatectomy, and Madi 2021 evaluated salvage rather than primary radical
# prostatectomy. For overlapping participant cohorts, only the user-selected
# representative publication is retained: Deng 2021a, Dall 2024, and Dalela
# 2017. The excluded companion publications cannot contribute secondary
# outcomes, which prevents publication-level double counting.
EXCLUDED_STUDY_IDS <- c(
  "Gaona_2022", "Madi_2021",
  "Barayan_2024", "Menon_2018", "Egan_2020", "Deng_a_2021"
)
EXCLUSION_REASONS <- c(
  Gaona_2022 = "Nonrobotic laparoscopic radical prostatectomy",
  Madi_2021 = "Salvage radical prostatectomy population",
  Barayan_2024 = "Duplicate publication of the Dalela randomized cohort; Dalela 2017 retained",
  Menon_2018 = "Duplicate publication of the Dalela randomized cohort; Dalela 2017 retained",
  Egan_2020 = "Overlapping publication of the Dall cohort; Dall 2024 retained",
  Deng_a_2021 = "Nested overlapping Nanchang publication; Deng 2021a retained"
)

# Multiple publications from the same participants are mapped to one study.
# The preferred report is selected separately for each outcome. Secondary
# reports may fill an outcome only when the representative report does not
# provide a usable comparative result; they never enter as independent studies.
STUDY_CLUSTER_MAP <- c(
  Barayan_2024 = "Dalela_Detroit_RCT_2015_2016",
  Dalela_2017 = "Dalela_Detroit_RCT_2015_2016",
  Menon_2018 = "Dalela_Detroit_RCT_2015_2016",
  Dall_2024 = "Kowalczyk_Georgetown_series",
  Egan_2020 = "Kowalczyk_Georgetown_series",
  Deng_2021 = "Nanchang_Deng_2021a",
  Deng_a_2021 = "Nanchang_Deng_2021a",
  Ota_2021 = "Nagoya_RS_anatomy_OAB"
)

# Used only when RUN_ALL_CONDITIONAL_OUTCOMES is FALSE.
PROTOCOL_OUTCOMES <- c(
  "continence_1m", "continence_3m", "continence_12m",
  "psm_overall", "operative_time", "blood_loss",
  "complications_overall", "clavien_1_2", "clavien_ge3"
)

get_script_dir <- function() {
  args <- commandArgs(trailingOnly = FALSE)
  file_arg <- grep("^--file=", args, value = TRUE)
  if (length(file_arg)) {
    return(dirname(normalizePath(sub("^--file=", "", file_arg[[1]]))))
  }
  frame_files <- vapply(sys.frames(), function(frame) {
    value <- frame$ofile
    if (is.null(value)) "" else as.character(value)
  }, character(1))
  frame_files <- frame_files[nzchar(frame_files)]
  if (length(frame_files)) return(dirname(normalizePath(tail(frame_files, 1))))
  getwd()
}

SCRIPT_DIR <- get_script_dir()
INPUT_FILE <- Sys.getenv(
  "NMA_INPUT_CSV",
  unset = file.path(SCRIPT_DIR, "Supplementary_Data_2_Validated_analysis_input.csv")
)
OUTPUT_DIR <- Sys.getenv(
  "NMA_OUTPUT_DIR",
  unset = file.path(SCRIPT_DIR, "results_nma_publication")
)
PRIMARY_FOREST_FILE <- Sys.getenv(
  "NMA_PRIMARY_FOREST",
  unset = file.path(SCRIPT_DIR, "Forest_plots_NMA_publication.pdf")
)
RCT_FOREST_FILE <- Sys.getenv(
  "NMA_RCT_FOREST",
  unset = file.path(SCRIPT_DIR, "Forest_plots_NMA_publication_RCT.pdf")
)
ROB_FILE <- Sys.getenv(
  "NMA_ROB_CSV",
  unset = file.path(SCRIPT_DIR, "risk_of_bias.csv")
)

required_packages <- c(
  "dplyr", "tidyr", "purrr", "stringr", "ggplot2", "meta", "netmeta"
)
missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages)) {
  stop(
    "Missing packages: ", paste(missing_packages, collapse = ", "),
    "\nInstall with:\ninstall.packages(c(",
    paste(sprintf('"%s"', missing_packages), collapse = ", "), "))"
  )
}

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(purrr)
  library(stringr)
  library(ggplot2)
  library(meta)
  library(netmeta)
})

if (!file.exists(INPUT_FILE)) stop("Input CSV not found: ", INPUT_FILE)
dir.create(OUTPUT_DIR, recursive = TRUE, showWarnings = FALSE)
dir.create(dirname(PRIMARY_FOREST_FILE), recursive = TRUE, showWarnings = FALSE)
dir.create(dirname(RCT_FOREST_FILE), recursive = TRUE, showWarnings = FALSE)
set.seed(RANDOM_SEED)

# -------------------------- OUTCOME DICTIONARY ------------------------------

outcomes <- tibble::tribble(
  ~outcome_id, ~outcome, ~protocol_role, ~type, ~timepoint, ~source_col, ~effect_measure, ~small_values, ~unit, ~status,
  "continence_immediate", "Immediate urinary continence", "Exploratory", "binary", "immediate", 41L, "OR", "undesirable", "events/total", "conditional",
  "continence_2w", "Urinary continence at 2 weeks", "Exploratory", "binary", "2 weeks", 42L, "OR", "undesirable", "events/total", "conditional",
  "continence_1m", "Urinary continence at 1 month", "Primary", "binary", "1 month", 43L, "OR", "undesirable", "events/total", "conditional",
  "continence_3m", "Urinary continence at 3 months", "Primary", "binary", "3 months", 44L, "OR", "undesirable", "events/total", "conditional",
  "continence_6m", "Urinary continence at 6 months", "Exploratory", "binary", "6 months", 45L, "OR", "undesirable", "events/total", "conditional",
  "continence_12m", "Urinary continence at 12 months", "Primary", "binary", "12 months", 46L, "OR", "undesirable", "events/total", "conditional",
  "bcr", "Biochemical recurrence", "Secondary", "binary", "not defined", 47L, "OR", "desirable", "events/total", "not_ready",
  "psm_overall", "Overall positive surgical margins", "Primary", "binary", "postoperative", 48L, "OR", "desirable", "events/total", "conditional",
  "potency_12m", "Potency/erectile function", "Not synthesized", "binary", "not harmonized", 51L, "OR", "undesirable", "evaluable events/total", "not_ready",
  "complications_overall", "Overall complications", "Secondary", "binary", "perioperative", 50L, "OR", "desirable", "events/total", "conditional",
  "transfusion", "Blood transfusion", "Exploratory", "binary", "perioperative", 53L, "OR", "desirable", "events/total", "conditional",
  "clavien_1_2", "Clavien-Dindo grade I-II complications", "Secondary", "binary", "perioperative", 71L, "OR", "desirable", "events/total", "conditional",
  "clavien_ge3", "Clavien-Dindo grade >=III complications", "Secondary", "binary", "perioperative", 72L, "OR", "desirable", "events/total", "conditional",
  "operative_time", "Operative time", "Secondary", "continuous", "intraoperative", 39L, "MD", "desirable", "minutes", "conditional",
  "console_time", "Console time", "Exploratory", "continuous", "intraoperative", 49L, "MD", "desirable", "minutes", "conditional",
  "blood_loss", "Estimated blood loss", "Secondary", "continuous", "intraoperative", 52L, "MD", "desirable", "mL (verify)", "conditional",
  "length_of_stay", "Length of hospital stay", "Exploratory", "continuous", "postoperative", 60L, "MD", "desirable", "days", "conditional",
  "time_to_continence", "Time to urinary continence", "Not synthesized", "continuous", "not harmonized", 40L, "MD", "desirable", "mixed reported units", "not_ready"
)

outcome_labels <- setNames(outcomes$outcome, outcomes$outcome_id)
primary_order <- c(
  "continence_1m", "continence_3m", "continence_12m", "psm_overall",
  "complications_overall", "clavien_1_2", "clavien_ge3",
  "operative_time", "blood_loss", "continence_immediate", "continence_2w",
  "continence_6m", "transfusion", "console_time", "length_of_stay"
)
rct_order <- c("continence_3m", "psm_overall", "complications_overall")

# ------------------------------- HELPERS ------------------------------------

clean_text <- function(x) {
  x <- gsub("\u00A0", " ", as.character(x), fixed = TRUE)
  x <- trimws(x)
  x[is.na(x)] <- ""
  x
}

is_missing_text <- function(x) {
  toupper(clean_text(x)) %in% c("", "NR", "NA", "N/A", "NONE", "-", "PENDING", "PENDIENTE")
}

parse_integer <- function(x) {
  x <- clean_text(x)
  if (is_missing_text(x) || !grepl("^[0-9]+$", x)) return(NA_real_)
  as.numeric(x)
}

safe_name <- function(x) {
  x |>
    iconv(to = "ASCII//TRANSLIT") |>
    tolower() |>
    gsub("[^a-z0-9]+", "_", x = _) |>
    gsub("^_|_$", "", x = _)
}

study_id_from <- function(author, year) {
  raw <- paste(clean_text(author), clean_text(year), sep = "_")
  gsub("^_|_$", "", gsub("[^A-Za-z0-9]+", "_", raw))
}

study_cluster_from <- function(study_id) {
  mapped <- unname(STUDY_CLUSTER_MAP[study_id])
  ifelse(is.na(mapped) | !nzchar(mapped), study_id, mapped)
}

report_priority <- function(study_cluster_id, outcome_id, source_report_id) {
  priority <- rep(1L, length(source_report_id))

  detroit <- study_cluster_id == "Dalela_Detroit_RCT_2015_2016"
  priority[detroit] <- match(
    source_report_id[detroit],
    c("Dalela_2017", "Menon_2018", "Barayan_2024"),
    nomatch = 99L
  )
  late_detroit <- detroit & outcome_id %in% c(
    "continence_6m", "continence_12m", "potency_12m"
  )
  priority[late_detroit] <- match(
    source_report_id[late_detroit],
    c("Menon_2018", "Dalela_2017", "Barayan_2024"),
    nomatch = 99L
  )
  long_term_detroit <- detroit & outcome_id == "bcr"
  priority[long_term_detroit] <- match(
    source_report_id[long_term_detroit],
    c("Barayan_2024", "Menon_2018", "Dalela_2017"),
    nomatch = 99L
  )

  georgetown <- study_cluster_id == "Kowalczyk_Georgetown_series"
  priority[georgetown] <- match(
    source_report_id[georgetown], c("Dall_2024", "Egan_2020"), nomatch = 99L
  )

  nanchang <- study_cluster_id == "Nanchang_Deng_2021a"
  priority[nanchang] <- match(
    source_report_id[nanchang], c("Deng_2021", "Deng_a_2021"), nomatch = 99L
  )

  priority
}

consolidate_linked_reports <- function(dat, data_type) {
  dat <- dat |>
    mutate(
      source_report_id = study_id,
      study_cluster_id = study_cluster_from(study_id),
      selection_priority = report_priority(
        study_cluster_id, outcome_id, source_report_id
      )
    )

  report_options <- dat |>
    group_by(study_cluster_id, outcome_id, source_report_id, selection_priority) |>
    summarise(
      usable_comparative = n_distinct(treatment[include_analysis]) >= 2,
      .groups = "drop"
    ) |>
    arrange(
      study_cluster_id, outcome_id,
      desc(usable_comparative), selection_priority, source_report_id
    ) |>
    group_by(study_cluster_id, outcome_id) |>
    slice_head(n = 1) |>
    ungroup() |>
    transmute(
      study_cluster_id,
      outcome_id,
      selected_source_report_id = source_report_id,
      selected_report_usable = usable_comparative,
      data_type = data_type
    )

  consolidated <- dat |>
    inner_join(report_options, by = c("study_cluster_id", "outcome_id")) |>
    filter(source_report_id == selected_source_report_id) |>
    mutate(study_id = study_cluster_id) |>
    select(-study_cluster_id, -selection_priority)

  list(data = consolidated, log = report_options)
}

normalize_treatment <- function(x) {
  x <- clean_text(x)
  if (x == "Retzius") return("Retzius")
  if (x %in% c("Ultra", "Ultrapreservation")) return("Ultra")
  if (startsWith(x, "Conventional")) return("Conventional")
  x
}

classify_design <- function(x) {
  value <- tolower(clean_text(x))
  if (value %in% c("randomized controlled trial", "randomized", "rct")) return("RCT")
  if (value == "prospective") return("Prospective nonrandomized")
  if (value == "retrospective") return("Retrospective")
  "Other/unclear"
}

parse_binary_value <- function(value, arm_n) {
  value <- clean_text(value)
  result <- list(
    events = NA_real_, total = NA_real_, parse_status = "unparsed",
    denominator_source = NA_character_, warning = ""
  )
  if (is_missing_text(value)) {
    result$parse_status <- "missing"
    return(result)
  }
  normalized <- gsub("[\u2013\u2014]", "-", value)
  fraction <- regexec("(?<![0-9.])([0-9]+)\\s*/\\s*([0-9]+)", normalized, perl = TRUE)
  fraction_match <- regmatches(normalized, fraction)[[1]]
  if (length(fraction_match) == 3) {
    result$events <- as.numeric(fraction_match[[2]])
    result$total <- as.numeric(fraction_match[[3]])
    result$parse_status <- "parsed"
    result$denominator_source <- "explicit"
  } else {
    count_match <- regexec("^\\s*([0-9]+)\\s*(?:\\([^)]*%?\\))?\\s*$", normalized, perl = TRUE)
    count_parts <- regmatches(normalized, count_match)[[1]]
    if (length(count_parts) == 2 && !is.na(arm_n)) {
      result$events <- as.numeric(count_parts[[2]])
      result$total <- arm_n
      result$parse_status <- "parsed"
      result$denominator_source <- "arm_population"
    }
  }
  if (result$parse_status == "parsed") {
    if (result$total <= 0 || result$events < 0 || result$events > result$total) {
      result$parse_status <- "invalid"
      result$warning <- "Events outside 0..total"
    } else if (!is.na(arm_n) && result$total > arm_n) {
      result$parse_status <- "invalid"
      result$warning <- "Denominator exceeds arm population"
    } else if (!is.na(arm_n) && result$total < arm_n) {
      result$warning <- "Evaluable denominator is smaller than arm population"
    }
  }
  result
}

parse_mean_sd <- function(value) {
  value <- clean_text(value)
  result <- list(mean = NA_real_, sd = NA_real_, parse_status = "unparsed")
  if (is_missing_text(value)) {
    result$parse_status <- "missing"
    return(result)
  }
  pattern <- "^\\s*(-?[0-9]+(?:[.,][0-9]+)?)\\s*(?:\u00B1|\\+\\s*/\\s*-|\\+/-)\\s*([0-9]+(?:[.,][0-9]+)?)\\s*$"
  match <- regexec(pattern, value, perl = TRUE)
  parts <- regmatches(value, match)[[1]]
  if (length(parts) == 3) {
    result$mean <- as.numeric(sub(",", ".", parts[[2]], fixed = TRUE))
    result$sd <- as.numeric(sub(",", ".", parts[[3]], fixed = TRUE))
    result$parse_status <- ifelse(is.na(result$sd) || result$sd < 0, "invalid", "parsed")
  }
  result
}

network_is_connected <- function(dat) {
  treatments <- sort(unique(dat$treatment))
  if (length(treatments) < 2) return(FALSE)
  study_sets <- split(dat$treatment, dat$study_id)
  adjacency <- setNames(lapply(treatments, function(x) character()), treatments)
  for (trts in study_sets) {
    trts <- unique(trts)
    if (length(trts) < 2) next
    pairs <- combn(trts, 2)
    for (i in seq_len(ncol(pairs))) {
      adjacency[[pairs[1, i]]] <- union(adjacency[[pairs[1, i]]], pairs[2, i])
      adjacency[[pairs[2, i]]] <- union(adjacency[[pairs[2, i]]], pairs[1, i])
    }
  }
  visited <- character()
  queue <- treatments[[1]]
  while (length(queue)) {
    node <- queue[[1]]
    queue <- queue[-1]
    if (node %in% visited) next
    visited <- union(visited, node)
    queue <- union(queue, setdiff(adjacency[[node]], visited))
  }
  setequal(visited, treatments)
}

keep_comparative_studies <- function(dat) {
  dat |>
    group_by(study_id) |>
    filter(n_distinct(treatment) >= 2) |>
    ungroup()
}

prepare_subset <- function(dat, study_ids = NULL) {
  if (!is.null(study_ids)) dat <- filter(dat, study_id %in% study_ids)
  dat <- keep_comparative_studies(dat)
  if (n_distinct(dat$study_id) < 2 || n_distinct(dat$treatment) < 2) return(NULL)
  if (!network_is_connected(dat)) return(NULL)
  dat
}

safe_pdf <- function(filename, draw, width = 9, height = 8) {
  grDevices::pdf(filename, width = width, height = height, onefile = TRUE)
  on.exit(grDevices::dev.off(), add = TRUE)
  tryCatch(draw(), error = function(e) {
    graphics::plot.new()
    graphics::text(0.5, 0.5, paste("Plot could not be generated:", conditionMessage(e)), cex = 0.9)
  })
}

# ----------------------- READ AND PREPARE THE CSV ---------------------------

raw_source <- read.csv(
  INPUT_FILE,
  header = FALSE,
  colClasses = "character",
  check.names = FALSE,
  fill = TRUE,
  quote = "\"",
  comment.char = "",
  na.strings = character(),
  fileEncoding = "UTF-8-BOM"
)

if (nrow(raw_source) < 3 || ncol(raw_source) != 89) {
  stop("Unexpected input shape. Expected 89 columns and at least 3 rows; found ",
       ncol(raw_source), " columns and ", nrow(raw_source), " rows.")
}

source_headers <- clean_text(unlist(raw_source[2, ], use.names = FALSE))
expected_headers <- c(
  `39` = "Operative time, min", `40` = "Time continence (weeks)",
  `41` = "Inmediate continence", `42` = "2 week constinence ",
  `43` = "1 month continence", `44` = "3 months continence",
  `45` = "6 months continence", `46` = "12 months continence",
  `47` = "BCR revisar tiempo", `48` = "surgical margins positive",
  `49` = "Console time, min", `50` = "Overall complications",
  `51` = "sexual potency erection at 12 month pendiente definir",
  `52` = "Blood loss", `53` = "Transfusion required, N ",
  `60` = "Lenght of stay, days", `71` = "Clavien 1-2",
  `72` = "Clavien >3"
)
header_check <- vapply(names(expected_headers), function(index) {
  clean_text(source_headers[[as.integer(index)]]) == clean_text(expected_headers[[index]])
}, logical(1))
if (!all(header_check)) {
  failed <- names(expected_headers)[!header_check]
  stop("Input columns have moved or were renamed at positions: ", paste(failed, collapse = ", "))
}

study_rows <- list()
arm_rows <- list()
current_study <- NULL

for (row_index in 3:nrow(raw_source)) {
  row <- clean_text(unlist(raw_source[row_index, ], use.names = FALSE))
  group_raw <- row[[7]]
  if (nzchar(group_raw)) {
    if (is.null(current_study)) stop("Arm row ", row_index, " has no preceding study header.")
    arm_rows[[length(arm_rows) + 1]] <- c(
      current_study,
      list(
        source_row = row_index,
        group_raw = group_raw,
        treatment = normalize_treatment(group_raw),
        arm_n = parse_integer(row[[8]])
      )
    )
    next
  }
  if (nzchar(row[[6]]) && any(nzchar(row[1:6]))) {
    current_study <- list(
      study_id = study_id_from(row[[6]], row[[3]]),
      author = row[[6]],
      year = suppressWarnings(as.numeric(row[[3]])),
      country = row[[4]],
      design_raw = row[[5]],
      design_class = classify_design(row[[5]]),
      doi = row[[1]],
      journal = row[[2]],
      study_header_row = row_index
    )
    study_rows[[length(study_rows) + 1]] <- current_study
  }
}

arms <- bind_rows(arm_rows)
studies <- bind_rows(study_rows) |>
  distinct(study_id, .keep_all = TRUE)
if (!nrow(arms)) stop("No treatment-arm rows were detected.")
if (any(!arms$treatment %in% TREATMENT_ORDER)) {
  stop("Unrecognized treatment labels: ",
       paste(sort(unique(arms$treatment[!arms$treatment %in% TREATMENT_ORDER])), collapse = ", "))
}

exclusion_log <- data.frame(
  study_id = EXCLUDED_STUDY_IDS,
  author = sub("_[0-9]{4}$", "", EXCLUDED_STUDY_IDS),
  year = as.numeric(sub("^.*_", "", EXCLUDED_STUDY_IDS)),
  reason = unname(EXCLUSION_REASONS[EXCLUDED_STUDY_IDS]),
  present_in_input = EXCLUDED_STUDY_IDS %in% studies$study_id,
  status = "Excluded before outcome preparation and meta-analysis",
  stringsAsFactors = FALSE
)
write.csv(exclusion_log, file.path(OUTPUT_DIR, "eligibility_exclusions.csv"),
          row.names = FALSE, na = "")
arms <- arms |>
  filter(!study_id %in% EXCLUDED_STUDY_IDS)
studies <- studies |>
  filter(!study_id %in% EXCLUDED_STUDY_IDS)

qc_rows <- list()
add_qc <- function(severity, scope, study_id, source_row, field, issue, action) {
  qc_rows[[length(qc_rows) + 1]] <<- data.frame(
    severity = severity, scope = scope, study_id = study_id,
    source_row = as.character(source_row), field = field,
    issue = issue, action_required = action,
    stringsAsFactors = FALSE
  )
}

binary_rows <- list()
continuous_rows <- list()

for (arm_index in seq_len(nrow(arms))) {
  arm <- arms[arm_index, ]
  source_row <- as.integer(arm$source_row[[1]])
  raw_row <- clean_text(unlist(raw_source[source_row, ], use.names = FALSE))
  if (is.na(arm$arm_n[[1]])) {
    add_qc("High", "Arm", arm$study_id[[1]], source_row, "Population size",
           "Arm population is not an integer", "Correct before analysis.")
  }
  for (outcome_index in seq_len(nrow(outcomes))) {
    spec <- outcomes[outcome_index, ]
    value <- raw_row[[spec$source_col[[1]]]]
    if (is_missing_text(value)) next
    common <- list(
      study_id = arm$study_id[[1]], author = arm$author[[1]], year = arm$year[[1]],
      design_class = arm$design_class[[1]], treatment = arm$treatment[[1]],
      group_raw = arm$group_raw[[1]], arm_n = arm$arm_n[[1]],
      outcome_id = spec$outcome_id[[1]], outcome = spec$outcome[[1]],
      protocol_role = spec$protocol_role[[1]], timepoint = spec$timepoint[[1]],
      effect_measure = spec$effect_measure[[1]], small_values = spec$small_values[[1]],
      unit = spec$unit[[1]], value_raw = value, source_col = spec$source_col[[1]],
      source_row = source_row
    )
    if (spec$type[[1]] == "binary") {
      parsed <- parse_binary_value(value, arm$arm_n[[1]])
      include <- identical(parsed$parse_status, "parsed") && spec$status[[1]] == "conditional"
      binary_rows[[length(binary_rows) + 1]] <- c(common, parsed, list(include_analysis = include))
      if (parsed$parse_status %in% c("unparsed", "invalid")) {
        add_qc("High", "Outcome", arm$study_id[[1]], source_row, spec$outcome_id[[1]],
               paste0("Unusable binary value: ", value, if (nzchar(parsed$warning)) paste0(". ", parsed$warning) else ""),
               "Verify events and denominator against the article.")
      } else if (nzchar(parsed$warning)) {
        add_qc("Moderate", "Outcome", arm$study_id[[1]], source_row, spec$outcome_id[[1]],
               paste0(value, ": ", parsed$warning),
               "Document evaluable population and missing-data handling.")
      }
    } else {
      parsed <- parse_mean_sd(value)
      include <- identical(parsed$parse_status, "parsed") && spec$status[[1]] == "conditional"
      continuous_rows[[length(continuous_rows) + 1]] <- c(common, parsed, list(include_analysis = include))
      if (parsed$parse_status == "unparsed") {
        add_qc("Low", "Outcome", arm$study_id[[1]], source_row, spec$outcome_id[[1]],
               paste0("Not an explicit mean +/- SD value: ", value),
               "Retain as reported or convert only under the prespecified conversion method.")
      }
    }
  }
}

binary <- bind_rows(binary_rows)
continuous <- bind_rows(continuous_rows)

# Exclude clinically impossible blood-loss entries without overwriting the
# source file. They remain visible in nma_continuous_prepared.csv and QC flags.
implausible_blood_loss <- continuous |>
  filter(include_analysis, outcome_id == "blood_loss", mean > 5000)
if (nrow(implausible_blood_loss)) {
  for (i in seq_len(nrow(implausible_blood_loss))) {
    flagged <- implausible_blood_loss[i, ]
    continuous <- continuous |>
      mutate(include_analysis = if_else(
        study_id == flagged$study_id[[1]] &
          treatment == flagged$treatment[[1]] &
          outcome_id == "blood_loss",
        FALSE,
        include_analysis
      ))
    add_qc(
      "Low", "Outcome", flagged$study_id[[1]], flagged$source_row[[1]],
      "blood_loss",
      paste0("Clinically impossible blood-loss mean: ", flagged$value_raw[[1]]),
      "Excluded from the validated analysis because the value could not be confirmed from the source report."
    )
  }
}

# A contrast with SD=0 in every reported arm has zero standard error and cannot
# enter an inverse-variance NMA. Keep it in the prepared table but exclude it.
zero_sd_groups <- continuous |>
  filter(include_analysis) |>
  group_by(study_id, outcome_id) |>
  summarise(
    comparative = n_distinct(treatment) >= 2,
    all_zero = all(sd == 0),
    rows = paste(source_row, collapse = ", "),
    .groups = "drop"
  ) |>
  filter(comparative, all_zero)
if (nrow(zero_sd_groups)) {
  for (i in seq_len(nrow(zero_sd_groups))) {
    group <- zero_sd_groups[i, ]
    continuous <- continuous |>
      mutate(include_analysis = if_else(
        study_id == group$study_id[[1]] & outcome_id == group$outcome_id[[1]],
        FALSE,
        include_analysis
      ))
    add_qc("Moderate", "Outcome", group$study_id[[1]], group$rows[[1]], group$outcome_id[[1]],
           "All reported arms have SD=0; the contrast has zero standard error",
           "Verify the SD in the article or exclude this study for this outcome.")
  }
}

# Consolidate companion publications after value-level QC so that report
# selection can fall back to a usable secondary report when the representative
# report lacks a comparative result for that specific outcome.
binary_consolidated <- consolidate_linked_reports(binary, "binary")
continuous_consolidated <- consolidate_linked_reports(continuous, "continuous")
binary <- binary_consolidated$data
continuous <- continuous_consolidated$data
consolidation_log <- bind_rows(
  binary_consolidated$log,
  continuous_consolidated$log
) |>
  filter(study_cluster_id %in% unique(unname(STUDY_CLUSTER_MAP))) |>
  arrange(study_cluster_id, outcome_id, data_type)
write.csv(
  consolidation_log,
  file.path(OUTPUT_DIR, "linked_report_consolidation.csv"),
  row.names = FALSE,
  na = ""
)

for (i in seq_len(nrow(studies))) {
  study <- studies[i, ]
  if (!nzchar(clean_text(study$doi[[1]]))) {
    add_qc("Moderate", "Study", study$study_id[[1]], study$study_header_row[[1]], "DOI",
           "DOI is missing", "Complete the identifier and check for duplicate reports.")
  }
}
add_qc("Low", "Analysis decision", "", "", "PSM pT2/pT3",
       "Stage-stratified positive-margin data were not available in a consistently extractable form.",
       "Not synthesized; overall positive surgical margins were retained.")
add_qc("Low", "Analysis decision", "", "", "BCR",
       "Biochemical-recurrence time horizons and event direction were not consistently reported.",
       "Not synthesized; raw study entries remain available in the source dataset.")
add_qc("Low", "Analysis decision", "", "", "Potency",
       "Definitions and baseline-potent denominators were not harmonized across reports.",
       "Not synthesized; raw study entries remain available in the source dataset.")
add_qc("Low", "Analysis decision", "", "", "Time to continence",
       "Reports used non-comparable time scales and summaries.",
       "Not synthesized; fixed-time continence outcomes were analyzed instead.")
add_qc("Low", "Unit", "", "", "Operative and console time",
       "Validated input headers are normalized to minutes; the original source-file labels are retained unchanged.",
       "Interpret mean differences in minutes.")
add_qc("Moderate", "Design", "", "", "Study design",
       "RCTs and nonrandomized studies are combined in the primary analysis",
       "Justify the combined analysis and prioritize design-specific sensitivity analyses.")

qc <- bind_rows(qc_rows)
write.csv(binary, file.path(OUTPUT_DIR, "nma_binary_prepared.csv"), row.names = FALSE, na = "")
write.csv(continuous, file.path(OUTPUT_DIR, "nma_continuous_prepared.csv"), row.names = FALSE, na = "")
write.csv(qc, file.path(OUTPUT_DIR, "qc_flags.csv"), row.names = FALSE, na = "")
write.csv(filter(qc, tolower(severity) == "high"),
          file.path(OUTPUT_DIR, "high_priority_qc_flags.csv"), row.names = FALSE, na = "")

# ----------------------------- NMA MODELS -----------------------------------

fit_binary_network <- function(dat) {
  pair <- pairwise(
    treat = dat$treatment,
    event = dat$events,
    n = dat$total,
    studlab = dat$study_id,
    sm = "OR",
    incr = 0.5,
    method.incr = "only0",
    allstudies = TRUE
  )
  netmeta(
    pair,
    common = FALSE,
    random = TRUE,
    method.tau = "REML",
    reference.group = REFERENCE_TREATMENT,
    small.values = unique(dat$small_values)[[1]],
    details.chkmultiarm = TRUE,
    prediction = FALSE
  )
}

fit_continuous_network <- function(dat) {
  pair <- pairwise(
    treat = dat$treatment,
    n = dat$arm_n,
    mean = dat$mean,
    sd = dat$sd,
    studlab = dat$study_id,
    sm = "MD"
  )
  netmeta(
    pair,
    common = FALSE,
    random = TRUE,
    method.tau = "REML",
    reference.group = REFERENCE_TREATMENT,
    small.values = unique(dat$small_values)[[1]],
    details.chkmultiarm = TRUE,
    prediction = FALSE
  )
}

write_netsplit <- function(split, filename) {
  if (inherits(split, "try-error") || is.null(split)) return(invisible(NULL))
  random <- split$random |> rename_with(~ paste0("nma_", .x), -comparison)
  direct <- split$direct.random |> rename_with(~ paste0("direct_", .x), -comparison)
  indirect <- split$indirect.random |> rename_with(~ paste0("indirect_", .x), -comparison)
  compare <- split$compare.random |> rename_with(~ paste0("difference_", .x), -comparison)
  out <- random |>
    left_join(direct, by = "comparison") |>
    left_join(indirect, by = "comparison") |>
    left_join(compare, by = "comparison")
  write.csv(out, filename, row.names = FALSE, na = "")
}

forest_direction <- function(outcome_id) {
  if (outcome_id %in% c(
    "continence_immediate", "continence_2w", "continence_1m",
    "continence_3m", "continence_6m", "continence_12m", "potency_12m"
  )) return("OR > 1 favours treatment versus Conventional")
  if (outcome_id %in% c(
    "bcr", "psm_overall", "psm_pt2", "psm_pt3",
    "complications_overall", "transfusion", "clavien_1_2", "clavien_ge3"
  )) return("OR < 1 favours treatment versus Conventional")
  "MD < 0 favours treatment versus Conventional"
}

favours_treatment_on_right <- function(outcome_id) {
  outcome_id %in% c(
    "continence_immediate", "continence_2w", "continence_1m",
    "continence_3m", "continence_6m", "continence_12m", "potency_12m"
  )
}

draw_forest <- function(nma, outcome_id, analysis_label) {
  label <- unname(outcome_labels[[outcome_id]])
  treatment_on_right <- favours_treatment_on_right(outcome_id)
  forest(
    nma,
    reference.group = REFERENCE_TREATMENT,
    pooled = "random",
    sortvar = nma$TE.random[, REFERENCE_TREATMENT],
    smlab = "",
    text.random = "Random-effects NMA",
    label.left = if (treatment_on_right) "Favours Conventional" else "Favours treatment",
    label.right = if (treatment_on_right) "Favours treatment" else "Favours Conventional"
  )
  try({
    grid::upViewport(0)
    grid::grid.text(
      paste0("Outcome: ", label),
      x = grid::unit(0.5, "npc"), y = grid::unit(0.035, "npc"),
      gp = grid::gpar(fontface = "bold", fontsize = 10, col = "#173B57")
    )
    grid::grid.text(
      paste(analysis_label, "|", forest_direction(outcome_id)),
      x = grid::unit(0.5, "npc"), y = grid::unit(0.015, "npc"),
      gp = grid::gpar(fontsize = 7.5, col = "#66727A")
    )
  }, silent = TRUE)
}

write_model_outputs <- function(nma, dat, outcome_dir, analysis_label, outcome_id) {
  dir.create(outcome_dir, recursive = TRUE, showWarnings = FALSE)
  order_present <- TREATMENT_ORDER[TREATMENT_ORDER %in% unique(dat$treatment)]
  if (!length(order_present)) order_present <- sort(unique(dat$treatment))
  
  ranking_p <- netrank(nma, method = "P-score")
  ranking_sucra <- netrank(nma, method = "SUCRA", nsim = SUCRA_SIMULATIONS)
  ranking <- data.frame(
    treatment = names(ranking_p$ranking.random),
    p_score = as.numeric(ranking_p$ranking.random),
    sucra = as.numeric(ranking_sucra$ranking.random)
  ) |>
    arrange(desc(p_score))
  write.csv(ranking, file.path(outcome_dir, paste0(analysis_label, "_ranking.csv")), row.names = FALSE)
  
  league <- netleague(nma, common = FALSE, random = TRUE, backtransf = TRUE)
  league_table <- as.data.frame(league$random, stringsAsFactors = FALSE)
  league_table <- cbind(treatment = rownames(league_table), league_table)
  rownames(league_table) <- NULL
  write.csv(league_table, file.path(outcome_dir, paste0(analysis_label, "_league_table.csv")),
            row.names = FALSE, na = "")
  
  heterogeneity <- data.frame(
    analysis = analysis_label,
    studies = n_distinct(dat$study_id),
    treatments = n_distinct(dat$treatment),
    tau = nma$tau,
    tau2 = nma$tau2,
    I2 = nma$I2,
    Q = nma$Q,
    df_Q = nma$df.Q,
    p_Q = nma$pval.Q
  )
  write.csv(heterogeneity, file.path(outcome_dir, paste0(analysis_label, "_heterogeneity.csv")),
            row.names = FALSE, na = "")
  
  design_decomp <- try(decomp.design(nma), silent = TRUE)
  if (!inherits(design_decomp, "try-error")) {
    write.csv(design_decomp$Q.decomp,
              file.path(outcome_dir, paste0(analysis_label, "_global_inconsistency.csv")),
              row.names = TRUE, na = "")
    write.csv(design_decomp$Q.het.design,
              file.path(outcome_dir, paste0(analysis_label, "_design_heterogeneity.csv")),
              row.names = TRUE, na = "")
  }
  
  split <- try(netsplit(nma, common = FALSE, random = TRUE, show = "all"), silent = TRUE)
  write_netsplit(split, file.path(outcome_dir, paste0(analysis_label, "_node_splitting.csv")))
  
  report_file <- file.path(outcome_dir, paste0(analysis_label, "_report.txt"))
  capture.output({
    cat("OUTCOME:", outcome_labels[[outcome_id]], "\n")
    cat("ANALYSIS:", analysis_label, "\n")
    cat("STUDIES:", n_distinct(dat$study_id), "\n")
    cat("TREATMENTS:", paste(sort(unique(dat$treatment)), collapse = ", "), "\n\n")
    print(summary(nma))
    cat("\n--- P-SCORE RANKING ---\n")
    print(ranking_p)
    cat("\n--- SUCRA RANKING ---\n")
    print(ranking_sucra)
    cat("\n--- GLOBAL INCONSISTENCY: DESIGN-BY-TREATMENT ---\n")
    print(design_decomp)
    cat("\n--- LOCAL INCONSISTENCY: NODE SPLITTING/SIDE ---\n")
    print(split)
  }, file = report_file)
  
  safe_pdf(
    file.path(outcome_dir, paste0(analysis_label, "_network.pdf")),
    function() netgraph(
      nma,
      thickness = "number.of.studies",
      number.of.studies = TRUE,
      points = TRUE,
      cex.points = 5,
      cex = 1.1,
      plastic = FALSE,
      multiarm = TRUE
    ),
    width = 9,
    height = 7
  )
  safe_pdf(
    file.path(outcome_dir, paste0(analysis_label, "_forest.pdf")),
    function() draw_forest(nma, outcome_id, analysis_label),
    width = 9,
    height = 8
  )
  if (!inherits(split, "try-error")) {
    safe_pdf(
      file.path(outcome_dir, paste0(analysis_label, "_node_splitting.pdf")),
      function() forest(split, common = FALSE, random = TRUE),
      width = 10,
      height = 8
    )
  }
  
  n_studies <- n_distinct(dat$study_id)
  if (RUN_SMALL_STUDY_TESTS && n_studies >= MIN_STUDIES_FOR_EGGER) {
    bias <- try(
      metabias(nma, order = order_present, pooled = "random", method.bias = "Egger"),
      silent = TRUE
    )
    capture.output(print(bias),
                   file = file.path(outcome_dir, paste0(analysis_label, "_egger_test.txt")))
    safe_pdf(
      file.path(outcome_dir, paste0(analysis_label, "_comparison_adjusted_funnel.pdf")),
      function() funnel(
        nma,
        order = order_present,
        pooled = "random",
        method.bias = "Egger",
        backtransf = FALSE
      ),
      width = 9,
      height = 7
    )
  } else {
    stale_small_study_outputs <- c(
      file.path(outcome_dir, paste0(analysis_label, "_egger_test.txt")),
      file.path(outcome_dir, paste0(analysis_label, "_comparison_adjusted_funnel.pdf"))
    )
    unlink(stale_small_study_outputs[file.exists(stale_small_study_outputs)])
    writeLines(
      paste0(
        "Egger/comparison-adjusted funnel not run: ", n_studies,
        " studies; prespecified threshold = ", MIN_STUDIES_FOR_EGGER, "."
      ),
      file.path(outcome_dir, paste0(analysis_label, "_small_study_tests_not_run.txt"))
    )
  }
  
  invisible(list(model = nma, data = dat, split = split))
}

risk_map <- data.frame(study_id = character(), overall = character())
if (file.exists(ROB_FILE)) {
  candidate <- read.csv(ROB_FILE, stringsAsFactors = FALSE, check.names = FALSE)
  names(candidate) <- safe_name(names(candidate))
  if (all(c("study_id", "overall") %in% names(candidate))) {
    risk_map <- candidate |>
      transmute(study_id = clean_text(study_id), overall = clean_text(overall)) |>
      filter(nzchar(study_id), nzchar(overall)) |>
      distinct(study_id, .keep_all = TRUE)
  } else {
    warning("risk_of_bias.csv was found but must contain study_id and overall columns.")
  }
}
if (!nrow(risk_map)) {
  writeLines(
    paste0(
      "Risk-of-bias sensitivity was not run because no completed risk_of_bias.csv was found.\n",
      "Expected columns: study_id, overall. Excluded judgements: High or Critical."
    ),
    file.path(OUTPUT_DIR, "risk_of_bias_not_completed.txt")
  )
}

run_outcome <- function(outcome_id, dat, type) {
  outcome_dir <- file.path(OUTPUT_DIR, safe_name(outcome_id))
  dir.create(outcome_dir, recursive = TRUE, showWarnings = FALSE)
  base <- dat |>
    filter(.data$outcome_id == .env$outcome_id, include_analysis)
  if (type == "binary") {
    base <- base |>
      filter(!is.na(events), !is.na(total), events >= 0, total > 0, events <= total)
  } else {
    base <- base |>
      filter(!is.na(arm_n), !is.na(mean), !is.na(sd), arm_n > 1, sd >= 0) |>
      group_by(study_id) |>
      filter(any(sd > 0)) |>
      ungroup()
  }
  base <- prepare_subset(base)
  if (is.null(base)) stop("Insufficient data or disconnected network")
  
  fitter <- if (type == "binary") fit_binary_network else fit_continuous_network
  primary <- fitter(base)
  primary_result <- write_model_outputs(
    primary, base, outcome_dir, "primary_all_designs", outcome_id
  )
  
  rct_result <- NULL
  sensitivity_notes <- character()
  if (RUN_RCT_SENSITIVITY) {
    rct <- prepare_subset(filter(base, design_class == "RCT"))
    if (!is.null(rct)) {
      rct_result <- tryCatch({
        rct_model <- fitter(rct)
        write_model_outputs(
          rct_model, rct, outcome_dir, "sensitivity_rct_only", outcome_id
        )
      }, error = function(e) {
        sensitivity_notes <<- c(
          sensitivity_notes,
          paste0("RCT-only analysis failed: ", conditionMessage(e))
        )
        NULL
      })
    } else {
      sensitivity_notes <- c(
        sensitivity_notes,
        "RCT-only analysis not run: insufficient data or disconnected network."
      )
    }
  }
  
  if (RUN_RISK_OF_BIAS_SENSITIVITY && nrow(risk_map)) {
    acceptable <- risk_map |>
      filter(!tolower(overall) %in% c("high", "critical")) |>
      pull(study_id)
    rob <- prepare_subset(filter(base, study_id %in% acceptable))
    if (!is.null(rob)) {
      tryCatch({
        rob_model <- fitter(rob)
        write_model_outputs(
          rob_model, rob, outcome_dir, "sensitivity_exclude_high_rob", outcome_id
        )
      }, error = function(e) {
        sensitivity_notes <<- c(
          sensitivity_notes,
          paste0("Risk-of-bias sensitivity failed: ", conditionMessage(e))
        )
        NULL
      })
    } else {
      sensitivity_notes <- c(
        sensitivity_notes,
        "Risk-of-bias sensitivity not run: insufficient data or disconnected network."
      )
    }
  }
  if (length(sensitivity_notes)) {
    writeLines(sensitivity_notes, file.path(outcome_dir, "sensitivity_notes.txt"))
  }
  
  list(
    manifest = data.frame(
      outcome_id = outcome_id,
      outcome = outcome_labels[[outcome_id]],
      type = type,
      status = "OK",
      studies = n_distinct(base$study_id),
      treatments = paste(sort(unique(base$treatment)), collapse = ", "),
      tau = primary$tau,
      I2 = primary$I2,
      message = "",
      stringsAsFactors = FALSE
    ),
    primary = primary_result,
    rct = rct_result
  )
}

# ------------------------------- RUN ----------------------------------------

binary_ids <- binary |>
  filter(include_analysis) |>
  distinct(outcome_id) |>
  pull(outcome_id)
continuous_ids <- continuous |>
  filter(include_analysis) |>
  distinct(outcome_id) |>
  pull(outcome_id)

if (!RUN_ALL_CONDITIONAL_OUTCOMES) {
  binary_ids <- intersect(binary_ids, PROTOCOL_OUTCOMES)
  continuous_ids <- intersect(continuous_ids, PROTOCOL_OUTCOMES)
}

manifest_rows <- list()
primary_models <- list()
rct_models <- list()

for (outcome_id in binary_ids) {
  message("Analyzing binary outcome: ", outcome_id)
  result <- tryCatch(
    run_outcome(outcome_id, binary, "binary"),
    error = function(e) list(
      manifest = data.frame(
        outcome_id = outcome_id,
        outcome = outcome_labels[[outcome_id]],
        type = "binary",
        status = "ERROR",
        studies = NA_integer_, treatments = NA_character_, tau = NA_real_, I2 = NA_real_,
        message = conditionMessage(e), stringsAsFactors = FALSE
      ),
      primary = NULL,
      rct = NULL
    )
  )
  manifest_rows[[length(manifest_rows) + 1]] <- result$manifest
  if (!is.null(result$primary)) primary_models[[outcome_id]] <- result$primary$model
  if (!is.null(result$rct)) rct_models[[outcome_id]] <- result$rct$model
}

for (outcome_id in continuous_ids) {
  message("Analyzing continuous outcome: ", outcome_id)
  result <- tryCatch(
    run_outcome(outcome_id, continuous, "continuous"),
    error = function(e) list(
      manifest = data.frame(
        outcome_id = outcome_id,
        outcome = outcome_labels[[outcome_id]],
        type = "continuous",
        status = "ERROR",
        studies = NA_integer_, treatments = NA_character_, tau = NA_real_, I2 = NA_real_,
        message = conditionMessage(e), stringsAsFactors = FALSE
      ),
      primary = NULL,
      rct = NULL
    )
  )
  manifest_rows[[length(manifest_rows) + 1]] <- result$manifest
  if (!is.null(result$primary)) primary_models[[outcome_id]] <- result$primary$model
  if (!is.null(result$rct)) rct_models[[outcome_id]] <- result$rct$model
}

manifest <- bind_rows(manifest_rows)
write.csv(manifest, file.path(OUTPUT_DIR, "analysis_manifest.csv"), row.names = FALSE, na = "")

write_combined_forests <- function(models, order, filename, analysis_label) {
  available <- order[order %in% names(models)]
  if (!length(available)) return(FALSE)
  grDevices::pdf(filename, width = 9, height = 8, onefile = TRUE)
  on.exit(grDevices::dev.off(), add = TRUE)
  for (outcome_id in available) {
    draw_forest(models[[outcome_id]], outcome_id, analysis_label)
  }
  TRUE
}

write_combined_forests(
  primary_models,
  primary_order,
  PRIMARY_FOREST_FILE,
  "Primary analysis: all study designs | random effects (REML)"
)
write_combined_forests(
  rct_models,
  rct_order,
  RCT_FOREST_FILE,
  "Sensitivity analysis: RCTs only | random effects (REML)"
)

capture.output({
  cat("Input:", normalizePath(INPUT_FILE), "\n")
  cat("Results:", normalizePath(OUTPUT_DIR), "\n")
  cat("Primary forest:", normalizePath(PRIMARY_FOREST_FILE, mustWork = FALSE), "\n")
  cat("RCT forest:", normalizePath(RCT_FOREST_FILE, mustWork = FALSE), "\n")
  cat("Date:", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z"), "\n")
  cat("Random seed:", RANDOM_SEED, "\n")
  cat("SUCRA simulations:", SUCRA_SIMULATIONS, "\n\n")
  print(sessionInfo())
}, file = file.path(OUTPUT_DIR, "session_info.txt"))

cat("\nAnalysis complete. Results:\n", normalizePath(OUTPUT_DIR), "\n")
print(manifest)