PISA: Interpretable Survival Analysis Pipeline

• The paper details an automated pipeline that constructs interpretable survival models via multi-objective symbolic regression and patient stratification flowcharts.
• PISA employs a dual approach with Cox regression and shallow survival trees to balance predictive performance with interpretability in clinical research.
• The framework integrates robust model pre-selection, multi-split validation, and Kaplan–Meier curves to ensure clinically useful risk stratifications.

Searching arXiv for the cited PISA paper and closely related interpretability/survival works. Pipeline for Interpretable Survival Analysis (PISA) is an interpretable-by-design survival analysis pipeline that provides multiple survival analysis models that trade off complexity and performance, using multiple-feature, multi-objective feature engineering to transform patient characteristics and time-to-event data into multiple survival analysis models, and converting every model into simple patient stratification flowcharts supported by Kaplan–Meier curves (Schlender et al., 13 Sep 2025). It is model-agnostic, but is illustrated through applications of Cox regression and shallow survival trees, the latter avoiding proportional hazards assumptions, and is intended to produce interpretable survival models, intuitive stratification flowcharts, and automated survival analysis workflows in clinical research (Schlender et al., 13 Sep 2025).

1. Definition and scope

PISA addresses a recurrent tension in survival analysis between predictive performance and interpretability. In the formulation used for PISA, survival analysis models the time until an event occurs while explicitly handling censoring, so that patients whose follow-up ends before the event is observed still provide partial information (Schlender et al., 13 Sep 2025). The pipeline is motivated by the observation that Cox regression assumes a hazard of the form

$h(t \mid x) = h_0(t)\,\exp(\beta^\top x),$

with an unspecified baseline hazard and a linear prognostic index, which imposes a proportional hazards assumption and linear effects in the original variables (Schlender et al., 13 Sep 2025). The same source states that deep survival models can capture more complex structure but are difficult to inspect, and that clinicians typically prefer simple, inspectable structures such as scores, rules, or flowcharts (Schlender et al., 13 Sep 2025).

Within this framework, PISA is not a single predictive model but an automated pipeline that starts from survival data and ends with multiple survival models spanning a complexity–accuracy spectrum, together with a simple patient stratification flowchart for each model (Schlender et al., 13 Sep 2025). Its central design principle is that interpretability is enforced at several design points: feature expressions are symbolic, only up to three constructed features per model are allowed, each feature is limited in size, survival models are kept simple, and final outputs are represented as flowcharts and Kaplan–Meier curves (Schlender et al., 13 Sep 2025).

A plausible implication is that PISA should be understood less as a replacement for a particular estimator than as a framework for generating a family of interpretable candidates from which a clinically acceptable model can be selected.

2. Survival modeling assumptions and representational basis

PISA assumes standard right-censored survival data consisting of covariates, time-to-event, and censoring indicators (Schlender et al., 13 Sep 2025). In the reported experiments, the datasets are already preprocessed clinical datasets, and categorical variables are already encoded, so the framework assumes that covariates are encoded as numerical, categorical, or Boolean variables, and that there are no special time-dependent covariates in the presented experiments (Schlender et al., 13 Sep 2025).

The pipeline is instantiated with two survival model families. The first is Cox proportional hazards regression, using constructed features $f=(f_1,\dots,f_k)$ and the hazard

$h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$

with ridge regularization parameter $1$ in the scikit-survival implementation (Schlender et al., 13 Sep 2025). The second is shallow survival trees, which recursively partition the patient population by selecting a split that maximizes separation according to the log-rank test statistic, with maximum depth $3$ and minimum samples per leaf $10$ (Schlender et al., 13 Sep 2025). The paper explicitly notes that these trees do not assume proportional hazards (Schlender et al., 13 Sep 2025).

This yields two distinct interpretive regimes. Cox-based PISA retains the proportional hazards assumption but applies it in a space of engineered non-linear features (Schlender et al., 13 Sep 2025). Tree-based PISA abandons proportional hazards but constrains model complexity through shallow depth and simple binary splits (Schlender et al., 13 Sep 2025). This suggests that PISA’s “complexity–accuracy spectrum” is partly a spectrum of survival-model assumptions.

3. Multi-objective symbolic feature engineering

At the core of PISA is multiple-feature, multi-objective symbolic regression using genetic programming, implemented with GP-GOMEA (Schlender et al., 13 Sep 2025). Each candidate solution is a feature set consisting of up to three symbolic expressions $(f_1,f_2,f_3)$, where each $f_j$ is represented as an expression tree with terminals corresponding to original variables and constants, and primitive operators

