Journal titlehttps://doi.org/10.xxxx/xxx Review

Zhao et al.

[$\ast$]Corresponding author<EMAIL_ADDRESS>
†These are joint last authors

# Tutorial: guidelines for survival analysis with omics data

Zhi Zhao∗1,2 John Zobolas1,2 Manuela Zucknick†1 Tero Aittokallio†1,2,3
${}^{\text{\sf 1}}$Oslo Centre for Biostatistics and Epidemiology (OCBE),
Department of Biostatistics, Faculty of Medicine, University of Oslo, 0372
Oslo, Norway
${}^{\text{\sf 2}}$Department of Cancer Genetics, Institute for Cancer
Research, Oslo University Hospital, 0310 Oslo, Norway
${}^{\text{\sf 3}}$Institute for Molecular Medicine Finland (FIMM), HiLIFE,
University of Helsinki, 00014, Finland

(2023; 2023)

###### Abstract

The identification of genomic, molecular and clinical markers predictive of
patient survival is important for developing personalized disease prevention,
diagnostic and treatment approaches. Modern omics technologies have made it
possible to investigate the prognostic impact of markers at multiple molecular
levels, including genomics, epigenomics (e.g. DNA methylation),
transcriptomics, proteomics and metabolomics, and how these potential risk
factors complement clinical characterization of patients for survival
prognosis. However, the massive sizes of the omics data sets pose challenges
for studying relationships between the molecular information and patients’
survival outcomes. We present a general workflow for survival analysis, with
emphasis on dealing with high-dimensional omics data as inputs when
identifying survival-associated omics features and validating survival models.
In particular, we focus on commonly used Cox-type penalized regressions and
hierarchical Bayesian models for feature selection in survival analysis, but
the framework and pipeline are applicable more generally. In cases where
multi-omics data are available for survival modelling, an extra caution is
needed to account for the underlying structure both within and between the
omics data sets and features. A step-by-step R tutorial using The Cancer
Genome Atlas survival and omics data for the execution and evaluation of
survival models has been made available at https://ocbe-
uio.github.io/survomics/survomics.html.

###### keywords:

Time-to-event data, omics, penalized regressions, sparse Bayesian models,
feature selection, survival prediction, model validation

##  1 Introduction

Personalized medicine improves diagnosis and treatment by making use of
patient-specific genomic and molecular markers that are indicative of disease
development. Survival time, i.e. time to an event of interest (for example
time to disease progression), is a widely-used end point and patient outcome
for many diseases, and therefore it has become popular to identify genomic,
molecular and clinical markers for survival prediction of patients suffering
from complex diseases such as cancer. In this tutorial, we will only consider
so-called right-censored survival data, where a patient has been followed for
a certain time period and the event of interest is either observed in this
time period or might occur at a later (as yet unobserved) time point. Right-
censored survival data include both a time and the status of each patient at
that time as joint outcomes, where time is a continuous variable and status is
a binary variable indicating whether the event of interest has been observed
up to the given time or not; in the latter case, we refer to this observation
as censored at the observed time. See Table 1 for an example illustration, and
more details of the survival data in Section “Data”. When using the status
label as an outcome in an ordinary logistic regression, the regression
coefficients will become increasingly uncertain and less reliable with
increasing follow-up time (Green and Symons,, 1983). When using the observed
time (or its transformation, e.g. logarithm of time) as an outcome in an
ordinary linear regression, the presence of censored observations (i.e.,
patients still alive by the end of follow-up period) causes considerable
difficulties for assessing the accuracy of predictions (Henderson,, 1995). For
example, Table 1 shows an example where the proportion of patients surviving
past 10 years is $1/5=20\%$ based on the observed data, but the (unobserved)
factual proportion surviving past 10 years is $3/5=60\%$.

Table 1: An example of right-censored time-to-event (survival) outcomes, illustrated in an example where the outcome of interest is the time from disease diagnosis to death. If a patient was alive at the last observed time point, the measurement is censored at that time point.  Patient | Observed time | Status at | Factual time to death,
---|---|---|---
ID | (years) | observed time | possibly unobserved (years)
1001 | 11 | censored | 20
1002 | 4 | dead | 4
1003 | 5 | censored | 12
1004 | 9 | dead | 9
1005 | 1 | censored | 11
|  |  |

The complexity of multiple types of genomic and molecular data poses further
challenges for building effective statistical models to identify associated
biomarkers that are prognostic for patient survival. For example, multi-omics
profiles, such as those from mRNA expression, DNA copy number, single-point
and other genetic mutations, may be available from the same patient, and these
high-dimensional data come with intra- and inter-dataset correlations,
heterogeneous measurement scales, missing values, technical variability and
other background noise (Hasin et al.,, 2021). Bøvelstad et al., (2007) and
Bøvelstad et al., (2009) compared various machine learning methods and showed
improved survival prediction performance by coefficient shrinkage methods that
combine several data sources, in this case clinical and gene expression data.
A recent large-scale benchmarking study for survival prediction using TCGA
multi-omics data demonstrated a slightly improved survival prediction
performance by taking into account the group structure of multi-omics data
sets (Herrmann et al.,, 2021). However, these studies did not evaluate the
survival model performance with respect to feature selection, i.e. selection
of most prognostic omics variables/features. Another recent study by Bommert
et al., (2022) performed a benchmark of filter methods for feature selection
in high-dimensional gene expression data for survival prediction, and
recommended using the simple variance filter before fitting a
$\ell_{2}$-regularized Cox regression model (Simon et al.,, 2011) for accurate
survival predictions and stable feature selection. Vinga, (2021) reviewed
structured penalized regressions for analyzing high-dimensional omics data for
survival prediction or evaluation. Multiple studies have shown that it is
possible to further improve the prediction accuracy and feature selection by
considering more complex structures, such as pathways, or by identifying
significant features among functional relationships between the omics features
(Wang et al.,, 2013; Chekouo et al.,, 2015; Kundu et al.,, 2018; Wang and
Liu,, 2020; Madjar et al.,, 2021; Li et al.,, 2022).

