# =============================================================================
# Reproducible statistical analysis for the final 461-patient BIA cohort
# =============================================================================
#
# This script uses base R only. It reads the supplied .xlsx workbook directly,
# validates its structure, reruns the reported fixed-horizon logistic analyses,
# and writes machine-readable analysis tables plus publication-quality TIFF plots.
#
# R version: 4.2 or newer recommended
#
# Command-line use:
#   Rscript Supplementary_Analysis_Code.R \
#     --data Supplementary_Data.xlsx \
#     --output Analysis_Results \
#     --bootstrap 1000 \
#     --seed 20260818
#
# Quick test without figures:
#   Rscript Supplementary_Analysis_Code.R --bootstrap 20 --no-plots
#
# The primary outcome is binary day-28 mortality. Accordingly, all prediction
# models below are binomial logistic regressions; no survival-time model is fit.
# The supplied workbook is the source-record-derived analysis dataset. All reported
# variables are complete for 461 patients; no statistical imputation is performed.
# =============================================================================

options(stringsAsFactors = FALSE, scipen = 999)

fail <- function(...) stop(paste0(...), call. = FALSE)

script_path <- function() {
  args <- commandArgs(trailingOnly = FALSE)
  hit <- grep("^--file=", args, value = TRUE)
  if (!length(hit)) return(NA_character_)
  # Retain the spelling supplied by Rscript. normalizePath() can recode
  # non-ASCII Windows user paths and prevent discovery of adjacent files.
  gsub("\\\\", "/", sub("^--file=", "", hit[1]))
}

parse_arguments <- function(args) {
  out <- list(data = NULL, output = NULL, bootstrap = 1000L,
              seed = 20260818L, plots = TRUE)
  i <- 1L
  while (i <= length(args)) {
    key <- args[i]
    if (key %in% c("--help", "-h")) {
      cat(paste(
        "Usage: Rscript Supplementary_Analysis_Code.R",
        "  [--data Supplementary_Data.xlsx]",
        "  [--output Analysis_Results]",
        "  [--bootstrap 1000] [--seed 20260818] [--no-plots]",
        sep = "\n"
      ), "\n")
      quit(save = "no", status = 0L)
    } else if (key == "--no-plots") {
      out$plots <- FALSE
      i <- i + 1L
      next
    } else if (key %in% c("--data", "--output", "--bootstrap", "--seed")) {
      if (i == length(args)) fail("Missing value after ", key)
      value <- args[i + 1L]
      if (key == "--data") out$data <- value
      if (key == "--output") out$output <- value
      if (key == "--bootstrap") out$bootstrap <- suppressWarnings(as.integer(value))
      if (key == "--seed") out$seed <- suppressWarnings(as.integer(value))
      i <- i + 2L
      next
    } else {
      fail("Unknown argument: ", key)
    }
  }
  if (is.na(out$bootstrap) || out$bootstrap < 1L) fail("--bootstrap must be a positive integer")
  if (is.na(out$seed)) fail("--seed must be an integer")
  out
}

args <- parse_arguments(commandArgs(trailingOnly = TRUE))
this_script <- script_path()
script_dir <- if (is.na(this_script)) getwd() else dirname(this_script)

if (is.null(args$data)) {
  candidates <- c(
    file.path(script_dir, "Supplementary_Data.xlsx"),
    file.path(script_dir, "Supplementary_Data.xlsx"),
    file.path(script_dir, "Supplementary_Data_461.xlsx"),
    file.path(getwd(), "Supplementary_Data.xlsx"),
    file.path(getwd(), "Supplementary_Data_461.xlsx")
  )
  candidates <- unique(candidates[file.exists(candidates)])
  if (!length(candidates)) {
    fail("Data file does not exist. Supply it with --data <path-to-Supplementary_Data.xlsx>.")
  }
  args$data <- candidates[1]
}
if (!file.exists(args$data)) fail("Data file does not exist: ", args$data)
# Preserve the user-supplied path spelling. This also keeps mapped-drive paths
# intact on Windows instead of resolving them to a locale-dependent long path.
data_file <- gsub("\\\\", "/", args$data)

if (is.null(args$output)) args$output <- file.path(script_dir, "Analysis_Results")
dir.create(args$output, recursive = TRUE, showWarnings = FALSE)
if (!dir.exists(args$output)) fail("Unable to create output directory: ", args$output)
output_dir <- gsub("\\\\", "/", args$output)

# -----------------------------------------------------------------------------
# Minimal named-sheet XLSX reader (base R; no external packages required)
# -----------------------------------------------------------------------------
xml_decode <- function(x) {
  x <- gsub("&lt;", "<", x, fixed = TRUE)
  x <- gsub("&gt;", ">", x, fixed = TRUE)
  x <- gsub("&quot;", "\"", x, fixed = TRUE)
  x <- gsub("&apos;", "'", x, fixed = TRUE)
  x <- gsub("&#39;", "'", x, fixed = TRUE)
  x <- gsub("&amp;", "&", x, fixed = TRUE)
  x
}

xml_matches <- function(text, pattern) {
  match_object <- gregexpr(pattern, text, perl = TRUE)
  result <- regmatches(text, match_object)[[1]]
  if (length(result) == 1L && identical(result, character(0))) character(0) else result
}

tag_text <- function(fragment) {
  pieces <- xml_matches(fragment, "<(?:[[:alnum:]_]+:)?t(?: [^>]*)?>.*?</(?:[[:alnum:]_]+:)?t>")
  if (!length(pieces)) return("")
  pieces <- gsub("^<(?:[[:alnum:]_]+:)?t(?: [^>]*)?>|</(?:[[:alnum:]_]+:)?t>$", "", pieces, perl = TRUE)
  xml_decode(paste0(pieces, collapse = ""))
}

column_number <- function(cell_reference) {
  letters <- sub("[0-9]+$", "", cell_reference)
  values <- utf8ToInt(toupper(letters)) - utf8ToInt("A") + 1L
  total <- 0L
  for (value in values) total <- total * 26L + value
  total
}

