§ INTRODUCTION

In “parametric” survival analysis, researchers commonly rely on the Weibull
distribution or alternative acceleration failure time models.  To achieve
greater flexibility and overcome the strict distributional assumptions
inherent in these models, researchers often need to turn to
non-/semi-parametric methods to analyse their survival data.  When it comes
to non-/semi-parametric approaches, however, overcoming issues tied to
interval-censored or truncated observations can prove challenging due
to their limited availability in standard software implementations.

Furthermore, when aiming to fit models that handle crossing or non-proportional
hazards, clustered observations, or dependent censoring, researchers often
find themselves navigating a complex landscape of diverse software
implementations.  Even the same models are sometimes hard to compare across
different implementations, because different parametrisations,
estimation strategies or optimisation procedures are used.
This becomes even more challenging when different
models computed from different packages shall be compared – underscoring
the usefulness of a unified framework
that facilitates smooth transitions between different models relevant for
the analysis of survival data and that is based on a common core infrastructure
for model parametrisation, estimation strategy, and optimisation procedure.

To tackle these challenges, researchers might consider exploring the
potential of the tram add-on package <cit.> in R
<cit.>, which offers a flexible framework for survival analysis.  The
tram package implements a user-friendly interface to smooth
transformation models <cit.>, which encompass
a range of classical survival models including acceleration failure models
and the Cox proportional hazards model, as well as many useful extensions
and novel model formulations.  The package can be easily installed from the
Comprehensive R Archive Network (CRAN):

R> install.packages("tram")
R> library("tram")

All models are implemented in a fully parametric fashion, allowing for
straightforward maximum likelihood inference.  The estimation of these
models is facilitated by the core infrastructure package mlt
<cit.> which provides a unified maximum-likelihood
framework for general transformation models
<cit.>.  We further leverage two add-on
packages from this family of packages for transformation modelling: The
tramME package <cit.> implementing mixed-effects and
non-linear additive effects for smooth transformation models, and the
trtf package <cit.> designed for estimating tree-based
survival models.

In this tutorial, we will explore a variety of models commonly utilised in
survival analysis.  The focus of this tutorial lies on the practical
implementation and interpretation of these models within the framework of
smooth transformation models, rather than on theoretical aspects.  Our
objective is to provide users with a practical understanding of how to apply
these models using R.  Through an application to data from a
randomised trial on rectal cancer therapy, we showcase how users can
seamlessly conduct their survival analysis tasks without the need to
navigate through a multitude of packages in R.

In Section <ref>, we discuss models for independent observations.
We start with the well-known Weibull model, and then, to introduce more
flexibility and overcome the log-linear log-cumulative hazard assumption
inherent to the Weibull model, we explore a fully parametric version of the
Cox model.  We further discuss the estimation of stratified log-cumulative
hazard functions to account for baseline risk variations across patient
strata.  Moving beyond the assumption of proportional hazards, we showcase
models that challenge this assumption.  We discuss a location-scale version
of the Cox model, accommodating scenarios with non-proportional/crossing
hazards and models estimating time-varying treatment effects.

Addressing scenarios where the assumption of independent censoring might not
be justified, we discuss a copula proportional hazards model, that
accommodates dependent censoring (Section <ref>).  For clustered
observations we employ mixed-effects Cox models and alternative models
featuring marginally interpretable hazard ratios in Section <ref>.
Our tutorial also explores the domain of personalised medicine, presenting
models that incorporate covariate-dependent treatment effects and survival
trees (Section <ref>).  In Section <ref>, we explore further
extensions, including topics like frailty models, model estimation using the
non-parametric likelihood, and the potential of using these models for
sample size estimation of new trials.

This tutorial consists of the main text introducing models and, very
briefly, demonstrating how to estimate those using smooth transformation
models in R. In addition,  we present numerical head-to-head
comparisons of user-interfaces and results obtained from alternative
packages available in the R universe in Appendix <ref>.
Both parts come with
much more detailed R code readers can consult for ways how to
interrogate fitted models (for example, how to plot model terms or how to
compute confidence intervals and to perform tests), which can be run with:

R> demo("survtram", package = "tram")

In our Appendix <ref>, we conduct a thorough comparison of a
subset of the models discussed here which can be estimated using alternative
implementations (in total 14 established CRAN packages) and corresponding results obtained with smooth transformation models from
tram and tramME.  This quality assurance task not only helped to validate the
unified implementation in tram but also lead to the identification of
problematic special cases and, in some instances, practically relevant
discrepancies between different implementations of the very same model.
Moreover, Appendix <ref> presents the different user interfaces
offered by various packages side-by-side, such that it becomes simple to
estimate and compare relatively complex models across independent
implementations.  For the analysis of future survival trials, an assessment
of the agreement of such estimates, standard errors, and possibly other
model aspects can help to increase faith in reported numerical results and
conclusions derived therefrom.

§ APPLICATION

In our tutorial, we will work with data from the CAO/ARO/AIO-04 two-arm
randomised controlled trial
<cit.>, a phase 3
study that aimed to compare an experimental regimen with the previously
established treatment regimen (control) for locally advanced rectal cancer.  In this
experimental regimen, Oxaliplatin was added to the control treatment of preoperative
Fluorouracil-based chemoradiotherapy and postoperative chemotherapy of
locally advanced rectal cancer.

The trial was conducted across
88 centers and included a cohort of
1'236 patients.  The patients were randomly allocated
to the two treatment arms $\rW \in \{0, 1\}$, receiving the experimental
treatment of Oxaliplatin added to Fluorouracil-based preoperative
chemoradiotherapy and postoperative chemotherapy (5-FU + Oxaliplatin, $\rW =
1$) or the control treatment using Fluorouracil only (5-FU, $\rW
= 0$).  Treatment allocation was performed using block-randomisation
stratified by study center $j = 1, \dots,
88$ and the stratum $\rs$, which is
defined by four categories consisting of a combination of clinical
N category,  lymph node involvement
(cN0 vs cN+), and clinical
T category  tumor grading (cT1-3 vs cT4).  The distribution of patients in the two treatment arms across
strata is shown in Table <ref>.

5-FU     5-FU + Oxaliplatin

cT1-34 : -cN0     163     156

cT41-3 : -cN0     8     7

cT1-34 : -cN+     411     417

cT41-3 : -cN+     41     33

Total     623     613

Number of patients in each treatment arm stratified by the combination of
clinical N and T category.

The primary endpoint is disease-free survival, defined
as the time $\rY \in \RR^{+}$ between randomisation and non-radical surgery
of the primary tumor (R2 resection), loco-regional recurrence after R0/1
resection, metastatic disease or progression, or death from any cause –
whichever occurred first.  The observed times encompass a mix of exact
observations $\ry$ for time to death or incomplete removal of the primary
tumor, interval-censored observations $\ry \in (\ubar{\ry}, \bar{\ry}]$ for the
time span from the previous follow-up $\ubar{\ry}$ to the time-point of detecting local or distant
metastases $\bar{\ry}$, and right-censored observations $\ry \in (\ry,
\infty)$ corresponding to the end of the follow-up period or instances where
patients were lost to follow-up.  The survivor curves of the primary
endpoint (disease-free survival) estimated by the non-parametric Turnbull
estimator <cit.> are shown for the two treatment arms in Figure <ref>.

Disease-free survival. The survivor functions of the two treatment arms estimated
by the non-parametric Turnbull method are shown together with the number at risk table.

A secondary endpoint considered in the study is overall survival, defined as
the time $\rY \in \RR^{+}$ from the point of randomisation to death
from any cause.  Notably, all observations $\ry$ for this endpoint are
exact or right-censored.  The corresponding survivor curves, estimated
non-parametrically by the Kaplan-Meier method <cit.>, are shown in Figure <ref>.

Overall survival. The survivor functions of the two treatment arms estimated
by the non-parametric Kaplan-Meier method are shown together with the number at risk table.

The primary data analysis for this trial was performed by
<cit.>.  In their analysis, the
treatment effect comparing the effect of the experimental treatment
to the effect of the control treatment on disease-free and overall survival was
assessed by adjusted log-rank tests and mixed-effects Cox models (both adjusting for the stratified randomisation process), treating both survival endpoints as exact.  We demonstrate
the process of fitting a fully parametric mixed-effects Cox model that
accounts for interval-censored event times of the primary endpoint in Section <ref>, also in terms
of a model featuring marginal interpretation of the estimated hazard ratio.

A subsequent post hoc analysis was carried out by
<cit.>.  In this analysis, the
possibility of an age-varying treatment effect on both the primary endpoint
of disease-free survival and the secondary endpoint of overall survival was
investigated.  In
Section <ref>, we demonstrate how such an analysis can be performed
within the discussed framework, taking into account interval-censoring.
We illustrate two approaches for estimating age-varying effects, using smoothly
age-varying hazard ratios and age-structured survival trees.

While the conducted analyses made serious efforts to address the research
questions effectively, they encountered limitations due to the lack of
software resources capable of appropriately handling interval-censored and
correlated observations for the analysis of the primary endpoint.  Notably,
the first R add-on package capable of estimating Cox models in
the presence of interval censoring was published in 2014
(). At the time of the
statistical analysis of the primary endpoint, it was impossible to fit
mixed-effects models with flexible baseline hazard to interval-censored outcomes.
This obstacle was one of the main
motivation to develop a comprehensive software package implementing a
general class of transformation models with applications in the domain of survival
analysis. The corresponding framework implementing smooth
transformation models <cit.> in R
helps to address such and other practically relevant limitations.
In this tutorial, we present analyses that the CAO/ARO/AIO-04 study investigators
would have liked to have been able to perform a decade ago.

