SurvivalNet

A Theano-era package for deep survival models on genomic profiles — Cox partial likelihood as a network loss, Bayesian architecture search, and per-patient gradient attribution — unrunnable on any current Python, carrying a confirmed defect in its Cox likelihood, and worth reading only for the attribution idea and as the cleanest worked example of pooled-cohort concordance inflation in this repo.

Purpose

Not one of my repositories, and the recommendation is emphatically not to adopt it. Nothing here has been run on Memorial data and there is no local_path because it has not been cloned into the user’s own code folder. It is recorded because the reading produced four things later work will want to cite: a confirmed mathematical defect in the published loss function, two statistical findings that independently reproduce criticisms already written up against MultiSurv, a measured lesson about what the phrase “c-index” means in code, and a feature attribution workflow that nothing in the current stack has.

Provenance, from the repository’s own metadata rather than from inference. setup.py declares author='Emory University', author_email='lee.cooper@emory.edu', and url='https://github.com/cooperlab/SurvivalNet' — so this is Lee Cooper’s group as it was constituted at Emory, and the code predates the move of the repository into the PathologyDataScience organisation where it now sits. That organisation is already a known quantity here: ScanTools came from it, and that page’s own survey of the org’s output named SurvivalNet as a sibling repository. This is the second repository from that group evaluated here.

The README points at two papers: an arXiv preprint (1609.08663) and Yousefi et al., Nature Scientific Reports 2017, doi 10.1038/s41598-017-11817-6. Neither has been read — everything below comes from the code and from the artifacts committed alongside it, and no claim here should be read as a claim about what the papers say. [unverified] as to whether the papers disclose any of the defects below.

Licence: Apache-2.0, and consistent. The LICENSE file is the stock Apache 2.0 text, and setup.py’s classifier agrees (License :: OSI Approved :: Apache Software License). No conflict, nothing to route to anyone. The one oddity is cosmetic: setup.py sets license=license_str, stuffing the entire 191-line licence text into a metadata field meant for a short identifier.

Currency, read 2026-07-26. 199 commits, spanning 2015 (3), 2016 (128) and 2017 (68). Last commit 83b99ff, 2017-10-19. Zero commits in the trailing year, and none in the eight years before it. No tags, so no released version ever. This is not “low activity”; the project stopped.

Data used

None locally. The repository ships one dataset, data/Brain_Integ.mat, which is a TCGA glioma cohort — see TCGA. Read directly from the file: 560 patients, 399 features, 64.1% censored, follow-up in days. The feature matrix is five pre-concatenated blocks — 176 protein, 126 mutation, 63 CNV gene, 25 CNV arm, 9 clinical. There are no missing values at all (zero NaN cells in either the raw or normalised matrix). The file also carries TCGA patient barcodes and molecular subtype labels; the barcodes are deliberately not reproduced on this page.

Two data problems visible in the shipped file. At least one patient has a negative survival time (the minimum is -1.0, and three patients have follow-up <= 0), which is not a meaningful follow-up duration and which calc_at_risk will silently sort to the front of the risk ordering. And 97 of the 560 follow-up times are tied, which matters for the c-index discussion below rather than being a defect in itself.

Methods

A feed-forward network over a flat feature vector, trained by maximising a Cox partial likelihood, with optional layer-wise unsupervised pre-training by sparse denoising autoencoders. The interesting parts are the three that the README advertises.

Cox partial likelihood as the loss (survivalnet/model/RiskLayer.py). One linear output per patient is the risk score; the loss is the partial log-likelihood computed with a reversed-cumulative-sum trick over patients sorted by increasing event time, with calc_at_risk supplying each patient’s risk-set start index. Ties are handled by giving every tied patient the first tied index, which is the Breslow convention. This is where the defect is — see below.

Bayesian hyperparameter search (examples/Bayesian_Optimization.py) over six dimensions: depth, width, dropout rate, choice of ReLU vs tanh, and L1 and L2 rates, using the bayesopt package.