In this tutorial, we describe a workflow for survival analysis with omics
data, review the commonly used statistical methods for feature selection and
survival prediction, and provide a step-by-step R tutorial using publicly
available multi-omics data from The Cancer Genome Atlas (TCGA) project
(http://cancergenome.nih.gov). The overall survival time, demographic and gene
expression data from primary invasive breast cancer patients in TCGA (TCGA-
BRCA) are retrieved from the Genomic Data Commons Data Portal data release
32.0-36.0. Compared to the previous reviews and benchmarks of survival models
in bioinformatic applications, this tutorial provides a complete workflow
ranging from data preparation to final calibrated models, with a particular
focus on building survival models using high-dimensional (multi-)omics data as
input, and as such covers both commonly used penalized regressions and
Bayesian models for survival analysis with high-dimensional data.

##  2 Data

We use multi-omics data and overall survival outcomes from cancer patients as
an example in this tutorial, but the methods are applicable also to other
diseases with similar data types in personalized medicine applications. The
ultimate goal of personalized medicine is to identify patient-specific risk
factors to guide disease prevention, diagnostic and treatment strategies. The
identification of potential risk factors for cancer patients often considers
clinical, demographic, genomic and molecular information, and their
associations with the patients’ time-to-event data (i.e. survival).

Figure 1: Pan-cancer survival data and multi-omics signatures. (a) Multiple
omics layers (modified from Haukaas et al., (2017) and Jendoubi, (2021)). The
network illustrates intra- and cross-layered biological features or molecules
(e.g. DNA methylation, mRNA, proteins, metabolites). (b) Different origins of
tumors represented by TCGA pan-cancer patient cohorts. (c) Overall survival
times of cancer patients from TCGA. Novel methods are needed to model the
high-dimensional multi-omics data and leverage information from heterogeneous
cancer cohorts for improved survival prognosis and biomarker identification.

###  2.1 Time-to-event data

Time-to-event or survival data contain the event of interest (e.g., death is
the event assumed in this section), together with the time from the beginning
of a study period either to the occurrence of the event, end of the study, or
patient loss to follow-up (i.e. right censoring; discussions of right, left
and interval censoring can be found in Leung et al., (1997)). Figure 1c shows
the survival or right-censored times since cancer diagnosis of patients. The
exact survival time of a patient may be not observed due to censoring.
Therefore, a patient has two outcome indices: censoring indicator $\delta$
(also called status) and observed time $\tilde{T}=\min\\{T^{*},C\\}$, where
$T^{*}$ is the exact survival time and $C$ is the censoring time. Indicator
$\delta$ can also be denoted as $\delta=\mathbbm{1}\\{T^{*}\leq C\\}$, where
$\mathbbm{1}\\{\cdot\\}$ is an indicator function. To characterize the
survival time of a patient, we can use survival function

$S(t)=\mathbbm{P}\\{T^{*}>t\\},$

which gives the probability of the patient’s survival beyond time $t$. Another
useful quantity is the hazard function

$h(t)=\lim_{\Delta t\to 0}\frac{\mathbbm{P}\\{t\leq T^{*}<t+\Delta t|T^{*}\geq
t\\}}{\Delta t},$

which is the instantaneous probability of the patient’s death at time $t$
conditional on the patient having survived up to that time point.

###  2.2 Clinical and demographic data

There are multiple sources of patient-level information that can be explored
to identify risk factors for cancer patients. One can start with routinely
collected and commonly used patient data, such as clinical and demographic
variables. For example, an older male patient with a low body mass index (BMI)
has a relatively high risk of gastric cancer (Nam et al.,, 2022). Table 2
illustrates selected clinical and demographic variables often available for
cancer patients. Clinical and demographic variables are considered important
sources of information for predicting survival and are often used to build
reference models for omics-based prognostic models (Herrmann et al.,, 2021;
Bommert et al.,, 2022).

Table 2: Examples of clinical and demographic variables Variable | Data type
---|---
Sex | binary
Body mass index (BMI, $kg/m^{2}$) | continuous
Ethnicity | nominal
Age at first diagnosis in years | integer
Pathological stage | ordinal
Therapy type (e.g. chemo-, hormone, immuno-) | nominal
|

###  2.3 Multi-omics data

Thanks to the rapid development of modern molecular biotechnology, large
amounts of human genomic and molecular data have become available from many
patient profiling projects. These projects often collect multiple levels of
molecular information such as genomics data for DNA variation, transcriptomics
data for mRNA expression, proteomics data for protein abundance and
metabolomics data for metabolite processes, as illustrated in Figure 1a. Among
the multiple omics levels, metabolomics is the closest to observable
phenotypes, such as tumor growth and proliferation (Carins et al.,, 2011). To
deeply understand the molecular biology of tumor development, integration of
multi-omics data may deliver new insights into the circuits of molecular
interactions that underlie the disease initiation and progression (Tarazona et
al.,, 2021).

DNA-level omics data often include single nucleotide polymorphisms (SNP), DNA
methylation, and somatic copy number variation, see illustration of these data
types in Table 3. Each SNP feature can be coded as $\\{0,1,2\\}$, according to
the number of minor alleles at the given locus in each individual, i.e., AA,
AG, GC or TT. DNA methylation reveals methyl groups added to the DNA molecule,
which can be quantified as a $\beta$-value $\beta=\frac{M}{M+U+a}\in[0,1]$,
where $M$ and $U$ are the fluorescence intensities of the methylated and
unmethylated DNA at a CpG locus, and the offset $a$ is often set to 100 to
stabilize the $\beta$-value when $M$ and $U$ are small. Somatic copy number
variation measures the number of repeated sections of the tumor genome. Table
3 shows the assumed distributions for the downstream statistical modelling.
For example, DNA methylation $\beta$-value is a proportion, often either close
to 1 or 0, which can be characterized by a mixture of two beta distributions.

Table 3: Illustration of multi-omics data and distribution assumptions for
data generated by commonly used high-throughput technologies.

Patient ID SNP Methylation Copy number variation Gene expression miRNA Protein
Metabolite ($\beta$-value) (number of copies) (reads) (reads)
(intensity/concentration) (intensity/concentration) 1001 1 0.2 0 5 5 0.07 0.07
1002 0 0.11 5 2 2 0.1 0.1 1003 1 0.95 10 0 0 9.6 9.6 1004 2 0.5 4 30 30 2.8
2.8 $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$
Technology SNP array Infinium BeadChip SNP array (Gistic2) RNA-Seq RNA-Seq
Mass spectrometry Mass spectrometry Distribution Ordinal Beta mixture Negative
binomial Negative binomial Negative binomial log-Gaussian log-Gaussian
assumptions log2 ratio: Gaussian log2 scale: Gaussian log2 scale: Gaussian
log2 scale: Gaussian log2 scale: Gaussian

RNA-level omics data usually include messenger RNA (mRNA) expression and
microRNA (miRNA) expression. Traditionally, DNA microarrays were used to
measure the expression levels of DNA sequences called probes, which acted as a
proxy for the amount of reads representing a genomic feature of interest. The
reads of microarray expression level can be characterized by a negative
binomial distribution directly, or a Gaussian distribution after log2
transformation (see Table 3). Nowadays, RNA sequencing (RNA-seq) has replaced
DNA microarrays, since it allows for the sequencing of the whole
transcriptome, while DNA microarrays only profile predefined transcripts or
gene probes. miRNAs are small noncoding regulatory RNAs that play an important
role in regulating gene expression and are highly evolutionary conserved in
mammals (Bartel,, 2018).

Protein-level omics data usually originates from mass spectrometry (MS)-based
proteomics profiling of a particular cell, tissue, organ, which can detect and
measure abundance levels of entire or phosphorylated proteins. In contrast to
global proteomics with MS, TCGA consortium has produced more targeted
proteomics profiles, involving a set of 180-250 protein features using
reverse-phase protein arrays (RPPA) (Akbani et al.,, 2014; Zhang et al.,,
2017). In contrast, The Clinical Proteomic Tumor Analysis Consortium (CPTAC)
profiling has produced more than 10000 protein features using MS technology
(Edwards et al.,, 2015). Protein expression data are often Gaussian
distributed after log2 transformation. In global proteomics, there are often
sample-specific missing values due to detection limits of protein
quantification.

Metabolite-level omics data has similar statistical properties to proteomics
data, and metabolomics profiling is usually done with MS-based technologies,
enabling the detection and quantification of many thousands of metabolic
features simultaneously. Nuclear magnetic resonance (NMR) spectroscopy is the
other main analytical technology to profile metabolic processes. An important
limitation of NMR spectroscopy is its relatively low sensitivity, which may
lead to relatively few detected metabolites. Metabolic concentrations often
follow logarithmic Gaussian distribution. Similar to large-scale proteomics,
missing values due to detection limits should be treated differently than
missing values due to measurement artefacts, which are more frequent in
metabolite-level omics data (Sun and Xia,, 2023).

Single-cell sequencing is becoming increasingly more prevalent in many
profiling studies. The modern single-cell multi-omics technologies can produce
multiple levels of omics data derived from the same samples, such as
transcriptomics and chromatin accessibility. The single-cell data types are
similar to the illustration in Table 3, but each measurement originates from
the level of individual cells. This article focuses mainly on survival
analysis either with single bulk omics or multi-omics data, where each omics
layer is treated independently. We briefly comment on single-cell data
analysis in Section “Towards single-cell data analysis”.

###  2.4 Missing data

Missing values are often observed in many types of high-dimensional omics data
due to various experimental reasons (Aittokallio,, 2009; Kong et al.,, 2022).
For example, mRNA transcriptomics data from microarrays have 1%-10% missing
values affecting up to 95% of genes due to corruption of image, scratches on
the slides, poor hybridization, inadequate resolution, fabrication errors (de
Brevern et al.,, 2004; Tuikkala et al.,, 2005); MS-based metabolomics data
have 10%–20% missing values affecting up to 80% of variables due to lack of
peak identification by chromatogram, limitation of computational detection,
measurement error, and deconvolution errors (Hrydziuszko and Viant,, 2012; Sun
and Xia,, 2023). The aforementioned technical reasons can lead to missing data
that is either missing at random (MAR) or missing not at random (MNAR). When
dealing with MNAR data, traditional imputation methods like multiple
imputation may introduce bias. It is thus recommended to remove omics features
which have large proportion (e.g. 50%) missingness over patients, and then
apply imputation methods (e.g. $k$-nearest neighbor imputation) for the rest
of the features with missing values before doing any statistical analysis or
modelling. Alternatively, imputation-free methods (e.g. mixture models) that
can deal with missing values can be applied directly (Taylor et al.,, 2022).
Single-cell RNA-seq (scRNA-seq) data has a vast number of zeros, so-called
gene dropout events, leading highly scarce data matrices. Jiang et al., (2022)
discussed the sources of biological and non-biological zeros in scRNA-seq data
and the existing approaches for handling them.

##  3 Survival analysis with low-dimensional input data

Let us assume we have data
$\mathcal{D}=\\{(\tilde{T}_{i},\delta_{i},\mathbf{X}_{i}):i=1,\cdots,n\\}$ for
$n$ patients, where $\tilde{T}_{i}$ is the observed survival time,
$\delta_{i}$ the censoring indicator and $\mathbf{X}_{i}$ contains $p$
covariates including clinical, demographic and multi-omics features. To
estimate a survival function $S(t)$ given the data $\mathcal{D}$, one needs to
keep track both of the number of patients at risk and those who left the study
at any time point (here we only consider the case of right censoring and
assume no delayed entry). At time $t$ there are $Y(t)=\sum_{i=1}^{n}Y_{i}(t)$
patients at risk, where $Y_{i}(t)=\mathbbm{1}\\{\tilde{T}_{i}\geq t\\}$ is the
indicator that patient $i$ is at risk. The non-parametric Kaplan-Meier (KM)
estimator (Kaplan and Meier,, 1958) uses the multiplication rule for
conditional probabilities to obtain an estimation of the survival function

$\hat{S}(t)=\prod_{\begin{subarray}{c}k=1\\\ T_{k}\leq
t\end{subarray}}^{K}\left\\{1-\frac{d_{k}}{Y(T_{k})}\right\\},$

where all events happen (e.g. patients die) at $K$ distinct times,
$T_{1}<T_{2}<\cdots<T_{K}$ $(K\leq n)$ and there are $d_{k}\geq 1$ events
happened at the time $T_{k}$. If no two events happen at the same time point,
$d_{k}=1$ and $k=1,\cdots,n$. The KM estimator gives an estimate of the
marginal survival function, i.e., when you disregard the information from the
covariates.

Figure 2a shows the KM curve for TCGA-BRCA primary breast cancer patients
data. Some basic statistics can be revealed from the survival curve. For
example, the estimated median survival time, i.e. the time when the survival
probability is 50%, of all the breast cancer patients is 10.8 years (dashed
line in Figure 2a), and 1-, 5- and 10-year survival probabilities are 0.988,
0.853 and 0.658, respectively. A log-rank test (Peto and Peto,, 1972) can be
used to test whether two groups of patients (e.g. with treatment
(pharmaceutical/radiation therapy) or nontreatment) have the same (null
hypothesis) or different survival functions (alternative hypothesis), and
provide a corresponding $p$-value (see Figure 2b). The log-rank test can also
be used to compare the survival probabilities of any subgroups of patients
based on other categorical variables.

Figure 2: Kaplan-Meier curves of TCGA-BRCA data. (a) KM curve of all the TCGA-
BRCA patients’ survival data. (b) KM curves of the TCGA-BRCA patients’
survival data grouped by treatment (i.e. pharmaceutical/radiation therapy) or
nontreatment. The log-rank test is used to compare the two survival
distributions corresponding to the two groups of patients.

In the case where multiple clinical, demographic or omics features are
available, one can first explore each variable’s association with survival
outcomes separately. For a categorical variable, KM curves and log-rank tests
can be used to investigate whether there is a difference between multiple
survival curves categorized by the variable. For a continuous variable $X$,
the semi-parametric Cox proportional hazards model (Cox model, Cox, (1972)) is
often used:

$h(t|X)=h_{0}(t)\exp\\{X\beta\\},$ (1)

where $h_{0}(t)$ is the baseline hazard function and is left unspecified. As
the name “proportional hazards” implies, the Cox-model estimated hazard
functions of two individuals (with different values of the covariate $X$) are
indeed proportional, because $h_{0}(t)$ does not depend on $X$ and is thus
assumed to be the same for all individuals. The functional form (1) describes
the log-linear relationship between variable $X$ and the hazard $h(t|X)$ at
any given timepoint $t$. It may be difficult to satisfy the log-linear
relationship based on the original scale of some omics features, e.g., gene
expression data from a DNA microarray study; in those cases, the use of
$\log_{2}$-transformation of the data (Table 3) can be helpful. The Cox model
can also be used to investigate the risk of a categorical variable if the
assumptions above are satisfied, which provides an effect estimate (hazard
ratio, HR) in addition to the $p$-value.

To reduce the computational burden of multivariable analysis (see next
Section), a heuristic is to consider preselection of a subset of omics
features. A simply approach is to include omics features at a specific
statistical significance level when fitting a univariate Cox model (1) for
each omics feature, which would usually be higher than the commonly used 0.05
threshold to avoid losing important features, e.g. 0.1 or 0.2. However, this
approach will focus on features which are independently associated with the
outcome and might miss variables that are important in combination. Another
quite simple approach is to preselect omics features by variance, since larger
variability usually implies higher biological information. For example, one
could preselect omics features explaining 50% of the total cumulative
variances (Zhao and Zucknick,, 2020). Feature preselection is a recommended
method to reduce dimensionality of hundreds of thousands of omics features and
improve the stability of final feature selection (Bommert et al.,, 2022).

When drawing final conclusions on survival differences solely from the
univariate Cox model (1) via $p$-values, it is important to adjust the
$p$-values for multiple comparisons to be able to control false positives
globally, e.g. by controlling the family-wise error rate or false discovery
rates. The problem of multiple testing is beyond the scope of this article,
but the interested readers are referred to a recent review article Korthauer
et al., (2019). Note also that univariate analysis does not consider any
relationships between multiple omics features, e.g. confounders for survival
(Clark et al.,, 2003). To avoid making seriously misleading conclusions in
such cases, it is necessary to perform multivariable survival analysis
(Bradburn et al.,, 2003).

##  4 Survival analysis with high-dimensional input data

Multi-omics data can have hundreds of thousands of variables measured at
various molecular levels, which challenges classical multivariable regression
models for time-to-event endpoints, since the number of variables is larger
than the number of patients (i.e. $p\gg n$). To proceed with survival
analysis, one option is to reduce the dimension of the multi-omics features
via unsupervised learning, and then investigate the association of the reduced
low-dimensional variables with survival outcomes. An alternative approach is
to directly use supervised learning methods, such as penalized regressions or
sparse Bayesian models (Table 4), which enable the modelling of high-
dimensional multi-omics features, and the selection of a few important
features associated with survival outcomes.

###  4.1 Unsupervised learning

Unsupervised learning methods aim to identify hidden patterns or data
groupings, and are for example useful when a phenotype (e.g. a cancer type) is
to be divided into several subtypes (e.g. to explain heterogeneity among
patients). For example, breast cancer has been traditionally categorized into
five conceptual molecular classes, originally using pairwise average-linkage
cluster analysis of DNA microarray data, to better understand tumor biology
and guide clinical decision making (Perou et al.,, 2000; Sørlie et al.,,
2001). Unsupervised methods learn underlying patterns from unlabeled data by
transforming high-dimensional omics features into a lower dimensional space.
Principal component analysis (PCA) is a classical multivariate technique that
represents high-dimensional features in a low-dimensional space by building
orthogonal (uncorrelated) linear combinations of the features that best
capture the variance in the high-dimensional data (Pearson,, 1901). Figure 3a
shows the first two principal components of PCA using the RNA-Seq data of
breast cancer patients from TCGA.

Different from the distance-based PCA with linear transformation, non-linear
techniques have recently emerged, such as $t$-stochastic neighbour embedding
($t$-SNE). $t$-SNE uses pairwise similarities of individuals based on
Kullback-Leibler divergence to project similarities into a lower dimensional
space, ensuring that individuals with similar omics features are close in the
generated embedding (van der Maaten and Hinton,, 2008; van der Maaten,, 2014),
see Figure 3b. The focus of $t$-SNE is on preserving local distances between
neighbouring data points. Another non-linear alternative to PCA is UMAP
(Uniform Manifold Approximation and Projection), which is a general dimension
reduction method built upon Riemannian geometry (McInnes et al.,, 2018).
Compared to other non-linear methods of dimension reduction such as $t$-SNE,
UMAP can sometimes provide better visualization quality in a shorter amount of
time, while preserving the global structure of the omics data better, see
Figure 3c.

Figure 3: Unsupervised projection and clustering of transcriptomic data from
TCGA breast cancer cohort. (a) PCA. (b) $t$-SNE. (c) UMAP.

Unsupervised methods only make use of the input data matrix $\mathbf{X}$ and
are agnostic to the survival information. To identify omics features as risk
factors associated with survival outcomes, a straightforward approach is to
use a few representative components from PCA, $t$-SNE or UMAP as covariates in
a multivariable Cox model. An alternative is to use semi-supervised methods
(Bair and Tibshirani,, 2004) that combine the clustering procedure and
survival modelling together. However, these approaches lose interpretability
of the individual omics features, since each component is a linear or
nonlinear combination of all omics features.

Table 4: Cox-type supervised learning methods Method | Feature selection via | Grouping | Uncertainty | Comment
---|---|---|---|---
|  | effects | quantification |
Penalized regressions: |  penalty |  |  |
Lasso Cox [Tibshirani, (1997)] |  $\ell_{1}$-norm | ✗ | ✗ |
Adaptive Lasso Cox [Zhang and Lu, (2007)] |  weighted $\ell_{1}$-norm | ✗ | ✗ | less false positives than Lasso
Elastic Net Cox [Simon et al., (2011)] |  $\ell_{1}/\ell_{2}$-norm | ✓ | ✗ |
Group-Lasso Cox [Kim et al., (2012)] |  $\ell_{2}$-norm | ✓ | ✗ | independent groups of features selected
Sparse Group-Lasso Cox [Simon et al., (2013)] |  $\ell_{1}/\ell_{2}$-norm | ✓ | ✗ |
SCAD Cox [Fan and Li, (2002)] | quadratic spline and symmetric penalty | ✓ | ✗ | selection of relatively large effects
SIS Cox [Fan et al., (2010)] | top ranked variables and any penalty | ✓ | ✗ | suited to ultra-high dimensions
Bayesian models: |  shrinkage prior |  |  |
[Lee et al., (2011)] | Lasso (Laplace) prior | ✗ | ✓ | selection of posterior mean with a cutoff
[Lee et al., (2015)] | Elastic Net, group/fused Lasso priors | ✓ | ✓ | selection of posterior mean with a cutoff
[Konrath et al., (2013)] | Spike-and-slab prior | ✗ | ✓ |
[Madjar et al., (2021)] | Spike-and-slab & MRF priors | ✓ | ✓ |
[Mu et al., (2021)] | Horseshoe prior | ✗ | ✓ | selection of posterior mean with a cutoff
|  |  |  |

###  4.2 Supervised learning via penalized regressions

For the purpose of personalized cancer medicine, one is typically interested
in identifying risk factor combinations from clinical and omics features.
These factors can be targeted (directly or indirectly) via therapeutic
strategies or used for diagnostics. Therefore, the objective is to identify a
parsimonious set of features linked to survival outcomes by utilizing the
wealth of information present in, for example, the vast amount of available
omics data. Penalized Lasso Cox regression (Tibshirani,, 1997) can be used to
select a few relevant omics features by estimating their coefficients as non-
zero (the non-relevant features’ coefficients are shrunk to zero) via
maximizing the penalized partial log-likelihood function of the regression
coefficients with $\ell_{1}$-norm penalty

$\frac{2}{n}\ell(\bm{\beta}|\mathcal{D})-\lambda\|\bm{\beta}\|_{\ell_{1}}$ (2)

Here, $2/n$ is a scaling factor for convenience,
$\mathcal{D}=\\{(\tilde{T}_{i},\delta_{i},\mathbf{X}_{i}):i=1,\cdots,n\\}$,
$\mathbf{X}_{i}$ includes $p$ (omics) features of the $i$-th patient,
$\lambda$ is a tuning parameter to control the overall penalty strength of the
coefficients, $\|\bm{\beta}\|_{\ell_{1}}=\sum_{j=1}^{p}|\beta_{j}|$, and the
partial log-likelihood is

$\ell(\bm{\beta}|\mathcal{D})=\log\prod_{i=1}^{n}\left\\{\frac{\exp(\mathbf{X}_{i}\bm{\beta})}{\sum_{l\in\mathcal{R}_{k}}\exp(\mathbf{X}_{l}\bm{\beta})}\right\\}^{\delta_{i}},$

where $\mathcal{R}_{k}=\\{l:Y_{l}(T_{k})=1\\}$ is the risk set at time
$T_{k}$. The Elastic Net Cox model (Simon et al.,, 2011) considers both the
Lasso feature selection and the grouping effect of correlated omics features
in ridge regression via a combination of the $\ell_{1}$\- and $\ell_{2}$-norm
(i.e. $\ell_{1}/\ell_{2}$-norm) penalty
$\lambda\left\\{\alpha\|\bm{\beta}\|_{\ell_{1}}+\frac{1}{2}(1-\alpha)\|\bm{\beta}\|_{\ell_{2}}^{2}\right\\}$
(where $\|\bm{\beta}\|_{\ell_{2}}^{2}=\sum_{j=1}^{p}|\beta_{j}|^{2}$) , which
can usually improve the prediction performance over the Lasso Cox model.
Figure 4 shows an example of Elastic Net Cox model feature selection from gene
expression features associated with breast cancer patients’ survival. Note
that often we wish to include a small set of well-established clinical risk
factors in a survival model. Since they are established as important
covariates, they can be included in the Lasso or Elastic Net Cox model as
mandatory covariates without penalization. Then the penalized partial log-
likelihood function becomes

$\frac{2}{n}\log\prod_{i=1}^{n}\left\\{\frac{\exp(\mathbf{X}_{0i}\bm{\beta}_{0}+\mathbf{X}_{i}\bm{\beta})}{\sum_{l\in\mathcal{R}_{k}}\exp(\mathbf{X}_{0l}\bm{\beta}_{0}+\mathbf{X}_{l}\bm{\beta})}\right\\}^{\delta_{i}}-\text{pen}(\bm{\beta}),$

where $\bm{\beta}_{0}$ are coefficients corresponding to the $i$-th
individual’s mandatory covariates $\mathbf{X}_{0i}$, pen($\bm{\beta}$) is a
$\ell_{1}$\- or $\ell_{1}/\ell_{2}$-norm penalty for feature selection of
omics features $\mathbf{X}_{i}$.

There are many alternative penalties that will achieve feature selection which
have been applied to Cox proportional hazards regression for survival
outcomes, such as the Adaptive Lasso Cox model that incorporates different
penalties for different coefficients to retain important variables (Zhang and
Lu,, 2007), Group-Lasso that performs group selection on (predefined) groups
of variables (Kim et al.,, 2012), Sparse Group-Lasso that introduces sparse
effects both on a group and within group level (Simon et al.,, 2013), the
smoothly clipped absolute deviation (SCAD) Cox model that overcomes
substantially biased estimates for large coefficients in ultra-sparse models
(Fan and Li,, 2002), sure independence screening (SIS) procedure in
combination with Cox model that speeds-up the feature selection dramatically
and can also improve the accuracy of estimation when dimensionality becomes
ultra-high, i.e., $\log(p)=\mathcal{O}(n^{\xi})$ for some $\xi>0$ (Fan et
al.,, 2010). However, all these penalized Cox models do not directly provide
uncertainty of feature selection or survival prediction. One empirical way for
uncertainty quantification is through additional resampling-based methods, see
Section “Survival model validation” for more details.

Figure 4: Coefficient trace plot of an Elastic Net Cox model for overall
survival prognosis of breast cancer patients from TCGA based on transcriptomic
data and mandatory demographic variables age and ethnicity. The y-axis shows
the magnitude of each feature’s coefficient given the strength of penalization
displayed on the x-axis (from left to right the penalization decreases). The
vertical gray line indicates the optimal $\lambda=0.032$ (maximizes the
partial likelihood via cross-validation) and its corresponding selected
features are marked with green (positive coeﬀicient) and red (negative
coeﬀicient) colors. Note that the demographic variables age and ethnicity were
not penalized, so that their coeﬀicient paths (black color) did not start from
zero in the figure.

###  4.3 Supervised learning via Bayesian priors

Bayesian inference is an appealing approach for survival analysis due to its
ability to provide straightforward uncertainty quantification (e.g. credible
intervals) of parameters and survival probabilities. For instance, Lee et al.,
(2011) proposed a Bayesian version of the Lasso Cox (Bayesian Lasso Cox) model
that provides credible intervals of coefficients fairly straightforward (see
Figure 5), but it is not easy to derive confidence intervals of coefficients
in a Lasso-type model. In general Bayesian regression models, one is
interested in the posterior density $f(\bm{\beta}|\mathcal{D},\lambda)$, which
is proportional to the product of the prior density and the likelihood
function. In the Bayesian Lasso Cox model, all the coefficients are assigned
to the same double exponential (DE, also known as Laplace, Figure 6a) prior

$\mathcal{DE}(\bm{\beta}|\lambda)=\prod_{j=1}^{p}\frac{\lambda}{2}\exp\\{-\lambda|\beta_{j}|\\}$

and a simple version of the derived log posterior density is

$\log
f(\bm{\beta}|\mathcal{D},\lambda)=\ell^{*}(\bm{\beta}|\mathcal{D})-\lambda\|\bm{\beta}\|_{\ell_{1}}+C,$
(3)

where $\ell^{*}(\bm{\beta}|\mathcal{D})$ is the full log-likelihood function
i.e. $\log(\int_{0}^{t}h_{0}(s)ds)+\ell(\bm{\beta}|\mathcal{D})$, $C$ is a
normalizing constant independent of $\bm{\beta}$ and the $\ell_{1}$-norm
penalty tends to choose only a few nonzero coefficients, which is similar to
(2) in the frequentist Lasso Cox model. Markov chain Monte Carlo (MCMC)
sampling can be performed for posterior inference of $\bm{\beta}$. Note that
instead of using the partial log-likelihood $\ell(\bm{\beta}|\mathcal{D})$ in
(3), a full log-likelihood function is used, which includes the baseline
hazard function. This can be achieved, for example, by assigning a stochastic
process prior, e.g. a gamma process, for the cumulative baseline hazard
function. More details about the prior setup and inference can be found in Lee
et al., (2011). Lee et al., (2015) extended the Laplace prior to Elastic Net
prior, fused Lasso prior and group Lasso prior, which are often more suitable
for correlated omics features in survival analysis. But Lee et al., (2011) and
Lee et al., (2015) assign the same shrinkage priors to all covariates
indiscriminately, see Zucknick et al., (2015) for an extended Bayesian Lasso
Cox model which permits the use of mandatory covariates.

Figure 5: Point estimation and uncertainty quantification of regression
coeﬀicients using Bayesian Lasso Cox model with Laplace prior for overall
survival prognosis of breast cancer patients from TCGA based on transcriptomic
data and mandatory demographic variables age and ethnicity. Solid dots
indicate the posterior mean over 20000 MCMC iterations excluding burn-in
period, and horizontal lines show the corresponding 95% credible intervals.

Note that the Bayesian version of the Lasso with Laplace-type priors in
practice do not result in automatic feature selection, because only the
posterior modes of the coefficients are equivalent to the frequentist Lasso
solution, while in Bayesian inference one usually focuses on the posterior
means as point estimates. As an alternative, a particular omics feature can be
excluded if the estimated credible interval of the corresponding coefficient
covers zero. For example, Figure 5 indicates the selected variables as age,
MELK, ORC6L, CDC20, KRT17, FOXC1, ESR1 and MMP11.

Figure 6: Density of shrinkage priors in Bayesian survival modelling. (a)
Density of the double exponential (Laplace) prior. (b) Density of the mixture
spike-and-slab prior. The spike component $\delta_{0}(\beta_{j})$ induces
$\beta_{j}=0$ and the slab component $\mathcal{N}(0,\tau_{j})$ induces
$\beta_{j}\neq 0$.

Stochastic search variable selection (SSVS) is an alternative approach to
identify important covariates (George and McCulloch,, 1993; Konrath et al.,,
2013). SSVS uses independent spike-and-slab priors for regression
coefficients, e.g.

$\displaystyle\beta_{j}|\gamma_{j},\tau_{j}^{2}$
$\displaystyle\sim\gamma_{j}\mathcal{N}(0,\tau_{j}^{2})+(1-\gamma_{j})\delta_{0}(\beta_{j}),$
(4) $\displaystyle\gamma_{j}$ $\displaystyle\sim\mathcal{B}ernoulli(\pi),$

where $\gamma_{j}$ ($j=1,\cdots,p$) is a latent variable (which can have a
Bernoulli prior with a fixed probability $\pi$) for feature selection
indicating $\beta_{j}\neq 0$ if $\gamma_{j}=1$ and $\beta_{j}=0$ if
$\gamma_{j}=0$, $\tau_{j}^{2}$ is an additional shrinkage parameter which can
be assigned with an additional prior (e.g. exponential or inverse gamma
prior), and $\delta_{0}(\cdot)$ is the Dirac delta function. Figure 6b shows
the two components of the spike-and-slab mixture distribution. Tang et al.,
(2017) assigned a spike-and-slab mixture double-exponential (i.e. Laplace)
prior for Cox model $\beta_{j}$, i.e.

$\beta_{j}|\gamma_{j},s_{0},s_{1}\sim\gamma_{j}\mathcal{DE}\left(\beta_{j}\Big{|}\frac{1}{s_{0}}\right)+(1-\gamma_{j})\mathcal{DE}\left(\beta_{j}\Big{|}\frac{1}{s_{1}}\right).$

A expectation maximization coordinate descent algorithm was implemented to
optimize the joint posterior distribution of parameters, which provides the
posterior modes of regression coefficients and corresponding frequentist
confidence intervals (Tang et al.,, 2017). Recently, Madjar et al., (2021)
proposed graph-structured feature selection priors for Cox model by assigning
a Markov random field prior on the latent variables, in which the graph helps
to identify pathways of functionally related genes or proteins that are
simultaneously prognostic in different patient cohorts. Formulation (4)
implies independence between the priors of the individual $\beta_{j}$ in the
slab component. In contrast, a variant of the spike-and-slab prior has a
$g$-prior slab (Zellner,, 1986; Held et al.,, 2016)

$\bm{\beta}_{\bm{\gamma}}|\bm{\gamma}\sim\mathcal{N}(\bm{0},g\mathcal{I}_{\bm{\beta}_{\bm{\gamma}},\bm{\beta}_{\bm{\gamma}}}^{-1}),$

where $\bm{\beta}_{\bm{\gamma}}=\\{\beta_{j}:\gamma_{j}=1,j=1,\cdots,p\\}$,
$\bm{\gamma}=\\{\gamma_{j}:j=1,\cdots,p\\}$, $g$ is either a scalar estimated
by Empirical Bayes or assigned with additional prior, and
$\mathcal{I}_{\bm{\beta}_{\bm{\gamma}},\bm{\beta}_{\bm{\gamma}}}^{-1}$ is the
expected Fisher information for $\bm{\beta}_{\bm{\gamma}}$.

Another popular shrinkage prior is the horseshoe prior, a continuous and
global-local shrinkage prior, in which the global parameters allow for sharing
information across omics features and the local parameters allow for
adjustments at the individual omics feature level (Carvalho et al.,, 2009; Mu
et al.,, 2021). In a similar setup to the Cox model with spike-and-slab priors
in (4), a horseshoe prior for the regression coefficient is

$\displaystyle\beta_{j}|\lambda^{2},\tau_{j}^{2}$
$\displaystyle\sim\mathcal{N}(0,\lambda_{j}^{2}\tau^{2}),$
$\displaystyle\lambda_{j}$ $\displaystyle\sim\text{C}^{+}(0,1),$
$\displaystyle\tau$ $\displaystyle\sim\text{C}^{+}(0,1),$

where the local parameter $\lambda_{j}$ and global parameter $\tau$ are both
half-Cauchy distributed $\text{C}^{+}(\cdot,\cdot)$. With the horseshoe prior,
the posterior mean of $\beta_{j}$ will be shrunk by a weight
$\kappa_{j}=\frac{1}{1+\lambda_{j}^{2}}\in(0,1)$ as in Figure 7, where
$\kappa_{j}\to 0$ induces $\beta_{j}\to 0$. Using a user-adjustable cutoff
value, many coefficients can be shrunk to zero, enabling the selection of only
a few associated omics features (with non-zero coefficients).

Figure 7: Density of the horseshoe prior in Bayesian survival modelling. The
shrinkage weight $\kappa_{j}$ close to $0$ shrinks $\beta_{j}$ towards zero,
and $\kappa_{j}$ close to $1$ allows $\beta_{j}$ to escape the shrinkage.

Although Bayesian models can quantify uncertainty of estimators more
straightforward than penalized regressions, most Bayesian Cox-type models for
high-dimensional covariates do not have user-friendly and standalone R
packages on CRAN or GitHub. The main reason is the high computational cost of
running a high-dimensional Bayesian Cox model. Advanced users with programming
capabilities can contact the corresponding authors for original scripts. Since
Bayesian priors are more flexible than frequentist Lasso-type penalties, it
can be easier to extend Bayesian models by changing shrinkage priors while
keeping almost the same algorithm framework. This means that the Bayesian
framework can be more suitable if one is interested in tailoring the shrinkage
effects, e.g. to include prior knowledge about the importance of omics
features, for example features corresponding to a molecular pathway that is
known to be affected in the disease under study (Zhao et al.,, 2023). For more
information on different Bayesian priors in cancer prognosis, we suggest a
recent review by Chu et al., (2022) which summarized other different shrinkage
priors on regression coefficients, such as Gaussian-gamma, Gaussian, Cauchy,
pMOM (product moment distribution), piMOM (product inverse moment
distribution) and peNMIG (parameter-expanded normal-mixture-of-inverse-gamma
distribution) priors.

##  5 Survival model validation

Model validation plays an important role in identifying potential issues, such
as model misspecification or overfitting. This is achieved by revisiting the
model’s specifications and assumptions following model estimation. For
example, the Cox model (1) requires proportional hazards and the logarithm of
the hazard to be linear with respect to the model covariates. The former
assumption can for example be checked by the cumulative Schoenfeld residuals
(Grambsch and Therneau,, 1994; Scheike and Zhang,, 2011), and the latter
assumption by plotting a nonlinear functional form (e.g. spline) for the
effect of a covariate. If the Cox model assumptions are not satisfied, one can
try certain transformations of covariates (e.g. Box-Cox power transformations,
Box and Cox, (1964)), allow time-varying coefficients or model interactions
among covariates (Lin et al.,, 2016; Ng et al.,, 2023), or investigate patient
stratification using unsupervised approaches (Garman et al.,, 2008). However,
the assumption checks are usually suitable only for low-dimensional models,
i.e., for a few clinical variables or a few factors projected from the high-
dimensional omics feature space. Johnson and Long, (2011) used heuristic
methods to investigate the Cox model assumptions by separately fitting
univariate Cox models one feature at one time and testing with the cumulative
Schoenfeld residuals correspondingly. An alternative approach is to loosen the
model assumptions for a more robust modelling approach. One example developed
specifically for feature selection in a high-dimensional space is concordance
regression (Dunkler et al.,, 2010; Pang et al.,, 2012).

###  5.1 Feature stability analysis

One important aspect in model validation when using omics or other high-
dimensional data is the potential instability of feature selection (Kalousis
et al.,, 2007). Feature selection using penalized regressions as described in
Section “Supervised learning via penalized regressions” heavily depends on the
values of the penalty parameters (e.g. for the $\lambda$ parameter in Lasso
Cox model (2)). The penalty parameters are often optimized by cross-validation
(CV) or other resampling methods, and the uncertainty associated with the
random selection of subsets may result in uncertainty in the feature
selection, e.g. different CV folds will typically result in different selected
features. A straightforward way to identify the most stable features is to
find the overlap of identified omics features among different data subsets
(e.g. CV folds or resamples) to avoid high false discovery rate (Zucknick et
al.,, 2008; Zhao and Zucknick,, 2020). One can also perform stability
selection (Meinshausen and Bühlmann,, 2010), which allows to select the most
stable features at a given Type-I error level for a Lasso or Elastic Net Cox
model (Sill et al.,, 2014).

For the Bayesian models in Section “Supervised learning via Bayesian priors”,
feature selection stability is naturally assessed by the uncertainty of
coefficients’ estimators, as reflected in the posterior variances of
$\beta_{j}$ or the posterior selection probabilities
$p(\gamma_{j}|\mathcal{D})$ (in SSVS), which is a natural benefit of utilizing
full Bayesian inference. Although the uncertainty in feature selection
introduces increased variability in the predicted survival probabilities, in
the Bayesian framework, this can be addressed quite naturally by averaging the
survival predictions over all models using Bayesian model averaging (Volinsky
et al.,, 1997). If one is interested in a single model, rather than model
averaging, the median probability model (Barbieri and Berger,, 2004) can be
used for uncertainty analyses in survival and high-dimensional omics data
(Madjar et al.,, 2021).

###  5.2 Survival prediction and calibration

To confirm that the identified clinical and omics features have prognostic
power with respect to the prediction of patients’ survival outcomes, a model
should be both accurate (low prediction error) and precise (low prediction
uncertainty). The simplest way to demonstrate the survival prognostic power is
to dichotomize the prognostic scores (i.e. linear predictor
$\mathbf{X}_{i}\bm{\beta}$ in the Cox model) by its median value, and then use
a log-rank test to compare the survival probabilities of the patients in the
two groups, see Figure 2b. Similarly, one can categorize the prognostic scores
by multiple quantiles (e.g., 25%, 50% median, and 75%) into multiple groups of
patients and perform a log-rank test.

To validate a prediction model systematically (Rahman et al.,, 2017), the
predictive performance of the model is commonly addressed by

  * •

discrimination: the ability of the model to distinguish between low and high
risk patients,

  * •

calibration: the agreement between the observed and predicted survival
probabilities, and

  * •

overall performance: the distance between the observed and predicted survival
probabilities.

####  5.2.1 Discrimination performance

If one focuses on survival prediction at a fixed time point (e.g. 5-year
survival probability), a receiver operating characteristic (ROC) curve can be
used to evaluate the prognostic (i.e. prediction or discrimination) ability of
the survival model, often summarized by its area under the ROC curve (AUC)
(Heagerty et al.,, 2000), see Figure 10a for an example. An AUC of 0.5 is
equivalent to the toss of a coin, and the closer the AUC is to 1, the more
predictive is the model. When making predictions at multiple time points, ROC
curves can be summarized as time-dependent AUC scores, i.e. AUC scores
calculated at prespecified time points. Alternatively, the concordance index
(C-index) provides a more global, time-independent assessment of the
discrimination ability of a prognostic model, such that a better model
predicts higher prognostic scores for patients with shorter survival times
(Harrell et al.,, 1982), i.e.

$C=\mathbb{P}\\{S(t|\mathbf{X}_{i}(\mathbf{t}))<S(t|\mathbf{X}_{j}(\mathbf{t}))\
|T_{i}<T_{j}\ \&\ \delta_{i}=0\\},$

which means in the absence of censoring, any pair of individuals $\\{i,j\\}$
with survival times $T_{i}<T_{j}$ is concordant if and only if
$S(t|\mathbf{X}_{i}(\mathbf{t}))<S(t|\mathbf{X}_{j}(\mathbf{t}))$ (equivalent
to ranking the prognostic scores
$\mathbf{X}_{i}\bm{\beta}>\mathbf{X}_{j}\bm{\beta}$ in a Cox model), where
$\mathbf{t}$ denotes the time instants where there are covariate variations.
The C-index can be expressed as a weighted average of the time-dependent AUC
over time (Heagerty and Zheng,, 2005). Therefore, its interpretation is
similar to the AUC, where a C-index of 0.5 indicates random predictions, while
a perfect prognostic model would have a C-index of 1. There are multiple types
of C-indices for survival modelling, in particular the most frequently used
Harrell’s (Harrell et al.,, 1982) and Uno’s C-index (Uno et al.,, 2011). Uno’s
C-index is more robust than Harrell’s C-index, in case there is dependence of
the concordance on the study-specific censoring distribution (Uno et al.,,
2011).

In the classical Cox modelling framework, both Harrell’s and Uno’s C-indices
only depend on the linear predictors $\mathbf{X}_{i}\bm{\beta}$, which is
independent of $t$. But if a model includes covariates with time-dependent
effects $\bm{\beta}(t)$ and/or time-dependent covariates $\mathbf{X}_{i}(t)$,
Harrell’s and Uno’s C-indices are difficult to be calculated, since they
require the calculation of survival functions for each individual over time.
In this context, Antolini et al., (2005) proposed a time-dependent C-index,
which assesses the concordance of a model’s survival distribution predictions
$\hat{S}(t)$. This means that Antolini’s C-index requires the full
specification of $\hat{S}(t)$, even though a C-index only compares the
survival probabilities between any pair of individuals, that is, it only
assesses whether the relative order of estimated survival probabilities is
concordant with observed survival outcomes (Blanche et al.,, 2019). Time-
dependent prediction indices can better evaluate a model including candidate
features with time-dependent effects and/or time-dependent features. To avoid
C-hacking among different C-indices in model comparison, Sonabend et al.,
(2022) recommended that if all models make survival distribution predictions,
then select a time-dependent C-index; otherwise choose a time-independent
measure (e.g. Uno’s C-index); if there is a combination of risk- and
distribution-predicting models, then choose a transformation method for
analysis (e.g. expected mortality).

####  5.2.2 Calibration performance

Calibration is to quantify the agreement between the observed and predicted
outcomes, which is useful for both internal and external model validation and
is recommended to report routinely. The calibration slope is commonly used
(van Houwelingen,, 2000), which is the slope of the regression of the
observed/actual survival probabilities on the model-predicted survival
probabilities. A survival model can be reported with the estimated $t$-year
survival probabilities in predefined subgroups, denoted as
$S_{\text{model}}(t|g)$ for subgroups $g=1,\cdots,G$. The observed/actual
$t$-year survival probabilities in the subgroups can be estimated by the KM
method, denoted as $S_{\text{KM}}(t|g)$. Using the $\ln(-\ln(\cdot))$-link,
the calibration model is

$\ln(-\ln(S_{\text{KM}}(t|g)))=\alpha+\beta\ln(-\ln(S_{\text{model}}(t|g)))+\epsilon,$
(5)

where $\epsilon$ is an error term. If the intercept $\alpha=0$ and the slope
$\beta=1$, it means that the survival prediction model is well calibrated. For
example, Figure 10b shows a calibration plot, visualizing the calibration of
the estimated 5-year survival probabilities (with 95% confidence interval by
bootstrapping) using the KM method for TCGA-BRCA patients grouped by the
quartiles of Cox-model predicted survival probabilities. Furthermore, one can
calibrate a Cox model in terms of the baseline cumulative hazard and
prognostic score. For non-proportional hazard models, calibration using the
model cumulative hazard function can be considered.

As an alternative to the calibration slope at a single time point, Andres et
al., (2018) and Haider et al., (2020) suggested the distributional
(D)-calibration for accounting survival probabilities across all time points.
This can be useful when assessing the entire post-treatment survival
prognosis, for example, assessing post liver transplantation survival utility
in Andres et al., (2018).

####  5.2.3 Overall performance

Scoring rules can evaluate the accuracy and confidence of probabilistic
predictions, and assess both discrimination and calibration (Gneiting and
Raftery,, 2007; Avati et al.,, 2020). The idea of scoring rules dates back to
Brier, (1950) which assigned a numerical score for verifying ensemble-based
probabilistic predictions of discrete outcomes.

Graf et al., (1999) proposed the time-dependent Brier score, which is the
expected mean-squared error of survival probability prediction at different
time points (Figure 8), i.e.

$\text{BS}(t)=\frac{1}{n}\sum_{i=1}^{n}\left\\{\frac{\hat{S}(t|t_{i})^{2}\mathbbm{1}\\{t_{i}\leq
t\ \&\ \delta_{i}=1\\}}{\hat{G}(t_{i})}\ +\
\frac{\\{1-\hat{S}(t|t_{i})\\}^{2}\mathbbm{1}\\{t_{i}>t\\}}{\hat{G}(t)}\right\\},$

where $t_{i}$ is the survival time of $i$-th individual, $\hat{S}(t|t_{i})$ is
the Cox-model predicted survival probability and $\hat{G}(t)$ is the KM
estimate of the censoring distribution. The benefit of the Brier score is that
it does not only measure discrimination, similar to evaluation measures like
the C-index, but also calibration performance of a survival model. The
integrated Brier score (IBS) is as a single measure of prediction accuracy
integrating BS($t$) over an entire follow-up time period. Hielscher et al.,
(2010) presented a comparison between the IBS and a $D$ measure (Schemper and
Henderson,, 2000), which is an integrated measure based on the mean absolute
deviation rather than the mean-squared error used in IBS. The $D$ measure is
more robust towards extreme observations and has a smaller variance than the
IBS.

To overcome potential overfitting when using feature selection and model
estimation, the survival predictions and model calibration should be evaluated
in an independent validation data set. As validation data are seldom
available, it is recommended to use resampling-based methods for estimating
the uncertainty of the survival model’s performance (Sill et al.,, 2014). This
can be done for example by repeatedly splitting the data to training and
validation sets, and evaluating a survival model’s performance on the
different validation sets using various discrimination or calibration indices.
The .632+ bootstrap estimator for a discrimination or calibration index
(Figure 8) can balance the apparent (training) error and the average out-of-
bag bootstrap error, and in addition accounts for the relative overfitting
based on a no-information error rate in high-dimensional settings (Schumacher
et al.,, 2007; Binder and Schumacher, 2008a, ).

Figure 8: Resampling-based prediction error curves (time-dependent Brier
score) for a Lasso Cox model applied to TCGA-BRCA data. The blue line shows
the prediction error curve of the null model (estimated by the Kaplan-Meier
estimator without covariate information), the red line shows the full apparent
error estimates (i.e., the errors as estimated when applying the model to the
entire training data set), and the black line corresponds to the .632+
bootstrap error estimates. The gray area indicates the pointwise 2.5% and
97.5% quantiles of the bootstrap error using 100 out-of-bag samples.

####  5.2.4 Graphical representation

After confirming that a model is valid (assumptions hold), accurate (low
prediction error), precise (uncertainty of performance measures properly
quantified) and its predictions generalizable beyond the training data set
(using independent validation data if available), a prognostic nomogram
(Kattan et al.,, 1999) can be used to summarize the prognostic effect of the
identified clinical and omics features on the risk of a specific year’s
overall survival (Figure 9), which may help the clinicians to enhance the
patient management and personalized treatment strategies. For example, the red
colored dots in Figure 9 show the information of the identified five variables
from an example patient and the corresponding scoring points. The summed
scoring points of 263 maps to the predicted 1-year, 3-year and 5-year survival
probabilities of this patient. Note that most nomograms treat the identified
variables independently in the risk calculation, even though there may be
significant interactions among the model features that were used in the
feature selection step. However, visualizing such interaction effects would
make the nomograms less accessible and interpretable, and so, there is still a
room for improvement in how to translate the multi-variate risk scores into
clinical practice.

Figure 9: Nomogram developed to estimate the overall survival probability for
TCGA-BRCA patients based on clinical (age) demographic (ethnicity) and three
selected mRNA features form a Lasso Cox model. The red coloured symbols
represent example patient’s information and predicted probabilities of 1-year,
3-year and 5-year survival.

When an independent validation data set is available, it is recommended to
report a calibration plot corresponding to the nomogram. Using independent
validation data to obtain $S_{\text{model}}(t|g)$ in the calibration model (5)
is for the generalization capacity of the model. Since we here do not have
independent validation data besides TCGA-BRCA data, Figure 10 shows an example
calibration plot at 5-year survival evaluation time point based on the built
Cox model in Figure 9 for a split 20% TCGA-BRCA data set.

Figure 10: Discrimination and calibration of survival prediction for 20% TCGA-
BRCA validation data. TCGA-BRCA data were split to a 80% training set and a
20% validation set. (a) Receiver operating characteristic (ROC) curve
estimated at 5-year survival evaluation time point based on the Cox model in
Figure 9. The AUC score is the area under the ROC curve. The 45-degree line
represents the performance of a random prediction of the outcome event with
AUC = 0.5. (b) Calibration plot estimating 5-year survival probabilities in
TCGA-BRCA patients grouped by the quartiles of the predicted probabilities
based on the Cox model in Figure 9. The actual 5-year survival probabilities
in each group (green colored dot and error bar) were estimated by the Kaplan-
Meier (KM) method with a 95% confidence interval by bootstrapping. The blue
colored cross indicates a bias-corrected KM estimate. The predicted survival
probability of each quartile group was the mean of the 5-year survival
probabilities based on the Cox model for the corresponding group.

##  6 Beyond penalized and Bayesian Cox models

In this tutorial, we mainly focused on penalized regressions and Bayesian
hierarchical models in the Cox proportional hazards framework. One can extend
this framework in several ways. For instance, stay in a likelihood-based
modelling framework, but replace the partial likelihood function of the semi-
parametric Cox model by alternative likelihood functions (which do not
necessarily need to imply proportional hazards), e.g. for parametric survival
models like exponential, Weibull, or accelerated failure time (AFT) models, or
for Aalen’s additive hazard model (Gorst-Rasmussen and Scheike,, 2012; Wang et
al.,, 2019). Alternatively, one can move to a more algorithmic machine
learning approach, such as tree-based boosting or bagging methods, e.g. random
survival forests (Hothorn et al.,, 2006; Ishwaran et al.,, 2008; Binder and
Schumacher, 2008b, ; Jaeger et al.,, 2019; Hornung and Wright,, 2019), and
(deep) neural networks (Wiegrebe et al.,, 2023).

Hothorn et al., (2006) and Ishwaran et al., (2008) introduced ensemble tree
methods for analyzing right-censored survival data, which construct ensembles
from base learners, e.g., binary survival trees for each omics feature.
Hothorn et al., (2006) also proposed a gradient boosting algorithm to predict
the survival time of patients with acute myeloid leukemia (AML), based on
clinical and omics features. Similarly, Binder and Schumacher, 2008b
developed a likelihood-based boosting method, which aims to maximize the Cox
partial likelihood function, for modelling time-to-event data based on high-
dimensional omics input data and which also allows the inclusion of a small
number of mandatory covariates. In general, one needs to be cautious if using
some machine learning methods that are not well-suited for high-dimensional
features. For example, Kvamme et al., (2019) proposed extensions of the Cox
model with neural networks, which are only valid if the number of covariates
is smaller than the number of samples, i.e. if $p<n$. A systematic review of
deep learning for survival analysis, which includes a survey of methods
suitable for high-dimensional data ($p>n$), is provided by Wiegrebe et al.,
(2023).

In the case of non-proportional hazards, many likelihood-based survival models
beyond Cox-type models have also been extended to account for high-dimensional
omics as input data. For example, Ma et al., (2006) combined Lin and Yin’s
additive hazard model (Lin and Ying,, 1994) with principal component
regression for dimension reduction of omics features, which was applied to the
study of gene expression-based survival prognosis for diffuse large B-cell
lymphoma. Engler and Li, (2009) added the Elastic Net penalty in an
accelerated ATF model, which assumes that the effect of a covariate
accelerates or decelerates the life course of patients. Schmid and Hothorn,
(2008) and Barnwal et al., (2022) used boosting algorithms to learn parametric
AFT models. Ha et al., (2014) considered the Lasso, SCAD and a penalized
h-likelihood for feature selection in frailty models which assume that
individuals have unobserved heterogeneity captured by a latent random term
$Z$, which adapts the Cox model (1) into $h(t|X)=Zh_{0}(t)\exp(X\beta)$.

####  6.0.1 Advanced survival models: cure models, competing risks and multi-
state models

In some situations, survival data may be different from Figure 1c (also
Section “Time-to-event data”), where it was presumed that all individuals will
eventually experience the event of interest. Liu et al., (2012) studied the
Lasso and SCAD feature selection for the proportional hazard mixture cure
model, in which a certain fraction of individuals will never experience the
event of interest. Tapak et al., (2015) investigated Lasso, Elastic Net and
likelihood-based boosting for microarray-based survival modelling with
competing risks, such as “progression” versus “death from non-cancer cause”,
i.e., the event of a patient can occur due to one of multiple distinct causes.
For a single individual who can experience several possible events, Dutta et
al., (2016) proposed a multi-state model to identify risk factors in different
stages of disease based on high-dimensional input data.

##  7 Towards single-cell data analysis

The cellular heterogeneity of complex sample mixtures pose challenges and also
opportunities for precision medicine and survival prediction. For example,
Zhou et al., (2019) showed that tumor microenvironment-related gene expression
signatures do not only accurately predict the survival among colon cancer
patients, but also serve as biomarkers for identifying patients who could
potentially benefit from adjuvant chemotherapy. scRNA-seq technologies provide
an unprecedented opportunity for dissecting the interplay between the cancer
cells and the associated tumor microenvironment, and the produced high-
dimensional data should also augment existing survival modelling approaches.
The emerging single-cell atlas will provide a detailed and quantitative
overview of tissue composition and organization, and advance both biomedical
research and clinical practice (Elmentaite et al.,, 2022).

Currently, the focus of statistical model development for single-cell data
analysis is to understand cell type composition and its impact on gene
regulation and transcriptional dynamics, usually based on only a small number
of samples/individuals. On the one hand, the underlying statistical models for
single-cell data analysis are still in development and continuously being re-
evaluated and challenged (Andrews et al.,, 2021; Kharchenko,, 2021). On the
other hand, survival analysis tackles disease or health problems on the
individual level, and usually requires a relatively large number of
individuals for sufficient statistical power. In particular, the development
of improved computational methods is urgently needed to enable the
consideration of multiple layers (e.g. individual-level, cellular-level,
molecular-level) of information when integrating groups of individuals and
omics data from a variety of molecular modalities. Therefore, current
approaches often map the findings from single-cell omics data to large-scale
bulk sequencing omics and survival data, rather than jointly analyzing single-
cell omics and survival data. For example, Guo et al., (2018) introduced a
generalizable approach to first study T-cells in 14 non-small cell lung cancer
(NSCLC) patients and to identify gene signatures from the tumor-enriched T
cell clusters, and then investigated these gene signatures using bulk RNA-seq
and survival data from larger TCGA-NSCLC patient cohort. Similarly, Zhang et
al., (2020) developed a scRNA-seq-based approach to reconstruct a multilayer
signaling network based on 16 glioma patients, and then investigated the
network genes using survival data from TCGA Chinese glioma genome atlas (TCGA-
CGGA) patients. However, direct joint analysis of survival and single-cell
omics data from multiple cellular hierarchies requires further methodological
developments and new statistical and machine learning methods.

##  8 Discussion

Although survival analysis faces many modelling challenges, mainly due to
censored outcomes, it represents a well-established methodology for finding
risk factors associated with patients’ survival. The identification of omics
biomarkers for survival prognosis may provide systematic means to guide
patient management and personalized treatment and diagnostic strategies. In
this tutorial, we provided a comprehensive workflow for survival analysis with
omics and clinical data, with a specific focus on feature selection of
survival-associated omics features and survival model validation. We covered
many penalized regressions and Bayesian models for feature selection and
survival prediction, accounting for their specific assumptions and
applications. Examples of real data and R scripts have been made available to
illustrate the use of different methods, which should help researchers to
choose and apply suitable methods for their survival analysis applications
(https://ocbe-uio.github.io/survomics/survomics.html). We note that this
review only considers methods for right-censored time-to-event data, i.e.
where all individuals are assumed to be followed continuously from the start
of the study, but where the follow-up period might end before the event (e.g.
death) was observed. Other types of censoring include interval censoring and
left truncation, and appropriate statistical methods dealing with these
censoring patterns should be chosen accordingly.

Most of the current methods for survival analysis do not take into account the
complex structures within and between multi-omics data, such as gene
regulation and DNA-protein interactions. Regulatory networks constructed
either based on prior biological knowledge or using data-driven, yet
biologically explainable approaches, may help establish useful methodologies
for survival analysis that are more effective for deriving biological insights
as well as enable improved clinical translation.

Another limitation of most of the reviewed methods is that they identify omics
features predictive of survival, but they cannot determine causal
relationships. Causality is a fundamental notion to understand omics features
causing disease progression, which will allow one to reliably intervene omics
features for targeted therapies. There are two popular models, Pearl’s
structural causal model (SCM) and Rubin’s causal model (RCM), both of which
introduce perturbations to draw causal inference. Farooq et al., (2023)
utilized SCM-based causal discovery approaches to unravel relationships
between omics features and survival of breast cancer patients. However, to
identify reliable causal relations for clinical applications, laboratory-based
experiments, e.g. clustered regularly interspaced short palindromic repeats
(CRISPR) techniques (Wang and Doudna,, 2023), are often necessary to verify
the functional relevance of the identified omics features. Causal mediation
analysis is an important tool, which considers the problem of decomposing the
causal effect of treatment/exposure into direct and indirect effects (Lange
and Hansen,, 2011; VanderWeele,, 2011). The direct effect corresponds to the
effect of a treatment directly on the survival outcome, while an indirect
effect corresponds to the effect of a treatment on the outcome that is due to
its effect on an intermediate variable (e.g. gene expression) that also has a
causal effect on the survival outcome. High-dimensional RCM-based mediation
analysis has been used to investigate the indirect effect transmitted by omics
features between an exposure and survival outcomes (Song et al.,, 2020, 2021).
Targeted learning also fills a much needed gap between statistical modelling
and causal inference (van der Laan and Rose,, 2011, 2018). Tuglus and van der
Laan, (2011) used targeted maximum likelihood estimation to provide an
interpretable causal measure of variable importance for the discovery of
biomarkers and omics features. Another way to formalize personalized medicine
is dynamic treatment regimes (Tsiatis et al.,, 2019; Deliu et al.,, 2023) that
encompasses causal inference and takes into account for the variability in
omics, environment and lifestyle factors for each individual to improve the
treatment of a particular patient.

## Data availability

Supplementary step-by-step R tutorial is available online at https://ocbe-
uio.github.io/survomics/survomics.html. TCGA data is publicly available at
https://portal.gdc.cancer.gov.

## Acknowledgements

The authors thank Professor Ørnulf Borgan for his valuable suggestions and
comments. This work was supported by grants from Helse Sør-Øst (grant 2020026
to TA), the Norwegian Cancer Society (216104 to TA), the Radium Hospital
Foundation, the Academy of Finland (grants 326238, 340141, 344698, 345803 to
TA), Cancer Society of Finland (to TA), the European Union’s Horizon 2020
research and innovation programme (grant ‘PANCAIM’ 101016851 to JZ and TA,
grant ‘RESCUER’ 847912 to MZ), the Innovative Medicines Initiative 2 Joint
Undertaking of the European Union’s Horizon 2020 research and innovation
programme and EFPIA and JDRF INTERNATIONAL (grant ‘imSAVAR’ 853988 to MZ), ERA
PerMed under the ERA-NET Cofund scheme of the European Union’s Horizon 2020
research and innovation framework programme (grant ERAPerMed2021_330_SYMMETRY
to MZ).

## References

  * Aittokallio, (2009) Aittokallio, T. (2009).  Dealing with missing values in large-scale studies: microarray data imputation and beyond.  Briefings in Bioinformatics, 11(2):253–264.
  * Akbani et al., (2014) Akbani, R., Ng, P. K. S., Werner, H. M. J., Shahmoradgoli, M., Zhang, F., Ju, Z., Liu, W., Yang, J.-Y., Yoshihara, K., Li, J., Ling, S., Seviour, E. G., Ram, P. T., Minna, J. D., Diao, L., Tong, P., Heymach, J. V., Hill, S. M., Dondelinger, F., Städler, N., Byers, L. A., Meric-Bernstam, F., Weinstein, J. N., Broom, B. M., Verhaak, R. G. W., Liang, H., Mukherjee, S., Lu, Y., and Mills, G. B. (2014).  A pan-cancer proteomic perspective on the cancer genome atlas.  Nature Communications, 5:3887.
  * Andres et al., (2018) Andres, A., Montano-Loza, A., Greiner, R., Uhlich, M., Jin, P., Hoehn, B., Bigam, D., Shapiro, J. A. M., and Kneteman, N. M. (2018).  A novel learning algorithm to predict individual survival after liver transplantation for primary sclerosing cholangitis.  PLOS ONE, 13(3):e0193523.
  * Andrews et al., (2021) Andrews, T. S., Kiselev, V. Y., McCarthy, D., and Hemberg, M. (2021).  Tutorial: guidelines for the computational analysis of single-cell rna sequencing data.  Nature Protocols, 16(1):1–9.
  * Antolini et al., (2005) Antolini, L., Boracchi, P., and Biganzoli, E. (2005).  A time-dependent discrimination index for survival data.  Statistics in Medicine, 24(24):3927–3944.
  * Avati et al., (2020) Avati, A., Duan, T., Zhou, S., Jung, K., Shah, N. H., and Ng, A. Y. (2020).  Countdown Regression: Sharp and Calibrated Survival Predictions.  In Adams, R. P. and Gogate, V., editors, Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pages 145–155. PMLR.
  * Bair and Tibshirani, (2004) Bair, E. and Tibshirani, R. (2004).  Semi-supervised methods to predict patient survival from gene expression data.  PLOS Biology, 2(4):e108.
  * Barbieri and Berger, (2004) Barbieri, M. M. and Berger, J. O. (2004).  Optimal predictive model selection.  The Annals of Statistics, 32(3):870–897.
  * Barnwal et al., (2022) Barnwal, A., Cho, H., and Hocking, T. (2022).  Survival regression with accelerated failure time model in XGBoost.  Journal of Computational and Graphical Statistics, 31(4):1292–1302.
  * Bartel, (2018) Bartel, D. P. (2018).  Metazoan MicroRNAs.  Cell, 173(1):20–51.
  * (11) Binder, H. and Schumacher, M. (2008a).  Adapting prediction error estimates for biased complexity selection in high-dimensional bootstrap samples.  Statistical Applications in Genetics and Molecular Biology, 7:1.
  * (12) Binder, H. and Schumacher, M. (2008b).  Allowing for mandatory covariates in boosting estimation of sparse high-dimensional survival models.  BMC Bioinformatics, 9:1.
  * Blanche et al., (2019) Blanche, P., Kattan, M. W., and Gerds, T. A. (2019).  The c-index is not proper for the evaluation of t-year predicted risks.  Biostatistics, 20(2):347–357.
  * Bommert et al., (2022) Bommert, A., Welchowski, T., Schmid, M., and Rahnenführer, J. (2022).  Benchmark of filter methods for feature selection in high-dimensional gene expression survival data.  Briefings in Bioinformatics, 23(1):bbab354.
  * Bøvelstad et al., (2007) Bøvelstad, H., Nygård, S., Størvold, H., Aldrin, M., Borgan, O., Frigessi, A., and Lingjærde, O. (2007).  Predicting survival from microarray data - a comparative study.  Bioinformatics, 23(16):2080–2087.
  * Bøvelstad et al., (2009) Bøvelstad, H. M., Nygård, S., and Borgan, O. (2009).  Survival prediction from clinico-genomic models - a comparative study.  BMC Bioinformatics, 10:413.
  * Box and Cox, (1964) Box, G. E. P. and Cox, D. R. (1964).  An analysis of transformations.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 26(2):211–252.
  * Bradburn et al., (2003) Bradburn, M. J., Clark, T. G., Love, S. B., and Altman, D. G. (2003).  Survival analysis part II: Multivariate data analysis – an introduction to concepts and methods.  British Journal of Cancer, 89(3):431–436.
  * Brier, (1950) Brier, G. W. (1950).  Verification of forecasts expressed in terms of probability.  Monthly Weather Review, 78(1):1–3.
  * Carins et al., (2011) Carins, R., Harris, I., and Mak, T. (2011).  Regulation of cancer cell metabolism.  Nature Reviews Cancer, 11(2):85–95.
  * Carvalho et al., (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In van Dyk, D. and Welling, M., editors, Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, volume 5 of Proceedings of Machine Learning Research, pages 73–80, Florida USA. PMLR.
  * Chekouo et al., (2015) Chekouo, T., Stingo, F. C., Doecke, J. D., and Do, K.-A. (2015).  mirna-target gene regulatory networks: A Bayesian integrative approach to biomarker selection with application to kidney cancer.  Biometrics, 71(2):428–438.
  * Chu et al., (2022) Chu, J., Sun, N. A., Hu, W., Chen, X., Yi, N., and Shen, Y. (2022).  The application of Bayesian methods in cancer prognosis and prediction.  Cancer Genomics & Proteomics, 19(1):1–11.
  * Clark et al., (2003) Clark, T. G., Bradburn, M., Love, S. B., and Altman, D. G. (2003).  Survival analysis part I: Basic concepts and first analyses.  British Journal of Cancer, 89:232–238.
  * Cox, (1972) Cox, D. R. (1972).  Regression models and life-tables.  Journal of the Royal Statistical Society: Series B (Methodological), 34(2):187–202.
  * de Brevern et al., (2004) de Brevern, A., Hazout, S., and Malpertuy, A. (2004).  Influence of microarrays experiments missing values on the stability of gene groups by hierarchical clustering.  BMC Bioinformatics, 5:114.
  * Deliu et al., (2023) Deliu, N., Williams, J. J., and Chakraborty, B. (2023).  Reinforcement learning in modern biostatistics: Constructing optimal adaptive interventions.  arXiv:2203.02605.
  * Dunkler et al., (2010) Dunkler, D., Schemper, M., and Heinze, G. (2010).  Gene selection in microarray survival studies under possibly non-proportional hazards.  Bioinformatics, 26(6):784–790.
  * Dutta et al., (2016) Dutta, S., Datta, S., and Datta, S. (2016).  Temporal prediction of future state occupation in a multistate model from high-dimensional baseline covariates via pseudo-value regression.  Journal of Statistical Computation and Simulation, 87(7):1363–1378.
  * Edwards et al., (2015) Edwards, N. J., Oberti, M., Thangudu, R. R., Cai, S., McGarvey, P. B., Jacob, S., Madhavan, S., and Ketchum, K. A. (2015).  The CPTAC data portal: A resource for cancer proteomics research.  Journal of Proteome Research, 14(6):2707–2713.
  * Elmentaite et al., (2022) Elmentaite, R., Conde, C. D., Yang, L., and Teichmann, S. A. (2022).  Single-cell atlases: shared and tissue-specific cell types across human organs.  Nature Reviews Genetics, 23(7):395–410.
  * Engler and Li, (2009) Engler, D. and Li, Y. (2009).  Survival analysis with high-dimensional covariates: An application in microarray studies.  Statistical Applications in Genetics and Molecular Biology, 8(1):1–22.
  * Fan et al., (2010) Fan, J., Feng, Y., and Wu, Y. (2010).  High-dimensional variable selection for Cox’s proportional hazards model.  IMS Collections, 6:70–86.
  * Fan and Li, (2002) Fan, J. and Li, R. (2002).  Variable Selection for Cox’s proportional Hazards Model and Frailty Model.  The Annals of Statistics, 30(1):74–99.
  * Farooq et al., (2023) Farooq, M., Hardan, S., Zhumbhayeva, A., Zheng, Y., Nakov, P., and Zhang, K. (2023).  Understanding breast cancer survival: Using causality and language models on multi-omics data.  arXiv:2305.18410.
  * Garman et al., (2008) Garman, K. S., Acharya, C. R., Edelman, E., Grade, M., Gaedcke, J., Sud, S., Barry, W., Diehl, A. M., Provenzale, D., Ginsburg, G. S., Ghadimi, B. M., Ried, T., Nevins, J. R., Mukherjee, S., Hsu, D., and Potti, A. (2008).  A genomic approach to colon cancer risk stratification yields biologic insights into therapeutic opportunities.  Proceedings of the National Academy of Sciences, 105(49):19432–19437.
  * George and McCulloch, (1993) George, E. and McCulloch, R. (1993).  Variable selection via Gibbs sampling.  Journal of the American Statistical Association, 88(423):881–889.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378.
  * Gorst-Rasmussen and Scheike, (2012) Gorst-Rasmussen, A. and Scheike, T. H. (2012).  Coordinate descent methods for the penalized semiparametric additive hazards model.  Journal of Statistical Software, 47(9):1–17.
  * Graf et al., (1999) Graf, E., Schmoor, C., Sauerbrei, W., and Schumacher, M. (1999).  Assessment and comparison of prognostic classification schemes for survival data.  Statistics in Medicine, 18(17-18):2529–2545.
  * Grambsch and Therneau, (1994) Grambsch, P. M. and Therneau, T. M. (1994).  Proportional hazards tests and diagnostics based on weighted residuals.  Biometrika, 81(3):515–526.
  * Green and Symons, (1983) Green, M. S. and Symons, M. J. (1983).  A comparison of the logistic risk function and the proportional hazards model in prospective epidemiologic studies.  Journal of Chronic Diseases, 36(10):715–723.
  * Guo et al., (2018) Guo, X., Zhang, Y., Zheng, L., Zheng, C., Song, J., Zhang, Q., Kang, B., Liu, Z., Jin, L., Xing, R., Gao, R., Zhang, L., Dong, M., Hu, X., Ren, X., Kirchhoff, D., Roider, H. G., Yan, T., and Zhang, Z. (2018).  Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing.  Nature Medicine, 24(7):978–985.
  * Ha et al., (2014) Ha, I. D., Pan, J., Oh, S., and Lee, Y. (2014).  Variable selection in general frailty models using penalized h-likelihood.  Journal of Computational and Graphical Statistics, 23(4):1044–1060.
  * Haider et al., (2020) Haider, H., Hoehn, B., Davis, S., and Greiner, R. (2020).  Effective ways to build and evaluate individual survival distributions.  Journal of Machine Learning Research, 21(85):1–63.
  * Harrell et al., (1982) Harrell, Frank E., J., Califf, R. M., Pryor, D. B., Lee, K. L., and Rosati, R. A. (1982).  Evaluating the yield of medical tests.  JAMA, 247(18):2543–2546.
  * Hasin et al., (2021) Hasin, Y., Seldin, M., and Lusis, A. (2021).  Multi-omics approaches to disease.  Genome Biology, 18:83.
  * Haukaas et al., (2017) Haukaas, T. H., Euceda, L. R., Giskeødegård, G. F., and Bathen, T. F. (2017).  Metabolic portraits of breast cancer by HR MAS MR spectroscopy of intact tissue samples.  Metabolites, 7(2):18.
  * Heagerty et al., (2000) Heagerty, P. J., Lumley, T., and Pepe, M. S. (2000).  Time-Dependent ROC Curves for Censored Survival Data and a Diagnostic Marker.  Biometrics, 56(2):337–344.
  * Heagerty and Zheng, (2005) Heagerty, P. J. and Zheng, Y. (2005).  Survival model predictive accuracy and ROC curves.  Biometrics, 61(1):92–105.
  * Held et al., (2016) Held, L., Gravestock, I., and Sabanés Bové, D. (2016).  Objective Bayesian model selection for Cox regression.  Statistics in Medicine, 35(29):5376–5390.
  * Henderson, (1995) Henderson, R. (1995).  Problems and prediction in survival-data analysis.  Statistics in Medicine, 14(2):161–184.
  * Herrmann et al., (2021) Herrmann, M., Probst, P., Hornung, R., Jurinovic, V., and Boulesteix, A.-L. (2021).  Large-scale benchmark study of survival prediction methods using multi-omics data.  Briefings in Bioinformatics, 22(3):bbaa167.
  * Hielscher et al., (2010) Hielscher, T., Zucknick, M., Werft, W., and Benner, A. (2010).  On the prognostic value of survival models with application to gene expression signatures.  Statistics in Medicine, 29(7-8):818–829.
  * Hornung and Wright, (2019) Hornung, R. and Wright, M. N. (2019).  Block forests: random forests for blocks of clinical and omics covariate data.  BMC Bioinformatics, 20:358.
  * Hothorn et al., (2006) Hothorn, T., Bühlmann, P., Dudoit, S., Molinaro, A., and van der Laan, M. J. (2006).  Survival ensembles.  Biostatistics, 7(3):355–373.
  * Hrydziuszko and Viant, (2012) Hrydziuszko, O. and Viant, M. R. (2012).  Missing values in mass spectrometry based metabolomics: an undervalued step in the data processing pipeline.  Metabolomics, 8(Suppl 1):S161–S174.
  * Ishwaran et al., (2008) Ishwaran, H., Kogalur, U. B., Blackstone, E. H., and Lauer, M. S. (2008).  Random survival forests.  The Annals of Applied Statistics, 2(3):841–860.
  * Jaeger et al., (2019) Jaeger, B. C., Long, D. L., Long, D. M., Sims, M., Szychowski, J. M., Min, Y.-I., Mcclure, L. A., Howard, G., and Simon, N. (2019).  Oblique random survival forests.  The Annals of Applied Statistics, 13(3):1847–1883.
  * Jendoubi, (2021) Jendoubi, T. (2021).  Approaches to integrating metabolomics and multi-omics data: A primer.  Metabolites, 11(3):184.
  * Jiang et al., (2022) Jiang, R., Sun, T., Song, D., and Li, J. J. (2022).  Statistics or biology: the zero-inflation controversy about scRNA-seq data.  Genome Biology, 23(1):31.
  * Johnson and Long, (2011) Johnson, B. A. and Long, Q. (2011).  Survival ensembles by the sum of pairwise differences with application to lung cancer microarray studies.  The Annals of Applied Statistics, 5(2A):1081–1101.
  * Kalousis et al., (2007) Kalousis, A., Prados, J., and Hilario, M. (2007).  Stability of feature selection algorithms: a study on high-dimensional spaces.  Knowledge and Information Systems, 12(1):95–116.
  * Kaplan and Meier, (1958) Kaplan, E. L. and Meier, P. (1958).  Nonparametric estimation from incomplete observations.  Journal of the American Statistical Association, 53(282):457–481.
  * Kattan et al., (1999) Kattan, M. W., Wheeler, T. M., and Scardino, P. T. (1999).  Postoperative nomogram for disease recurrence after radical prostatectomy for prostate cancer.  Journal of Clinical Oncology, 17(5):1499–1499.
  * Kharchenko, (2021) Kharchenko, P. V. (2021).  The triumphs and limitations of computational methods for scRNA-seq.  Nature Methods, 18(7):723–732.
  * Kim et al., (2012) Kim, J., Sohn, I., Jung, S.-H., Kim, S., and Park, C. (2012).  Analysis of survival data with group lasso.  Communications in Statistics - Simulation and Computation, 41(9):1593–1605.
  * Kong et al., (2022) Kong, W., Hui, H. W. H., Peng, H., and Goh, W. W. B. (2022).  Dealing with missing values in proteomics data.  Protenomics, 22(23-24):2200092.
  * Konrath et al., (2013) Konrath, S., Fahrmeir, L., and Kneib, T. (2013).  Bayesian smoothing, shrinkage and variable selection in hazard regression.  In Becker, C., Fried, R., and Kuhnt, S., editors, Robustness and Complex Data Structures, pages 149–170. Springer, Berlin.
  * Korthauer et al., (2019) Korthauer, K., Kimes, P. K., Duvallet, C., Reyes, A., Subramanian, A., Teng, M., Shukla, C., Alm, E. J., and Hicks, S. C. (2019).  A practical guide to methods controlling false discoveries in computational biology.  Genome Biology, 20:118.
  * Kundu et al., (2018) Kundu, S., Cheng, Y., Shin, M., Manyam, G., Mallick, B. K., and Baladandayuthapani, V. (2018).  Bayesian variable selection with graphical structure learning: Applications in integrative genomics.  PLOS ONE, 13(7):e0195070.
  * Kvamme et al., (2019) Kvamme, H., Ørnulf Borgan, and Scheel, I. (2019).  Time-to-event prediction with neural networks and Cox regression.  Journal of Machine Learning Research, 20(129):1–30.
  * Lange and Hansen, (2011) Lange, T. and Hansen, J. V. (2011).  Direct and indirect effects in a survival context.  Epidemiology, 22(4):575–581.
  * Lee et al., (2011) Lee, K. H., Chakraborty, S., and Sun, J. (2011).  Bayesian variable selection in semiparametric proportional hazards model for high dimensional survival data.  The International Journal of Biostatistics, 7(1):Article 21.
  * Lee et al., (2015) Lee, K. H., Chakraborty, S., and Sun, J. (2015).  Survival prediction and variable selection with simultaneous shrinkage and grouping priors.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 8(2):114–127.
  * Leung et al., (1997) Leung, K., Elashoff, R., and Afifi, A. (1997).  Censoring issues in survival analysis.  Annual Review of Public Health, 18:83–104.
  * Li et al., (2022) Li, X., Li, M., Xiang, J., Zhao, Z., and Shang, X. (2022).  SEPA: signaling entropy-based algorithm to evaluate personalized pathway activation for survival analysis on pan-cancer data.  Bioinformatics, 38(9):2536–2543.
  * Lin and Ying, (1994) Lin, D. Y. and Ying, Z. (1994).  Semiparametric analysis of the additive risk model.  Biometrika, 81(1):61–71.
  * Lin et al., (2016) Lin, H., Tan, M. T., and Li, Y. (2016).  A semiparametrically efficient estimator of single-index varying coefficient Cox proportional hazards models.  Statistica Sinica, 26(2):779–807.
  * Liu et al., (2012) Liu, X., Peng, Y., Tu, D., and Liang, H. (2012).  Variable selection in semiparametric cure models based on penalized likelihood, with application to breast cancer clinical trials.  Statistics in Medicine, 31(24):2882–2891.
  * Ma et al., (2006) Ma, S., Kosorok, M. R., and Fine, J. P. (2006).  Additive risk models for survival data with high-dimensional covariates.  Biometrics, 62(1):202–210.
  * Madjar et al., (2021) Madjar, K., Zucknick, M., Ickstadt, K., and Rahnenführer, J. (2021).  Combining heterogeneous subgroups with graph-structured variable selection priors for Cox regression.  BMC Bioinformatics, 22:586.
  * McInnes et al., (2018) McInnes, L., Healy, J., and Melville, J. (2018).  UMAP: Uniform manifold approximation and projection for dimension reduction.  arXiv:1802.03426.
  * Meinshausen and Bühlmann, (2010) Meinshausen, N. and Bühlmann, P. (2010).  Stability selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72(4):417–473.
  * Mu et al., (2021) Mu, J., Liu, Q., Kuo, L., and Hu, G. (2021).  Bayesian variable selection for the Cox regression model with spatially varying coefficients with applications to louisiana respiratory cancer data.  Biometrical Journal, 63(8):1607–1622.
  * Nam et al., (2022) Nam, S. Y., Jeong, J., Lee, W. K., and Jeon, S. W. (2022).  Sex-specific effect of body mass index and fasting glucose on gastric cancer risk and all causes mortality; a cohort study of 5.17 million.  International Journal of Obesity, 46(9):1644–1651.
  * Ng et al., (2023) Ng, H. M., Jiang, B., and Wong, K. Y. (2023).  Penalized estimation of a class of single-index varying-coefficient models for integrative genomic analysis.  Biometrical Journal, 65(1):2100139.
  * Pang et al., (2012) Pang, H., George, S. L., Hui, K., and Tong, T. (2012).  Gene selection using iterative feature elimination random forests for survival outcomes.  IEEE/ACM Transactions on Computational Biology and Bioinformatics, 9(5):1422–1431.
  * Pearson, (1901) Pearson, K. (1901).  LIII. on lines and planes of closest fit to systems of points in space.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11):559–572.
  * Perou et al., (2000) Perou, C., Sørlie, T., Eisen, M., van de Rijn, M., Jeffrey, S., Rees, C., Pollack, J., Ross, D., Johnsen, H., Akslen, L., Fluge, O., Pergamenschikov, A., Williams, C., Zhu, S., Lonning, P., Børresen-Dale, A., Brown, P., and Botstein, D. (2000).  Molecular portraits of human breast tumours.  Nature, 406(6797):747–752.
  * Peto and Peto, (1972) Peto, R. and Peto, J. (1972).  Asymptotically efficient rank invariant test procedures.  Journal of the Royal Statistical Society. Series A (General), 135(2):185.
  * Rahman et al., (2017) Rahman, M. S., Ambler, G., Choodari-Oskooei, B., and Omar, R. Z. (2017).  Review and evaluation of performance measures for survival prediction models in external validation settings.  BMC Medical Research Methodology, 17:60.
  * Scheike and Zhang, (2011) Scheike, T. H. and Zhang, M.-J. (2011).  Analyzing competing risk data using the R timereg package.  Journal of Statistical Software, 38(2):1–15.
  * Schemper and Henderson, (2000) Schemper, M. and Henderson, R. (2000).  Predictive accuracy and explained variation in Cox regression.  Biometrics, 56(1):249–255.
  * Schmid and Hothorn, (2008) Schmid, M. and Hothorn, T. (2008).  Flexible boosting of accelerated failure time models.  BMC Bioinformatics, 9:269.
  * Schumacher et al., (2007) Schumacher, M., Binder, H., and Gerds, T. (2007).  Assessment of survival prediction models based on microarray data.  Bioinformatics, 23(14):1768–1774.
  * Sill et al., (2014) Sill, M., Hielscher, T., Becker, N., and Zucknick, M. (2014).  c060: Extended inference with lasso and elastic-net regularized Cox and generalized linear models.  Journal of Statistical Software, 62(5):1–22.
  * Simon et al., (2011) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011).  Regularization paths for Cox’s proportional hazards model via coordinate descent.  Journal of Statistical Software, 39(5):1–13.
  * Simon et al., (2013) Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2013).  A sparse-group lasso.  Journal of Computational and Graphical Statistics, 22(2):231–245.
  * Sonabend et al., (2022) Sonabend, R., Bender, A., and Vollmer, S. (2022).  Avoiding C-hacking when evaluating survival distribution predictions with discrimination measures.  Bioinformatics, 38(17):4178–4184.
  * Song et al., (2021) Song, Y., Zhou, X., Kang, J., Aung, M. T., Zhang, M., Zhao, W., Needham, B. L., Kardia, S. L. R., Liu, Y., Meeker, J. D., Smith, J. A., and Mukherjee, B. (2021).  Bayesian hierarchical models for high-dimensional mediation analysis with coordinated selection of correlated mediators.  Statistics in Medicine, 40(27):6038–6056.
  * Song et al., (2020) Song, Y., Zhou, X., Zhang, M., Zhao, W., Liu, Y., Kardia, S. L. R., Roux, A. V. D., Needham, B. L., Smith, J. A., and Mukherjee, B. (2020).  Bayesian shrinkage estimation of high dimensional causal mediation effects in omics studies.  Biometrics, 76(3):700–710.
  * Sørlie et al., (2001) Sørlie, T., Perou, C. M., Tibshirani, R., Aas, T., Geisler, S., Johnsen, H., Hastie, T., Eisen, M. B., van de Rijn, M., Jeffrey, S. S., Thorsen, T., Quist, H., Matese, J. C., Brown, P. O., Botstein, D., Lønning, P. E., and Børresen-Dale, A.-L. (2001).  Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications.  Proceedings of the National Academy of Sciences, 98(19):10869–10874.
  * Sun and Xia, (2023) Sun, J. and Xia, Y. (2023).  Pretreating and normalizing metabolomics data for statistical analysis.  Genes & Diseases.
  * Tang et al., (2017) Tang, Z., Shen, Y., Zhang, X., and Yi, N. (2017).  The spike-and-slab lasso Cox model for survival prediction and associated genes detection.  Bioinformatics, 33(18):2799–2807.
  * Tapak et al., (2015) Tapak, L., Saidijam, M., Sadeghifar, M., Poorolajal, J., and Mahjub, H. (2015).  Competing risks data analysis with high-dimensional covariates: An application in bladder cancer.  Genomics, Proteomics & Bioinformatics, 13(3):169–176.
  * Tarazona et al., (2021) Tarazona, S., A., A.-L., and Conesa, A. (2021).  Undisclosed, unmet and neglected challenges in multi-omics studies.  Nature Computational Science, 1:395–402.
  * Taylor et al., (2022) Taylor, S., Ponzini, M., Wilson, M., and Kim, K. (2022).  Comparison of imputation and imputation-free methods for statistical analysis of mass spectrometry data with missing data.  Briefings in Bioinformatics, 23(1):bbab353.
  * Tibshirani, (1997) Tibshirani, R. (1997).  The lasso method for variable selection in the Cox model.  Statistics in Medicine, 16(4):385–395.
  * Tsiatis et al., (2019) Tsiatis, A. A., Davidian, M., Holloway, S. T., and Laber, E. B. (2019).  Dynamic Treatment Regimes: Statistical Methods for Precision Medicine.  Chapman and Hall/CRC, Boca Raton, FL.
  * Tuglus and van der Laan, (2011) Tuglus, C. and van der Laan, M. J. (2011).  Targeted methods for biomarker discovery.  In Targeted Learning: Causal Inference for Observational and Experimental Data, pages 367–382. Springer, New York.
  * Tuikkala et al., (2005) Tuikkala, J., Elo, L., Nevalainen, O. S., and Aittokallio, T. (2005).  Improving missing value estimation in microarray data with gene ontology.  Bioinformatics, 22(5):566–572.
  * Uno et al., (2011) Uno, H., Cai, T., Pencina, M. J., D’Agostino, R. B., and Wei, L. J. (2011).  On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data.  Statistics in Medicine, 30(10):1105–1117.
  * van der Laan and Rose, (2011) van der Laan, M. J. and Rose, S. (2011).  Targeted Learning: Causal Inference for Observational and Experimental Data.  Springer, New York.
  * van der Laan and Rose, (2018) van der Laan, M. J. and Rose, S. (2018).  Targeted Learning in Data Science: Causal Inference for Complex Longitudinal Studies.  Springer, New York.
  * van der Maaten, (2014) van der Maaten, L. (2014).  Accelerating t-SNE using tree-based algorithms.  Journal of Machine Learning Research, 15(93):3221–3245.
  * van der Maaten and Hinton, (2008) van der Maaten, L. and Hinton, G. (2008).  Visualizing data using t-SNE.  Journal of Machine Learning Research, 9(86):2579–2605.
  * van Houwelingen, (2000) van Houwelingen, H. C. (2000).  Validation, calibration, revision and combination of prognostic survival models.  Statistics in Medicine, 19(24):3401–3415.
  * VanderWeele, (2011) VanderWeele, T. J. (2011).  Causal mediation analysis with survival data.  Epidemiology, 22(4):582–585.
  * Vinga, (2021) Vinga, S. (2021).  Structured sparsity regularization for analyzing high-dimensional omics data.  Briefings in Bioinformatics, 22(1):77–87.
  * Volinsky et al., (1997) Volinsky, C. T., Madigan, D., Raftery, A. E., and Kronmal, R. A. (1997).  Bayesian model averaging in proportional hazard models: Assessing the risk of a stroke.  Journal of the Royal Statistical Society Series C: Applied Statistics, 46(4):433–448.
  * Wang et al., (2019) Wang, C., Li, Q., Song, X., and Dong, X. (2019).  Bayesian adaptive lasso for additive hazard regression with current status data.  Statistics in Medicine, 38(20):3703–3718.
  * Wang and Doudna, (2023) Wang, J. Y. and Doudna, J. A. (2023).  CRISPR technology: A decade of genome editing is only the beginning.  Science, 379(6629):eadd8643.
  * Wang et al., (2013) Wang, W., Baladandayuthapani, V., Holmes, C. C., and Do, K.-A. (2013).  Integrative network-based Bayesian analysis of diverse genomics data.  BMC Bioinformatics, 14(13):S8.
  * Wang and Liu, (2020) Wang, W. and Liu, W. (2020).  Integration of gene interaction information into a reweighted Lasso-Cox model for accurate survival prediction.  Bioinformatics, 36(22-23):5405–5414.
  * Wiegrebe et al., (2023) Wiegrebe, S., Kopper, P., Sonabend, R., and Bender, A. (2023).  Deep learning for survival analysis: A review.  arXiv:2305.14961.
  * Zellner, (1986) Zellner, A. (1986).  On assessing prior distributions and Bayesian regression analysis with $g$-prior distributions.  In Goel, R. K. and Zellner, A., editors, Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti. Studies in Bayesian Econometrics and Statistics, Vol. 6, chapter 15, pages 223–243. North-Holland.
  * Zhang and Lu, (2007) Zhang, H. H. and Lu, W. (2007).  Adaptive Lasso for Cox’s proportional hazards model.  Biometrika, 94(3):691–703.
  * Zhang et al., (2020) Zhang, J., Guan, M., Wang, Q., Zhang, J., Zhou, T., and Sun, X. (2020).  Single-cell transcriptome-based multilayer network biomarker for predicting prognosis and therapeutic response of gliomas.  Briefings in Bioinformatics, 21(3):1080–1097.
  * Zhang et al., (2017) Zhang, Y., Kwok-Shing Ng, P., Kucherlapati, M., Chen, F., Liu, Y., Tsang, Y. H., de Velasco, G., Jeong, K. J., Akbani, R., Hadjipanayis, A., Pantazi, A., Bristow, C. A., Lee, E., Mahadeshwar, H. S., Tang, J., Zhang, J., Yang, L., Seth, S., Lee, S., Ren, X., Song, X., Sun, H., Seidman, J., Luquette, L. J., Xi, R., Chin, L., Protopopov, A., Westbrook, T. F., Shelley, C. S., Choueiri, T. K., Ittmann, M., Van Waes, C., Weinstein, J. N., Liang, H., Henske, E. P., Godwin, A. K., Park, P. J., Kucherlapati, R., Scott, K. L., Mills, G. B., Kwiatkowski, D. J., and Creighton, C. J. (2017).  A pan-cancer proteogenomic atlas of PI3K/AKT/mTOR pathway alterations.  Cancer Cell, 31(6):820–832.e3.
  * Zhao et al., (2023) Zhao, Z., Banterle, M., Lewin, A., and Zucknick, M. (2023).  Multivariate bayesian structured variable selection for pharmacogenomic studies.  arXiv:2101.05899.
  * Zhao and Zucknick, (2020) Zhao, Z. and Zucknick, M. (2020).  Structured penalized regression for drug sensitivity prediction.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 69(3):525–545.
  * Zhou et al., (2019) Zhou, R., Zeng, D., Zhang, J., Sun, H., Wu, J., Li, N., Liang, L., Shi, M., Bin, J., Liao, Y., Huang, N., and Liao, W. (2019).  A robust panel based on tumour microenvironment genes for prognostic prediction and tailoring therapies in stage I-III colon cancer.  EBioMedicine, 42:420–430.
  * Zucknick et al., (2008) Zucknick, M., Richardson, S., and Stronach, E. A. (2008).  Comparing the characteristics of gene expression profiles derived by univariate and multivariate classification methods.  Statistical Applications in Genetics and Molecular Biology, 7(1):Article 7.
  * Zucknick et al., (2015) Zucknick, M., Saadati, M., and Benner, A. (2015).  Nonidentical twins: Comparison of frequentist and Bayesian lasso for Cox models.  Biometrical Journal, 57(6):959–981.