# I-SPY2 real-data analysis
# =============================================================================
# Produces the real-data exhibits for the harm-aware priority-score paper:
#   1. a LaTeX data-description table (cell counts and raw pCR rates),
#   2. thresholded posterior priority layers at a single level alpha,
#   3. a posterior priority-rank interval plot. 
#
# Estimation convention (single flat Beta(1,1) prior on each cell rate):
#   * Point estimates use the posterior mean rate_pm = (y + 1) / (n + 2).
#   * Uncertainty is propagated through the full Beta(y + 1, n - y + 1)
#     posterior (certification, rank intervals).
#   * The raw rate y / n is reported in the data table for transparency only.
# No minimum cell-size filter is applied; the full trial is used.
#
# Requires in scope: Fun.R  (priorityscore(), get_best_treatment()).
source("./Paper/Fun.R")
# =============================================================================

library(readxl); library(dplyr); library(ggplot2);library(ggh4x)

# -----------------------------------------------------------------------------
# GLOBAL CONFIG
# -----------------------------------------------------------------------------
DATA_PATH <- "./Paper/ISPY2-JCCELL-TableS2.xlsx"
OUT       <- "output"
N_DRAW    <- 10000
SEED_REAL      <- 20250520
ALPHA_MAIN <- 0.30            # single alpha for the main priority-tier figure

# One shared, colourblind-safe subtype palette used by EVERY figure so a given
# subtype has the same colour in the priority-layer plot and the rank plot.
SUBTYPE_BASE <- c("#0072B2", "#009E73", "#D55E00", "#CC79A7",
                  "#E69F00", "#56B4E9", "#333333", "#999999")
make_subtype_palette <- function(subtypes) {
  subs <- sort(unique(subtypes))
  setNames(grDevices::colorRampPalette(SUBTYPE_BASE)(length(subs)), subs)
}


# =============================================================================
# 0. Data loading
# =============================================================================
# Same function from Fun.Simulations 
load_ispy2_cells <- function(path, sheet = 1, skip = 1,
                             subtype = "Receptor Subtype") {
  raw <- read_excel(path, sheet = sheet, skip = skip)
  
  raw %>%
    rename(subtype_clean = !!sym(subtype),
           arm = `Arm (short name)`) %>%
    mutate(subtype_clean = na_if(subtype_clean, "NA")) %>%
    filter(!is.na(subtype_clean)) %>%
    group_by(subtype_clean, arm) %>%
    summarise(n_cell = n(), n_pcr = sum(pCR), .groups = "drop") %>%
    mutate(rate_raw = n_pcr / n_cell,
           rate_pm  = (n_pcr + 1) / (n_cell + 2)) %>%
    rename(subtype = subtype_clean)
}