§ INDEPENDENT OBSERVATIONS

§.§ Survival models

§.§.§ Weibull proportional hazards model

Probably one of the most renowned parametric models in survival analysis
is the Weibull model
<cit.>, where the response $\rY$ conditional on treatment
assignment $\rW = \rw$ is assumed to follow a Weibull
distribution.  We consider the Weibull proportional hazards model with
survivor functions conditional on the treatment assignment parametrised as,

\begin{eqnarray} \label{mod:WEI}
\snul(\ry) &=& \Prob(\rY > \ry \mid \rW = 0) =
\sMEV{\eparm_1 + \eparm_2 \log(\ry)},
\qquad \qquad \eparm_2 > 0,\nonumber\\
\sone(\ry) &=& \Prob(\rY > \ry \mid \rW = 1) =
\sMEV{\eparm_1 + \eparm_2 \log(\ry) - \eshiftparm},\nonumber\\
\sw(\ry) &=& \Prob(\rY > \ry \mid \rW = \rw) =
\sMEV{\eparm_1 + \eparm_2 \log(\ry) - \eshiftparm \rw}.
\end{eqnarray}

The log-cumulative baseline hazard $\log(-\log(\snul(\ry))) =
\log(\Haznul(\ry))$ here is given by $\h(\ry) = \eparm_1 + \eparm_2
\log(\ry)$, assuming a linear shift $\eshiftparm$ on the scale of log-time $\log(\ry)$.  The model not only assumes proportional
hazards, with hazard ratio $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)} =
\nicefrac{\hazone(\ry)}{\haznul(\ry)} = \exp(-\eshiftparm)$, but is also an
accelerated failure time model

\begin{eqnarray*} \label{mod:WEIAFT}
\log(\rY) &=& \frac{- \eparm_1 + \tilde\eshiftparm \rw + \rZ}{\eparm_2}, \qquad \rZ \sim \MEV\\
&=& - \frac{\eparm_1}{\eparm_2} + \frac{\tilde\eshiftparm}{\eparm_2} \rw + \frac{1}{\eparm_2}\rZ,
\end{eqnarray*}

with the errors $\rZ$ following a minimum extreme value distribution (MEV).
Consequently, $\rY \mid \rW = \rw$ follows a Weibull distribution
(, Chapter 2) with shape parameter
$\eparm_1 \eparm_2^{-1}$, scale parameter $\eparm_2$ and log-acceleration
factor $\tilde\eshiftparm = \eshiftparm \eparm_2^{-1}$.  This implies that,
for the experimental arm the time $\rY$ is accelerated by
$\exp(\tilde\eshiftparm)$, that is $\rY_1 = \exp(\tilde\eshiftparm) \rY_0$,
thus the probability of survival of the experimental arm is given by
$\sone(\ry) = \snul(\exp(-\tilde\eshiftparm) \ry)$.

Alternatively, different distributions for $\rZ$, such as the normal or
logistic distribution, can be specified, leading to the formulation of
log-normal or log-logistic models.

Parameter estimation of the Weibull model is straightforward using maximum
likelihood, because the distribution function can be directly evaluated and
thus allows to effectively handle interval-censored or truncated data, as
will be discussed in Section <ref>.

§.§.§ Flexible proportional hazards model

The assumption of a log-linear log-cumulative baseline hazard function $\h$, implied by the
Weibull model, can be relaxed by replacing
$\log(\Haznul(\ry)) = \h(\ry) = \eparm_1 + \eparm_2
\log(\ry)$ with a more flexible function
$\h(\ry) = \basisy(\ry)^\top \parm$ defined in terms of smooth spline basis
functions $\basisy$ and corresponding parameters $\parm$.  This yields the following model

\begin{eqnarray} \label{mod:COX}
\snul(\ry) &=& \Prob(\rY > \ry \mid \rW = 0) = \sMEV{\h(\ry)},\nonumber\\
\sone(\ry) &=& \Prob(\rY > \ry \mid \rW = 1) = \sMEV{\h(\ry) + \eshiftparm},\nonumber\\
\sw(\ry) &=& \Prob(\rY > \ry \mid \rW = \rw) = \sMEV{\h(\ry) + \eshiftparm \rw},
\end{eqnarray}

where the hazard ratio is given by $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)} =
\nicefrac{\hazone(\ry)}{\haznul(\ry)} = \exp(\eshiftparm)$.
<cit.> and <cit.> proposed
parametrising $\h(\ry) = \basisy(\ry)^\top \parm$ as polynomials in
Bernstein form $\h(\ry) = \bern{\dimparm - 1}(\ry)^\top \parm$ of order
$\dimparm - 1$.  In the context of skewed survival data, it might be
suitable to apply Bernstein polynomials to the logarithm of the survival
time $\ry$, leading to $\h(\ry) = \bern{\dimparm - 1}(\log(\ry))^\top
\parm$.  This model is a fully parametric version of the otherwise semi-parametric Cox
proportional hazards model <cit.>.  The latter treats $\h$ as an
infinite dimensional nuisance parameter which is profiled out from the
likelihood.  This gives rise to the partial likelihood allowing to estimate
the log-hazard ratio $\eshiftparm$ <cit.>.  In contrast, all
parameters of model (<ref>) are estimated simultaneously by maximum
likelihood.  The parametrisation of $\h$ in terms of basis functions
and corresponding parameters allows to specify of a flexible, yet fully
parametric, monotonically increasing log-cumulative hazard function.
This is achieved under appropriate constraints $\eparm_p \le \eparm_{p + 1}$
for $p \in 1, \dots, \dimparm - 1$ <cit.>.
Adopting this specific parametrisation for the log-cumulative baseline
hazard function $\log(\Haznul(\ry)) = \h(\ry)$ facilitates the computation
of the associated density $\dnul(\ry)$ and distribution function
$\pnul(\ry)$, thus allowing for straightforward parameter estimation using
maximum likelihood.  This holds true even when dealing with
interval-censored or truncated observations.

§.§.§ Stratified proportional hazards model

Accounting for the variations in baseline risks among different patient
strata identified by variable $\rs$, one can employ stratified models that
incorporate stratum-specific log-cumulative baseline hazard functions
$\h(\ry \mid \rs)$.  These models can be defined by

\begin{eqnarray*} \label{mod:STRAT}
\snul(\ry \mid \rs) &=& \Prob(\rY > \ry \mid \rS = \rs, \rW = 0) =
\sMEV{\h(\ry \mid \rs)},\nonumber\\
\sone(\ry \mid \rs) &=& \Prob(\rY > \ry \mid \rS = \rs, \rW = 1) =
\sMEV{\h(\ry \mid \rs) + \eshiftparm},\nonumber\\
\sw(\ry \mid \rs) &=& \Prob(\rY > \ry \mid \rS = \rs, \rW = \rw) =
\sMEV{\h(\ry \mid \rs) + \eshiftparm \rw},
\end{eqnarray*}

with $\h(\ry \mid \rs) = \basisy(\ry)^\top \parm(\rs)$ and global hazard
ratio $\nicefrac{\Hazone(\ry \mid \rs)}{\Haznul(\ry \mid \rs)} =
\nicefrac{\hazone(\ry \mid \rs)}{\haznul(\ry \mid \rs)} = \exp(\eshiftparm)$,
assuming that the treatment effect is equivalent
across all patient strata $\rs$.

One could further relax the restriction of a global treatment effect, allowing
for an interaction of the treatment assignment $\rW = \rw$ and the stratum
$\rs$ by formulating the log-cumulative hazard as $\h(\ry \mid \rs) + \rw
\eshiftparm + \gammavec^\top(\rw \times \rs)$.

§.§.§ Non-proportional hazards model

Extending beyond the proportional hazards assumption, an additional
covariate-dependent model term can be estimated.  <cit.>
introduced the multi-parameter extension to the Weibull model (<ref>) in the context
of survival analysis, specifically outlining its use for interval-censored observations <cit.>.

A similar extension can be made to the more flexible, parametric, Cox
model (<ref>),
by additionally estimating a scale term $\escaleparm$ for the experimental arm

\begin{eqnarray*} \label{mod:SCOX}
\snul(\ry) &=& \Prob(\rY > \ry \mid \rW = 0) =
\sMEV{\h(\ry)},\nonumber\\
\sone(\ry) &=& \Prob(\rY > \ry \mid \rW = 1) =
\sMEV{\sqrt{\exp(\escaleparm)}\,\h(\ry) + \eshiftparm},\nonumber\\
\sw(\ry) &=& \Prob(\rY > \ry \mid \rW = \rw) =
\sMEV{\sqrt{\exp(\escaleparm\rw)}\,\h(\ry) + \eshiftparm \rw}.
\end{eqnarray*}

In this case the ratio of the cumulative hazards, $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)}$, is a
non-proportional function of $\ry$.
The corresponding bivariate score test
(Section 3.1.2 of ) further allows to test the null
hypothesis of equal survival,  $\eshiftparm = \escaleparm = 0$,
without relying on the assumption of proportional hazards.