Per-patient gradient attribution (survivalnet/analysis/). This is the genuinely distinctive capability and the reason the page exists. RiskCohort backpropagates from the risk output to the input features for each patient, giving a per-patient vector of ∂risk/∂feature. FeatureAnalysis then unit-normalises each patient’s gradient vector, ranks features by the absolute value of the mean gradient across patients, and emits ranked box plots, a risk-gradient clustering heatmap, per-feature Kaplan-Meier plots, and — the part worth stealing — GSEA-Preranked .rnk and .gct files, so that feature-level attributions can be carried up to pathway level. PathwayAnalysis does that aggregation directly, though its docstring restricts it to models trained on pure gene expression. The idea is written up on Gradient-Based Feature Attribution.

Three optimizers, and this is not a cosmetic choice — GD, GDLS (gradient descent with a Wolfe line search) and BFGS. examples/Run.py sets opt = 'GDLS'.

Current state / open questions

Nothing here runs on a current Python, and that alone closes the adoption question. Of 28 .py files, 11 cannot even be parsed by Python 3 — including survivalnet/train.py, the central training entry point, and all four examples/ scripts. Ten fail on Python 2 print statements; Optimization.py fails separately on inconsistent tabs and spaces. import survivalnet therefore raises SyntaxError immediately, before any real code runs, because the package __init__ chain reaches BFGS.py on the second line. setup.py declares only Programming Language :: Python :: 2, so the authors agree.

Worth stating precisely, because it cuts against the obvious conclusion: the syntax is the shallow problem. 17 of the 28 files parse cleanly, and SurvivalAnalysis.py was loaded standalone and executed correctly on Python 3 — so a syntax port is roughly eleven files of mostly mechanical print fixes, not a rewrite. The real blocker is underneath: Theano==0.8.2, for a framework whose own development stopped, and bayesopt==0.3, a C++ library with Python bindings. The README’s Docker image exists precisely so users can “avoid installation of the /bayesopt/ package and other dependencies”, which is the authors telling you the install was already hard in 2017. Porting the framework is the rewrite, and there is no reason to do it when jsurvival and MultiSurv between them cover what this does.

The Cox partial likelihood is not the Cox partial likelihood. This is the finding that matters most, and it was verified independently before being written here. RiskLayer.py:42 is

partial_sum = Te.cumsum(exp)[::-1] + 1

Because exp is computed on prediction - prediction.max() and the maximum is added back two lines later, the arithmetic that actually runs is

log( Σ_{j in risk set} exp(pred_j)  +  exp(max_k pred_k) )

in place of the correct log( Σ_{j in risk set} exp(pred_j) ) — as though a phantom extra patient carrying the cohort-maximum risk score sat in every risk set. An independent check reproduced the arithmetic in NumPy against a brute-force partial likelihood and found the per-patient discrepancy matches log((S+1)/S) to floating-point precision, where S is the shifted risk-set sum; on a small synthetic cohort the objective differed by −18.93 versus −13.81. It is not a numerical-stability guard: the risk-set sum always contains the patient’s own strictly positive term and so can never be zero. It is not a Breslow or Efron tie correction, because it applies uniformly to untied singleton risk sets too. And it does not cancel out of the gradient — a finite-difference check on a four-parameter model gave [-10.57, 7.79, 6.16, -0.52] against a correct [-4.89, -1.27, 0.55, 0.92], differing in sign on two coordinates. Git history rules out the innocent reading: the +1 was present in the file’s first commit (April 2016) when there was no max-subtraction at all, and the later “Fixed NaNs” commit that introduced the max-subtraction left the +1 untouched as context. The distortion is proportionally worst for the smallest, latest-time risk sets, approaching log 2 for a singleton.

What is not established: the effect on final trained performance. That would need the model trained, which was not done. The claim is about the correctness of the objective and its gradient, not about how much the published numbers move.

The three optimizers optimise three different objectives. Read side by side:

Path	Objective actually built
GD (Model.build_finetune_functions)	cost - (lambda1*L1 + lambda2*L2_sqr)
GDLS (GDLS.py:22,27)	cost - L1 - L2_sqr — coefficients dropped, penalty weight pinned at 1.0
BFGS (BFGS.py:27,32)	cost — no regularisation at all