# =============================================================================
# 1. DATA-DESCRIPTION TABLE  ->  LaTeX (booktabs)
# =============================================================================
# Goal      : a descriptive table of the analysed cells (counts + raw rates).
# Arguments :
#   cells    : from load_ispy2_cells()
#   file     : output .tex path
#   caption  : LaTeX caption text
#   label    : LaTeX \label key
#   digits   : decimals for the raw rate
# Returns   : (invisibly) the wide summary data.frame; writes `file`.
# Layout    : wide table -- rows = arms, columns = subtypes,
#             each cell = "n_pcr/n_cell (rate)". A long-format alternative is
#             given (commented) at the bottom of this function.
# -----------------------------------------------------------------------------
write_data_table_latex <- function(cells, file,
                                   caption = "I-SPY2 cells by receptor subtype and treatment arm: number of pathologic complete responses over cell size, with raw pCR rate in parentheses.",
                                   label   = "tab:ispy2_data",
                                   digits  = 2) {
  
  subs <- sort(unique(cells$subtype))
  arms <- sort(unique(cells$arm))
  
  # entry string per (arm, subtype); "--" where the cell is empty
  entry <- function(a, s) {
    r <- cells[cells$arm == a & cells$subtype == s, ]
    if (nrow(r) == 0) return("--")
    sprintf("%d/%d\\,(%.*f)", r$n_pcr[1], r$n_cell[1], digits, r$rate_raw[1])
  }
  
  # column totals (patients and pooled pCR rate per subtype)
  tot_n <- sapply(subs, function(s) sum(cells$n_cell[cells$subtype == s]))
  tot_y <- sapply(subs, function(s) sum(cells$n_pcr [cells$subtype == s]))
  tot_r <- tot_y / tot_n
  
  colspec <- paste0("l", paste(rep("c", length(subs)), collapse = ""))
  esc <- function(x) gsub("([&%_#])", "\\\\\\1", x)   # escape LaTeX specials
  
  lines <- c(
    "\\begin{table}[t]",
    "\\centering",
    sprintf("\\caption{%s}", esc(caption)),
    sprintf("\\label{%s}", label),
    sprintf("\\begin{tabular}{%s}", colspec),
    "\\toprule",
    paste0("Arm & ", paste(esc(subs), collapse = " & "), " \\\\"),
    "\\midrule",
    vapply(arms, function(a)
      paste0(esc(a), " & ",
             paste(vapply(subs, function(s) entry(a, s), character(1)),
                   collapse = " & "), " \\\\"),
      character(1)),
    "\\midrule",
    paste0("Total ($n$) & ",
           paste(sprintf("%d", tot_n), collapse = " & "), " \\\\"),
    paste0("Pooled rate & ",
           paste(sprintf("%.*f", digits, tot_r), collapse = " & "), " \\\\"),
    "\\bottomrule",
    "\\end{tabular}",
    "\\par\\vspace{2pt}",
    "{\\footnotesize\\raggedright Cell entries are $n_{\\mathrm{pCR}}/n_{\\mathrm{cell}}$ with the raw pCR rate in parentheses; ``--'' marks arm$\\times$subtype combinations with no patients.\\par}",
    "\\end{table}"
  )
  writeLines(lines, file)
  message(sprintf("Wrote LaTeX data table -> %s  (%d arms x %d subtypes)",
                  file, length(arms), length(subs)))
  
  # ---- Long-format alternative (uncomment to emit instead) -----------------
  # cells %>% arrange(subtype, arm) %>%
  #   transmute(Subtype = subtype, Arm = arm, `$n_{cell}$` = n_cell,
  #             `$n_{pCR}$` = n_pcr, Rate = round(rate_raw, digits))
  invisible(cells)
}


# =============================================================================
# 2. Posterior scores
# =============================================================================
# For each candidate switch (subtype, arm != z*), draw Beta(y+1, n-y+1)
# posteriors for the cell rates and compute per-draw priority scores via Fun.R.
#
# Arguments:
#   cells       : load_ispy2_cells() output.
#   method      : "CATE", "Frechet", or "Indep".
#   n_draw,seed : posterior sample size and RNG seed. The seed is fixed so that
#                 all methods share the same Beta draws, coupling ranks across
#                 methods.
# Returns:
#   list(post = n_draw x M matrix of draws over ALL arm x subtype nodes, 
#   prob = M x M matrix of P(score_i < score_j), labels = "subtype::arm").
#   The draw-specific optimal arm carries score 0 in that draw.
# -----------------------------------------------------------------------------
posterior_scores <- function(cells, method, n_draw = N_DRAW, seed = SEED_REAL) {
  set.seed(seed)
  ptype <- c(CATE = "CATE", Frechet = "Frechet", Indep = "Indep")[[method]]
  
  subtypes <- unique(cells$subtype)
  post_list <- list(); label_list <- character(0)
  
  for (s in subtypes) {
    sub <- cells %>% filter(subtype == s) %>% arrange(arm)
    K   <- nrow(sub)
    if (K < 2) next                                # need >= 2 arms to compare
    
    gmat <- matrix(NA_real_, nrow = n_draw, ncol = K,
                   dimnames = list(NULL, sub$arm))
    for (k in seq_len(K))
      gmat[, k] <- rbeta(n_draw, sub$n_pcr[k] + 1,
                         sub$n_cell[k] - sub$n_pcr[k] + 1)
    
    # Draw-specific destination: z*(x) is a functional of the parameters, so it
    # is recomputed within every posterior draw. 
    # Exact ties in which.max have probability zero under continuous Beta draws.
    pol_vec <- rep(1 / K, K)
    sc <- matrix(NA_real_, nrow = n_draw, ncol = K,
                 dimnames = list(NULL, sub$arm))
    for (d in seq_len(n_draw)) {
      z_d      <- get_best_treatment(gmat[d, ])
      sc[d, ]  <- priorityscore(gmat[d, ], pol_vec, z_d, type = ptype)$score
      sc[d, z_d] <- 0                              
    }
    
    colnames(sc) <- paste0(s, "::", sub$arm)       # ALL arms retained
    post_list[[s]] <- sc
    label_list <- c(label_list, colnames(sc))
  }
  
  post <- do.call(cbind, post_list)
  M <- ncol(post)
  prob <- matrix(NA_real_, M, M, dimnames = list(label_list, label_list))
  for (i in seq_len(M)) for (j in seq_len(M))
    prob[i, j] <- mean(post[, i] < post[, j], na.rm = TRUE)
  
  list(post = post, prob = prob, labels = label_list)
}