§.§.§ Time-varying hazards model

Accounting for changing effects of the treatment over time, we can further
extend beyond the proportional hazards assumptions and estimate a model
incorporating a time-varying treatment effect,

\begin{eqnarray*} \label{mod:TCOX}
\snul(\ry) &=& \Prob(\rY > \ry \mid \rW = 0) =
\sMEV{\h(\ry)},\\
\sone(\ry) &=& \Prob(\rY > \ry \mid \rW = 1) =
\sMEV{\h(\ry) + \eshiftparm(\ry)},\\
\sw(\ry) &=& \Prob(\rY > \ry \mid \rW = \rw) =
\sMEV{\h(\ry) + \eshiftparm(\ry)\rw}.
\end{eqnarray*}

Here, the model introduces a time-varying shift $\eshiftparm(\ry)$ in the log-cumulative hazard
function $\log(\Hazone(\ry)) = \log(\Haznul(\ry)) + \eshiftparm(\ry)$, avoiding
the assumption of a constant log-hazard ratio $\eshiftparm$. The shift
$\eshiftparm(\ry)$ is also parameterised in terms of a
polynomial in Bernstein form, thus allowing to estimate a time-varying ratio
of the cumulative hazards $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)} = \exp(\eshiftparm(\ry))$.

§.§ Likelihood

From the above models, the log-likelihoods for exact or independently right-,
left- or interval-censored and truncated observations are easily deducible.
We here consider the most general case where the log-cumulative hazard
function is given by $\h(\ry \mid \rw, \rs) = \sqrt{\exp(\escaleparm \rw)}\,\h(\ry
\mid \rs) + \eshiftparm \rw$.

For exact continuous observations $(\ry, \rw, \rs)$, the corresponding likelihood
contributions are given by

\begin{eqnarray*} \label{eq:cll}
\ell(\parm(\rs), \eshiftparm, \escaleparm \mid \rY = \ry) =
\log\left\{\dZ\left[\h(\ry \mid \rw, \rs) \right]\right\} +
\log\left\{\h^\prime(\ry \mid \rw, \rs)\right\},
\end{eqnarray*}

with the standardised minimum extreme value density $\dZ(\rz) = \exp(\rz -
\exp(\rz))$ and the derivative of the log-cumulative hazard function with
respect to $\ry$, $\h^\prime(\ry \mid \rw, \rs) = \sqrt{\exp(
\escaleparm \rw)}\, \basisy^\prime(\ry)^\top \parm(\rs)$.

Because the transformation function $\h$, defining the log-cumulative
baseline hazard function, is parametrised in terms of polynomials in
Bernstein form, where the basis functions $\bern{\dimparm - 1}(\ry) \in
\RR^\dimparm$ are specific beta densities <cit.>, it is
straightforward to obtain the derivatives of the basis functions with
respect to $\ry$, leading to $\h^\prime(\ry \mid \rs) =
\bern{\dimparm - 1}^\prime(\ry)^\top \parm(\rs)$.

Under independent left-, right- or interval-censored event times
$(\ubar{\ry}, \bar{\ry}]$ the exact log-likelihood contribution is

\begin{eqnarray*} \label{eq:ill}
\ell(\parm(\rs), \eshiftparm, \escaleparm \mid \rY \in (\ubar{\ry}, \bar{\ry}])  =
\log\left\{\Prob(\rY \in (\ubar{\ry}, \bar{\ry}] \mid \rw, \rs)\right\} =
\log\left\{\sw(\ubar\ry \mid \rw, \rs) -
\sw(\bar\ry \mid \rw, \rs)\right\}.
\end{eqnarray*}

For observations that are right-censored at time $\ry$ the datum is given by
$(\ubar{\ry}, \bar{\ry}] = (\ry, \infty)$ and for left-censored observations
it is $(\ubar{\ry}, \bar{\ry}] = (0, \ry]$.

In cases where event times are subject to random left-, right-, or
interval-truncation $(\ry_{\text{l}}, \ry_{\text{r}}] \subset \RR^{+}$, the
above log-likelihood contributions change to

\begin{eqnarray*} \label{eq:tll}
\ell(\parm(\rs), \eshiftparm, \escaleparm \mid \rY \in (\ubar{\ry}, \bar{\ry}]) -
\ell(\parm(\rs), \eshiftparm, \escaleparm \mid \rY \in (\ry_\text{l}, \ry_\text{r}])
\end{eqnarray*}

with $\ry_\text{l} = 0$ for right-truncated and $\ry_\text{r} = \infty$ for
left-truncated observations.  Such considerations are relevant in scenarios
involving a late entry approach, for instance, resulting in left-truncated
observations, where one is interested in modelling $\Prob(\ry > \rY \mid \rY \in
(\ry_{\text{l}}, \infty))$, or for modelling time-varying covariates.

§.§ Application