So lambda1 and lambda2 are honoured only on the GD path. examples/Run.py selects GDLS, which ignores them and applies an unscaled penalty. Compounding it, Bayesian_Optimization.py’s search bounds are lb = [1, 10, 0., 0., 0., 0.] and ub = [10, 500, 1., 1., 0., 0.] — the last two entries have lower bound equal to upper bound equal to zero, so the two regularisation dimensions cannot vary and the advertised six-dimensional search is effectively four-dimensional. Two independent mechanisms, both pointing the same way: the regularisation hyperparameters do not do what the README implies.

The penalty that does get applied covers one layer. In Model.py the self.L1 += and self.L2_sqr += statements sit at two tabs of indentation, outside the four-tab body of the for i in xrange(self.n_layers) loop that builds the hidden layers. hidden_layer is therefore the leaked loop variable, and only the last hidden layer’s weight matrix enters the penalty, alongside the risk layer’s. On a 10-layer network — and the search space allows 10 — nine of the ten hidden layers are unpenalised. Two honest qualifications: every layer is still correctly included in the gradient update (self.params is extended inside the loop), so this affects only the penalty term; and the shipped non-search default is n_layers = 1, where “last layer” and “only layer” coincide and the bug is inert. It bites when the Bayesian search picks depth > 1. Biases are never penalised on any layer, which is conventional and fine.

How much the unscaled GDLS penalty actually distorts training was not measured — it would need the model trained. A rough scale check suggests the penalty is comparable to rather than dominant over the partial log-likelihood, so the defensible claim is the mechanism (the coefficients are ignored, and the nominal lambda = 0 setting still gets a full-weight penalty), not a magnitude. [unverified] as to effect size.

Minibatch pre-training silently does nothing. train.py:71 computes n_batches = len(train_set) / (pretrain_config['pt_batchsize'] or len(train_set)), but train_set is a dict of four keys, so len(train_set) is 4 rather than the patient count. With pt_batchsize=50, Python 2 integer division gives 4/50 == 0, xrange(0) iterates zero times, and the entire layer-wise pre-training loop is skipped with no error and no warning. It appears to work only when pt_batchsize is None, where 4/4 == 1 and the data happens not to be sliced. examples/Run.py ships with pre-training disabled anyway, so the advertised “layer-wise unsupervised pre-training” is off by default and quietly broken when switched on the obvious way.

The attribution code recompiles per patient and is capped at one hidden layer. RiskCohort loops over patients calling _RiskBackpropagate, which calls theano.function(...) inside the per-patient call — a fresh Theano compile for every patient, where one compile and 560 evaluations would do. And its givens supplies only Model.masks[0], so a model with more than one hidden layer leaves the remaining dropout masks neither as inputs nor as givens. This was checked empirically rather than argued: a two-layer ifelse-based dropout chain mirroring DropoutHiddenLayer was built against Theano 1.0.5 and, supplied with only layer 0’s mask, raised MissingInputError on the layer-1 mask — ifelse’s runtime laziness does not exempt an unresolved input from the compile-time completeness check. So any architecture the Bayesian search picks with depth > 1 cannot be interpreted by the shipped analysis code at all, and depth is a first-class search dimension. (A separate oddity in the same function swaps the Model.o and Model.at_risk variables relative to its argument names; it was traced and is genuinely inert, because the gradient taken is of risk_layer.output, which depends on neither.)

Two statistical findings that reproduce MultiSurv’s, independently

Both are traceable to data/Brain_Integ.mat and to the authors’ own committed output in results/, and both are more sharply evidenced here than in the MultiSurv note, because this repository ships the attribution ranking that MultiSurv does not.

1. The demo cohort is three diseases wearing one label, and the model is told which. The Subtypes field splits the 560 patients into the standard molecular classes, and their outcomes are not comparable:

Subtype	n	Median follow-up (days)	Event rate
IDHmut-codel	145	656	12.4%
IDHmut-non-codel	209	785	20.6%
IDHwt	203	360	68.0%