# =============================================================================
# 3. THRESHOLDED POSTERIOR PRIORITY LAYERS
# =============================================================================
# peel_layers : tiers from groups that dominate no other group.
#   j dominates i iff prob[i,j] > 1 - alpha; columns count the groups
#   dominated by each group.
# -----------------------------------------------------------------------------
peel_layers <- function(prob, alpha) {
  dom <- prob > (1 - alpha); diag(dom) <- FALSE
  labels <- rownames(dom); remaining <- seq_len(nrow(dom))
  layers <- list(); ln <- 0L
  while (length(remaining) > 0) {
    sub <- dom[remaining, remaining, drop = FALSE]
    # First tier: groups that dominate no other remaining group.
    top <- which(colSums(sub) == 0)
    if (length(top) == 0) {
      # The thresholded pairwise relation need not be transitive and can have
      # cycles. Find strongly connected components and remove every source
      # component (no dominance edge to another remaining component) together.
      n_sub <- nrow(sub)
      seen <- rep(FALSE, n_sub); finish <- integer(0)
      visit_forward <- function(v) {
        seen[v] <<- TRUE
        for (w in which(sub[v, ])) if (!seen[w]) visit_forward(w)
        finish <<- c(finish, v)
      }
      for (v in seq_len(n_sub)) if (!seen[v]) visit_forward(v)
      seen <- rep(FALSE, n_sub); components <- list()
      visit_reverse <- function(v, component) {
        seen[v] <<- TRUE; component <- c(component, v)
        for (w in which(sub[, v])) if (!seen[w])
          component <- visit_reverse(w, component)
        component
      }
      for (v in rev(finish)) if (!seen[v])
        components[[length(components) + 1L]] <- visit_reverse(v, integer(0))
      has_outgoing <- vapply(components, function(component) {
        other <- setdiff(seq_len(n_sub), component)
        if (!length(other)) FALSE else any(sub[other, component, drop = FALSE])
      }, logical(1))
      top <- sort(unlist(components[!has_outgoing], use.names = FALSE))
    }
    ln <- ln + 1L
    layers[[ln]] <- sort(labels[remaining[top]])
    remaining <- setdiff(remaining, remaining[top])
  }
  layers
}


