TCGA-BRCA results: benchmark and stability study

The data: ER status is visible before any model is fit

TCGA-BRCA, 1,043 primary tumors with RNA-seq (806 ER+, 237 ER- by immunohistochemistry, an independent clinical label). Before any feature selection runs, the ER program is already visible in the raw data. That is the point of the positive controls below: a method that cannot find ESR1 in this dataset is not being conservative, it is broken.

PCA of the 1,500 most variable genes, colored by ER status. PC1 alone separates ER+ from ER- with AUC 0.914. The biology is not subtle.

Expression of the four ER positive controls by ER status. ESR1 and PGR show near-complete separation; these are the genes any honest method must recover.

Top 50 most variable genes (z-scored), samples ordered ER+ then ER-. The ER program appears as a coherent block, which is why group-wise selection matches the biology.

1. Before the fix: the thesis method as written

Running the thesis pipeline as written on TCGA-BRCA (RNA only, 1,043 samples, 1,500 genes) gave a warning sign: held-out AUC was 0.974, but the selection missed every known ER driver (ESR1, PGR, FOXA1, GATA3 all at selection frequency 0), with no PAM50 or Oncotype DX enrichment. The diagnosis: the SKAT-based adaptive weights were invalid for continuous data (see Methods), and the lambda rule implied a false discovery rate around 0.65.

	Genes selected	Held-out AUC	ESR1/PGR/FOXA1/GATA3	PAM50	Oncotype DX
V1 (thesis as written)	302	0.974	all missed (freq 0)	none	none
V3 (fixed weights + FDR lambda)	224	0.954	all at freq 1.0	3.8x, p = 0.0003	5.0x, p = 0.012
V3 standard weights (thesis correlation-rank)	75	0.926	all missed (freq < 0.5)	none	none

The fixed adaptive method recovers the ER program: top genes include PGR, GFRA1, STAC2, FABP7, ELF5, all ER-associated, and all four positive controls stay selected even at the stringent 0.9 cutoff (152 genes, PAM50 p = 0.0015, Oncotype DX p = 0.0038). The thesis's standard correlation-rank weights latch onto a tightly correlated but ER-irrelevant block (S100 genes) and miss the biology. Outcome-blind weighting is the wrong prior here.

2. Head-to-head benchmark: five arms

Adding methylation (n = 737 with both modalities), I compared five arms on the ER-IHC endpoint, 5 fresh train/test splits each. Hover the bars for stability and biology numbers.

Held-out AUC by method. Error bars are one SD over 5 splits.

Selection stability for the same arms: mean pairwise Jaccard of the selected sets across the 5 splits.

Arm	Selected (SD)	AUC (SD)	Jaccard	PAM50 sig runs	Oncotype sig runs	Controls (of 5)
Baseline LASSO	30 (12)	0.931 (0.017)	0.201	2/5	4/5	ESR1 5/5; PGR 0/5; FOXA1 0/5; GATA3 2/5
EF-Std	23 (51)	0.899 (0.000)	0.600	0/5	0/5	all missed
EF-Adap	236 (67)	0.930 (0.017)	0.229	5/5	4/5	ESR1 5/5; PGR 4/5; FOXA1 3/5; GATA3 3/5
LF-Std	107 (47)	0.917 (0.018)	0.363	0/5	0/5	nearly all missed
LF-Adap	498 (148)	0.924 (0.015)	0.327	5/5	2/5	ESR1 5/5; PGR 1/5; FOXA1 5/5; GATA3 5/5

Read: the adaptive arms win on biology. Standard group lasso (EF-Std) is the most stable arm (Jaccard 0.600) yet selects nothing meaningful: all four positive controls missed, no signature enrichment. That is the trap this whole study is designed to avoid.

Signature enrichment by arm: -log10 of the median hypergeometric p-value over 5 splits, for the PAM50 and Oncotype DX signatures. Dashed line marks p = 0.05. Only the adaptive arms clear it consistently.