$$\{+, -, \times, \div, x^2, \le, =, \text{NOT}, \text{AND}, \text{OR}\}$$

(Schlender et al., 13 Sep 2025). Feature expressions are limited to a maximum of 15 nodes in the expression tree (Schlender et al., 13 Sep 2025).

The symbolic feature search is explicitly multi-objective. For each feature set $F$, PISA fits an elementary survival model using only the constructed features and evaluates performance by the concordance index on 25 stratified train–test splits of the internal data (Schlender et al., 13 Sep 2025). Rather than average all 25 values, it summarizes performance by the interquartile mean of the 25 C-indices, discarding the lowest and highest quartiles and averaging the middle 50%; this becomes the performance objective $f=(f_1,\dots,f_k)$0, to be maximized (Schlender et al., 13 Sep 2025). The second objective is complexity, defined as the total number of nodes in the expressions actually used by the survival model,

$f=(f_1,\dots,f_k)$1

to be minimized (Schlender et al., 13 Sep 2025).

The output of one GP-GOMEA run is therefore a Pareto front of non-dominated models, where no model is simultaneously better in both complexity and performance (Schlender et al., 13 Sep 2025). The search uses a maximum of 50 generations or stops after 5 generations without improvement in hypervolume (Schlender et al., 13 Sep 2025).

A concise summary of the feature-engineering stage is given below.

Component	Specification	Source
Feature set size	Up to three constructed features	(Schlender et al., 13 Sep 2025)
Expression size	Max 15 nodes per feature	(Schlender et al., 13 Sep 2025)
Primitive operators	$f=(f_1,\dots,f_k)$2	(Schlender et al., 13 Sep 2025)
Objectives	Interquartile-mean C-index and expression complexity	(Schlender et al., 13 Sep 2025)
Search engine	GP-GOMEA	(Schlender et al., 13 Sep 2025)

This feature-engineering design is the mechanism by which PISA embeds interpretability before model fitting rather than after it. The symbolic expressions are not explanatory overlays on a separate black-box model; they are the actual covariates used downstream (Schlender et al., 13 Sep 2025).

4. Model generation, pre-selection, and robustness

PISA does not select a single candidate from one stochastic optimization run. Instead, it performs repeated internal splitting and pre-selection. The internal dataset is split 30 times into 75% training and 25% internal validation, stratified by event status (Schlender et al., 13 Sep 2025). For each split, the feature engineering process is run to obtain a Pareto set of models (Schlender et al., 13 Sep 2025). Each candidate model is then evaluated using 1000 bootstrap resamples, jointly on the training and validation sets, ensuring that patients used to train are not used for testing in that bootstrap sample (Schlender et al., 13 Sep 2025). From these bootstraps, 95% confidence intervals for the mean C-index on training and validation data are computed (Schlender et al., 13 Sep 2025).

A model is retained in the pre-selected set if no other model or baseline model has a 95% confidence interval that strictly dominates it in performance (Schlender et al., 13 Sep 2025). The resulting subset contains robust, non-dominated models spanning the complexity–performance trade-off (Schlender et al., 13 Sep 2025). The paper also states that users can enforce additional constraints, such as a minimum acceptable C-index, inclusion of specific original features, or maximum complexity thresholds (Schlender et al., 13 Sep 2025).

This stage is important because the symbolic search is stochastic and the elementary survival models depend on data splitting (Schlender et al., 13 Sep 2025). A plausible implication is that the pre-selection stage serves as a stability filter over both feature discovery and survival-model fitting.

5. Conversion to patient stratification flowcharts

A distinguishing feature of PISA is that every model is converted into a patient stratification flowchart supported by Kaplan–Meier curves (Schlender et al., 13 Sep 2025). The conversion proceeds differently depending on the survival model family.

For survival-tree models, the leaves of the trained survival tree define the initial patient groups directly (Schlender et al., 13 Sep 2025). For Cox-based models, patients are first grouped by quantiles of the prognostic index $f=(f_1,\dots,f_k)$3, using

$f=(f_1,\dots,f_k)$4

which yields 6 initial risk groups with extra resolution at the extremes (Schlender et al., 13 Sep 2025). Because these groups are defined on the Cox score rather than explicit rules, PISA then trains a shallow decision tree classifier to predict the risk-group label from the engineered features, producing explicit rule-based flowchart conditions (Schlender et al., 13 Sep 2025).

For model-agnostic use, the paper describes a fallback in which each constructed feature is discretized, using every value if it has fewer than 7 distinct values or 4 quantiles otherwise, and initial groups are formed as combinations of feature bins (Schlender et al., 13 Sep 2025).