read_xlsx_analysis_sheet <- function(path, preferred_sheet = "Data") {
  extraction_dir <- tempfile("xlsx-")
  dir.create(extraction_dir)
  on.exit(unlink(extraction_dir, recursive = TRUE, force = TRUE), add = TRUE)
  utils::unzip(path, exdir = extraction_dir)

  workbook_path <- file.path(extraction_dir, "xl", "workbook.xml")
  relationships_path <- file.path(extraction_dir, "xl", "_rels", "workbook.xml.rels")
  if (!file.exists(workbook_path) || !file.exists(relationships_path)) {
    fail("The workbook is missing its worksheet index")
  }
  workbook_xml <- paste(readLines(workbook_path, warn = FALSE, encoding = "UTF-8"), collapse = "")
  relationships_xml <- paste(readLines(relationships_path, warn = FALSE, encoding = "UTF-8"), collapse = "")
  sheet_entries <- xml_matches(workbook_xml, "<(?:[[:alnum:]_]+:)?sheet(?: [^>]*)?/>")
  relationship_entries <- xml_matches(relationships_xml, "<Relationship(?: [^>]*)?/>")
  sheet_names <- vapply(sheet_entries, function(entry) {
    parts <- regmatches(entry, regexec(" name=\"([^\"]+)\"", entry, perl = TRUE))[[1]]
    if (length(parts) >= 2L) xml_decode(parts[2]) else ""
  }, character(1))
  selected <- which(tolower(sheet_names) == tolower(preferred_sheet))
  if (!length(selected)) {
    fail("Worksheet '", preferred_sheet, "' was not found; available sheets: ", paste(sheet_names, collapse = ", "))
  }
  sheet_entry <- sheet_entries[selected[1]]
  id_parts <- regmatches(sheet_entry, regexec(" r:id=\"([^\"]+)\"", sheet_entry, perl = TRUE))[[1]]
  if (length(id_parts) < 2L) fail("Worksheet '", preferred_sheet, "' has no relationship identifier")
  relationship_id <- id_parts[2]
  relationship_entry <- relationship_entries[vapply(relationship_entries, function(entry) {
    parts <- regmatches(entry, regexec(" Id=\"([^\"]+)\"", entry, perl = TRUE))[[1]]
    length(parts) >= 2L && identical(parts[2], relationship_id)
  }, logical(1))]
  if (!length(relationship_entry)) fail("Worksheet relationship was not found for '", preferred_sheet, "'")
  target_parts <- regmatches(relationship_entry[1], regexec(" Target=\"([^\"]+)\"", relationship_entry[1], perl = TRUE))[[1]]
  if (length(target_parts) < 2L) fail("Worksheet target was not found for '", preferred_sheet, "'")
  target <- sub("^/xl/", "", xml_decode(target_parts[2]))
  sheet_path <- file.path(extraction_dir, "xl", target)
  if (!file.exists(sheet_path)) fail("Worksheet file is missing for '", preferred_sheet, "': ", target)

  shared <- character(0)
  shared_path <- file.path(extraction_dir, "xl", "sharedStrings.xml")
  if (file.exists(shared_path)) {
    shared_xml <- paste(readLines(shared_path, warn = FALSE, encoding = "UTF-8"), collapse = "")
    shared_items <- xml_matches(shared_xml, "<(?:[[:alnum:]_]+:)?si(?: [^>]*)?>.*?</(?:[[:alnum:]_]+:)?si>")
    shared <- vapply(shared_items, tag_text, character(1))
  }

  sheet_xml <- paste(readLines(sheet_path, warn = FALSE, encoding = "UTF-8"), collapse = "")
  rows <- xml_matches(sheet_xml, "<(?:[[:alnum:]_]+:)?row(?: [^>]*)?>.*?</(?:[[:alnum:]_]+:)?row>")
  if (length(rows) < 2L) fail("Worksheet '", preferred_sheet, "' has no analyzable data rows")

  parsed_rows <- vector("list", length(rows))
  maximum_column <- 0L
  for (i in seq_along(rows)) {
    cells <- xml_matches(rows[i], "<(?:[[:alnum:]_]+:)?c(?: [^>]*)?(?:>.*?</(?:[[:alnum:]_]+:)?c>|/>)")
    row_values <- list()
    for (cell in cells) {
      reference_match <- regexec(" r=\"([^\"]+)\"", cell, perl = TRUE)
      reference_parts <- regmatches(cell, reference_match)[[1]]
      if (length(reference_parts) < 2L) next
      column <- column_number(reference_parts[2])
      maximum_column <- max(maximum_column, column)

      type_match <- regexec(" t=\"([^\"]+)\"", cell, perl = TRUE)
      type_parts <- regmatches(cell, type_match)[[1]]
      cell_type <- if (length(type_parts) >= 2L) type_parts[2] else "n"
      value_match <- regexec("<(?:[[:alnum:]_]+:)?v>(.*?)</(?:[[:alnum:]_]+:)?v>", cell, perl = TRUE)
      value_parts <- regmatches(cell, value_match)[[1]]
      raw_value <- if (length(value_parts) >= 2L) value_parts[2] else ""

      value <- if (cell_type == "s") {
        index <- suppressWarnings(as.integer(raw_value)) + 1L
        if (is.na(index) || index < 1L || index > length(shared)) "" else shared[index]
      } else if (cell_type == "inlineStr") {
        tag_text(cell)
      } else if (cell_type == "b") {
        ifelse(raw_value == "1", "1", "0")
      } else {
        xml_decode(raw_value)
      }
      row_values[[as.character(column)]] <- value
    }
    parsed_rows[[i]] <- row_values
  }

  matrix_values <- matrix("", nrow = length(parsed_rows), ncol = maximum_column)
  for (i in seq_along(parsed_rows)) {
    if (!length(parsed_rows[[i]])) next
    indices <- as.integer(names(parsed_rows[[i]]))
    matrix_values[i, indices] <- unlist(parsed_rows[[i]], use.names = FALSE)
  }

  headers <- trimws(matrix_values[1, ])
  keep <- nzchar(headers)
  if (anyDuplicated(headers[keep])) fail("Duplicate column names were found in the workbook")
  result <- as.data.frame(matrix_values[-1, keep, drop = FALSE], stringsAsFactors = FALSE,
                          check.names = FALSE)
  names(result) <- headers[keep]

  for (name in names(result)) {
    values <- result[[name]]
    numeric_values <- suppressWarnings(as.numeric(values))
    nonempty <- nzchar(values)
    if (all(!nonempty | !is.na(numeric_values))) result[[name]] <- numeric_values
  }
  result
}

d <- read_xlsx_analysis_sheet(data_file, preferred_sheet = "Data")

