---
title: "3. Naive analysis vs rnaSentry: what the guardrails prevent"
author: "Abdelkarim Hani Ghrieb"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{3. Naive analysis vs rnaSentry: what the guardrails prevent}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r impact-setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# The three-tier structure

The case studies are written as three acts. Tier 1 runs the analysis the way
a standard, unguarded pipeline would, and prints the headline a researcher
would have believed. Tier 2 runs the *same data* through `rnaSentry` and
shows which guardrail fired, with the statistic that justified it. Tier 3
states the counterfactual: the error the guardrails prevented, and the true
signal the guarded analysis still recovered.

Four failure modes -- the ways a prognostic-signature study most often goes
wrong -- are covered by the case studies:

| # | Failure mode | Where it is demonstrated |
|---|---|---|
| 1 | **Wrong direction**: an inverted concordance convention turns a strong signature into a "failure" | `vignette("case-study-brca")` |
| 2 | **Confounded design**: batch-driven genes masquerade as prognostic; collinear terms are double-counted | `vignette("case-study-confounder-audit")` |
| 3 | **Underpowered discovery**: a mean C-index hides coin-flip folds | `vignette("case-study-small-cohort")` |
| 4 | **False transfer**: screening-internal CV optimism is mistaken for cross-cohort validity | README / main vignette (LUAD honest negative, below) |

# The master impact table