Every split, not just the mean: held-out AUC for each of the 5 fresh train/test splits per arm. The headline means hide real split-to-split variation, especially for the unstable arms.

3. The stability study: S1 to S3

The adaptive models were accurate but unstable (EF-Adap Jaccard 0.229, LF-Adap 0.327). Within one fixed split the method is perfectly reproducible (Jaccard 1.00 across seeds), so the instability was genuine sampling variation. Three stabilization stages, consensus groups (S1), shrunk weights (S2), stability selection (S3), attack exactly that. Hover each bubble for the full numbers.

Prediction vs reproducibility for the adaptive fusion models, TCGA-BRCA ER-IHC endpoint (RNA + methylation, n = 737), 5 fresh train/test splits per arm. Bubble size is the mean number of selected features. Dashed lines trace each family's path: current to S1 (consensus groups) to S2 (shrunk weights) to S3 (stability selection). EF-Adap S3 was adopted: Jaccard 0.68, 123 features, AUC 0.948, controls all 5/5.

Variant	Selected (SD)	AUC (SD)	Jaccard	PAM50	Oncotype	Controls
EF-Adap current	236 (67)	0.930 (0.017)	0.229	5/5	4/5	ESR1 5/5; PGR 4/5; FOXA1 3/5; GATA3 3/5
EF-Adap S1	206.6 (35.2)	0.940 (0.008)	0.482	5/5	5/5	all 5/5
EF-Adap S2	206.2 (35.7)	0.940 (0.008)	0.483	5/5	5/5	all 5/5
EF-Adap S3	122.8 (17.9)	0.948 (0.006)	0.680	5/5	5/5	all 5/5
LF-Adap current	498 (148)	0.924 (0.015)	0.327	5/5	2/5	ESR1 5/5; PGR 1/5; FOXA1 5/5; GATA3 5/5
LF-Adap S1	437.2 (40.5)	0.946 (0.012)	0.535	5/5	2/5	ESR1 5/5; PGR 5/5; FOXA1 4/5; GATA3 5/5
LF-Adap S2	427.8 (46.2)	0.945 (0.013)	0.557	5/5	2/5	ESR1 5/5; PGR 5/5; FOXA1 4/5; GATA3 5/5
LF-Adap S3	206.8 (26.1)	0.944 (0.016)	0.588	5/5	1/5	ESR1 2/5; PGR 1/5; FOXA1 2/5; GATA3 2/5

The late-fusion caution. LF-Adap S3 gains stability (Jaccard 0.588) but loses the positive controls in most splits. If I had optimized for stability alone, I would have picked a model that throws away the ER drivers. That is why every stage is judged on prediction, stability, and biology together.

One more calibration: simple gene baselines predict ER status nearly as well (ESR1 alone: AUC 0.943; four-gene set: 0.939, over 10 splits). The paper's claim is therefore selection quality, stability plus biology, not prediction. Framing it that way is the honest read.

5. What the adopted model actually selects

EF-Adap S3 is the adopted model, so here is what it picks. Seventy genes are selected in all 5 splits; the four ER positive controls (ESR1, PGR, FOXA1, GATA3, in green) are all among them. The rest of the top set is dominated by ER-associated genes (TFF1, TFF3, AGR3, SCGB2A2), which is exactly what a method claiming to recover the ER program should produce.

Selection frequency for the adopted EF-Adap S3 model: top genes by number of splits selecting them, over 5 fresh train/test splits. Green bars are the four ER positive controls.

4. Secondary endpoint: PAM50 Basal-like vs rest

A consistency check only, since the labels come from the same RNA. The fixed pipeline selects 164 genes, held-out AUC 0.994, Jaccard 0.73. It recovers EGFR (frequency 1.0), keratins KRT6A/KRT16, and the luminal program (SCGB2A2, PIP, TFF1) as negative predictors. KRT5/KRT14 were correctly not selected: their univariate AUCs for these calls are only 0.77 and 0.66 in this dataset. The method follows the data, not the textbook.