A 12.4%-versus-68.0% spread is the entire prognostic range of glioma. And the variables that define those classes are input features: IDH1_Mut and IDH2_Mut are among the 126 mutation features, 1p_CNVArm and 19q_CNVArm among the CNV arm features, and three of the nine clinical features are histological_type-Is-… glioblastoma indicators. The authors’ own results/Gradients.rnk — 399 features ranked by signed mean risk gradient — shows the model leaning hard on that split: the seven strongest protective features are CDKN2A (−0.188), 10q (−0.185), 10p (−0.167), SMARCA4, PTEN, IDH1 (−0.101) and IDH2 (−0.098), while the two strongest risk-increasing ones are age_at_initial_pathologic_diagnosis (+0.179) and histological_type-Is-untreated primary (de novo) gbm (+0.151).

But it gets there through the glioblastoma markers, not through the codeletion definition — worth stating precisely, because the obvious version of this claim is wrong. 1p and 19q, the arms whose codeletion defines the best-prognosis class, rank 280th and 373rd of 399 with small gradients (+0.021 and +0.063). So the attribution recovers the large contrast — IDHwt at 68% events against everything else — via IDH mutation plus the chromosome-10 loss, CDKN2A deletion and PTEN alterations that travel with glioblastoma, and largely ignores the finer 12.4%-versus-20.6% codel/non-codel distinction. (The signs are biologically coherent throughout under a deletion-negative coding: chromosome 10 loss raises risk, 1p/19q codeletion lowers it.)

Either way the consequence for the metric is the same. A concordance index pooled across these three groups is measuring the group separation, which is exactly the mechanism Concordance Index describes — reproduced four years before MultiSurv, in molecular rather than pan-cancer data.

2. A treatment is used as a predictor, and the model learns its indication. radiation_therapy-Is-yes_Clinical is one of the nine clinical inputs, and it ranks fifth from the top of the risk-increasing end of the authors’ own gradient list (+0.094). The confounding is measurable in the shipped file: 42.8% of IDHmut-codel patients were irradiated, against 68.4% of IDHmut-non-codel and 81.8% of IDHwt; the event rate is 44.1% among irradiated patients versus 19.4% among the rest. So the model has learned “received radiotherapy → higher risk”, which inverts the causal direction of the treatment and reflects who gets treated. This is criticism 3 of sources/papers/vale-silva-2021-multisurv.md appearing independently in a different repository, and here the attribution output shows the model actively using it rather than merely having access to it.

The c-index implementation, measured rather than assumed

Worth recording precisely, because the first reading looked like a bug and measurement said otherwise. SurvivalAnalysis.c_index is a pure-Python O(n²) double loop, and it is the one file in the package that parses on Python 3 — so it can be imported and tested directly, which was done against a reference Harrell’s C on the shipped cohort.

• For continuous risk scores, which is what the model emits, it agrees with the reference to within a constant −0.0006 across five random seeds. The published numbers are not materially affected by it.
• The offset traces to a real but tiny deviation: when a censoring time exactly equals an event time the pair enters the denominator, but no branch can ever award it credit.
• The larger deviation is that tied predictions get no half credit. A constant-risk null model scores 0.0003, not 0.5. Artificially discretising the risk scores to two levels drops the value 0.071 below the reference.

So: fine as used, but not a drop-in c-index, and the “chance level is 0.5” intuition does not hold for it. One genuine point in its favour — this model produces a single time-invariant risk score per patient, so Harrell’s C is the correct choice here, unlike the non-proportional setting where Concordance Index shows Antolini’s Ctd is required.

The committed results artifact was not produced by the shipped configuration. Run.py declares N_SHUFFLES = 20 and saves c_index_list.mat after the loop, so a completed run would store 20 values. The committed results/c_index_list.mat holds one value, 0.8968. Whatever produced that file was not the script as configured, so the number should not be quoted as the package’s performance, and it is deliberately not connected to the heterogeneity finding above — that finding rests on the cohort composition and the gradient ranking, which are solid on their own.

