Long-term cancer survival prediction using multimodal deep learning
MultiSurv — a pan-cancer deep learning model that fuses clinical, four omics, and whole-slide image data into 30 years of predicted survival probabilities, trained and tested entirely within TCGA.
The study
| Field | Value |
|---|---|
| Cohort | TCGA, 33 cancer types, 11,081 patients (234 dropped for incomplete follow-up labels) |
| Source | NCI Genomic Data Commons |
| Split | 80/10/10, stratified by cancer type, at patient level — verified in data/labels.tsv: 8,880 train / 1,109 val / 1,092 test |
| External validation | None. TCGA only |
| Scanner/magnification | Diagnostic FFPE slides; 40× only — 20× slides and slides with unreadable aperio.AppMag were dropped |
| Metrics | Antolini time-dependent concordance (Ctd), integrated Brier score (IBS), integrated negative binomial log-likelihood (INBLL); 1,000-resample bootstrap CIs |
Modality availability is uneven, which matters for reading Table 1 (see criticism 2):
| Modality | Patients | Features fed to the model |
|---|---|---|
| Clinical | 11,081 | 10 (9 categorical + age) |
| DNAm | 10,257 | 5,000 highest-variance probes (beta values) |
| CNV | 10,325 | 2,000 highest-variance genes, coded −1/0/1 |
| miRNA | 9,616 | 1,881 (all, RPM) |
| mRNA | 9,605 | 1,000 highest-variance genes (FPKM-UQ) |
| WSI | 8,376 | 5 patches of 299 px, sampled per forward pass |
Architecture
Each modality gets its own submodel emitting a 512-d vector; the vectors are fused; a 4-layer fully-connected block and a 30-unit sigmoid layer emit one conditional survival probability per 1-year interval, and the cumulative product is the survival curve.
- Clinical — entity embeddings per categorical feature plus one batch-normed continuous feature (age).
- Omics — plain fully-connected stacks (3–5 layers).
- WSI — ImageNet-pretrained ResNeXt-50, frozen up to
layer3, then max-pooling over the patch features. 50 patches were generated per slide offline at OpenSlide level 0 (512 px, downsampled to 299); 5 are drawn at random per patient per forward pass. - Fusion — six methods are implemented (
cat,max,sum,prod,embrace,attention); the paper uses element-wise maximum. - Loss — discrete-time negative log-likelihood in the Gensheimer Nnet-survival style. Censored patients are credited for every interval they survived past its midpoint. This is a genuinely non-proportional, fully parametric formulation — the paper’s main methodological claim, and it is properly implemented in
src/loss.py.
The authors test and reject proportional hazards for 6 of 10 clinical features and for 12–18 of 50 principal components in each omics block, which justifies the non-PH design.
Headline results (MultiSurv, test set)
| Input | Ctd | IBS |
|---|---|---|
| Clinical alone | 0.809 (0.793–0.825) | 0.143 |
| Clinical + mRNA (best) | 0.822 (0.805–0.837) | 0.138 |
| Clinical + DNAm | 0.808 | 0.134 |
| All six, incl. WSI | 0.787 | 0.152 |
| WSI alone | 0.569 (0.543–0.597) | 0.220 |
| CNV alone | 0.617 | 0.210 |
Baselines beaten on unimodal data: CPH, RSF, DeepSurv, CoxTime, DeepHit, MTLR, Nnet-survival. Per-cancer Ctd ranges from 0.554 (lung squamous) to 0.988 (thyroid).
Criticisms
Ordered by how much they should change your reading. Everything here is traceable to the article or to a named file in the repository; where I could not settle something, it says so.
1. The pan-cancer concordance number is mostly cancer-type stratification, not patient-level risk discrimination. data/preprocess_clinical.ipynb cell 105 lists project_id — the TCGA cancer type — as the first categorical input, and src/sub_models.py gives it a 33→17 embedding. So the model is told which cancer the patient has. A pan-cancer concordance index pools comparisons across patients with different diseases, and thyroid carcinoma versus pancreatic adenocarcinoma is an easy comparison that any model knowing the tumour type will get right. That inflates the headline 0.82 well above what the model achieves within a disease, which is the clinically useful quantity. The authors do report per-cancer results, and those are the honest numbers: a median well below the pan-cancer figure, and two entities (LUSC 0.554, SARC 0.589) barely above chance.
2. Table 1’s rows are not computed on the same patients. MultimodalDataset drops patients missing all selected modalities (src/dataset.py:33-40). For a unimodal model that means every patient lacking that one modality is dropped, so unimodal WSI is scored on the WSI-available subset (~8,376/11,081 overall, so on the order of 800 of the 1,092 test patients — the per-split counts are not derivable from the repo) while unimodal clinical is scored on all 1,092. Since WSI availability is not uniform across TCGA projects, part of the gap between 0.569 and 0.809 is which cancers each was scored on. Table 2 does not have this problem: clinical is present for everyone, so every clinical-inclusive combination is evaluated on the same 1,092 patients.
3. Two of the ten clinical features are treatments, used as predictors. treatments_pharmaceutical_treatment_or_therapy and treatments_radiation_treatment_or_therapy are inputs. Whether a patient received chemotherapy or radiotherapy is not available at diagnosis, is recorded during follow-up, and is assigned by indication — sicker patients get treated differently, and patients who die early may never be recorded as treated at all. So this is not a baseline prognostic model, and some of the clinical submodel’s strong standalone performance is information that a real prediction-at-diagnosis setting would not have. The paper does not discuss this.
4. A missing modality is not a no-op — it contributes a constant vector to the fusion. src/dataset.py:188,209 substitute an all-zero tensor for a missing modality. But the submodel still runs, and its linear layers have biases, so an all-zero input produces a non-zero, patient-independent 512-d representation. I confirmed this by instantiating the repo’s own model: the missing-modality vector had L2 norm 0.646 and was identical for every patient in the batch. Because every submodel output is post-ReLU and therefore non-negative, element-wise maximum fusion means that constant can only ever raise the fused value — it is never inert. At random initialisation it dominated 28.8% of the 512 dimensions; the trained weights are not published so the trained fraction is unknown, but the structural point does not depend on the weights. For the ~13% of test patients with no mRNA, clinical + mRNA is therefore “clinical maxed against a constant”, not “clinical alone”. This weakens the 0.822-vs-0.809 comparison in a direction that is hard to sign.
Related and cleaner: EmbraceNet._get_selection_probabilities tries to detect missing data by testing whether the representation is all-zero, which by the above it never is. I verified it assigns a fully missing modality a selection probability of 0.5 — i.e. treats it as present. The embrace fusion path does not do what it says. It is not the fusion the paper used, so this does not touch the published numbers, but the same broken test also gates the “non-missing features” filter in MultiSurv.forward that feeds the learned-representation figure (Fig. 3).
5. The best combination is a maximum over ten correlated estimates. 0.822 is the winner of ten modality combinations evaluated on one test set. Even with the test set legitimately blinded until a single final evaluation pass — which the paper states and I have no artifact contradicting — the maximum of ten correlated estimates is optimistically biased as a prediction of new-data performance, and its bootstrap CI does not account for having been selected. It also overlaps clinical-alone’s CI (0.805–0.837 vs 0.793–0.825). The paper’s conclusion that multimodal beats unimodal is directionally plausible but not established by this comparison.
6. Feature selection used all three splits. The top-variance filters that reduce mRNA to 1,000 genes, DNAm to 5,000 probes and CNV to 2,000 genes were computed across every patient: src/scripts/compute_gene_variance.py (“variance across all samples”, no split filter), compute_DNAm_variance.py (filters to patients in labels.tsv, which is all 11,081), and preprocess_omics.ipynb cell 209 (var = cnv.var() on the full matrix). So the test set was not fully held out while the feature space was being defined. Worth stating precisely: this is unsupervised selection — it uses X and never touches survival — so the practical optimism is usually small, and it is a much milder problem than criticism 3. To the authors’ credit, min-max scaling and the label encoders were correctly fit on train only.
7. Test-set predictions are stochastic and not averaged. _get_patches does random.sample(patch_files, 5) on every access, including at evaluation, with no seed and no averaging over draws. Any WSI-inclusive result is therefore one random draw of 5 patches out of ~50 per patient, and re-running the evaluation notebook gives different numbers. No run-to-run variability is reported.
8. Five random 299-px patches cannot represent a gigapixel slide, and the results show it. WSI alone reaches Ctd 0.569 — barely above chance — and adding WSI to the full model lowers Ctd from 0.822 to 0.787. The authors name patch sampling as their main limitation, which is fair, but the framing matters: this paper is not evidence that histology adds prognostic information over clinical and molecular data. It is evidence that this particular sampling scheme fails to extract it.
9. No stain normalisation anywhere. src/transforms.py implements flips, 90° rotations and tensor conversion — no colour handling beyond dataset mean/std. TCGA slides come from dozens of tissue source sites with visibly different H&E, and site correlates with both cancer type and patient demographics, so site-specific colour is learnable signal that is unrelated to biology. Choosing 40×-only slides is a real methodological strength that removes one confounder, but it does not touch this one. I have not cited a specific study for the TCGA site-signature effect here — the mechanism above is what I verified from the code and the cohort. [unverified] as to published magnitude.
10. No external validation, and the CIs understate uncertainty. Everything is TCGA. The bootstrap resamples the test set, so the intervals capture test-set sampling only — not training variability. The model was trained once per configuration; no seed repeats are reported.
What is genuinely good
Worth saying plainly, because the criticisms above are long:
- The discrete-time non-PH formulation is correct, well-motivated by their own PH tests, and cleanly implemented.
- The split is patient-level and stratified by cancer type. No slide-level leakage.
- Dropping 20× slides so that all image data is at one magnification is a control most WSI papers skip.
- Scaling and label encoding are fit on train only.
- Baselines are numerous, current for 2021, and beaten fairly on unimodal data.
- Calibration is reported (IBS and INBLL), not just discrimination — still uncommon.
- The code is MIT-licensed, complete enough to follow end to end, and the README is honest that it is research code rather than a library.
Open questions
- Trained weights are not published, so criticism 4’s trained-model magnitude cannot be settled from the repository.
- Per-modality test-set counts are not derivable — the per-patient data directories are not in the repo — so criticism 2 is stated as a mechanism with an approximate n.
- What within-cancer Ctd looks like averaged across entities, which is the number that should have been the headline. Reported per cancer in the paper’s figures but never summarised.