Once initial groups are formed, PISA iteratively merges them until all remaining groups are statistically well separated. For all pairs of groups, a log-rank test is computed; if all pairwise tests are significant at $f=(f_1,\dots,f_k)$5, merging stops; otherwise, the pair with highest similarity is merged and the procedure repeats (Schlender et al., 13 Sep 2025). The final groups are assigned a group-level risk value defined by the median overall survival of that group on the internal data (Schlender et al., 13 Sep 2025).

Kaplan–Meier curves are then computed for each terminal group both on the internal data and on the external validation data using the standard estimator

$f=(f_1,\dots,f_k)$6

where $f=(f_1,\dots,f_k)$7 is the number of events at time $f=(f_1,\dots,f_k)$8 and $f=(f_1,\dots,f_k)$9 is the number at risk just before $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$0 (Schlender et al., 13 Sep 2025).

This flowchart conversion is central to PISA’s clinical presentation layer. A plausible implication is that it translates model complexity into a bounded set of visually inspectable decision paths, allowing clinicians to use the system as a stratification instrument rather than only a risk-scoring device.

6. Complexity–accuracy trade-offs and empirical evaluation

PISA’s explicit purpose is to return multiple models trading off complexity and performance (Schlender et al., 13 Sep 2025). Complexity is measured by expression size, while performance is evaluated using the concordance index during search and external validation with 1000 bootstrap samples and 95% confidence intervals (Schlender et al., 13 Sep 2025). Stratifications derived from the fitted models are also evaluated with the C-index, although they usually incur a small reduction relative to the underlying full models because all patients in a group receive the same static risk (Schlender et al., 13 Sep 2025).

The empirical demonstration covers the GBSG/Rotterdam breast cancer dataset, the SUPPORT dataset, and a symptomatic spinal bone metastases cohort with external validation on the Dutch Bone Metastases Study (Schlender et al., 13 Sep 2025). The paper reports that Cox-based PISA models on GBSG reach external C-index values around 0.688 for top models and are comparable to DeepSurv, with some PISA models matching or exceeding CoxKAN variants (Schlender et al., 13 Sep 2025). Survival-tree-based PISA models on GBSG reach external C-index values around 0.671 for top models and are comparable to or better than baseline Cox on original features, with some tree-based PISA models matching DeepSurv performance despite depth $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$1 and no proportional hazards assumption (Schlender et al., 13 Sep 2025).

For the spinal metastases dataset, both Cox- and tree-based PISA models outperform their respective baselines on external validation, and PISA’s stratifications outperform the manually designed 4-group stratification from the original departmental study, though often using 5–8 groups (Schlender et al., 13 Sep 2025). On SUPPORT, the paper states that when operators are restricted to simple arithmetic and logical forms, PISA does not always reach deep state-of-the-art models that use richer operators, illustrating the performance–interpretability trade-off (Schlender et al., 13 Sep 2025).

The following table summarizes the reported evaluation dimensions.

Dataset / setting	Reported finding	Source
GBSG, Cox-based PISA	External C-index around 0.688 for top models	(Schlender et al., 13 Sep 2025)
GBSG, tree-based PISA	External C-index around 0.671 for top models	(Schlender et al., 13 Sep 2025)
Spinal metastases	PISA stratifications outperform the manually designed 4-group stratification	(Schlender et al., 13 Sep 2025)
SUPPORT	Restricting operators can reduce performance relative to richer deep models	(Schlender et al., 13 Sep 2025)

This evaluation structure supports the claim that PISA is not optimized solely for maximum discrimination; it is optimized for the Pareto frontier between discrimination and symbolic simplicity (Schlender et al., 13 Sep 2025).

7. Relation to adjacent interpretable survival methodologies

PISA is model-agnostic, but its design places it in a broader line of research on interpretable survival analysis. One closely related direction is symbolic or structured feature engineering for survival models. The paper introducing PISA explicitly compares against Cox regression, survival trees, DeepSurv, CoxKAN, and a symbolic-regression approach of Rovito et al. where applicable (Schlender et al., 13 Sep 2025). In that comparison, PISA uses symbolic regression as a front-end rather than as a post-hoc approximation (Schlender et al., 13 Sep 2025).