# =============================================================================
# 4. PRIORITY-LAYER PLOT — single alpha, methods side by side
# =============================================================================
# Goal      : thresholded posterior priority layers, ONE panel per method
#             (CATE / Indep / Fréchet) at a SINGLE alpha.
# Arguments :
#   cells    : load_ispy2_cells() output
#   methods  : panel order, left -> right
#   alpha    : single certification level
#   dest_threshold : nodes with Pr(arm is draw optimum | data) above this are
#                    ringed with a dashed border to mark them as destinations
#   max_per_row    : nodes per sub-row before a tier wraps
# Returns   : (invisibly) list(nodes, layer_counts, p_best); writes the PDF.
# Reading   : y = priority tier (top = first tier); nodes coloured by subtype,
#             labelled with ArmShort. Dashed ring = arm is frequently the
#             draw-specific destination z*(x), so it carries score 0 in those
#             draws and is not a switching candidate there.
# -----------------------------------------------------------------------------
plot_priority_layers <- function(cells,
                                 methods = c("CATE", "Indep", "Frechet"),
                                 alpha   = ALPHA_MAIN,
                                 n_draw  = N_DRAW, seed = SEED_REAL,
                                 dest_threshold = 0.5,
                                 max_per_row    = 7,
                                 output_dir = ".",
                                 filename   = "ispy2_priority_tiers.pdf") {
  
  RING_RX <- 0.105     # ring radii, data units (x spans ~[-1, 1])
  RING_RY <- 0.135     # (y: tiers are 1 apart, sub-row offsets +/- 0.20)
  
  method_lab <- c(CATE = "CATE", Indep = "Non-neg. dep.", Frechet = "Fr\u00e9chet")
  post_by_method <- lapply(methods, function(m)
    posterior_scores(cells, m, n_draw, seed))
  names(post_by_method) <- methods
  
  # Pr(arm is the draw-specific optimum | data). Method-invariant: z*(x) depends
  # only on the drawn rates, not on the scoring rule, so one method suffices.
  p_best <- colMeans(post_by_method[[1]]$post == 0)
  
  node_rows <- list(); layer_counts <- integer(0)
  
  for (m in methods) {
    pr     <- post_by_method[[m]]$prob
    layers <- peel_layers(pr, alpha)
    layer_counts[m] <- length(layers)
    
    npos <- do.call(rbind, lapply(seq_along(layers), function(li) {
      items  <- layers[[li]]
      n      <- length(items)
      tier_y <- -(length(layers) - li + 1L)
      nr     <- ceiling(n / max_per_row)
      row_id <- rep(seq_len(nr), each = ceiling(n / nr), length.out = n)
      off <- if (nr == 1) rep(0, n)
      else scales::rescale(row_id, to = c(0.20, -0.20), from = c(1, nr))
      xs <- unlist(lapply(split(seq_len(n), row_id), function(ix)
        if (length(ix) == 1) 0 else
          scales::rescale(seq_along(ix), to = c(-1, 1),
                          from = c(0.5, length(ix) + 0.5))), use.names = FALSE)
      data.frame(item_id = items,
                 layer   = length(layers) - li + 1L,
                 y       = tier_y + off,
                 y_tier  = tier_y,
                 x       = xs,
                 stringsAsFactors = FALSE)
    }))
    
    npos$subtype <- sub("::.*$", "", npos$item_id)
    npos$arm     <- sub("^.*::", "", npos$item_id)
    npos$p_best  <- unname(p_best[npos$item_id])
    npos$is_dest <- npos$p_best > dest_threshold
    npos$facet   <- factor(method_lab[[m]], levels = unname(method_lab[methods]))
    
    node_rows[[m]] <- npos
  }
  
  nodes <- bind_rows(node_rows)
  y_scales <- setNames(
    lapply(levels(nodes$facet), function(f) {
      br <- sort(unique(nodes$y_tier[nodes$facet == f]))
      scale_y_continuous(breaks = br, labels = paste0("Tier ", -br),
                         name = "Priority tier",
                         expand = expansion(mult = 0.22))
    }),
    levels(nodes$facet)
  )
  pal   <- make_subtype_palette(cells$subtype)
  
  # Dashed rings for destination-heavy nodes (geom_point has no linetype
  # aesthetic, so the ring is drawn as an explicit closed path).
  ang   <- seq(0, 2 * pi, length.out = 61)
  rings <- nodes %>% filter(is_dest) %>%
    rowwise() %>%
    do(data.frame(item_id = .$item_id, facet = .$facet,
                  rx = .$x + RING_RX * cos(ang),
                  ry = .$y + RING_RY * sin(ang),
                  stringsAsFactors = FALSE)) %>%
    ungroup() %>%
    mutate(grp = paste0(item_id, "|", facet))
  
  pub_theme <- theme_bw(base_size = 10) +
    theme(legend.position = "bottom",
          legend.title    = element_text(size = 9, face = "bold"),
          legend.text     = element_text(size = 8),
          strip.text      = element_text(size = 11, face = "bold"),
          strip.background = element_rect(fill = "grey94", colour = NA),
          axis.title.x    = element_blank(), axis.text.x = element_blank(),
          axis.ticks      = element_blank(), panel.grid = element_blank(),
          panel.border    = element_rect(colour = "grey40", fill = NA, linewidth = 0.3),
          panel.spacing   = unit(0.8, "lines"),
          axis.text.y     = element_text(size = 8),
          axis.title.y    = element_text(size = 9),
          plot.caption    = element_text(size = 8, hjust = 0))
  
  p <- ggplot() +
    geom_path(data = rings, aes(x = rx, y = ry, group = grp),
              linetype = "22", linewidth = 0.45, colour = "grey25") +
    geom_point(data = nodes, aes(x = x, y = y, fill = subtype),
               shape = 21, size = 5, stroke = 0.4, colour = "grey20") +
    ggrepel::geom_text_repel(
      data = nodes, aes(x = x, y = y, label = arm),
      size = 3.1, fontface = "bold", colour = "grey10",
      box.padding = 0.35, point.padding = 0.45, min.segment.length = 0.15,
      segment.size = 0.2, segment.colour = "grey65",
      max.overlaps = Inf, seed = 1) +
    facet_wrap(~ facet, ncol = length(methods), scales = "free") +
    facetted_pos_scales(y = y_scales) +
    scale_fill_manual(values = pal, name = "Subtype") +
    scale_x_continuous(expand = expansion(mult = 0.30)) +
    pub_theme +
    guides(fill = guide_legend(nrow = 1, override.aes = list(size = 4)))
  
  n_rows_total <- max(vapply(split(nodes, nodes$facet), function(d)
    nrow(unique(d[, c("y_tier", "y")])), numeric(1)))
  
  ggsave(file.path(output_dir, filename), plot = p,
         width  = 3.7 * length(methods) + 0.8,
         height = max(5.2, 0.55 * n_rows_total + 2.2), device = DEVICE)
  
  message("Priority layer counts: ",
          paste(sprintf("%s=%d", names(layer_counts), layer_counts), collapse = ", "))
  message(sprintf("Destination-marked nodes (p_best > %.2f): %s",
                  dest_threshold,
                  paste(sort(unique(nodes$item_id[nodes$is_dest])), collapse = ", ")))
  
  invisible(list(nodes = nodes, layer_counts = layer_counts, p_best = p_best))
}