| Failure mode | Naive headline | Guarded headline | Guardrails that fired | What a guarded analysis can claim |
|---|---|---|---|---|
| Wrong direction (GSE20685) | "CV C = 0.217, signature fails to generalize (winner's curse)" | held-out C = **0.783**, direction correct, far above the 0.582 random-gene control | `reverse = TRUE` convention; `C + C_rev = 1` parity tests; direction-sensitive regression tests | the signature is not a 0.22 collapse but a directionally correct held-out signal (CV alone does not exceed selection-noise controls; see BRCA Limitations) |
| Confounded design (synthetic batch cohort) | "six batch-driven genes are prognostic; batch and region are two independent risk factors" | `batch`/`region` flagged redundant (Cramer's V = 0.82); `batch` surrogate eta-squared 0.97; minimal design `~ batch`; adjusted screening drops the batch-artifact genes | `redundant_variable`, surrogate-association flags, `adjusted_screening` | the discovered signature is not a batch echo |
| Underpowered discovery (30 samples, 23 events) | "cross-validated C = 0.74 -- decent" | fold range **0.50-1.00**, sd 0.14, 4.6 events/parameter < 5 | `events_per_parameter`, fold-level reporting | the concordance is an average over coin-flip folds; hypothesis-generating only |
| False transfer (GSE31210 -> GSE50081, LUAD) | "CV C = 0.86, the signature is validated" | CV C = 0.862 vs permutation null 0.836 (delta +0.026); external C = **0.540**, log-rank p = **0.266** -- transfer *not* demonstrated | permutation-null calibration; honest-negative external validation; power analysis | no false claim of cross-cohort validity; the negative is a power finding (35 events), not a tool failure |

# Why the first failure mode is worth dwelling on

A sign-convention bug is the one error class that flips a paper's conclusion
without any red flag in the data. `survival::concordance()` defaults to
*larger score means longer survival*; a Cox risk score means the opposite.
Calling the default on a risk score returns `1 - Harrell's C`. The identity
below is what the parity test suite now asserts on every fold:

```{r sign-convention}
library(rnaSentry)
library(SummarizedExperiment)
library(survival)

# Synthetic cohort for the sign-convention demo (real GEO via GEOquery
# shown in case-study-brca; synthetic used here for offline builds).
# NOTE: seed deliberately differs from any GEO accession so future
# `git log -S "set.seed(20685)"` audits stay unambiguous.
set.seed(42)
n_genes <- 3000
n_samples <- 327
genes <- paste0("GENE", seq_len(n_genes))
samples <- paste0("Sample", seq_len(n_samples))
mat <- matrix(rnorm(n_genes * n_samples, mean = 6, sd = 1.5),
              nrow = n_genes, ncol = n_samples,
              dimnames = list(genes, samples))
risk0 <- colMeans(mat[1:5, , drop = FALSE])
event_time <- rexp(n_samples, rate = 0.08 * exp(0.6 * scale(risk0)[, 1]))
censor_time <- rexp(n_samples, rate = 0.04)
time <- pmin(event_time, censor_time)
event <- as.integer(event_time < censor_time)
coldata <- S4Vectors::DataFrame(time = time, event = event,
                                row.names = samples)
se <- SummarizedExperiment(assays = list(logcounts = mat), colData = coldata)
score <- colMeans(assays(se)$logcounts[1:5, , drop = FALSE])
c_default <- as.numeric(concordance(Surv(se$time, se$event) ~
                                    score)$concordance[1])
c_reverse <- as.numeric(concordance(Surv(se$time, se$event) ~ score,
                                    reverse = TRUE)$concordance[1])
c(c_default = c_default, c_reverse = c_reverse, sum = c_default + c_reverse)
```

A naive reading of the default concordance convention reports 0.349 and
would rationalize it as winner's curse; `1 - 0.349 = 0.651` is the correct,
biologically expected value all along. The same class of error has appeared
in published survival literature; rnaSentry pins the convention in code, in
parity tests, and in the rendered report.

## Pipeline execution on the same cohort

Running `rnaSentry` on the synthetic cohort above demonstrates the guardrail
flags in action:

```{r pipeline-demo, warning = TRUE}
run <- run_rnaSentry(se, time_col = "time", event_col = "event",
                     top_n = 5, repeats = 2, folds = 3, seed = 42,
                     render_report = FALSE)
run$stages$build_signature$flags
```

Every flag in the ledger is backed by a statistic. The `events_per_parameter`
warning, the `reverse = TRUE` convention, and the fold-level concordance
distribution are the guardrails that prevent the three failure modes
demonstrated in the case studies above.

```{r unlock-impact}
run$stages$build_signature <-
    lock_signature(run$stages$build_signature, lock = FALSE)
```

# The LUAD honest negative (GSE31210 -> GSE50081)

The fourth failure mode is demonstrated on independent GEO cohorts in the
README and main vignette; the numbers are reproduced via
`validation/repro_live_geo.R` in the GitHub repository. Because the data are
downloaded live, the walkthrough is narrative here rather than a
build-time vignette:

- **Discovery** GSE31210 (LUAD, 35 deaths), **external validation** GSE50081
  (n = 128, 52 deaths).
- Screening-internal CV C = 0.862 against a permutation null of 0.836
  (delta +0.026) -- even survival-permuted data yields a high null CV (about
  0.8 on ~21k genes), so the CV must be calibrated against a matched null.
- External validation: **C = 0.540, log-rank p = 0.266** -- transfer was
  *not* demonstrated, and the tool reports an honest negative rather than a
  "validated signature."
- A power analysis in the simulation study (`validation/simulate_study.R`)
  quantifies why: transfer power < 0.30 at ~35 discovery events, and >= 0.80
  only near ~300 events for a C ~ 0.65 signature.

## Demonstrating the honest-negative concept with synthetic data

The same pattern -- high screening-internal CV but poor external transfer --
can be reproduced on a small synthetic cohort:

```{r luad-null-demo}
set.seed(7)
n_genes <- 120
n_samples <- 80
cnt <- matrix(rpois(n_genes * n_samples, lambda = 500),
              nrow = n_genes, ncol = n_samples,
              dimnames = list(paste0("gene", seq_len(n_genes)),
                              paste0("S", seq_len(n_samples))))
sig_expr <- colMeans(cnt[1:5, , drop = FALSE])
risk <- scale(sig_expr)[, 1] * 0.4
event_time <- rexp(n_samples, rate = 0.03 * exp(0.8 * risk))
censor_time <- rexp(n_samples, rate = 0.02)
cd <- S4Vectors::DataFrame(time = pmin(event_time, censor_time),
                           event = as.integer(event_time < censor_time),
                           row.names = colnames(cnt))
se_null <- SummarizedExperiment(assays = list(counts = cnt), colData = cd)
discovery <- se_null[, seq_len(60)]
external <- se_null[, 61:80]

sig_null <- build_signature(discovery, "time", "event", top_n = 5,
                            repeats = 2, folds = 3, seed = 7)
cat("Screening-internal CV C:", mean(sig_null$cv_results$c_index), "\n")

km_null <- km_curve(sig_null, discovery)
val_null <- validate_external(sig_null, external,
                              cutpoint = km_null$cutpoint)
cat("External C:", val_null$concordance, "\n")
cat("External log-rank p:", val_null$log_rank_p, "\n")
```

The screening-internal CV is optimistic (high C) but the external validation
on the held-out half is poor -- exactly the honest-negative pattern the
tool is designed to surface rather than hide.

# Reading this package's claims

- **Every number above is reproducible.** The breast-cancer case study downloads
  the live GSE20685 series from GEO via `GEOquery` (which caches the series
  matrix locally), with an evaluated synthetic fallback of identical
  structure if the download is unavailable, and the live-GEO LUAD numbers
  are reproduced by `validation/repro_live_geo.R` in the GitHub repository.
- **The guardrails are testable.** Each flag in the report's audit trail is a
  `(check, severity, detail, stage)` record backed by a statistic, and each
  has a regression test asserting the condition that fires it.
- **This is a demonstration, not a clinical claim.** All signatures are
  provisional.

# Session info

```{r session-info}
sessionInfo()
```