Other recent survival-specific interpretability approaches emphasize different components of the pipeline. “Interpretable Prediction and Feature Selection for Survival Analysis” proposes DyS, a feature-sparse generalized additive model that combines feature selection and interpretable prediction into one model, particularly for large survival datasets (Ness et al., 2024). “Functional Decomposition and Shapley Interactions for Interpreting Survival Models” introduces SurvFD and SurvSHAP-IQ to analyze time-dependent and interaction effects in hazard and survival functions (Langbein et al., 18 Feb 2026). “Automated and Interpretable Survival Analysis from Multimodal Data” describes a multimodal pipeline where imaging features are interpreted via Grad-CAM, clinical variables are modeled as symbolic expressions through genetic programming, and risk estimation uses a transparent Cox regression (Malafaia et al., 25 Sep 2025). These works suggest complementary interpretations of “interpretable survival analysis”: sparse generalized additive structure (Ness et al., 2024), interaction-aware time-indexed decomposition (Langbein et al., 18 Feb 2026), and multimodal symbolic-plus-visual explanation (Malafaia et al., 25 Sep 2025).

PISA differs in that its defining output is not only an interpretable predictor but a set of flowchart-based patient stratifications derived from a multi-objective symbolic search (Schlender et al., 13 Sep 2025). This suggests that PISA occupies a particular niche within interpretable survival modeling: automated generation of clinically usable rule-based representations from symbolic engineered features and simple survival learners.

8. Example models and clinical interpretation

The paper illustrates PISA using explicit symbolic features. In GBSG, a top tree-based model uses

$h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$2

and

$h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$3

where $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$4 acts as a risk score and $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$5 represents progesterone receptor status (Schlender et al., 13 Sep 2025). The paper states that $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$6 makes the tumor-size contribution grow faster for the largest size group, $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$7 increases risk with more positive nodes, postmenopausal status increases risk, and hormonal therapy reduces risk through the $h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$8 term (Schlender et al., 13 Sep 2025). PISA then converts such features into a flowchart whose thresholds split patients into groups with distinct Kaplan–Meier curves on internal and external data (Schlender et al., 13 Sep 2025).

In the spinal metastases cohort, a tree-based PISA model without brain metastases uses

$h(t \mid f) = h_0(t)\,\exp(\beta^\top f),$9

and

$1$0

capturing poor performance status, visceral metastases, and unfavorable or moderate profile in the symbolic score (Schlender et al., 13 Sep 2025). Another model includes brain metastases through

$1$1

yielding a flowchart in which the absence of visceral and brain metastases is protective when combined with performance score and profile (Schlender et al., 13 Sep 2025).

These examples show how PISA’s engineered features are intended to be read as explicit risk formulas rather than opaque latent dimensions (Schlender et al., 13 Sep 2025). A plausible implication is that the symbolic-expression layer can function as a bridge between conventional clinical reasoning and automated nonlinear feature construction.

9. Limitations and prospective extensions

The PISA paper identifies several limitations. First, performance and interpretability remain in tension: limiting tree depth, number of engineered features, expression length, and operator set may prevent the pipeline from matching richer state-of-the-art models on some problems (Schlender et al., 13 Sep 2025). Second, the complexity measure is a crude proxy for interpretability because it counts expression length without distinguishing operator type or number of distinct original variables (Schlender et al., 13 Sep 2025). Third, symbolic regression can overfit, and the current approach relies mainly on multi-split and bootstrap validation rather than a dedicated anti-overfitting mechanism (Schlender et al., 13 Sep 2025). Fourth, multiple-feature multi-objective GP remains computationally expensive, and the authors identify improved C++ and GPU implementations as future work (Schlender et al., 13 Sep 2025). Finally, the presented experiments use standard right-censored survival data with no time-varying covariates, and future extensions are proposed for time-varying covariates and joint models (Schlender et al., 13 Sep 2025).

Related work suggests several plausible extensions. Multi-objective symbolic regression paired with clinically transparent downstream models also appears in multimodal settings, where genetic programming is used to derive clinical symbolic features and Cox regression is used for final risk prediction (Malafaia et al., 25 Sep 2025). Interaction-aware survival explanation methods indicate that additive explanations can fail on hazard or survival scales and that explicit interaction analysis may be needed when PISA models are extended beyond simple Cox or shallow tree settings (Langbein et al., 18 Feb 2026). Sparse generalized additive survival models such as DyS suggest another possible backend for a PISA-like workflow when direct interpretable prediction and built-in feature selection are desired (Ness et al., 2024).

Taken together, these observations suggest that PISA is best seen as a configurable architecture for interpretable survival model generation whose current published form emphasizes symbolic feature engineering, simple survival learners, and flowchart-based stratification (Schlender et al., 13 Sep 2025).