# =============================================================================
# 5. POSTERIOR RANK INTERVALS  (helper)
# =============================================================================
# Goal      : per candidate item, the posterior distribution of its PRIORITY
#             RANK (1 = highest priority to eliminate) under each method.
# Arguments :
#   cells       : load_ispy2_cells() output
#   methods     : methods to rank (default all three)
#   rank_scope  : "global" (rank across ALL candidates, all subtypes) or
#                 "within_subtype". Default "global".
#   n_draw,seed: passed through to posterior_scores()
# Returns   : tidy data.frame  item_id, subtype, arm, method,p_best,
#             rmed, rlo50, rhi50, rlo95, rhi95   (rank summaries)
# Coupling  :all methods share the seed, hence the same Beta draws AND the
#             same per-draw destinations (z*(x) depends only on the drawn
#             rates, not on the scoring rule).
# -----------------------------------------------------------------------------
posterior_rank_intervals <- function(cells,
                                     methods = c("CATE", "Indep", "Frechet"),
                                     rank_scope = "global",
                                     n_draw = N_DRAW, seed = SEED_REAL) {
  
  ps    <- lapply(methods, function(m) posterior_scores(cells, m, n_draw, seed))
  names(ps) <- methods
  posts <- lapply(ps, `[[`, "post")
  
  common <- Reduce(intersect, lapply(posts, colnames))
  posts  <- lapply(posts, function(P) P[, common, drop = FALSE])
  subtype_of <- sub("::.*$", "", common)
  
  # Pr(node is the draw-specific destination | data). Method-invariant.
  p_best <- colMeans(posts[[1]] == 0)
  
  # per-draw ranks (higher score = rank 1); ties averaged
  rank_matrix <- function(P) {
    R <- matrix(NA_real_, nrow(P), ncol(P), dimnames = dimnames(P))
    if (rank_scope == "global") {
      for (d in seq_len(nrow(P)))
        R[d, ] <- rank(-P[d, ], ties.method = "average")
    } else {
      for (s in unique(subtype_of)) {
        cols <- which(subtype_of == s)
        for (d in seq_len(nrow(P)))
          R[d, cols] <- rank(-P[d, cols], ties.method = "average")
      }
    }
    R
  }
  
  out <- lapply(methods, function(m) {
    P  <- posts[[m]]
    R  <- rank_matrix(P)
    qs <- apply(R, 2, quantile, c(.025, .25, .5, .75, .975), na.rm = TRUE)
    data.frame(item_id = common,
               subtype = subtype_of,
               arm     = sub("^.*::", "", common),
               method  = m,
               p_best  = unname(p_best[common]),
               rlo95 = qs[1, ], rlo50 = qs[2, ], rmed = qs[3, ],
               rhi50 = qs[4, ], rhi95 = qs[5, ],
               row.names = NULL, stringsAsFactors = FALSE)
  })
  bind_rows(out)
}