Now, turning our attention to the CAO/ARO/AIO-04 two-arm
randomised controlled trial, we aim to estimate the previously discussed
models using the unified maximum likelihood framework provided by the
R add-on package tram <cit.>.  We fit the models
to the primary endpoint of disease-free survival $\rY$ (stored in the object
iDFS of class `Surv')
estimating the treatment effect corresponding to the assignment $\rW =
\var{randarm}$. It should be noted that iDFS contains a mix of exact,
interval-, and right-censored observations.

We start by fitting the Weibull model (Section <ref>) using the Survreg

R> Survreg(iDFS   randarm, data = CAOsurv, dist = "weibull")

Coefficient     Estimate     Std. Error     95%-Wald CI

$\eparm_1$     $ -6.231 $     0.265     $-6.752$ to $-5.711$

$\eparm_2$     $ 0.733 $     0.036     $\phantom{-}0.663$ to $\phantom{-}0.803$

$\eshiftparm$     $ 0.229 $     0.106     $\phantom{-}0.020$ to $\phantom{-}0.438$

Log-Likelihood     Likelihood Ratio Test     Score Test     Permutation Score Test

$-$2'281.17     $p$ = 0.031     $p$ = 0.031     $p$ = 0.035

The model quantifies the treatment effect through a hazard ratio
$\exp(-\hat\eshiftparm) = 0.795$, comparing the hazard of the
treatment arm to the hazard of the control arm.  The results indicate
a reduction in hazard for patients receiving the
experimental treatment compared to the control treatment.
This reduction in hazard translates to a
prolonged disease-free survival time in the experimental arm.
Since the model is a proportional hazards counterpart of the Weibull accelerated
failure time model fitted by survreg from the survival package
the estimate can also be translated into a log-acceleration factor
$\hat{\tilde\eshiftparm} = \hat\eshiftparm \hat\eparm_2^{-1} =
0.312$.  This implies that
the disease-free survival time $\rY$ is prolonged by
$\exp(\tilde\eshiftparm) = 1.367$ in the
experimental arm, compared to the control arm.

Next, we fit the flexible proportional hazards model (Section <ref>) using the
Coxph function from the tram package.

R> Coxph(iDFS   randarm, data = CAOsurv, log_first = TRUE)

Coefficient     Estimate     Std. Error     95%-Score CI

$\eshiftparm$     $ -0.234 $     0.106     $-0.442$ to $-0.025$

Log-Likelihood     Likelihood Ratio Test     Score Test     Permutation Score Test

$-$2'255.85     $p$ = 0.028     $p$ = 0.028     $p$ = 0.028

The fitted model is a fully parametric version of the famous Cox model,
otherwise estimated semi-parametrically using the partial likelihood (as
implemented in the survival package in the coxph function).
Here the log-cumulative hazard function is specified in terms of
polynomials in Bernstein form, by default of order $\dimparm - 1 = 6$, which
is defined on the scale of $\log(\ry)$ (setting log_first = TRUE),
thus specifying the transformation function $\h(\ry)$ in terms of
$\bern{6}(\log(\ry))^\top\parm$.  The fully parametric approach enables the
straightforward incorporation of interval-censored disease-free survival times.
Figure <ref> illustrates the estimated log-cumulative baseline hazard
function $\Haznul(\ry)$ along with the 95%-confidence band,
revealing a non-linear function of log-time ($\log(\ry)$).
The estimated hazard ratio is $\exp(\hat\eshiftparm) = 0.792$,
indicating reduced hazards in the experimental arm.
The software implementation further allows
to compute a corresponding 95%-permutation score confidence interval which
ranges from $0.642$ to $0.976$.

Flexible proportional hazards model.
Log-cumulative baseline hazard function and corresponding 95%-confidence
band estimated by the model.

To further accommodate for varying log-cumulative baseline hazard functions $\Haznul(\ry \mid \rs)$
across patient strata $\rs$ (here identified by $\code{strat}$), we can fit
a stratified model (Section <ref>).

R> Coxph(iDFS | strat   randarm, data = CAOsurv, log_first = TRUE)

Coefficient     Estimate     Std. Error     95%-Score CI

$\eshiftparm$     $ -0.232 $     0.107     $-0.440$ to $-0.024$

Log-Likelihood     Likelihood Ratio Test     Score Test     Permutation Score Test

$-$2'231.85     $p$ = 0.029     $p$ = 0.029     $p$ = 0.030

The model estimates separate smooth log-cumulative baseline hazards for each
stratum $\rs$, as illustrated in Figure <ref>, but provides an estimate
of the global hazard ratio $\exp(\hat\eshiftparm) = 0.793$,
indicating a reduction of the hazard in the experimental arm by
0.793 relative to the hazard in the control arm across all

Stratified proportional hazards model.
Log-cumulative baseline hazard functions $\Haznul(\ry \mid \rs)$ estimated
by the model are shown separately for each stratum $\rs$.

Moving away from the proportional hazards assumption, we can fit a non-proportional hazards model (a location-scale version of the Cox model, Section <ref>) using the same function.

R> Coxph(iDFS   randarm | randarm, data = CAOsurv, log_first = TRUE)

Coefficient     Estimate     Std. Error     95%-Wald CI

$\eshiftparm$     $ -0.076 $     0.164     $-0.398$ to $0.247$

$\escaleparm$     $ 0.289 $     0.202     $-0.106$ to $0.685$

Log-Likelihood     Likelihood Ratio Test     Biv. Permutation Score Test

$-$2'254.83     $p$ = 0.032     $p$ = 0.014

The ratio of the cumulative hazards $\nicefrac{\widehat\Hazone(\ry)}{\widehat\Haznul(\ry)}$,
shown in Figure <ref>, no longer remains proportional but varies over time.
The curve indicates a pronounced initial reduction in cumulative hazards for the
experimental arm compared to the control arm, which
gradually decreases as time progresses.  This suggests that the treatment
effect is stronger in the beginning.  The corresponding bivariate score
test, which tests the null hypothesis of equal survival without assuming
proportional hazards, further indicates evidence for non-equal disease-free
survival times.

Non-proportional hazards model. $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)}$ (solid line) estimated by the model is shown alongside the constant hazard
ratio estimated from the proportional hazards model (dashed line) over time

Finally, we fit the model featuring a time-varying treatment effect (Section <ref>).

R> Coxph(iDFS | randarm   1, data = CAOsurv, log_first = TRUE)

The treatment effect $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)} =
\exp(\eshiftparm(\ry))$ is a function of time, as shown in
Figure <ref>.  The curve again demonstrates a reduction in hazards
for the experimental arm compared to the control arm, which is more
substantial in the beginning and gradually becomes less prominent as time

Time-varying hazards model. $\nicefrac{\Hazone(\ry)}{\Haznul(\ry)}$
and corresponding 95%-confidence
bands over time
$\ry$ (solid line) estimated by the model is shown alongside the constant hazard
ratio estimated from the proportional hazards model (dashed line). The log-likelihood of the model is

§ DEPENDENT CENSORING

§.§ Copula proportional hazards model

Until now, the models, we have discussed, have been constructed under the
assumption of independent censoring, which implies that the censoring times
$\rC$ are independent of the event times $\rY$ given the treatment
assignment $\rW = \rw$, that is $\rY \indep \rC \mid \rW = \rw$.  We can
however move beyond relying on this assumption and allow the model to
capture potential dependence between the censoring times $\rC$ and event times $\rY$.

We explore models discussed in recent work of <cit.>,
which involve a joint model for $\rY$ and $\rC$ employing a bivariate
Gaussian copula of
the marginal survivor functions $\sYw(\ry)$ and
$\sCw(\rc)$ of $\rY$ and $\rC$ respectively,

\begin{eqnarray*} \label{mod:DEPC}
\Prob(\rY \leq \ry, \rC \leq \rc \mid \rW = \rw) = \Phi_{0, \mR(\lparm)}
\left\{
\Phi^{-1}\left[1 - \sYw(\ry)\right],
\Phi^{-1}\left[1 - \sCw(\rc)\right]
\right\}
\end{eqnarray*}

with correlation matrix

\begin{eqnarray*} \label{mod:DEPC:CORR}
\mR(\lparm) =
\left[\begin{array}{cc}
1 & \nicefrac{-\lparm}{\sqrt{1 + \lparm^2}}\\
\nicefrac{-\lparm}{\sqrt{1 + \lparm^2}} & 1
\end{array} \right], \qquad \lparm \in (-\infty, \infty)
\end{eqnarray*}

to account for the association between $\rY$ and $\rC$.
<cit.> recently demonstrated that the above model
maintains identifiability when the marginal survivor functions $\sYw(\ry)$ and
$\sCw(\rc)$ are described by a
flexible proportional hazards model (<ref>) and a model that assumes
a Weibull distribution (<ref>),
respectively.  This allows to estimate the dependence parameter $\lparm$ and
other model parameters from the observed data.
A dependence parameter $\lparm$ of zero corresponds to uncorrelated event
times $\rY$ and censoring times $\rC$, thus indicating lack of evidence for dependent censoring.

§.§ Application

Returning to our application, where we previously assumed independent
censoring of disease-free survival times,
we now aim to address the potential scenario where loss of
follow-up time $\rC \in \RR^{+}$ in the two arms is not independent of the
overall survival time  $\rY \in \RR^{+}$ (secondary endpoint).

The observed times can be categorised into the following event types (specified in DepCevent):
The event of interest (overall survival), loss of
follow-up (dependent censoring), and end of follow-up period
(administrative/independent censoring).

5-FU     5-FU + Oxaliplatin

Administrative censoring     469     466

Event of interest     106      96

Loss of follow-up      48      51

The model accommodating dependent censoring can also be fitted using the
Coxph function by appropriately specifying the event in the
`Surv' object.

R> Coxph(Surv(OStime, event = DepCevent)   randarm, data = CAOsurv)

Coefficient     Estimate     Std. Error     95%-Wald CI

$\eshiftparm_\rY$     $ -0.029 $     0.139     $-0.302$ to $0.243$

$\lparm$     $ -0.481 $     0.805     $-2.058$ to $1.097$

The joint model is fitted based on a Gaussian copula, estimating a
flexible proportional hazards marginal model (<ref>) for the time of overall survival $\rY$
and a Weibull proportional hazards marginal model (<ref>)
for the time of loss of follow-up $\rC$, while
accounting for independent right-censoring at the end of the follow-up

The estimated hazard ratio assessing the treatment effect on overall survival,
is $\exp(\hat\eshiftparm_\rY) = 0.971$ with a
95%-confidence interval from $0.739$ to $1.276$.  This indicates
no evidence for prolonged overall survival in the experimental arm.  The dependence
parameter $\lparm$ is estimated to be $ -0.481 $,
corresponding to a Kendall's $\tau = 0.285$.
The corresponding 95%-confidence interval from $-2.058$ to $1.097$
for $\lparm$ does include zero, thus providing no evidence for a dependence between time
of overall survival $\rY$ and loss of follow-up $\rC$ given the treatment assignment $\rW = \rw$.

§ DEPENDENT OBSERVATIONS

§.§ Survival models

Up to this point, the models we have discussed have been built upon the
assumption of independent observations.  However, this
assumption may not hold in situations where observations are clustered, such
as for multi-center trials where observations from the same center are

§.§.§ Mixed-effects proportional hazards model

In order to address this challenge, we can leverage a flexible mixed-effects
proportional hazards model as proposed by <cit.>.
This approach extends the previously discussed smooth transformation models by conditioning on an
unobserved cluster-specific random effect $R = r$ that accounts for the dependence within

\begin{eqnarray*} \label{mod:COXME}
\snul(\ry \mid R = r) &=&  \Prob(\rY > \ry \mid \rW = 0, R = r) =
\sMEV{\h(\ry) + r}, \qquad R \sim \ND(0, \tau^2),\\
\sone(\ry \mid R = r) &=& \Prob(\rY > \ry \mid \rW = 1, R = r) =
\sMEV{\h(\ry) + \eshiftparm + r},\\
\sw(\ry \mid R = r) &=&  \Prob(\rY > \ry \mid \rW = \rw, R = r) =
\sMEV{\h(\ry) + \eshiftparm \rw + r}.
\end{eqnarray*}

This formulation provides a fully parametric version of the Cox proportional
hazards model (<ref>), incorporating multivariate normally distributed random
effects with a zero mean and variance $\tau^2$.  The treatment
effect $\eshiftparm$ is interpreted as a log-hazard ratio conditional on
unobserved random effects. For more in-depth information on likelihood-based inference,
see <cit.> and <cit.>.

§.§.§ Marginalised proportional hazards model

Furthermore, we can explore the model proposed by
<cit.>, where the marginal survivor functions are
characterised by models (<ref>), while the correlations within
clusters are captured by a joint multivariate normal distribution.  This
joint modeling approach yields a marginalised formulation for the random
intercept model,

\begin{eqnarray*} \label{mod:MCOX}
\snul(\ry) &=&  \Prob(\rY > \ry \mid \rW = 0) =
\sMEV{\frac{\h(\ry)}{\sqrt{\mparm^2 + 1}}},\\
\sone(\ry) &=&  \Prob(\rY > \ry \mid \rW = 1) =
\sMEV{\frac{\h(\ry) + \eshiftparm}{\sqrt{\mparm^2 + 1}}},\\
\sw(\ry) &=&  \Prob(\rY > \ry \mid \rW = \rw) =
\sMEV{\frac{\h(\ry) + \eshiftparm \rw}{\sqrt{\mparm^2 + 1}}}.\\
\end{eqnarray*}

Here, $\mparm^2$ represents the variance of a cluster-specific intercept.
Using this framework, it becomes possible to quantify the treatment effect
using the marginal hazard ratio computed by $\exp\left(\nicefrac{\eshiftparm}{
\sqrt{\mparm^2 + 1}}\right)$.

Further details on the models, including likelihood-based inference, can be
found in <cit.>.

§.§ Application

To estimate mixed-effects smooth transformation models (Section <ref>) we can use the
tramME package <cit.>, available from

R> install.packages("tramME")
R> library("tramME")

Including a random-intercept for the block used in the randomisation, which
is a combination of the centers $j = 1, \dots
88$ and the stratum $\rs$ ($j \times
\rs = \code{Block}$) in the model, we can account for potential correlation
between patients from the same block.  The corresponding mixed-effects
proportional hazards model can be estimated using the CoxphME

R> CoxphME(iDFS   randarm + (1 | Block), data = CAOsurv, log_first = TRUE)

Coefficient     Estimate     Std. Error     95%-Wald CI

$\eshiftparm$     $ -0.239 $     0.107     $-0.449$ to $-0.029$

$\tau^2$     0.078

The model provides an estimate of the log-hazard ratio $\eshiftparm$, which
is conditional on the unobserved random effects.  The estimated hazard ratio
of $\exp(\hat\eshiftparm) = 0.787$ indicates a prolonged
disease-free survival time in the experimental arm. The estimated variance of the random effect
$R$ is relatively small, with
$\hat{\tau}^2 = 0.078$.
We can further examine
the corresponding marginal estimates of the survivor curves or related
measures by integrating out the random effects (for more details,
see <cit.>).  The corresponding marginal survivor curves for
patients from all blocks are depicted in Figure <ref>.

Mixed-effects proportional hazards model.  Marginal survivor curves
obtained from integrating out the random effects from the model, along with
the non-parametric Turnbull estimates.

The estimated mixed-effects proportional hazards model using CoxphME
translates the analysis conducted in <cit.> into
the smooth transformation framework.  The aim of the primary analysis of
<cit.> was to fit a Cox model for clustered
observations accurately estimating the treatment effect and corresponding
standard errors.  However, at the time of the primary analysis it was not
feasible to estimate the mixed-effects Cox model while accounting for
interval-censoring.  Fortunately, here the discrepancies
between the reported results from the model ignoring interval-censoring and
the fitted one, accounting for it, are negligible.

To obtain a marginal hazard ratio we can estimate the model that capture the
dependence within clusters using a joint multivariate normal distribution
(Section <ref>), which can be fitted using the mtram
function from the tram package.  Estimation is straightforward for
completely exact or interval-censored outcomes within a cluster.  Since
iDFS comprises a mix of different types of outcomes, we handle exact
event times by treating them as interval-censored, accomplished by adding a
4-day window (stored in the object iDFS2 of class `Surv', see Section 5 of <cit.> for details).

R> mtram(Coxph(iDFS2   randarm, data = CAOsurv, log_first = TRUE),
+    formula =   (1 | Block), data = CAOsurv)

Coefficient     Estimate     Std. Error     95%-Wald CI

$\eshiftparm$     $ -0.238 $     0.107     $-0.448$ to $-0.027$

$\mparm$     $ 0.174 $     0.121     $-0.062$ to $\phantom{-}0.411$

The corresponding estimate of the marginal hazard ratio is
$\exp(\nicefrac{\hat\eshiftparm}{\sqrt{\hat\mparm^2 + 1}}) =
0.791$ with empirical 95%-confidence
intervals from $0.647$ to $0.978$.  The results indicate that the
hazard for patients in the experimental arm is reduced by
$ 0.791 $ compared to the hazard in the control arm,
regardless of the block.

§ PERSONALISED MEDICINE

In the context of personalised medicine, our focus now turns towards
modeling heterogeneous treatment effects, that is, log-hazard ratios which
might depend on patient's age.  This involves a
shift towards capturing the individualised responses to treatment,
acknowledging that a global treatment effect may not adequately represent
how different age groups respond to treatment.

§.§ Survival models

§.§.§ Age-varying hazards model

To detect varying treatment effects based on age we can employ models which
estimate an age-varying hazard ratio $\exp(\eshiftparm(\ra))$ <cit.>,

\begin{eqnarray*} \label{mod:HTECOX}
\snul(\ry) &=& \Prob(\rY > \ry \mid \rW = 0, \rA = \ra) =
\sMEV{\h(\ry)},\\
\sone(\ry) &=& \Prob(\rY > \ry \mid \rW = 1, \rA = \ra) =
\sMEV{\h(\ry) + \eshiftparm(\ra)},\\
\sw(\ry) &=& \Prob(\rY > \ry \mid \rW = \rw,  \rA = \ra) =
\sMEV{\h(\ry) + \eshiftparm(\ra) \rw}.
\end{eqnarray*}

This formulation aligns with the model estimated in the analysis of

Such models could be further extended to additionally capture variations in
baseline risks across age by including an age-dependent
log-cumulative baseline hazard function: $\log(\Haznul(\ry \mid \ra)) =
\h(\ry \mid \ra) = \basisy(\ry)^\top\parm + \beta_0(\ra)$.

§.§.§ Tree-based age-varying hazards model

Furthermore, for estimating heterogeneous treatment effects, tree-based Cox
models can also be employed

\begin{eqnarray*} \label{mod:TRT}
\snul(\ry \mid \rA = \ra) &=& \Prob(\rY > \ry \mid \rW = 0, \rA = \ra) =
\sMEV{\h(\ry \mid \ra)},\\
\sone(\ry \mid \rA = \ra) &=& \Prob(\rY > \ry \mid \rW = 1, \rA = \ra) =
\sMEV{\h(\ry \mid \ra) + \eshiftparm(\ra)},\\
\sw(\ry \mid \rA = \ra) &=& \Prob(\rY > \ry \mid \rW = \rw,  \rA = \ra) =
\sMEV{\h(\ry \mid \ra) + \eshiftparm(\ra) \rw},
\end{eqnarray*}

allowing to partition both the log-cumulative baseline hazard
$\log(\Haznul(\ry \mid \ra)) = \h(\ry \mid \ra) =
\basisy(\ry)^\top\parm(\ra)$ and the treatment effect $\eshiftparm(\ra)$
with respect to the different age groups.  In contrast to the model in
Section <ref>, both the log-cumulative baseline hazard and the
log-hazard ratio $\eshiftparm$ depend on age, in this case via an
age-structured tree.

§.§ Application

The hazards model featuring an age-varying treatment effect
(Section <ref>) can be fitted using the tramME package

R> CoxphME(iDFS   randarm + s(age, by = as.ordered(randarm), fx = TRUE, k = 6),
+    data = CAOsurv, log_first = TRUE)

The model estimates a global treatment effect and an
additional smooth effect for age in the experimental arm, allowing the
calculation of an age-varying hazard ratio $\exp(\eshiftparm(\ra))$.

The estimated age-varying hazard ratio curve, shown in Figure <ref>,
indicates that the experimental treatment is more effective than the control
treatment for patients aged $40-70$ years, while for older patients the control
treatment reduces the hazard compared to the experimental treatment.
However, the corresponding confidence interval is notably wide and mostly
overlaps with a hazard ratio of 1.

Age-varying hazards model. Hazard ratio and corresponding 95%-confidence interval estimated by the model shown
along age. The log-likelihood of the corresponding model is $-$2'251.25.

Fitting a hazards model partitioning the log-cumulative baseline hazards and
treatment effect by age, a survival tree (Section <ref>) can be estimated using the trtf package <cit.>.

R> install.packages("trtf")
R> library("trtf")

R> trafotree(Coxph(iDFS   randarm, data = CAOsurv, log_first = TRUE),
+    formula = iDFS   randarm | age, data = CAOsurv,
+    control = ctree_control(teststat = "maximum", minbucket = 40))

The survivor functions corresponding to the terminal nodes of the estimated
tree are shown in Figure <ref>.  The results again indicate that the
experimental treatment is more effective for younger patients, while the
control treatment is more effective for older patients. This result is also
in line with the one previously obtained by <cit.>.

Tree-based age-varying hazards model.
Survival tree depicting the estimated survivor curves of the age-groups
corresponding to the terminal nodes of the partitioned hazards model. The corresponding
in-sample log-likelihood is $-$2'250.42.

§ OTHER EXTENSIONS

§.§ Frailty proportional hazards model

In cases where the assumption of a homogeneous study population falls short,
frailty models offer a valuable alternative.  These models address
situations where the population under consideration comprises heterogeneous
individuals with varying risk levels, reflecting the concept of frailty

To handle such scenarios in the framework of smooth transformation models,
the approach discussed by <cit.> can be employed.  The
corresponding frailty proportional hazards model introduces an unobservable
multiplicative frailty effect $\rF = \rf$ on the hazard,

\begin{eqnarray*} \label{mod:frailty}
\snul(\ry \mid \rF = \rf) &=&
\Prob(\rY > \ry \mid \rW = 0, \rF = \rf) = \sMEV{\h(\ry) + \log(\rf)},\\
\sone(\ry \mid \rF = \rf) &=&
\Prob(\rY > \ry \mid \rW = 1, \rF = \rf) = \sMEV{\h(\ry) + \log(\rf) + \eshiftparm},\\
\sw(\ry \mid \rF = \rf) &=&
\Prob(\rY > \ry \mid \rW = \rw, \rF = \rf) = \sMEV{\h(\ry) + \log(\rf) + \eshiftparm \rw}.
\end{eqnarray*}

The model implies that the proportional hazards assumption, that is
$\nicefrac{\Hazone(\ry \mid \rf)}{\Haznul(\ry \mid \rf)} =
\nicefrac{\hazone(\ry \mid  \rf)}{\haznul(\ry \mid\rf)} = \eshiftparm$,
holds conditional
on frailty $\rF = \rf$.  The frailty $\rF$ specifies a latent random
term, assumed to have a certain non-negative distribution, such as the
gamma, inverse Gaussian or positive stable distribution
<cit.>, which, for identifiability, is scaled such that
$\Ex(\rF) = 1$.

The proportional hazards model with gamma frailty can be fitted in
tram, using the Coxph function specifying the frailty
distribution, for instance frailty = "Gamma".

§.§ Non-parametric likelihood

In this tutorial, we have primarily focused on the implementation of smooth
parametrisation for the log-cumulative baseline hazard function $\h$.
Nevertheless, it is important to highlight that researchers also have the
option to utilize the discussed models that incorporate a non-parametric
version of the transformation function $\h$, should they wish to do so.
The corresponding non-parametric transformation function $\h$ is specified
in terms of $\basisy(\ry_k)^\top \parm = \eparm_k$, where a parameter $\eparm_k$
is assigned to each distinct event time $\ry_k$ with $k = 1, \dots, K - 1$.
A head-to-head comparison of the smooth parametrisation
and the non-parametric version can be found in <cit.>.

§.§ Link function

Undoubtedly, the proportional hazards model stands as a cornerstone in
survival analysis, prominently emerging from specifying the complementary
log-log link (the cumulative distribution function of the standard minimum
extreme value), wherein $\h$ parametrises the log-cumulative baseline
hazards function.  Nevertheless, it is worth noting that researchers have
the opportunity to explore other link functions for all the models shown
above, such as the logit link (as also discussed in detail in
<cit.>), or the probit or log-log link.

For example, specifying a flexible proportional odds model “only” requires to change the link
function from complementary log-log to logit; such a model can be estimated via

R> Colr(iDFS   randarm, data = CAOsurv, log_first = TRUE)

Coefficient     Estimate     Std. Error     95%-Score CI

$\eshiftparm$     $ -0.297 $     0.125     $-0.541$ to $-0.053$

Log-Likelihood     Likelihood Ratio Test     Score Test     Permutation Score Test

$-$2'254.67     $p$ = 0.017     $p$ = 0.017     $p$ = 0.018

This inherent versatility of link functions facilitates the construction
of alternative models, encompassing mean or odds models, by specifying a probit or
logit link respectively.  These alternative
models are well known from the class of acceleration failure time models
(with log-linear $\h$), but extend seamlessly to the more flexible
framework of smooth transformation models.  Moreover, by selecting the
log-log link, it is possible to define a reverse time hazards model.  For a
comprehensive overview, see Table 1 of <cit.>.

§.§ Sample size estimation

The framework of smooth transformation models can also be valuable for
researchers involved in designing new studies.  Simulating from the
illustrated models (using simulate) offers a straightforward approach
for tasks such as sample size estimation. Because the transformation
function $\h$ is smooth, it is relatively simple to invert this function
numerically, such that is becomes possible to apply probability integral
transforms for generating new event times from a fitted model analogously to

§ DISCUSSION

Motivated by the complexities researchers face when navigating various
software implementations for their survival analysis tasks, we highlight the
potential of leveraging the unified framework of smooth transformation
models <cit.> implemented in the tram
R add-on package.  Together with other add-on packages such as
the tramME <cit.> and trtf <cit.>
packages, the tram package <cit.> provides a unified
maximum-likelihood approach to survival analysis, effortlessly extending
classical models to accommodate complex scenarios, thereby offering a
versatile toolkit for the analysis of survival data.

Throughout this tutorial, we present practical strategies for effectively
addressing prominent challenges inherent to survival analysis in
R.  These challenges include interval-censored outcomes,
non-proportional or crossing hazards, dependent censoring, clustered observations
or heterogeneous treatment effects.

The comparative overview of implementations in
Appendix <ref> serves as a crucial validation tool, facilitating
comparisons across multiple package implementations and allowing
researchers to make informed decisions about the most suitable
package for their specific analysis needs.

While we have demonstrated models with a single covariate $\rw$, it is
important to note that the models allow for multivariable linear predictors
$\rx^\top \shiftparm$.  Further, conditional transformation models
<cit.> featuring complex, non-linear covariate
effects can be estimated using the tbm package <cit.>.  The $L_1$ and $L_2$ penalised versions of the models can
also be estimated using the tramnet package
<cit.>.  Additionally, when dealing with multiple event times per subject, multivariate survival models can be estimated
with package tram using the mmlt function

The framework also provides the basis for many interesting extensions.  For
instance, to address scenarios involving cured patients, a noteworthy
extension could involve cure mixture models.  The implementations of such
models in the smooth transformation model framework could, for example,
provide a fully parametric version of the semi-parametric cure mixture model
proposed by <cit.>, with a survival model similar
to the model presented in Section <ref>.
Future work could also explore the implementation of
collapsible models similar to the smooth acceleration failure time model
discussed by <cit.>.  Additionally,
alternative strategies such as the marginalisation of hazard ratios, as
suggested by <cit.>, could be explored further.

In this tutorial we have presented different configurations of smooth
transformation models for survival analysis.  The modular architecture of
this framework not only provides different options, but also
facilitates the combination of individual components.  For example, one can
easily incorporate covariate dependent hazard ratios coupled with random
effects into a smooth transformation model in tramME.  This
flexibility provides users with a unified toolkit to seamlessly address a spectrum of
complex survival analysis tasks.

§ ACKNOWLEDGMENTS

The authors acknowledge financial support by the Swiss National Science
Foundation (grant number 200021_219384).

§ COMPARATIVE OVERVIEW OF R IMPLEMENTATIONS

In the following, we compare the implementations of the models from the tram
<cit.> and tramME <cit.> packages shown in the tutorial, with
alternative methods available in various established R packages
from CRAN.

We fit the  models to the primary endpoint of disease-free survival, which
comprises a mixture of exact times and right- and interval-censored event times (encoded in
iDFS, an object of class `Surv').  In order to further compare the models
with other implementations that cannot handle
interval-censored outcomes, we treat the interval-censored observations as
if they were observed exactly (encoded in DFS, also a `Surv'

We contrast treatment effect estimates (Estimate) and corresponding standard
errors (Std. Error) estimated
by the fitted models.  It is important to note the difference in
interpretation of the estimates (Interpretation).  Additionally, we provide
the in-sample log-likelihood (logLik) of the fitted models, with penalised or
semi-parametric/partial likelihoods highlighted in grey.

§.§ Weibull models

The survival package <cit.>, the icenReg
package <cit.> and the flexsurv package <cit.>
provide alternative implementations of these
models.  Both the survival and icenReg package (specifying model =
"aft") implement acceleration failure time Weibull models, where the effect
can be interpreted as log-acceleration factor (log-AF).  An alternative parametrisation of such models is in
terms of proportional hazards Weibull models, estimating
log-hazard ratios (log-HRs), instead.  This is how the
model, fitted by Survreg, is implemented in the tram package
<cit.>, which can be directly compared to the analogous
parametrisation in the icenReg (specifying model = "ph") and the flexsurv package
(with dist = "weibullPH").  All Weibull
models can handle interval-censoring, owing to the parametric nature of the

The Weibull models can be fitted to the interval-censored outcomes as follows:

R> tram::Survreg(iDFS   randarm, data = CAOsurv, dist = "weibull")
R> icenReg::ic_par(iDFS   randarm, data = CAOsurv, dist = "weibull",
+    model = "ph")
R> flexsurv::flexsurvreg(iDFS   randarm, data = CAOsurv, dist = "weibullPH")
R> survival::survreg(iDFS   randarm, data = CAOsurv, dist = "weibull")
R> icenReg::ic_par(iDFS   randarm, data = CAOsurv, dist = "weibull",
+    model = "aft")

Function     Package     Interpretation     Estimate     Std. Error     logLik

Survreg     tram     log-HR     $ -0.229 $     $ 0.106 $     $-$2'281.17

ic_par     icenReg     log-HR     $ -0.229 $     $ 0.106 $     $-$2'281.17

flexsurvreg     flexsurv     log-HR     $ -0.229 $     $ 0.106 $     $-$2'281.17

survreg     survival     log-AF     $ 0.312 $     $ 0.146 $     $-$2'281.17

ic_par     icenReg     log-AF     $ 0.312 $     $ 0.146 $     $-$2'281.17

As expected, all packages provide equivalent models fits.

§.§ Flexible proportional hazards models

Flexible versions of the proportional hazards model are implemented in
several packages, of which the following accommodate interval-censored
The rstpm2 package
<cit.> and the flexsurv package
provide parametric versions of the model by using splines (analogously to the
approach discussed by <cit.>). We set k = 3 for the number
of knots in the spline for flexsurvspline from the flexsurv package
Alternatively, the icensReg package <cit.> provides a
semi-parametric implementation of the model that can handle interval censoring.

The corresponding models can be fitted to the interval-censored outcomes as follows:

R> tram::Coxph(iDFS   randarm, data = CAOsurv, log_first = TRUE)
R> rstpm2::stpm2(Surv(time = iDFStime, time2 = iDFStime2, event = iDFSevent,
+      type = "interval")   randarm, data = CAOsurv)
R> flexsurv::flexsurvspline(iDFS   randarm, data = CAOsurv, k = 3)
R> icenReg::ic_sp(iDFS   randarm, data = CAOsurv, model = "ph")

Function     Package     Interpretation     Estimate     Std. Error     logLik

Coxph     tram     log-HR     $ -0.234 $     $ 0.106 $     $-$2'255.85

stpm2     rstpm2     log-HR     $ -0.232 $     $ 0.107 $     $-$2'250.48

flexsurvspline     flexsurv     log-HR     $ -0.236 $     $ 0.097 $     $-$2'318.18

ic_sp     icenReg     log-HR     $ -0.230 $     $-$     $-$1'977.29

The models fit similarly across all four packages.
Due to the fact that the computations of the standard errors of ic_sp from the
icenReg package rely on computationally expensive bootstrap sampling,
we did not report any standard errors for this approach.  Also the log-likelihood (in
grey) of the semi-parametric model from the icenReg is not comparable
to the otherwise fully parametric implementations.

The ICsurv package <cit.>
could also potentially handle interval-censored event times. The
TransModel package <cit.>,
featuring an alternative implementation of linear transformation model,
could also serve as an interesting comparator. However, we encountered
difficulties and errors when trying to fit the model using these two packages.

In scenarios where interval-censoring is not taken into account, there are
several other implementations available for fitting corresponding models.
The coxph function from the survival package provides a
semi-parametric approach for exact or right-censored observations
<cit.>. (Note, that again the likelihood of the fitted model is not
comparable to the other fully parametric-models and thus marked in grey).
The rms package <cit.> implements a semi-parametric model, in line with
the model from package survival.

R> tram::Coxph(DFS   randarm, data = CAOsurv, log_first = TRUE)
R> survival::coxph(DFS   randarm, data = CAOsurv)
R> rms::cph(DFS   randarm, data = CAOsurv)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Coxph     tram     log-HR     $ -0.231 $     $ 0.106 $     $-$3'277.59

coxph     survival     log-HR     $ -0.228 $     $ 0.106 $     $-$2'430.66

cph     rms     log-HR     $ -0.228 $     $ 0.106 $     $-$2'430.66

§.§ Stratified proportional hazards models

For comparing stratified flexible proportional hazards models we can again
utilize the model from the rstpm2 and flexsurv packages,
which employ stratified spline-basis functions. The two models can be fitted to the
interval-censored event times as follows

R> tram::Coxph(iDFS | strat   randarm, data = CAOsurv, log_first = TRUE)
R> rstpm2::stpm2(Surv(time = iDFStime, time2 = iDFStime2, event = iDFSevent,
+      type = "interval")   randarm + strata(strat), data = CAOsurv)
R> flexsurv::flexsurvspline(iDFS   randarm + gamma1(strat) + gamma2(strat),
+    data = CAOsurv, k = 3)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Coxph     tram     log-HR     $ -0.232 $     $ 0.107 $     $-$2'231.85

stpm2     rstpm2     log-HR     $ -0.220 $     $ 0.107 $     $-$2'242.88

flexsurvspline     flexsurv     log-HR     $ -0.263 $     $ 0.104 $     $-$2'307.64

The results from the three models are very similar

Now, ignoring interval-censoring, we can, once again, contrast the implementation of
the semi-parametric models from the
survival package and the rms package:

R> tram::Coxph(DFS | strat   randarm, data = CAOsurv, log_first = TRUE)
R> survival::coxph(DFS   randarm + strata(strat), data = CAOsurv)
R> rms::cph(DFS   randarm + strat(strat), data = CAOsurv)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Coxph     tram     log-HR     $ -0.229 $     $ 0.107 $     $-$3'253.85

coxph     survival     log-HR     $ -0.222 $     $ 0.107 $     $-$2'089.54

cph     rms     log-HR     $ -0.222 $     $ 0.107 $     $-$2'089.54

The three model fits are practically equivalent.

We can proceed to compare the stratified version of the Weibull model, for which
we also will ignore interval-censoring due to the
fact that the utilised eha package <cit.> lacks support for interval-censored
data.  Additionally, it is worth highlighting that there is a distinction
from the model fitted using survreg from the survival package
<cit.>.  This model only stratifies the scale parameter of the
Weibull distribution, whereas the models from the eha package and the tram package
estimate both strata-dependent scale and shape parameters.
The survreg function from the survival package fits an acceleration-failure time
Weibull model, while the eha package implements a proportional-hazards Weibull model,
analogously to the Survreg implementation from the tram package.
The model can be fitted to the exact event times, ignoring interval censoring,
as follows

R> tram::Survreg(DFS | strat   randarm, data = CAOsurv)
R> eha::phreg(DFS   randarm + strata(strat), data = CAOsurv)
R> survival::survreg(DFS   randarm + strata(strat), data = CAOsurv)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Survreg     tram     log-HR     $ -0.219 $     $ 0.107 $     $-$3'277.35

phreg     eha     log-HR     $ -0.219 $     $ 0.107 $     $-$3'277.35

survreg     survival     log-AF     $ 0.274 $     $ 0.133 $     $-$3'280.87

The fit of the survreg function from survival package is slightly different.
In contrast, the parametrisation and fits of the model from the eha and
the tram are equivalent.

§.§ Non-proportional hazards models

To the best of our knowledge, there is currently no implementation  available that
estimates an analogous model to the flexible
non-proportional (location-scale) hazards model implemented in tram.

However we can contrast the implementation of the less-flexible Weibull model
with the gamlss package <cit.>
using the WEI2 distribution and the gamlss.cens package <cit.> to account for

R> tram::Survreg(iDFS   randarm | randarm, data = CAOsurv,
+    remove_intercept = FALSE)
R> gamlss::gamlss(formula = iDFS   randarm, sigma.fo =   randarm,
+    family = gamlss.cens::cens(family = "WEI2", type = "interval"),
+    data = CAOsurv[, c("iDFS", "randarm")],
+    control = gamlss.control(n.cyc = 300, trace = FALSE))

Since the scale term in the tram package and the
gamlss package are parameterised differently, we only show the estimates and standard errors of the
location parameter below.

Function     Package     Estimate     Std. Error     logLik

Survreg     tram     $ -0.849 $     $ 0.536 $     $-$2'280.47

gamlss     gamlss     $ -0.948 $     $ 0.542 $     $-$2'280.53

The two implementations provide almost equivalent models fits.

The mpr package <cit.> also offers an implementation for a
non-proportional Weibull model, it, however, does not support interval-censored data. Thus we fit the
models ignoring interval-censoring.

R> tram::Survreg(DFS   randarm | randarm, data = CAOsurv,
+    remove_intercept = FALSE)
R> mpr::mpr(DFS   list(  randarm,   randarm), data = CAOsurv)

Again, we only show the estimates and standard errors of the location parameter below.

Function     Package     Estimate     Std. Error     logLik

Survreg     tram     $ -0.976 $     $ 0.568 $     $-$3'290.43

mpr     mpr     $ -0.975 $     $ 0.567 $     $-$3'290.43

The two implementations also provide equivalent models fits.

§.§ Time-varying hazards model

We can compare the time-varying hazards model from
tram to the flexsurv package <cit.> which allows to
estimate time-varying covariate effects. We start by examining the models
for the interval-censored event times.

R> tram::Coxph(iDFS | randarm   1, data = CAOsurv, log_first = TRUE)
R> flexsurv::flexsurvspline(iDFS   randarm + gamma1(randarm) + gamma2(randarm),
+    data = CAOsurv, k = 3)

The in-sample log-likelihood is $-$2'299.35 for the flexsurv model
and $-$2'321.58 for the tram model. The estimated
time-varying ratios of the cumulative hazards are shown in the plot below.

We will now explore the same models  ignoring interval-censoring.

R> tram::Coxph(DFS | randarm   1, data = CAOsurv, log_first = TRUE)
R> flexsurv::flexsurvspline(DFS   randarm + gamma1(randarm) + gamma2(randarm),
+    data = CAOsurv, k = 3)

The in-sample log-likelihood is $-$2'299.35 for the flexsurv model
and $-$2'321.58 for the tram model, with the computed time-varying ratios of the cumulative hazards shown below.

The time-varying effects estimated from DFS show good agreement, the ratios
differ rather substantially when the models are estimated on the interval-censored data (iDFS).

§.§ Mixed-effects proportional hazards models

The implementation of a mixed-effects proportional
hazards model with flexible log-cumulative baseline hazards for
interval-censored event times is unique in the tramME package <cit.>.
While the rstpm2 package
also accommodates interval-censored event times,
we were not able to fit the corresponding mixed-effects model to our data.

Thus, to contrast the models with other implementations we, again,
need to ignore interval-censoring. For exact event times, we can then
compare the fitted model with the
fully parametric spline-based version from
the rstpm2 package <cit.> and
the semi-parametric model estimated by the
coxme package <cit.>, employing Gaussian random effects using a Laplace
approximation <cit.>.

R> tramME::CoxphME(DFS   randarm + (1 | Block), data = CAOsurv,
+    log_first = TRUE)
R> rstpm2::stpm2(Surv(DFStime, DFSevent)   randarm, data = CAOsurv,
+    cluster = "Block", RandDist = "LogN")
R> coxme::coxme(DFS   randarm + (1 | Block), data = CAOsurv)

Function     Package     Interpretation     Estimate     Std. Error     logLik

CoxphME     tramME     log-HR     $ -0.237 $     $ 0.107 $     $-$3'277.23

stpm2     rstpm2     log-HR     $ -0.234 $     $ 0.107 $     $-$3'272.86

coxme     coxme     log-HR     $ -0.231 $     $ 0.107 $     $-$2'414.48

The fitted models from the three packages agree very well.

§.§ Age-varying hazards model

We can compare the age-varying hazards model from
package tramME to the implementation in the mgcv package
<cit.> which estimates a smooth Cox model via the partial likelihood maximisation.
As the model from the mgcv package only accommodates right-censored observations
we again fit the models ignoring interval-censoring.

R> tramME::CoxphME(DFS   randarm + s(age, by = as.ordered(randarm),
+      fx = TRUE, k = 6), data = CAOsurv, log_first = TRUE)
R> mgcv::gam(DFStime   randarm + s(age, by = as.ordered(randarm),
+      fx = TRUE, k = 6), data = CAOsurv, family = cox.ph(), weights = DFSevent)

The in-sample log-likelihood of the model from the package mgcv is $-$2'426.04 (partial log-likelihood)
and $-$3'272.88 for the tramME model.
The estimated age-varying hazard ratios and corresponding 95%-confidence intervals
are shown in the plot below.

The hazard ratio curves and confidence intervals estimated by the two packages are equivalent.

§.§ Frailty proportional hazards models

For models featuring a gamma frailty, we can contrast implementations
using a semi-parametric approach or the spline-based approach from the
rstpm2 package <cit.>.  The coxph model from the survival package uses a semi-parametric approach and
estimates the frailty term using penalised regression <cit.>.
The frailtyEM <cit.> and the frailtypack package <cit.> also feature models with
semi-parametric baseline hazards. Again we fit the models ignoring interval-censoring.

R> tram::Coxph(DFS   randarm, data = CAOsurv, log_first = TRUE,
+    frailty = "Gamma")
R> rstpm2::stpm2(Surv(DFStime, DFSevent)   randarm, data = CAOsurv,
+    cluster = "id", RandDist = "Gamma")
R> survival::coxph(DFS   randarm + frailty(id, distribution = "gamma"),
+    data = CAOsurv)
R> frailtyEM::emfrail(DFS   randarm + cluster(id), data = CAOsurv)
R> frailtypack::frailtyPenal(DFS   randarm + cluster(id), data = CAOsurv,
+    RandDist = "Gamma", n.knots = 10, kappa = 1)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Coxph     tram     log-HR     $ -0.632 $     $ 0.467 $     $-$3'263.81

stpm2     rstpm2     log-HR     $ -0.685 $     $ 0.670 $     $-$3'264.88

coxph     survival     log-HR     $ -0.406 $     $ 0.159 $     $-$1'944.22

emfrail     frailtyEM     log-HR     $ -0.384 $     $ 0.153 $     $-$2'430.45

frailtyPenal     frailtypack     log-HR     $ -0.660 $     $ 0.248 $     $-$3'259.82

The fitted models vary considerably across packages.

§.§ Flexible proportional odds models

We can compare the fit of the flexible proportional odds model with the implementation in
the rstpm2 package <cit.> and package flexsurv <cit.>.
The Gprop.odds function from
package timereg <cit.> can
also estimate a flexible proportional odds model using the partial likelihood, thus
we again compare the models ignoring interval-censoring.

R> tram::Colr(DFS   randarm, data = CAOsurv, log_first = TRUE)
R> rstpm2::stpm2(Surv(DFStime, DFSevent)   randarm, data = CAOsurv,
+    link.type = "PO")
R> flexsurv::flexsurvspline(iDFS   randarm, data = CAOsurv, k = 3,
+    scale = "odds")
R> timereg::Gprop.odds(DFS   prop(randarm), data = CAOsurv)

Function     Package     Interpretation     Estimate     Std. Error     logLik

Colr     tram     log-OR     $ -0.297 $     $ 0.125 $     $-$3'276.28

stpm2     rstpm2     log-OR     $ -0.294 $     $ 0.125 $     $-$3'272.44

flexsurvspline     flexsurv     $-$log-OR     $ 0.265 $     $ 0.122 $     $-$2'323.93

Gprop.odds     timereg     log-OR     $ -0.268 $     $ 0.125 $

The fitted models practically agree among the four packages.
Note, that we were not able to find the in-sample log-likelihood from the model object of
the timereg package and thus do not report it here.

§.§ Computational details

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;
LAPACK version 3.10.0

[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=de_CH.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=de_CH.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=de_CH.UTF-8       LC_NAME=C
[9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=de_CH.UTF-8 LC_IDENTIFICATION=C

time zone: Europe/Zurich
tzcode source: system (glibc)

attached base packages:
[1] parallel  splines   grid      stats     graphics  grDevices utils
[8] datasets  methods   base

other attached packages:
[1] timereg_2.0.5     mgcv_1.9-0        frailtypack_3.5.1
[4] doBy_4.6.20       survC1_1.0-3      boot_1.3-28.1
[7] frailtyEM_1.0.1   parfm_2.7.7       optimx_2023-10.21
[10] coxme_2.2-18.1    bdsmatrix_1.3-6   gamlss.cens_5.0-7
[13] gamlss_5.4-20     nlme_3.1-164      gamlss.dist_6.1-1
[16] gamlss.data_6.0-2 mpr_1.0.6         flexsurv_2.2.2
[19] rstpm2_1.6.3      eha_2.11.2        ICsurv_1.0.1
[22] rms_6.7-1         Hmisc_5.1-1       TransModel_2.3
[25] icenReg_2.0.16    coda_0.19-4.1     fastGHQuad_1.0.1
[28] Rcpp_1.0.12       xtable_1.8-4      colorspace_2.1-0
[31] coin_1.4-3        multcomp_1.4-25   TH.data_1.1-2
[34] MASS_7.3-60.0.1   survival_3.5-7    tramME_1.0.5
[37] ATR_0.1-1         SparseGrid_0.8.2  trtf_0.4-2
[40] partykit_1.2-20   libcoin_1.0-10    tram_1.0-2
[43] mvtnorm_1.2-4     mlt_1.5-0         basefun_1.1-4
[46] variables_1.1-1   knitr_1.45

loaded via a namespace (and not attached):
[1] rstudioapi_0.15.0     magrittr_2.0.3        modeltools_0.2-23
[4] nloptr_2.0.3          rmarkdown_2.25        vctrs_0.6.5
[7] minqa_1.2.6           base64enc_0.1-3       htmltools_0.5.7
[10] polspline_1.1.24      polynom_1.4-1         broom_1.0.5
[13] deSolve_1.40          Formula_1.2-5         parallelly_1.36.0
[16] pracma_2.4.4          htmlwidgets_1.6.4     sandwich_3.1-0
[19] rootSolve_1.8.2.4     zoo_1.8-12            TMB_1.9.10
[22] alabama_2023.1.0      lifecycle_1.0.4       iterators_1.0.14
[25] pkgconfig_2.0.3       Matrix_1.6-4          R6_2.5.1
[28] fastmap_1.1.1         future_1.33.1         digest_0.6.34
[31] numDeriv_2016.8-1.1   fansi_1.0.6           compiler_4.3.2
[34] microbenchmark_1.4.10 withr_3.0.0           htmlTable_2.4.2
[37] backports_1.4.1       BB_2019.10-1          lava_1.7.3
[40] quantreg_5.97         tools_4.3.2           foreign_0.8-86
[43] mstate_0.3.2          future.apply_1.11.1   nnet_7.3-19
[46] glue_1.7.0            quadprog_1.5-8        inum_1.0-5
[49] checkmate_2.3.1       cluster_2.1.6         generics_0.1.3
[52] orthopolynom_1.0-6.1  gtable_0.3.4          tidyr_1.3.1
[55] sn_2.1.1              data.table_1.15.0     Deriv_4.1.3
[58] utf8_1.2.4            foreach_1.5.2         pillar_1.9.0
[61] stringr_1.5.1         dplyr_1.1.4           lattice_0.22-5
[64] SparseM_1.81          tidyselect_1.2.0      muhaz_1.2.6.4
[67] expint_0.1-8          gridExtra_2.3         stats4_4.3.2
[70] xfun_0.41             expm_0.999-8          statmod_1.5.0
[73] matrixStats_1.2.0     stringi_1.8.3         evaluate_0.23
[76] codetools_0.2-19      msm_1.7.1             bbmle_1.0.25.1
[79] coneproj_1.17         tibble_3.2.1          cli_3.6.2
[82] rpart_4.1.23          munsell_0.5.0         globals_0.16.2
[85] MatrixModels_0.5-3    assertthat_0.2.1      ggplot2_3.4.4
[88] lme4_1.1-35.1         listenv_0.9.1         matrixcalc_1.0-6
[91] scales_1.3.0          purrr_1.0.2           rlang_1.1.3
[94] mnormt_2.1.1