required_columns <- c(
  "study_id", "age", "gender", "height", "weight", "bmi", "followup_days",
  "mortality", "smm", "pbf", "ecw_tbw", "bmr", "pha", "smi",
  "hypertension", "diabetes", "ckd", "chf", "mv", "vasopressor", "rrt",
  "temperature", "map", "heart_rate", "respiratory_rate", "fio2", "pao2",
  "aado2", "ph", "sodium", "potassium", "creatinine", "acute_renal_failure",
  "hematocrit", "wbc", "gcs", "aps_score", "age_score",
  "chronic_health_score", "apache", "daily_fluid_intake"
)
missing_columns <- setdiff(required_columns, names(d))
if (length(missing_columns)) fail("Required columns are missing: ", paste(missing_columns, collapse = ", "))
if (ncol(d) != 41L) fail("Expected 41 variables but found ", ncol(d))
if (nrow(d) != 461L) fail("Expected 461 patients but found ", nrow(d))
if (anyDuplicated(d$study_id)) fail("Duplicate study_id values were found")
if (anyNA(d)) fail("Missing cells were found in the analysis dataset")
if (!all(d$mortality %in% c(0, 1))) fail("mortality must be coded 0/1")
if (!all(d$gender %in% c(0, 1))) fail("gender must be coded 0/1")
if (sum(d$mortality == 1) != 92L) fail("Expected 92 deaths but found ", sum(d$mortality == 1))
temperature_points <- function(x) ifelse(x >= 41, 4, ifelse(x >= 39, 3, ifelse(x >= 38.5, 1, ifelse(x >= 36, 0, ifelse(x >= 34, 1, ifelse(x >= 32, 2, ifelse(x >= 30, 3, 4)))))))
map_points <- function(x) ifelse(x >= 160, 4, ifelse(x >= 130, 3, ifelse(x >= 110, 2, ifelse(x >= 70, 0, ifelse(x >= 50, 2, 4)))))
heart_rate_points <- function(x) ifelse(x >= 180, 4, ifelse(x >= 140, 3, ifelse(x >= 110, 2, ifelse(x >= 70, 0, ifelse(x >= 55, 2, ifelse(x >= 40, 3, 4))))))
respiratory_rate_points <- function(x) ifelse(x >= 50, 4, ifelse(x >= 35, 3, ifelse(x >= 25, 1, ifelse(x >= 12, 0, ifelse(x >= 10, 1, ifelse(x >= 6, 2, 4))))))
oxygenation_points <- function(fio2, pao2, aado2) ifelse(fio2 >= 0.5, ifelse(aado2 >= 500, 4, ifelse(aado2 >= 350, 3, ifelse(aado2 >= 200, 2, 0))), ifelse(pao2 > 70, 0, ifelse(pao2 >= 61, 1, ifelse(pao2 >= 55, 3, 4))))
ph_points <- function(x) ifelse(x >= 7.7, 4, ifelse(x >= 7.6, 3, ifelse(x >= 7.5, 1, ifelse(x >= 7.33, 0, ifelse(x >= 7.25, 2, ifelse(x >= 7.15, 3, 4))))))
sodium_points <- function(x) ifelse(x >= 180, 4, ifelse(x >= 160, 3, ifelse(x >= 155, 2, ifelse(x >= 150, 1, ifelse(x >= 130, 0, ifelse(x >= 120, 2, ifelse(x >= 111, 3, 4)))))))
potassium_points <- function(x) ifelse(x >= 7, 4, ifelse(x >= 6, 3, ifelse(x >= 5.5, 1, ifelse(x >= 3.5, 0, ifelse(x >= 3, 1, ifelse(x >= 2.5, 2, 4))))))
creatinine_points <- function(x, acute_renal_failure) {
  base <- ifelse(x >= 3.5, 4, ifelse(x >= 2, 3, ifelse(x >= 1.5, 2, ifelse(x >= 0.6, 0, 2))))
  ifelse(acute_renal_failure == 1, 2 * base, base)
}
hematocrit_points <- function(x) ifelse(x >= 60, 4, ifelse(x >= 50, 2, ifelse(x >= 46, 1, ifelse(x >= 30, 0, ifelse(x >= 20, 2, 4)))))
wbc_points <- function(x) ifelse(x >= 40, 4, ifelse(x >= 20, 2, ifelse(x >= 15, 1, ifelse(x >= 3, 0, ifelse(x >= 1, 2, 4)))))
age_points <- function(x) ifelse(x >= 75, 6, ifelse(x >= 65, 5, ifelse(x >= 55, 3, ifelse(x >= 45, 2, 0))))

aps_recalculated <- temperature_points(d$temperature) + map_points(d$map) +
  heart_rate_points(d$heart_rate) + respiratory_rate_points(d$respiratory_rate) +
  oxygenation_points(d$fio2, d$pao2, d$aado2) + ph_points(d$ph) +
  sodium_points(d$sodium) + potassium_points(d$potassium) +
  creatinine_points(d$creatinine, d$acute_renal_failure) +
  hematocrit_points(d$hematocrit) + wbc_points(d$wbc) + (15 - d$gcs)
age_recalculated <- age_points(d$age)
bmi_recalculated <- d$weight / (d$height / 100)^2
aps_mismatch <- sum(d$aps_score != aps_recalculated)
age_score_mismatch <- sum(d$age_score != age_recalculated)
bmi_mismatch <- sum(abs(d$bmi - bmi_recalculated) > 0.025)
apache_mismatch <- sum(d$apache != d$aps_score + d$age_score + d$chronic_health_score)
if (aps_mismatch != 0L) fail("APACHE II acute physiology score mismatches: ", aps_mismatch)
if (age_score_mismatch != 0L) fail("APACHE II age-point mismatches: ", age_score_mismatch)
if (bmi_mismatch != 0L) fail("BMI formula mismatches: ", bmi_mismatch)
if (apache_mismatch != 0L) fail("APACHE II component arithmetic mismatches: ", apache_mismatch)
if (any(d$apache != round(d$apache))) fail("APACHE II total scores must be integers")