# =============================================================================
# 6. RANK-INTERVAL PLOT 
# =============================================================================
# Goal      : horizontal plot of posterior priority-RANK
#             credible intervals for CATE, Independence and Fréchet.
# Arguments :
#   cells      : load_ispy2_cells() output
#   methods    : which methods to show (order sets the dodge order)
#   order_by   : method whose median rank orders the y-axis AND sets the panel
#                split. Default "CATE"
#   rank_scope : "global" | "within_subtype"  (passed to helper)
#   dest_threshold : nodes with Pr(destination) above this get a second label
#                    line under the arm name,
#   n_col      : number of side-by-side panels (default 3)
#   dodge      : offset between methods within an item's row
# Returns   : (invisibly) the plotted data.frame
# -----------------------------------------------------------------------------
plot_rank_intervals <- function(cells,
                                methods    = c("CATE", "Indep", "Frechet"),
                                order_by   = "CATE",
                                rank_scope = "global",
                                n_draw = N_DRAW, seed = SEED_REAL,
                                dest_threshold = 0.5,
                                n_col  = 3,
                                dodge  = 0.5,
                                output_dir = ".",
                                filename = "ispy2_rank_intervals.pdf") {
  
  df <- posterior_rank_intervals(cells, methods, rank_scope, n_draw, seed)
  df$is_dest <- df$p_best > dest_threshold
  
  # global order: by `order_by` median rank, best (rank 1) first
  anchor <- if (order_by %in% methods) order_by else methods[1]
  ord <- df %>% filter(method == anchor) %>% arrange(rmed) %>% pull(item_id)
  
  block_id <- cut(seq_along(ord), breaks = n_col, labels = FALSE)
  panel_of <- setNames(block_id, ord)
  df$panel <- factor(panel_of[as.character(df$item_id)],
                     levels = seq_len(n_col),
                     labels = paste0("Rank ",
                                     sapply(seq_len(n_col), function(k)
                                       sprintf("%d\u2013%d",
                                               min(which(block_id == k)),
                                               max(which(block_id == k))))))
  
  df$item_id <- factor(as.character(df$item_id), levels = rev(ord))
  df$method  <- factor(df$method, levels = methods)
  
  method_lab <- c(CATE = "CATE", Indep = "Non-neg. dep.", Frechet = "Fr\u00e9chet")
  
  # y-axis labels: destination-marked nodes get a second line with the
  # probability the arm is the draw-specific destination
  lab_lookup <- df %>% distinct(item_id, arm, subtype, p_best, is_dest) %>%
    mutate(label = ifelse(
      is_dest,
      sprintf("%s (%s)\nz*: %.2f", arm, subtype, p_best),
      paste0(arm, " (", subtype, ")")))
  lab_vec <- setNames(lab_lookup$label, as.character(lab_lookup$item_id))
  
  pal    <- make_subtype_palette(cells$subtype)
  shapes <- setNames(c(16, 17, 15, 18)[seq_along(methods)], methods)
  rank_dodge <- position_dodge(width = dodge, reverse = TRUE)
  max_r  <- max(df$rhi95, na.rm = TRUE)
  max_per_panel <- max(table(block_id))
  
  pub_theme <- theme_bw(base_size = 10) +
    theme(legend.position   = "bottom",
          legend.title      = element_text(size = 9, face = "bold"),
          legend.text       = element_text(size = 8),
          legend.box        = "vertical",
          strip.text        = element_text(size = 9, face = "bold"),
          strip.background  = element_rect(fill = "grey94", colour = NA),
          axis.title.y      = element_blank(),
          axis.title.x      = element_text(size = 10, face = "bold"),
          axis.text.y       = element_text(size = 7.5, lineheight = 0.85),
          panel.grid.major.y = element_blank(),
          panel.grid.minor  = element_blank(),
          panel.spacing     = unit(1.0, "lines"),
          panel.border      = element_rect(colour = "grey40", fill = NA, linewidth = 0.3))
  
  p <- ggplot(df, aes(y = item_id, colour = subtype, shape = method)) +
    geom_linerange(aes(xmin = rlo95, xmax = rhi95), linewidth = 0.4,
                   position = rank_dodge) +
    geom_linerange(aes(xmin = rlo50, xmax = rhi50), linewidth = 1.0,
                   position = rank_dodge) +
    geom_point(aes(x = rmed), size = 1.9, fill = "white", stroke = 0.5,
               position = rank_dodge) +
    facet_wrap(~ panel, ncol = n_col, scales = "free_y") +
    scale_colour_manual(values = pal, name = "Subtype") +
    scale_shape_manual(values = shapes,
                       labels = method_lab[methods], name = "Method") +
    scale_y_discrete(labels = lab_vec) +
    scale_x_continuous(breaks = scales::breaks_width(2),
                       limits = c(0.5, max_r + 0.5)) +
    labs(x = "Posterior priority rank  (1 = highest priority to eliminate; median, 50/95% CI)") +
    guides(colour = guide_legend(order = 1, override.aes = list(shape = 15, size = 3), nrow = 1),
           shape  = guide_legend(order = 2, override.aes = list(colour = "grey20"), nrow = 1)) +
    pub_theme
  
  ggsave(file.path(output_dir, filename), plot = p,
         width = 3.6 * n_col + 0.6,
         height = 0.42 * max_per_panel + 2.0,   # taller rows for the 2-line labels
         device = DEVICE)
  invisible(df)
}

# =============================================================================
# RUNNER
# =============================================================================
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)

cells <- load_ispy2_cells(DATA_PATH, subtype="Receptor Subtype")
message(sprintf("Loaded %d cells across %d subtypes / %d arms; sizes %d–%d.",
                nrow(cells), length(unique(cells$subtype)),
                length(unique(cells$arm)), min(cells$n_cell), max(cells$n_cell)))

# 1. data-description table
write_data_table_latex(cells, file.path(OUT, "ispy2_data_table.tex"))

# 2. Priority layers (single alpha, three methods side by side)
plot_priority_layers(cells, methods = c("CATE", "Indep", "Frechet"),
           alpha = ALPHA_MAIN, output_dir = OUT,
           dest_threshold = 0.5,
           filename = "ispy2_priority_tiers.pdf")

# 3. posterior priority-rank interval plot
plot_rank_intervals(cells, methods = c("CATE", "Indep", "Frechet"),
                    order_by = "CATE", rank_scope = "global",
                    dest_threshold = 0.9, # not showing destination probability
                    output_dir = OUT)