Both line-search optimizers step the wrong way when the line search fails. GDLS computes rho_t = -gf_t as its search direction and uses it on the normal path, but the except _LineSearchError fallback is theta + gf_t * 1e-5 — the same vector with the opposite sign, so it moves against the direction the code itself just computed as correct. BFGS.py:73 contains the identical sign flip at 1e-4, which makes a deliberate “small perturbation to escape a stall” reading hard to sustain, since the correctly-signed rho_t is in scope one variable name away. GDLS also sets self.stop = True on that failure, and a repository-wide search finds no code anywhere that ever reads .stop — so training simply continues, and because Run.py passes earlystp=False, nothing else halts it either. Whether this measurably moved any published number is not established: the failure rate is unknown without running the code, and each epoch recomputes a fresh correctly-signed gradient rather than compounding the bad step. [unverified] as to impact.

Early stopping is unsafe in two ways, and it is on by default in the library API. train()’s signature defaults to earlystp=True, and train() takes no validation-set parameter at all — only train_set and test_set — so isOverfitting is fed the test set’s c-index history. Model selection and final reporting would then use the same data. On top of that, when it fires it breaks and reports the best epoch, but returns the model holding the weights from the epoch where the break happened: a repository-wide search finds no deepcopy, no best_* variable, and the three reset_weight methods that exist are never called from anywhere. So max_iter is a number, not a restored model. Both problems are dormant in the shipped script — Run.py passes earlystp=False — which means the published results are not affected, but anyone calling the library directly gets them by default.

Smaller things, noted without ceremony. Model.__init__ carries a branch commented as linear Cox regression with no hidden layers, and it is unreachable: construction dies earlier, at self.n_hidden = hidden_layers_sizes[0], on an empty list. (The README never advertises a linear mode, so this is dead in-code intent rather than a broken feature.) train.py’s pre-training progress line has three format arguments for two placeholders, so the cost it claims to print never appears. Run.py passes data between processes by pickling train_set and val_set into the current working directory. Feature selection is not at issue here because there is none — the 399 features are fixed in the shipped file.

What to use instead

Need	Use instead
Cox model on tabular covariates, regularised	jsurvival — Cox with Schoenfeld diagnostics, plus LASSO/Ridge/Elastic Net, already the group’s tool
A survival loss for a neural network	MultiSurv’s src/loss.py — ~50 lines, PyTorch, discrete-time and non-proportional, and correct
Hyperparameter search	Any current framework; bayesopt==0.3 is a C++ build for a dead Python 2 package
Per-feature attribution	The idea is worth keeping, the implementation is not — see Gradient-Based Feature Attribution

Open questions a later reader might close: whether the 2017 paper reports the +1 in its stated likelihood or a correct one, which would say whether the code diverged from the method or the method itself was described this way [unverified]; and whether the pathway-level .rnk handoff is worth reimplementing for the group’s own work, which is the only capability here with no current equivalent in the stack.

Related: Concordance Index — this repository is the second independent worked example of the pooled-cohort inflation that page is built on, and its c_index supplies the separate, implementation-level lesson about what the metric means in code.

Related: MultiSurv — the same evaluation reflex four years later, from a different group; read together they show treatments-as-predictors and cohort-heterogeneity inflation are habits of the field rather than one team’s slip.

Related: Multimodal Fusion — Brain_Integ.mat is a clean example of the early fusion that page describes as usually wrong: 9 clinical variables flattened into one vector with 176 protein and 126 mutation features.

Related: Gradient-Based Feature Attribution — the one durable idea in the package, written up separately because the decision it raises outlives this code.

Related: ScanTools — the other repository evaluated here from the same organisation, and the page whose survey of that org first flagged this one.

Derived from: repository source read 2026-07-26 at commit 83b99ff — survivalnet/{train.py,model/*.py,optimization/*.py,analysis/*.py}, examples/*.py, setup.py, requirements.txt, LICENSE, README.md, plus the committed artifacts data/Brain_Integ.mat, results/Gradients.rnk and results/c_index_list.mat, and git log for commit dates and the RiskLayer.py history.