integrity <- data.frame(
  check = c("patients", "deaths", "survivors", "variables", "missing_cells",
            "duplicate_study_ids", "aps_recalculation_mismatches", "age_score_mismatches",
            "apache_component_mismatches", "bmi_formula_mismatches", "noninteger_apache_totals"),
  value = c(nrow(d), sum(d$mortality == 1), sum(d$mortality == 0), ncol(d),
            sum(is.na(d)), anyDuplicated(d$study_id), aps_mismatch, age_score_mismatch,
            apache_mismatch, bmi_mismatch, sum(d$apache != round(d$apache))),
  status = "PASS"
)
write.csv(integrity, file.path(output_dir, "data_integrity.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Statistical helper functions
# -----------------------------------------------------------------------------
z975 <- qnorm(0.975)

delong_components <- function(y, score) {
  positive <- score[y == 1]
  negative <- score[y == 0]
  comparison <- outer(positive, negative, FUN = function(a, b) (a > b) + 0.5 * (a == b))
  list(positive = rowMeans(comparison), negative = colMeans(comparison), auc = mean(comparison))
}

auc_summary <- function(y, score) {
  parts <- delong_components(y, score)
  variance <- var(parts$positive) / length(parts$positive) +
    var(parts$negative) / length(parts$negative)
  standard_error <- sqrt(max(variance, 0))
  list(auc = parts$auc,
       lower = max(0, parts$auc - z975 * standard_error),
       upper = min(1, parts$auc + z975 * standard_error),
       variance = variance)
}

paired_delong <- function(y, first, second) {
  a <- delong_components(y, first)
  b <- delong_components(y, second)
  covariance_positive <- cov(cbind(a$positive, b$positive))
  covariance_negative <- cov(cbind(a$negative, b$negative))
  covariance <- covariance_positive / length(a$positive) +
    covariance_negative / length(a$negative)
  variance <- covariance[1, 1] + covariance[2, 2] - 2 * covariance[1, 2]
  standard_error <- sqrt(max(variance, 0))
  difference <- a$auc - b$auc
  z_value <- if (standard_error > 0) difference / standard_error else 0
  c(difference = difference, ci_low = difference - z975 * standard_error,
    ci_high = difference + z975 * standard_error, z = z_value,
    p_value = if (standard_error > 0) 2 * pnorm(-abs(z_value)) else 1)
}

youden_probability <- function(y, probability) {
  best <- NULL
  for (threshold in sort(unique(probability))) {
    predicted <- probability >= threshold
    tp <- sum(predicted & y == 1)
    fp <- sum(predicted & y == 0)
    tn <- sum(!predicted & y == 0)
    fn <- sum(!predicted & y == 1)
    sensitivity <- tp / (tp + fn)
    specificity <- tn / (tn + fp)
    index <- sensitivity + specificity - 1
    replace <- is.null(best) || index > best$youden_index + 1e-14 ||
      (abs(index - best$youden_index) <= 1e-14 && sensitivity > best$sensitivity + 1e-14) ||
      (abs(index - best$youden_index) <= 1e-14 &&
         abs(sensitivity - best$sensitivity) <= 1e-14 && threshold < best$threshold)
    if (replace) {
      best <- list(threshold = threshold, sensitivity = sensitivity, specificity = specificity,
                   youden_index = index, tp = tp, fp = fp, tn = tn, fn = fn)
    }
  }
  best
}

youden_continuous <- function(y, x, lower_is_risk) {
  oriented_score <- if (lower_is_risk) -x else x
  oriented <- auc_summary(y, oriented_score)
  best <- NULL
  for (threshold in sort(unique(x))) {
    predicted <- if (lower_is_risk) x <= threshold else x >= threshold
    tp <- sum(predicted & y == 1)
    fp <- sum(predicted & y == 0)
    tn <- sum(!predicted & y == 0)
    fn <- sum(!predicted & y == 1)
    sensitivity <- tp / (tp + fn)
    specificity <- tn / (tn + fp)
    index <- sensitivity + specificity - 1
    replace <- is.null(best) || index > best$youden_index + 1e-14 ||
      (abs(index - best$youden_index) <= 1e-14 && sensitivity > best$sensitivity + 1e-14) ||
      (abs(index - best$youden_index) <= 1e-14 &&
         abs(sensitivity - best$sensitivity) <= 1e-14 && threshold < best$threshold)
    if (replace) {
      best <- list(threshold = threshold, sensitivity = sensitivity, specificity = specificity,
                   youden_index = index, tp = tp, fp = fp, tn = tn, fn = fn,
                   auc = oriented$auc, auc_low = oriented$lower, auc_high = oriented$upper)
    }
  }
  best
}

effect_table <- function(model, scales) {
  coefficients <- summary(model)$coefficients
  predictors <- names(scales)
  rows <- lapply(predictors, function(predictor) {
    beta <- coefficients[predictor, "Estimate"]
    standard_error <- coefficients[predictor, "Std. Error"]
    scale_value <- unname(scales[predictor])
    data.frame(
      predictor = predictor, scale = scale_value, beta = beta, standard_error = standard_error,
      odds_ratio = exp(beta * scale_value),
      ci_low = exp((beta - z975 * standard_error) * scale_value),
      ci_high = exp((beta + z975 * standard_error) * scale_value),
      z = coefficients[predictor, "z value"], p_value = coefficients[predictor, "Pr(>|z|)"]
    )
  })
  do.call(rbind, rows)
}

format_mean_sd <- function(x) sprintf("%.1f (%.1f)", mean(x), sd(x))
format_median_iqr <- function(x, digits = 1L) {
  q <- quantile(x, c(0.25, 0.5, 0.75), names = FALSE)
  sprintf(paste0("%.", digits, "f (%.", digits, "f- %.", digits, "f)"), q[2], q[1], q[3])
}
format_binary <- function(x) sprintf("%d (%.1f%%)", sum(x == 1), 100 * mean(x == 1))

continuous_table_row <- function(variable, label, method, digits = 1L) {
  overall <- d[[variable]]
  survivors <- d[d$mortality == 0, variable]
  deaths <- d[d$mortality == 1, variable]
  if (method == "mean_sd") {
    p <- t.test(survivors, deaths)$p.value
    values <- c(format_mean_sd(overall), format_mean_sd(survivors), format_mean_sd(deaths))
    test <- "Welch t-test"
  } else {
    p <- wilcox.test(survivors, deaths, exact = FALSE)$p.value
    values <- c(format_median_iqr(overall, digits), format_median_iqr(survivors, digits),
                format_median_iqr(deaths, digits))
    test <- "Wilcoxon rank-sum"
  }
  data.frame(variable = label, overall = values[1], survivors = values[2],
             non_survivors = values[3], test = test, p_value = p)
}

binary_table_row <- function(variable, label) {
  table_values <- table(d[[variable]], d$mortality)
  expected <- suppressWarnings(chisq.test(table_values, correct = FALSE)$expected)
  if (any(expected < 5)) {
    p <- fisher.test(table_values)$p.value
    test <- "Fisher exact"
  } else {
    p <- chisq.test(table_values, correct = FALSE)$p.value
    test <- "Pearson chi-square"
  }
  data.frame(variable = label, overall = format_binary(d[[variable]]),
             survivors = format_binary(d[d$mortality == 0, variable]),
             non_survivors = format_binary(d[d$mortality == 1, variable]),
             test = test, p_value = p)
}

# -----------------------------------------------------------------------------
# Table 1
# -----------------------------------------------------------------------------
table1 <- do.call(rbind, list(
  continuous_table_row("age", "Age, years", "median_iqr", 1),
  binary_table_row("gender", "Male sex"),
  continuous_table_row("bmi", "BMI, kg/m^2", "mean_sd"),
  continuous_table_row("apache", "APACHE II score", "mean_sd"),
  binary_table_row("hypertension", "Hypertension"),
  binary_table_row("diabetes", "Type 2 diabetes"),
  binary_table_row("ckd", "Chronic kidney disease"),
  binary_table_row("chf", "Chronic heart failure"),
  binary_table_row("mv", "Mechanical ventilation"),
  binary_table_row("vasopressor", "Vasopressor use"),
  binary_table_row("rrt", "Renal replacement therapy"),
  continuous_table_row("daily_fluid_intake", "Daily fluid intake, mL", "median_iqr", 0),
  continuous_table_row("smi", "SMI, kg/m^2", "median_iqr", 2),
  continuous_table_row("pha", "Phase angle, degrees", "mean_sd"),
  continuous_table_row("ecw_tbw", "ECW/TBW", "median_iqr", 3)
))
write.csv(table1, file.path(output_dir, "table1.csv"), row.names = FALSE, fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Logistic models and discrimination
# -----------------------------------------------------------------------------
y <- as.integer(d$mortality)
model_apache <- glm(mortality ~ apache, data = d, family = binomial())
model_bia <- glm(mortality ~ smi + pha + ecw_tbw, data = d, family = binomial())
model_core <- glm(mortality ~ apache + smi + pha + ecw_tbw, data = d, family = binomial())
model_sensitivity <- glm(
  mortality ~ apache + smi + pha + ecw_tbw + age + gender + bmi + ckd + chf + daily_fluid_intake,
  data = d, family = binomial()
)

probability_apache <- fitted(model_apache)
probability_bia <- fitted(model_bia)
probability_core <- fitted(model_core)
probability_sensitivity <- fitted(model_sensitivity)

core_scales <- c(apache = 5, smi = 1, pha = 1, ecw_tbw = 0.05)
sensitivity_scales <- c(apache = 5, smi = 1, pha = 1, ecw_tbw = 0.05,
                        age = 10, gender = 1, bmi = 1, ckd = 1, chf = 1,
                        daily_fluid_intake = 1000)
core_table <- effect_table(model_core, core_scales)
sensitivity_table <- effect_table(model_sensitivity, sensitivity_scales)
write.csv(core_table, file.path(output_dir, "core_model.csv"), row.names = FALSE, fileEncoding = "UTF-8")
write.csv(sensitivity_table, file.path(output_dir, "sensitivity_model.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

performance_row <- function(name, probability) {
  auc <- auc_summary(y, probability)
  data.frame(model = name, auc = auc$auc, ci_low = auc$lower, ci_high = auc$upper,
             brier = mean((probability - y)^2))
}
performance <- do.call(rbind, list(
  performance_row("APACHE II alone", probability_apache),
  performance_row("BIA-only", probability_bia),
  performance_row("Extended APACHE II + BIA", probability_core),
  performance_row("Clinically adjusted sensitivity", probability_sensitivity)
))
write.csv(performance, file.path(output_dir, "model_performance.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

delong_row <- function(first_name, second_name, first_probability, second_probability) {
  result <- paired_delong(y, first_probability, second_probability)
  data.frame(comparison = paste(first_name, "vs", second_name), model_first = first_name,
             model_second = second_name, difference_first_minus_second = result["difference"],
             ci_low = result["ci_low"], ci_high = result["ci_high"], z = result["z"],
             p_value = result["p_value"])
}
delong <- do.call(rbind, list(
  delong_row("Extended APACHE II + BIA", "APACHE II alone", probability_core, probability_apache),
  delong_row("Extended APACHE II + BIA", "BIA-only", probability_core, probability_bia),
  delong_row("BIA-only", "APACHE II alone", probability_bia, probability_apache)
))
write.csv(delong, file.path(output_dir, "delong_comparisons.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Exploratory operating thresholds
# -----------------------------------------------------------------------------
probability_threshold_row <- function(name, probability) {
  threshold <- youden_probability(y, probability)
  auc <- auc_summary(y, probability)
  data.frame(model = name, threshold = threshold$threshold,
             rule = "predicted probability >= threshold",
             sensitivity = threshold$sensitivity, specificity = threshold$specificity,
             youden_index = threshold$youden_index, tp = threshold$tp, fp = threshold$fp,
             tn = threshold$tn, fn = threshold$fn, auc = auc$auc,
             auc_ci_low = auc$lower, auc_ci_high = auc$upper)
}
continuous_threshold_row <- function(name, x, lower_is_risk) {
  threshold <- youden_continuous(y, x, lower_is_risk)
  data.frame(model = name, threshold = threshold$threshold,
             rule = if (lower_is_risk) "value <= threshold" else "value >= threshold",
             sensitivity = threshold$sensitivity, specificity = threshold$specificity,
             youden_index = threshold$youden_index, tp = threshold$tp, fp = threshold$fp,
             tn = threshold$tn, fn = threshold$fn, auc = threshold$auc,
             auc_ci_low = threshold$auc_low, auc_ci_high = threshold$auc_high)
}
thresholds <- do.call(rbind, list(
  continuous_threshold_row("SMI", d$smi, TRUE),
  continuous_threshold_row("Phase angle", d$pha, TRUE),
  continuous_threshold_row("ECW/TBW", d$ecw_tbw, FALSE),
  probability_threshold_row("APACHE II alone", probability_apache),
  probability_threshold_row("BIA-only", probability_bia),
  probability_threshold_row("Extended APACHE II + BIA", probability_core)
))

for (sex_code in c(0, 1)) {
  subset_index <- d$gender == sex_code
  sex_result <- youden_continuous(y[subset_index], d$smi[subset_index], TRUE)
  thresholds <- rbind(thresholds, data.frame(
    model = if (sex_code == 0) "SMI female exploratory" else "SMI male exploratory",
    threshold = sex_result$threshold, rule = "value <= threshold",
    sensitivity = sex_result$sensitivity, specificity = sex_result$specificity,
    youden_index = sex_result$youden_index, tp = sex_result$tp, fp = sex_result$fp,
    tn = sex_result$tn, fn = sex_result$fn, auc = sex_result$auc,
    auc_ci_low = sex_result$auc_low, auc_ci_high = sex_result$auc_high
  ))
}
write.csv(thresholds, file.path(output_dir, "thresholds.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Calibration helpers and nonlinear ECW/TBW sensitivity specification
# -----------------------------------------------------------------------------
calibration_statistics <- function(outcome, probability) {
  clipped <- pmin(pmax(probability, 1e-12), 1 - 1e-12)
  predictor <- qlogis(clipped)
  slope_model <- suppressWarnings(glm(outcome ~ predictor, family = binomial()))
  intercept_model <- suppressWarnings(glm(outcome ~ 1, offset = predictor, family = binomial()))
  c(intercept = unname(coef(intercept_model)[1]), slope = unname(coef(slope_model)[2]))
}

rcs_basis <- function(x, knots) {
  penultimate <- knots[length(knots) - 1L]
  last <- knots[length(knots)]
  span <- last - knots[1]
  columns <- lapply(knots[seq_len(length(knots) - 2L)], function(knot) {
    first <- pmax(x - knot, 0)^3
    second <- pmax(x - penultimate, 0)^3 * (last - knot) / (last - penultimate)
    third <- pmax(x - last, 0)^3 * (penultimate - knot) / (last - penultimate)
    (first - second + third) / span^2
  })
  do.call(cbind, columns)
}

ecw_knots <- as.numeric(quantile(d$ecw_tbw, c(0.05, 0.35, 0.65, 0.95), names = FALSE))
ecw_basis_original <- rcs_basis(d$ecw_tbw, ecw_knots)
ecw_basis_center <- colMeans(ecw_basis_original)
ecw_basis_scale <- apply(ecw_basis_original, 2, sd)
augment_ecw_spline <- function(data) {
  basis <- rcs_basis(data$ecw_tbw, ecw_knots)
  basis <- sweep(basis, 2, ecw_basis_center, "-")
  basis <- sweep(basis, 2, ecw_basis_scale, "/")
  data$rcs_non_linear_1 <- basis[, 1]
  data$rcs_non_linear_2 <- basis[, 2]
  data
}
ecw_augmented <- augment_ecw_spline(d)
model_ecw_nonlinear <- glm(
  mortality ~ apache + smi + pha + ecw_tbw + rcs_non_linear_1 + rcs_non_linear_2,
  data = ecw_augmented, family = binomial()
)
probability_ecw_nonlinear <- fitted(model_ecw_nonlinear)
ecw_lr_statistic <- max(
  0,
  2 * (as.numeric(logLik(model_ecw_nonlinear)) - as.numeric(logLik(model_core)))
)
ecw_lr_p_value <- pchisq(ecw_lr_statistic, df = 2, lower.tail = FALSE)

# -----------------------------------------------------------------------------
# Bootstrap internal validation (full refit in each resample)
# -----------------------------------------------------------------------------
set.seed(args$seed, kind = "Mersenne-Twister")
n <- nrow(d)
auc_optimism <- numeric(args$bootstrap)
brier_optimism <- numeric(args$bootstrap)
intercept_optimism <- numeric(args$bootstrap)
slope_optimism <- numeric(args$bootstrap)
nonlinear_auc_optimism <- numeric(args$bootstrap)
nonlinear_brier_optimism <- numeric(args$bootstrap)
nonlinear_intercept_optimism <- numeric(args$bootstrap)
nonlinear_slope_optimism <- numeric(args$bootstrap)
successful <- 0L
failed <- 0L
nonlinear_successful <- 0L
nonlinear_failed <- 0L
for (iteration in seq_len(args$bootstrap)) {
  indices <- sample.int(n, n, replace = TRUE)
  bootstrap_data <- d[indices, , drop = FALSE]
  if (length(unique(bootstrap_data$mortality)) < 2L) {
    failed <- failed + 1L
    next
  }
  bootstrap_model <- tryCatch(
    suppressWarnings(glm(mortality ~ apache + smi + pha + ecw_tbw,
                         data = bootstrap_data, family = binomial())),
    error = function(e) NULL
  )
  if (is.null(bootstrap_model) || !isTRUE(bootstrap_model$converged)) {
    failed <- failed + 1L
    next
  }
  probability_bootstrap <- fitted(bootstrap_model)
  probability_original <- predict(bootstrap_model, newdata = d, type = "response")
  bootstrap_y <- as.integer(bootstrap_data$mortality)
  successful <- successful + 1L
  auc_optimism[successful] <- auc_summary(bootstrap_y, probability_bootstrap)$auc -
    auc_summary(y, probability_original)$auc
  brier_optimism[successful] <- mean((probability_bootstrap - bootstrap_y)^2) -
    mean((probability_original - y)^2)
  calibration_bootstrap <- calibration_statistics(bootstrap_y, probability_bootstrap)
  calibration_original <- calibration_statistics(y, probability_original)
  intercept_optimism[successful] <- calibration_bootstrap["intercept"] -
    calibration_original["intercept"]
  slope_optimism[successful] <- calibration_bootstrap["slope"] -
    calibration_original["slope"]

  nonlinear_bootstrap_data <- augment_ecw_spline(bootstrap_data)
  nonlinear_bootstrap_model <- tryCatch(
    suppressWarnings(glm(
      mortality ~ apache + smi + pha + ecw_tbw + rcs_non_linear_1 + rcs_non_linear_2,
      data = nonlinear_bootstrap_data, family = binomial()
    )),
    error = function(e) NULL
  )
  if (is.null(nonlinear_bootstrap_model) || !isTRUE(nonlinear_bootstrap_model$converged)) {
    nonlinear_failed <- nonlinear_failed + 1L
    next
  }
  nonlinear_probability_bootstrap <- fitted(nonlinear_bootstrap_model)
  nonlinear_probability_original <- predict(
    nonlinear_bootstrap_model, newdata = ecw_augmented, type = "response"
  )
  nonlinear_successful <- nonlinear_successful + 1L
  nonlinear_auc_optimism[nonlinear_successful] <-
    auc_summary(bootstrap_y, nonlinear_probability_bootstrap)$auc -
    auc_summary(y, nonlinear_probability_original)$auc
  nonlinear_brier_optimism[nonlinear_successful] <-
    mean((nonlinear_probability_bootstrap - bootstrap_y)^2) -
    mean((nonlinear_probability_original - y)^2)
  nonlinear_calibration_bootstrap <- calibration_statistics(
    bootstrap_y, nonlinear_probability_bootstrap
  )
  nonlinear_calibration_original <- calibration_statistics(
    y, nonlinear_probability_original
  )
  nonlinear_intercept_optimism[nonlinear_successful] <-
    nonlinear_calibration_bootstrap["intercept"] -
    nonlinear_calibration_original["intercept"]
  nonlinear_slope_optimism[nonlinear_successful] <-
    nonlinear_calibration_bootstrap["slope"] - nonlinear_calibration_original["slope"]
}
if (successful == 0L) fail("No bootstrap refits succeeded")
if (nonlinear_successful == 0L) fail("No nonlinear bootstrap refits succeeded")
auc_optimism <- auc_optimism[seq_len(successful)]
brier_optimism <- brier_optimism[seq_len(successful)]
intercept_optimism <- intercept_optimism[seq_len(successful)]
slope_optimism <- slope_optimism[seq_len(successful)]
nonlinear_auc_optimism <- nonlinear_auc_optimism[seq_len(nonlinear_successful)]
nonlinear_brier_optimism <- nonlinear_brier_optimism[seq_len(nonlinear_successful)]
nonlinear_intercept_optimism <- nonlinear_intercept_optimism[seq_len(nonlinear_successful)]
nonlinear_slope_optimism <- nonlinear_slope_optimism[seq_len(nonlinear_successful)]
apparent_auc <- performance$auc[performance$model == "Extended APACHE II + BIA"]
apparent_brier <- performance$brier[performance$model == "Extended APACHE II + BIA"]
apparent_calibration <- calibration_statistics(y, probability_core)
corrected_calibration <- c(
  intercept = unname(apparent_calibration["intercept"] - mean(intercept_optimism)),
  slope = unname(apparent_calibration["slope"] - mean(slope_optimism))
)
bootstrap_validation <- data.frame(
  model = "Extended APACHE II + BIA", seed = args$seed,
  requested_refits = args$bootstrap, successful_refits = successful, failed_refits = failed,
  apparent_auc = apparent_auc, mean_auc_optimism = mean(auc_optimism),
  optimism_corrected_auc = apparent_auc - mean(auc_optimism),
  apparent_brier = apparent_brier, mean_brier_optimism = mean(brier_optimism),
  optimism_corrected_brier = apparent_brier - mean(brier_optimism)
)
write.csv(bootstrap_validation, file.path(output_dir, "bootstrap_validation.csv"),
          row.names = FALSE, fileEncoding = "UTF-8")

calibration_validation <- data.frame(
  assessment = c("apparent", "optimism_corrected"),
  intercept = c(apparent_calibration["intercept"], corrected_calibration["intercept"]),
  slope = c(apparent_calibration["slope"], corrected_calibration["slope"]),
  brier = c(apparent_brier, apparent_brier - mean(brier_optimism)),
  bootstrap_refits = c(NA_integer_, successful)
)
write.csv(
  calibration_validation,
  file.path(output_dir, "calibration_validation.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

nonlinear_apparent_auc <- auc_summary(y, probability_ecw_nonlinear)$auc
nonlinear_apparent_brier <- mean((probability_ecw_nonlinear - y)^2)
nonlinear_apparent_calibration <- calibration_statistics(y, probability_ecw_nonlinear)
nonlinear_corrected_calibration <- c(
  intercept = unname(nonlinear_apparent_calibration["intercept"] -
    mean(nonlinear_intercept_optimism)),
  slope = unname(nonlinear_apparent_calibration["slope"] - mean(nonlinear_slope_optimism))
)
nonlinear_ecw_comparison <- data.frame(
  specification = c("linear", "ECW/TBW restricted cubic spline"),
  nonlinear_df = c(0L, 2L),
  lr_chisq = c(NA_real_, ecw_lr_statistic),
  lr_p_value = c(NA_real_, ecw_lr_p_value),
  apparent_auc = c(apparent_auc, nonlinear_apparent_auc),
  optimism_corrected_auc = c(
    apparent_auc - mean(auc_optimism),
    nonlinear_apparent_auc - mean(nonlinear_auc_optimism)
  ),
  apparent_brier = c(apparent_brier, nonlinear_apparent_brier),
  optimism_corrected_brier = c(
    apparent_brier - mean(brier_optimism),
    nonlinear_apparent_brier - mean(nonlinear_brier_optimism)
  ),
  optimism_corrected_calibration_intercept = c(
    corrected_calibration["intercept"], nonlinear_corrected_calibration["intercept"]
  ),
  optimism_corrected_calibration_slope = c(
    corrected_calibration["slope"], nonlinear_corrected_calibration["slope"]
  ),
  bootstrap_refits = c(successful, nonlinear_successful)
)
write.csv(
  nonlinear_ecw_comparison,
  file.path(output_dir, "nonlinear_ecw_comparison.csv"),
  row.names = FALSE,
  fileEncoding = "UTF-8"
)

# -----------------------------------------------------------------------------
# Apparent calibration groups for plotting
# -----------------------------------------------------------------------------
clipped_probability <- pmin(pmax(probability_core, 1e-12), 1 - 1e-12)
linear_predictor <- qlogis(clipped_probability)
calibration_model <- glm(y ~ linear_predictor, family = binomial())
calibration_in_large <- glm(y ~ 1, offset = linear_predictor, family = binomial())
rank_groups <- pmin(10L, pmax(1L, ceiling(rank(clipped_probability, ties.method = "average") / n * 10)))
calibration <- data.frame(
  section = c("metric", "metric", "metric"),
  item = c("calibration_intercept", "calibration_slope", "apparent_brier"),
  group = NA_integer_, n = NA_integer_, mean_predicted = NA_real_, observed_event_rate = NA_real_,
  value = c(coef(calibration_in_large)[1], coef(calibration_model)[2], apparent_brier)
)
for (group in sort(unique(rank_groups))) {
  selected <- rank_groups == group
  calibration <- rbind(calibration, data.frame(
    section = "risk_group", item = "apparent_calibration", group = group,
    n = sum(selected), mean_predicted = mean(clipped_probability[selected]),
    observed_event_rate = mean(y[selected]), value = NA_real_
  ))
}
write.csv(calibration, file.path(output_dir, "calibration.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Correlations and variance inflation factors
# -----------------------------------------------------------------------------
correlation_variables <- c("age", "bmi", "apache", "smi", "pha", "ecw_tbw", "daily_fluid_intake")
correlation <- cor(d[, correlation_variables], method = "pearson")
correlation_output <- data.frame(variable = rownames(correlation), correlation,
                                 row.names = NULL, check.names = FALSE)
write.csv(correlation_output, file.path(output_dir, "correlation_matrix.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

vif_variables <- c("apache", "smi", "pha", "ecw_tbw", "age", "gender", "bmi", "ckd", "chf",
                   "daily_fluid_intake")
inverse_correlation <- solve(cor(d[, vif_variables]))
vif_output <- data.frame(predictor = vif_variables, vif = diag(inverse_correlation))
write.csv(vif_output, file.path(output_dir, "vif.csv"), row.names = FALSE, fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Restricted cubic spline sensitivity tests: four fixed empirical knots
# -----------------------------------------------------------------------------
nonlinearity_specification <- list(
  apache = list(model = model_core), smi = list(model = model_core),
  pha = list(model = model_core), ecw_tbw = list(model = model_core),
  age = list(model = model_sensitivity), bmi = list(model = model_sensitivity),
  daily_fluid_intake = list(model = model_sensitivity)
)
nonlinearity_rows <- lapply(names(nonlinearity_specification), function(variable) {
  base_model <- nonlinearity_specification[[variable]]$model
  knots <- as.numeric(quantile(d[[variable]], c(0.05, 0.35, 0.65, 0.95), names = FALSE))
  if (any(diff(knots) <= 0)) fail("Non-distinct spline knots for ", variable)
  basis <- scale(rcs_basis(d[[variable]], knots))
  augmented <- d
  augmented$rcs_non_linear_1 <- basis[, 1]
  augmented$rcs_non_linear_2 <- basis[, 2]
  full_formula <- update(formula(base_model), . ~ . + rcs_non_linear_1 + rcs_non_linear_2)
  nonlinear_model <- glm(full_formula, data = augmented, family = binomial())
  statistic <- max(0, 2 * (as.numeric(logLik(nonlinear_model)) - as.numeric(logLik(base_model))))
  data.frame(variable = variable, knot_05 = knots[1], knot_35 = knots[2],
             knot_65 = knots[3], knot_95 = knots[4], nonlinear_df = 2,
             likelihood_ratio_chisq = statistic, p_value = pchisq(statistic, df = 2, lower.tail = FALSE),
             converged = nonlinear_model$converged)
})
nonlinearity <- do.call(rbind, nonlinearity_rows)
write.csv(nonlinearity, file.path(output_dir, "nonlinearity_tests.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Exploratory apparent decision-curve analysis in the development cohort
# -----------------------------------------------------------------------------
decision_thresholds <- seq(0.01, 0.50, by = 0.01)
decision_curve_for_model <- function(model_name, probability) {
  do.call(rbind, lapply(decision_thresholds, function(threshold) {
    predicted <- probability >= threshold
    tp <- sum(predicted & y == 1)
    fp <- sum(predicted & y == 0)
    odds <- threshold / (1 - threshold)
    data.frame(
      model = model_name,
      threshold_probability = threshold,
      true_positive = tp,
      false_positive = fp,
      model_net_benefit = tp / n - fp / n * odds,
      treat_all_net_benefit = sum(y == 1) / n - sum(y == 0) / n * odds,
      treat_none_net_benefit = 0
    )
  }))
}
decision_curve <- rbind(
  decision_curve_for_model("APACHE II alone", probability_apache),
  decision_curve_for_model("BIA-only", probability_bia),
  decision_curve_for_model("Extended APACHE II + BIA", probability_core)
)
write.csv(decision_curve, file.path(output_dir, "decision_curve.csv"), row.names = FALSE,
          fileEncoding = "UTF-8")

# -----------------------------------------------------------------------------
# Figures
# -----------------------------------------------------------------------------
roc_coordinates <- function(y, score) {
  thresholds <- c(Inf, sort(unique(score), decreasing = TRUE), -Inf)
  rows <- lapply(thresholds, function(threshold) {
    predicted <- score >= threshold
    data.frame(false_positive_rate = sum(predicted & y == 0) / sum(y == 0),
               true_positive_rate = sum(predicted & y == 1) / sum(y == 1), threshold = threshold)
  })
  do.call(rbind, rows)
}

if (isTRUE(args$plots)) {
  roc_apache <- roc_coordinates(y, probability_apache)
  roc_bia <- roc_coordinates(y, probability_bia)
  roc_core <- roc_coordinates(y, probability_core)
  extended_threshold <- thresholds[thresholds$model == "Extended APACHE II + BIA", ]

  tiff(file.path(output_dir, "Figure_1_ROC.tif"), width = 7, height = 6.5, units = "in",
       res = 600, compression = "lzw")
  par(mar = c(5, 5, 2, 1) + 0.1, las = 1)
  plot(roc_apache$false_positive_rate, roc_apache$true_positive_rate, type = "l", lwd = 2.4,
       col = "#0072B2", xlab = "1 - Specificity", ylab = "Sensitivity", xlim = c(0, 1), ylim = c(0, 1),
       panel.first = grid(col = "#E6E6E6", lty = 1))
  lines(roc_bia$false_positive_rate, roc_bia$true_positive_rate, lwd = 2.4, col = "#E69F00")
  lines(roc_core$false_positive_rate, roc_core$true_positive_rate, lwd = 2.8, col = "#009E73")
  abline(0, 1, lty = 2, col = "#777777")
  points(1 - extended_threshold$specificity, extended_threshold$sensitivity,
         pch = 19, col = "#009E73", cex = 1.1)
  legend("bottomright",
         legend = sprintf(c("APACHE II alone (AUC %.3f)", "BIA-only (AUC %.3f)",
                            "Extended APACHE II + BIA (AUC %.3f)"),
                          performance$auc[1:3]),
         col = c("#0072B2", "#E69F00", "#009E73"), lwd = c(2.4, 2.4, 2.8),
         bty = "n", cex = 0.86)
  dev.off()

  tiff(file.path(output_dir, "Figure_2_Correlation_Heatmap.tif"), width = 7.5, height = 7,
       units = "in", res = 600, compression = "lzw")
  par(mar = c(8, 8, 2, 2))
  palette_values <- colorRampPalette(c("#2166AC", "#F7F7F7", "#B2182B"))(101)
  correlation_labels <- c("Age", "BMI", "APACHE II", "SMI", "PhA", "ECW/TBW", "24-h fluid intake")
  image(seq_len(ncol(correlation)), seq_len(nrow(correlation)), t(correlation[nrow(correlation):1, ]),
        col = palette_values, zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "")
  axis(1, at = seq_len(ncol(correlation)), labels = correlation_labels, las = 2, cex.axis = 0.85)
  axis(2, at = seq_len(nrow(correlation)), labels = rev(correlation_labels), las = 2, cex.axis = 0.85)
  for (row in seq_len(nrow(correlation))) for (column in seq_len(ncol(correlation))) {
    text(column, nrow(correlation) - row + 1, sprintf("%.2f", correlation[row, column]), cex = 0.72)
  }
  box()
  dev.off()

  tiff(file.path(output_dir, "Figure_3_Logistic_Forest_Plot.tif"), width = 8.5, height = 5.5,
       units = "in", res = 600, compression = "lzw")
  labels <- c("APACHE II (per 5 points)", "SMI (per 1 kg/m2)", "Phase angle (per 1 degree)",
              "ECW/TBW (per 0.05)")
  positions <- rev(seq_len(nrow(core_table)))
  par(mar = c(5, 13, 2, 2), las = 1)
  plot(NA, NA, log = "x", xlim = c(0.3, 7), ylim = c(0.5, 4.5),
       yaxt = "n", ylab = "", xlab = "Odds ratio (log scale)")
  grid(nx = NA, ny = NULL, col = "#E6E6E6")
  segments(core_table$ci_low, positions, core_table$ci_high, positions, lwd = 2, col = "#0072B2")
  points(core_table$odds_ratio, positions, pch = 19, cex = 1.1, col = "#0072B2")
  axis(2, at = positions, labels = labels, las = 1)
  abline(v = 1, lty = 2, col = "#777777")
  dev.off()

  calibration_groups <- calibration[calibration$section == "risk_group", ]
  tiff(file.path(output_dir, "Supplementary_Figure_S1_Calibration.tif"), width = 6.5, height = 6.5,
       units = "in", res = 600, compression = "lzw")
  par(mar = c(5, 5, 2, 1) + 0.1, las = 1)
  plot(calibration_groups$mean_predicted, calibration_groups$observed_event_rate,
       type = "p", pch = 19, col = "#009E73", xlim = c(0, 1), ylim = c(0, 1),
       xlab = "Mean predicted probability", ylab = "Observed event rate",
       panel.first = grid(col = "#E6E6E6"))
  corrected_x <- seq(0.001, 0.999, length.out = 300)
  corrected_y <- plogis(
    corrected_calibration["intercept"] + corrected_calibration["slope"] * qlogis(corrected_x)
  )
  lines(corrected_x, corrected_y, lwd = 2.2, col = "#0072B2")
  abline(0, 1, lty = 2, col = "#777777")
  legend(
    "topleft",
    legend = c("Apparent decile groups", "Optimism-corrected logistic calibration", "Ideal"),
    col = c("#009E73", "#0072B2", "#777777"),
    pch = c(19, NA, NA), lty = c(NA, 1, 2), lwd = c(NA, 2.2, 1), bty = "n", cex = 0.78
  )
  dev.off()

  tiff(file.path(output_dir, "Supplementary_Figure_S2_Decision_Curve.tif"), width = 7, height = 6,
       units = "in", res = 600, compression = "lzw")
  par(mar = c(5, 5, 2, 1) + 0.1, las = 1)
  dca_apache <- decision_curve[decision_curve$model == "APACHE II alone", ]
  dca_bia <- decision_curve[decision_curve$model == "BIA-only", ]
  dca_extended <- decision_curve[decision_curve$model == "Extended APACHE II + BIA", ]
  plot(dca_extended$threshold_probability, dca_extended$model_net_benefit, type = "l",
       lwd = 2.7, col = "#009E73", xlim = c(0.05, 0.50),
       ylim = range(c(decision_curve$model_net_benefit, decision_curve$treat_all_net_benefit, 0)),
       xlab = "Threshold probability", ylab = "Net benefit", panel.first = grid(col = "#E6E6E6"))
  lines(dca_bia$threshold_probability, dca_bia$model_net_benefit,
        lwd = 2.3, col = "#E69F00", lty = 1)
  lines(dca_apache$threshold_probability, dca_apache$model_net_benefit,
        lwd = 2.3, col = "#0072B2", lty = 1)
  lines(dca_extended$threshold_probability, dca_extended$treat_all_net_benefit,
        lwd = 2, col = "#D55E00", lty = 2)
  abline(h = 0, lwd = 2, col = "#555555", lty = 3)
  legend("topright",
         legend = c("Extended model", "BIA-only", "APACHE II alone", "Treat all", "Treat none"),
         col = c("#009E73", "#E69F00", "#0072B2", "#D55E00", "#555555"),
         lty = c(1, 1, 1, 2, 3), lwd = c(2.7, 2.3, 2.3, 2, 2),
         bty = "o", bg = "white", box.col = "white", cex = 0.88)
  dev.off()
}

# -----------------------------------------------------------------------------
# Human-readable audit report
# -----------------------------------------------------------------------------
report_lines <- capture.output({
  cat("VERIFIED R ANALYSIS REPORT\n")
  cat("==========================\n\n")
  cat("Data file:", data_file, "\n")
  cat("Patients:", nrow(d), "| Deaths:", sum(y == 1), "| Survivors:", sum(y == 0), "\n")
  cat("Variables:", ncol(d), "| Missing cells:", sum(is.na(d)), "\n")
  cat("APACHE II acute-physiology mismatches:", aps_mismatch, "\n")
  cat("APACHE II age-point mismatches:", age_score_mismatch, "\n")
  cat("APACHE II component arithmetic mismatches:", apache_mismatch, "\n")
  cat("BMI formula mismatches:", bmi_mismatch, "\n\n")

  cat("Core model equation\n")
  coefficients <- coef(model_core)
  cat(sprintf("logit(p) = %.6f %+.6f*APACHEII %+.6f*SMI %+.6f*PhA %+.6f*ECW_TBW\n",
              coefficients[1], coefficients["apache"], coefficients["smi"],
              coefficients["pha"], coefficients["ecw_tbw"]))
  cat("p = 1 / (1 + exp(-logit(p)))\n\n")

  cat("Core model effects\n")
  print(core_table, row.names = FALSE, digits = 5)
  cat("\nModel performance\n")
  print(performance, row.names = FALSE, digits = 5)
  cat("\nPaired DeLong comparisons\n")
  print(delong, row.names = FALSE, digits = 5)
  cat("\nExploratory thresholds\n")
  print(thresholds, row.names = FALSE, digits = 5)
  cat("\nBootstrap internal validation\n")
  print(bootstrap_validation, row.names = FALSE, digits = 6)
  cat("\nCalibration\n")
  print(calibration_validation, row.names = FALSE, digits = 6)
  cat("\nMaximum VIF:", sprintf("%.3f", max(vif_output$vif)), "\n")
  cat("\nNonlinearity sensitivity tests\n")
  print(nonlinearity[, c("variable", "likelihood_ratio_chisq", "p_value")],
        row.names = FALSE, digits = 5)
  cat("\nLinear versus nonlinear ECW/TBW sensitivity model\n")
  print(nonlinear_ecw_comparison, row.names = FALSE, digits = 6)
  cat("\nExploratory apparent decision-curve analysis\n")
  cat("Net benefit was calculated in the development cohort without optimism correction.\n")
  cat("\nR session\n")
  cat(R.version.string, "\n")
  cat("Platform:", R.version$platform, "\n")
  cat("Seed:", args$seed, "| Bootstrap refits requested:", args$bootstrap, "\n")
})
writeLines(report_lines, file.path(output_dir, "analysis_report.txt"), useBytes = TRUE)

cat("ANALYSIS COMPLETE\n")
cat("Output directory:", output_dir, "\n")
cat("Successful bootstrap refits:", successful, "of", args$bootstrap, "\n")