Bayesian Survival Analysis Using the \pkgrstanarm \proglangR Package

We assume that a true event time for individual $i$ ($i=1,...,N$) exists and can be denoted $T_{i}^{*}$. However, in practice $T_{i}^{*}$ may not be observed due to left, right, or interval censoring. We therefore observe outcome data $\mathcal{D}_{i}=\{T_{i},T_{i}^{U},T_{i}^{E},d_{i}\}$ for individual $i$ where:

• •

  $T_{i}$ denotes the observed event or censoring time;

• •

  $T_{i}^{U}$ denotes the observed upper limit for interval censored individuals;

• •

  $T_{i}^{E}$ denotes the observed entry time (the time at which an individual became at risk of experiencing the event); and

• •

  $d_{i}\in\{0,1,2,3\}$ denotes an event indicator taking value 0 if individual $i$ was right censored (i.e. $T_{i}^{*}>T_{i}$), value 1 if individual $i$ was uncensored (i.e. $T_{i}^{*}=T_{i}$), value 2 if individual $i$ was left censored (i.e. $T_{i}^{*}<T_{i}$), or value 3 if individual $i$ was interval censored (i.e. $T_{i}<T_{i}^{*}<T_{i}^{U}$).

Our modelling approaches are twofold. First, we define a class of models on the hazard scale. This includes both proportional and non-proportional hazard regression models. Second, we define a class of models on the scale of the survival time. These are often known as accelerated failure time (AFT) models and can include both time-fixed and time-varying acceleration factors.

These two classes of models and their respective features are described in the following sections.

Under all of the previous model formulations our linear predictor can be defined as:

where $\boldsymbol{X}_{i}(t)=[1,x_{i1}(t),...,x_{iP}(t)]$ denotes a vector of covariates with $x_{ip}(t)$ denoting the observed value of $p^{th}$ $(p=1,...,P)$ covariate for the $i^{th}$ $(i=1,...,N)$ individual at time $t$, and $\boldsymbol{\beta}(t)=[\beta_{0},\beta_{1}(t),...,\beta_{P}(t)]$ denotes a vector of parameters with $\beta_{0}$ denoting an intercept parameter and $\beta_{p}(t)$ denoting the possibly time-varying coefficient for the $p^{th}$ covariate.

The definition of the linear predictor in Equation 20 can be extended to allow for shared frailty or other clustering effects.

Suppose that the individuals in our sample belong to a series of clusters. The clusters may represent for instance hospitals, families, or GP clinics. We denote the $i^{th}$ individual ($i=1,...,N_{j}$) as a member of the $j^{th}$ cluster ($j=1,...,J$). Moreover, to indicate the fact that individual $i$ is now a member of cluster $j$ we index the observed data (i.e. event times, event indicator, and covariates) with a subscript $j$, that is $T_{ij}^{*}$, $\mathcal{D}_{ij}=\{T_{ij},T_{ij}^{U},T_{ij}^{E},d_{ij}\}$ and $X_{ij}(t)$, as well as estimated quantities such as the hazard rate, cumulative hazard, survival probability, and linear predictor, that is $h_{ij}(t)$, $H_{ij}(t)$, $S_{ij}(t)$, and $\eta_{ij}(t)$.

To allow for intra-cluster correlation in the event times we include cluster-specific random effects in the linear predictor as follows:

where $\boldsymbol{Z}_{ij}$ denotes a vector of covariates for the $i^{th}$ individual in the $j^{th}$ cluster, with an associated vector of cluster-specific parameters $\boldsymbol{b}_{j}$. We assume that the cluster-specific parameters are normally distributed such that $\boldsymbol{b}_{j}\sim N(0,\boldsymbol{\Sigma}_{b})$ for some variance-covariance matrix $\boldsymbol{\Sigma}_{b}$. We assume that $\boldsymbol{\Sigma}_{b}$ is unstructured, that is each variance and covariance term is allowed to be different.

In most cases $\boldsymbol{b}_{j}$ will correspond to just a cluster-specific random intercept (often known as a "shared frailty" term) but more complex random effects structures are possible.

For simplicitly of notation Equation 26 also assumes just one clustering factor in the model (indexed by $j=1,...,J$). However it is possible to extend the model to multiple clustering factors. For example, suppose that the $i^{th}$ individual was clustered within the $j^{th}$ hospital that was clustered within the $k^{th}$ geographical region. Then we would have hospital-specific random effects $\boldsymbol{b}_{j}\sim N(0,\boldsymbol{\Sigma}_{b})$ and region-specific random effects $\boldsymbol{u}_{k}\sim N(0,\boldsymbol{\Sigma}_{u})$ and assume $\boldsymbol{b}_{j}$ and $\boldsymbol{u}_{k}$ are independent for all $(j,k)$. Multiple clustering factors are accommodated as part of the survival modelling functionality in \pkgrstanarm.

The log posterior for the $i^{th}$ individual in the $j^{th}$ cluster can be specified as:

where $\log p(\mathcal{D}_{ij}\mid\boldsymbol{\theta},\boldsymbol{b}_{j})$ is the log likelihood for the outcome data, $\log p(\boldsymbol{b}_{j}\mid\boldsymbol{\theta})$ is the log likelihood for the distribution of any cluster-specific parameters (i.e. random effects) when relevant, and $\log p(\boldsymbol{\theta})$ represents the log likelihood for the joint prior distribution across all remaining unknown parameters.

Allowing for the three forms of censoring (left, right, and interval censoring) and potential delayed entry (i.e. left truncation) the log likelihood for the survival model takes the form:

where $I(x)$ is an indicator function taking value 1 if $x$ is true and 0 otherwise. That is, each individual’s contribution to the likelihood depends on the type of censoring for their event time.

The last term on the right hand side of Equation 28 accounts for delayed entry. When an individual is at risk from time zero (i.e. no delayed entry) then $T_{ij}^{E}=0$ and $S_{ij}(0)=1$ meaning that the last term disappears from the likelihood.

For each of the parameters a number of prior distributions are available. Default choices exist, but the user can explicitly specify the priors if they wish.

Estimation in \pkgrstanarm is based on either full Bayesian inference (Hamiltonian Monte Carlo) or approximate Bayesian inference (either mean-field or full-rank variational inference). The default is full Bayesian inference, but the user can change this if they wish. The approximate Bayesian inference algorithms are much faster, but they only provide approximations for the joint posterior distribution and are therefore not recommended for final inference.

Hamiltonian Monte Carlo is a form of Markov chain Monte Carlo (MCMC) in which information about the gradient of the log posterior is used to more efficiently sample from the posterior space. \proglangStan uses a specific implementation of Hamiltonian Monte Carlo known as the No-U-Turn Sampler (NUTS) (Hoffman:2014). A benefit of NUTS is that the tuning parameters are handled automatically during a "warm-up" phase of the estimation. However the \pkgrstanarm modelling functions provide arguments that allow the user to retain control over aspects such as the number of MCMC chains, number of warm-up and sampling iterations, and number of computing cores used.

If our survival model does not contain any clustering effects (i.e. it is not a multilevel survival model) then our prediction framework is more straightforward. Let $\mathcal{D}=\{\mathcal{D}_{i};i=1,...,N\}$ denote the entire collection of outcome data in our sample and let $T_{\max}=\max\{T_{i},T_{i}^{U},T_{i}^{E};i=1,...,N\}$ denote the maximum event or censoring time across all individuals in our sample.

Suppose that for some individual $i^{*}$ (who may or may not have been in our sample) we have covariate vector $\boldsymbol{x}_{i^{*}}$. Note that the covariate data must be time-fixed. The predicted probability of being event-free at time $0<t\leq T_{\max}$, denoted $\hat{S}_{i^{*}}(t)$, can be generated from the posterior predictive distribution:

We approximate this posterior predictive distribution by drawing from $p(\hat{S}_{i^{*}}(t)\mid\boldsymbol{x}_{i^{*}},\boldsymbol{\theta}^{(l)})$ where $\boldsymbol{\theta}^{(l)}$ is the $l^{th}$ $(l=1,...,L)$ MCMC draw from the posterior distribution $p(\boldsymbol{\theta}\mid\mathcal{D})$.

When there are clustering effects in the model (i.e. multilevel survival models) then our prediction framework requires conditioning on the cluster-specific parameters. Let $\mathcal{D}=\{\mathcal{D}_{ij};i=1,...,N_{j},j=1,...,J\}$ denote the entire collection of outcome data in our sample and let $T_{\max}=\max\{T_{ij},T_{ij}^{U},T_{ij}^{E};i=1,...,N_{j},j=1,...,J\}$ denote the maximum event or censoring time across all individuals in our sample.

Suppose that for some individual $i^{*}$ (who may or may not have been in our sample) and who is known to come from cluster $j^{*}$ (which may or may not have been in our sample) we have covariate vectors $\boldsymbol{x}_{i^{*}j^{*}}$ and $\boldsymbol{z}_{i^{*}j^{*}}$. Note again that the covariate data is assumed to be time-fixed.

If individual $i^{*}$ does in fact come from a cluster $j^{*}=j$ (for some $j\in\{1,...,J\}$) in our sample then the predicted probability of being event-free at time $0<t\leq T_{\max}$, denoted $S_{i^{*}j}(t)$, can be generated from the posterior predictive distribution:

Since cluster $j$ was included in our sample data it is easy for us to approximate this posterior predictive distribution by drawing from $p(\hat{S}_{i^{*}j}(t)\mid\boldsymbol{x}_{i^{*}j},\boldsymbol{z}_{i^{*}j},\boldsymbol{\theta}^{(l)},\boldsymbol{b}_{j}^{(l)})$ where $\boldsymbol{\theta}^{(l)}$ and $\boldsymbol{b}_{j}^{(l)}$ are the $l^{th}$ $(l=1,...,L)$ MCMC draws from the joint posterior distribution $p(\boldsymbol{\theta},\boldsymbol{b}_{j}\mid\mathcal{D})$.

Alternatively, individual $i^{*}$ may come from a new cluster $j^{*}\neq j$ (for all $j\in\{1,...,J\}$) that was not in our sample. The predicted probability of being event-free at time $0<t\leq T_{\max}$ is therefore denoted $\hat{S}_{i^{*}j^{*}}(t)$ and can be generated from the posterior predictive distribution:

where $\boldsymbol{\tilde{b}}_{j^{*}}$ denotes the cluster-specific parameters for the new cluster. We can obtain draws for $\boldsymbol{\tilde{b}}_{j^{*}}$ during estimation of the model (in a similar manner as for $\boldsymbol{b}_{j}$). At the $l^{th}$ iteration of the MCMC sampler we obtain $\boldsymbol{\tilde{b}}_{j^{*}}^{(l)}$ as a random draw from the posterior distribution of the cluster-specific parameters and store it for later use in predictions. The set of random draws $\boldsymbol{\tilde{b}}_{j^{*}}^{(l)}$ for $l=1,...,L$ then allow us to essentially marginalise over the distribution of the cluster-specific parameters. This is the approach used in \pkgrstanarm to generate survival predictions for individuals in new clusters that were not part of the original sample.

In some instances we want to evaluate the predicted survival probability conditional on a last known survival time. This is known as a conditional survival probability.

Suppose that individual $i^{*}$ is known to be event-free up until $C_{i^{*}}$ and we wish to predict the survival probability at some time $t>C_{i^{*}}$. To do this we draw from the conditional posterior predictive distribution:

or – equivalently – for multilevel survival models we have individual $i^{*}$ in cluster $j^{*}$ who is known to be event-free up until $C_{i^{*}j^{*}}$:

All of the previously discussed predictions require conditioning on some covariate values $\boldsymbol{x}_{ij}$ and $\boldsymbol{z}_{ij}$. Even if we have a multilevel survival model and choose to marginalise over the distribution of the cluster-specific parameters, we are still obtaining predictions at some known unique values of the covariates.

However sometimes we wish to generate an "average" survival probability. One possible approach is to predict at the mean value of all covariates (Cupples:1995). However this doesn’t always make sense, especially not in the presence of categorical covariates. For instance, suppose our covariates are gender and a treatment indicator. Then predicting for an individual at the mean of all covariates might correspond to a 50% male who was 50% treated. That does not make sense and is not what we wish to do.

A better alternative is to average over the individual survival probabilties. This essentially provides an approximation to marginalising over the joint distribution of the covariates. At any time $t$ it is possible to obtain a so-called standardised survival probability, denoted $\hat{S}^{*}(t)$, by averaging the individual-specific survival probabilities:

where $\hat{S}_{i}(t)$ is the predicted survival probability for individual $i$ ($i=1,...,N^{P}$) at time $t$, and $N^{P}$ is the number of individuals included in the predictions. For multilevel survival models the calculation is similar and follows quite naturally (details not shown).

Note however that if $N^{P}$ is not sufficiently large (for example we predict individual survival probabilities using covariate data for just $N^{P}=2$ individuals) then averaging over their covariate distribution may not be meaningful. Similarly, if we estimated a multilevel survival model and then predicted standardised survival probabilities based on just $N^{P}=2$ individuals from our sample, the joint distribution of their cluster-specific parameters would likely be a poor representation of the distribution of cluster-specific parameters for the entire sample and population.

It is therefore better to calculate standardised survival probabilities by setting $N^{P}$ equal to the total number of individuals in the original sample (i.e. $N^{P}=N$. This approach can then also be used for assessing the fit of the survival model in \pkgrstanarm (see the \codeps_check() function described in Section 5). Posterior predictive draws of the standardised survival probability are evaluated at a series of time points between 0 and $T_{\max}$ using all individuals in the estimation sample and the predicted standardised survival curve is overlaid with the observed Kaplan-Meier survival curve.

The \pkgrstanarm package is built on top of the \pkgrstan \proglangR package (Stan:2019), which is the \proglangR interface for \proglangStan. Models in \pkgrstanarm are written in the \proglangStan programming language, translated into \proglangC++ code, and then compiled at the time the package is built. This means that for most users – who install a binary version of \pkgrstanarm from the Comprehensive R Archive Network (CRAN) – the models in \pkgrstanarm will be pre-compiled. This is beneficial for users because there is no compilation time either during installation or when they estimate a model.

Survival models in \pkgrstanarm are implemented around the \codestan_surv() modelling function.

The function signature for \codestan_surv() is: {Schunk} {Sinput} R> stan_surv(formula, data, basehaz = "ms", basehaz_ops, qnodes = 15, + prior = normal(), prior_intercept = normal(), prior_aux, + prior_smooth = exponential(autoscale = FALSE), + prior_covariance = decov(), prior_PD = FALSE, + algorithm = c("sampling", "meanfield", "fullrank"), + adapt_delta = 0.95, …)

The following provides a brief description of the main features of each of these arguments:

• •

  The \codeformula argument accepts objects built around the standard \proglangR formula syntax (see \codestats::formula()). The left hand side of the formula should be an object returned by the \codeSurv() function in the \pkgsurvival package (Therneau:2019). Any random effects structure (for multilevel survival models) can be specified on the right hand side of the formula using the same syntax as the \pkglme4 \proglangR package (Bates:2015) as shown in the example in Section LABEL:sec:multilevelmodel.

  By default, any covariate effects specified in the fixed-effect part of the model \codeformula are included under a proportional hazards assumption (for models estimated using a hazard scale formulation) or under the assumption of time-fixed acceleration factors (for models estimated using an AFT formulation). Time-varying effects are specified in the model \codeformula by wrapping the covariate name in the \codetve() function. For example, if we wanted to estimate a time-varying effect for the covariate \codesex then we could specify \codetve(sex) in the model formula, e.g. \codeformula = Surv(time, status)   tve(sex) + age. The \codetve() function is a special function that only has meaning when used in the \codeformula of a model estimated using \codestan_surv(). Its functionality is demonstrated in the worked examples in Sections LABEL:sec:tvebs and LABEL:sec:tvepw.

• •

  The \codedata argument accepts an object inheriting the class ‘\codedata.frame’, in other words the usual \proglangR data frame.

• •

  The choice of parametric baseline hazard (or baseline survival distribution for AFT models) is specified via the \codebasehaz argument. For the M-spline (\code"ms") and B-spline (\code"bs") models additional options related to the spline degree $\delta$, knot locations $\boldsymbol{k}$, or degrees of freedom $L$ can be specified as a list and passed to the \codebasehaz_ops argument. For example, specifying \codebasehaz = "ms" and \codebasehaz_ops = list(degree = 2, knots = c(10,20)) would request a baseline hazard modelled using quadratic M-splines with two internal knots located at $t=10$ and $t=20$.

• •

  The argument \codeqnodes is a control argument that allows the user to specify the number of quadrature nodes when quadrature is required (as described in Section 3.2).

• •

  The \codeprior family of arguments allow the user to specify the prior distributions for each of the parameters, as follows:

  • –
    \code

    prior relates to the time-fixed regression coefficients;

  • –
    \code

    prior_intercept relates to the intercept in the linear predictor;

  • –
    \code

    prior_aux relates to the so-called "auxiliary" parameters in the baseline hazard ($\gamma$ for the Weibull and Gompertz models or $\boldsymbol{\gamma}$ for the M-spline and B-spline models);

  • –
    \code

    prior_smooth relates to the hyperparameter $\tau_{p}$ when the $p^{th}$ covariate has a time-varying effect; and

  • –
    \code

    prior_covariance relates to the covariance matrix for the random effects when a multilevel survival model is being estimated.

• •

  The remaining arguments (\codeprior_PD, \codealgorithm, and \codeadapt_delta) are optional control arguments related to estimation in Stan:

  • –

    Setting \codeprior_PD = TRUE states that the user only wants to draw from the prior predictive distribution and not condition on the data.

  • –

    The \codealgorithm argument specifies the estimation routine to use. This includes either Hamiltonian Monte Carlo (\code"sampling") or one of the variational Bayes algorithms (\code"meanfield" or \code"fullrank"). The model specification is agnostic to the chosen \codealgorithm. That is, the user can choose from any of the available algorithms regardless of the specified model.

  • –

    The \codeadapt_delta argument controls the target average acceptance probability. It is only relevant when \codealgorithm = "sampling" in which case \codeadapt_delta should be between 0 and 1, with higher values leading to smaller step sizes and therefore a more robust sampler but longer estimation times.

The model returned by \codestan_surv() is an object of class ‘\codestansurv’ and inheriting the ‘\codestanreg’ class. It is effectively a list with a number of important attributes. There are a range of post-estimation functions that can be called on ‘\codestansurv’ (and ‘\codestanreg’) objects – some of the most important ones are described in Section 5.3.

The \pkgrstanarm package provides a range of post-estimation functions that can be used after fitting the survival model. This includes functions for inference (e.g. reporting parameter estimates), diagnostics (e.g. assessing model fit), and generating predictions. We highlight the most important ones here:

• •

  The \codeprint() and \codesummary() functions provide reports of parameter estimates and some summary information on the data (e.g. number of observations, number of events, etc). They each provide varying levels of detail. For example, the \codesummary() method provides diagnostic measures such as Gelman and Rubin’s Rhat statistic (Gelman:1992) for assessing convergence of the MCMC chains and the number of effective MCMC samples. On the other hand, the \codeprint() method is more concise and does not provide this level of additional detail.

• •

  The \codefixef() and \coderanef() functions report the fixed effect and random effect parameter estimates, respectively.

• •

  The \codeposterior_survfit() function is the primary function for generating survival predictions. The type of prediction is specified via the \codetype arguments and can currently be any of the following:

  • –
    \code

    "surv": the estimated survival probability;

  • –
    \code

    "cumhaz": the estimated cumulative hazard;

  • –
    \code

    "haz": the estimated hazard rate;

  • –
    \code

    "cdf": the estimated failure probability;

  • –
    \code

    "logsurv": the estimated log survival probability;

  • –
    \code

    "logcumhaz": the estimated log cumulative hazard;

  • –
    \code

    "loghaz": the estimated log hazard rate; or

  • –
    \code

    "logcdf": the estimated log failure probability.

  There are additional arguments to \codeposterior_survfit() that control the time at which the predictions are generated (\codetimes), whether they are generated across a time range (referred to as extrapolation, see \codeextrapolate), whether they are conditional on a last known survival time (\codecondition), and whether they are averaged across individuals (referred to as standardised predictions, see \codestandardise). The returned predictions are a data frame with a special class called ‘\codesurvfit.stansurv’. The ‘\codesurvfit.stansurv’ class has both \codeprint() and \codeplot() methods that can be called on it. These will be demonstrated as part of the examples in Section 6.

• •

  The \codeloo() and \codewaic() functions report model fit statistics. The former is based on approximate leave-one-out cross validation (Vehtari:2017) and is recommended. The latter is a less preferable alternative that reports the Widely Applicable Information Criterion (WAIC) criterion (Watanabe:2010). Both of these functions are built on top of the \pkgloo \proglangR package (Vehtari:2019). The values (objects) returned by either \codeloo() or \codewaic() can also be passed to the \codeloo_compare() function to compare different models estimated on the same dataset. This will be demonstrated as part of the examples in Section 6.

• •

  The \codelog_lik() function generates a pointwise log likelihood matrix. That is, it calculates the log likelihood for each observation (either in the original dataset or some new dataset) using each MCMC draw of the model parameters.

• •

  The \codeplot() function allows for a variety of plots depending on the input to the \codeplotfun argument. The default is to plot the estimated baseline hazard (\codeplotfun = "basehaz"), but alternatives include a plot of the estimated time-varying hazard ratio for models with time-varying effects (\codeplotfun = "tve"), plots summarising the parameter estimates (e.g. posterior densities or posterior intervals), and plots providing diagnostics (e.g. MCMC trace plots).

• •

  The \codeps_check() function provides a quick diagnostic check for the fitted survival function. It is based on the estimation sample and compares the predicted standardised survival curve to the observed Kaplan-Meier survival curve.

In this example we fit a proportional hazards model with a flexible parametric baseline hazard modelled using cubic M-splines. This is an example of the default survival model estimated by the \codestan_surv() function in \pkgrstanarm.

We will use the German Breast Cancer Study Group dataset (Royston:2002). The data consist of $N=686$ patients with primary node positive breast cancer recruited between 1984-1989. The primary response is time to recurrence or death. Median follow-up time was 1084 days. Overall, there were 299 (44%) events and the remaining 387 (56%) individuals were right censored. We concern our analysis here with a 3-category baseline covariate for cancer prognosis (good/medium/poor).

Let us fit the proportional hazards model: {Schunk} {Sinput} R> mod1 <- stan_surv(formula = Surv(recyrs, status)   group, + data = bcancer, + chains = CHAINS, + cores = CORES, + seed = SEED, + iter = ITER)

Since there are no time-varying effects in the model (i.e. we did not wrap any covariates in the \codetve() function) there is a closed form expression for the cumulative hazard and survival function and so the model is relatively fast to fit. Specifically, the model takes  3.5 seconds for each MCMC chain based on the default 2000 MCMC iterations (1000 warm up, 1000 sampling) on a standard desktop.

We can easily obtain the estimated hazard ratios for the 3-category group covariate using the generic \codeprint() method for ‘\codestansurv’ objects, as follows: {Schunk} {Sinput} R> print(mod1, digits = 2) {Soutput} stan_surv baseline hazard: M-splines on hazard scale formula: Surv(recyrs, status)   group observations: 686 events: 299 (43.6right censored: 387 (56.4delayed entry: no —— Median MAD_SD exp(Median) (Intercept) -0.65 0.18 NA groupMedium 0.82 0.17 2.28 groupPoor 1.60 0.15 4.95 m-splines-coef1 0.00 0.00 NA m-splines-coef2 0.02 0.01 NA m-splines-coef3 0.40 0.07 NA m-splines-coef4 0.06 0.05 NA m-splines-coef5 0.21 0.12 NA m-splines-coef6 0.30 0.16 NA

—— * For help interpreting the printed output see ?print.stanreg * For info on the priors used see ?prior_summary.stanreg

From this output we see that individuals in the groups with \codePoor or \codeMedium prognosis have much higher rates of death relative to the group with \codeGood prognosis. The hazard of death in the \codePoor prognosis group is approximately 5-fold higher than the hazard of death in the \codeGood prognosis group. Similarly, the hazard of death in the \codeMedium prognosis group is approximately 2-fold higher than the hazard of death in the \codeGood prognosis group.

It may also be of interest to compare the different types of parametric baseline hazards we could have used for this model (some more restricting, others more flexible). Let us fit several models, each with a different baseline hazard: {Schunk} {Sinput} R> mod1_exp <- update(mod1, basehaz = "exp") R> mod1_weibull <- update(mod1, basehaz = "weibull") R> mod1_gompertz <- update(mod1, basehaz = "gompertz") R> mod1_bspline <- update(mod1, basehaz = "bs") R> mod1_mspline1 <- update(mod1, basehaz = "ms") R> mod1_mspline2 <- update(mod1, basehaz = "ms", basehaz_ops = list(df = 9))

The default action of the \codeplot() method for ‘\codestansurv’ objects is to plot the estimated baseline hazard. We will use this method to plot the baseline hazard for each of the competing models. First, we write a little helper function to adjust the y-axis limits and add a centered title on each plot, as follows: {Schunk} {Sinput} R> plotfun <- function(model, title) + plot(model, plotfun = "basehaz") + + coord_cartesian(ylim = c(0,0.4)) + + labs(title = title) + + theme(plot.title = element_text(hjust = 0.5)) + and then we generate each of the plots: {Schunk} {Sinput} R> p_exp <- plotfun(mod1_exp, "Exponential") R> p_weibull <- plotfun(mod1_weibull, "Weibull") R> p_gompertz <- plotfun(mod1_gompertz, "Gompertz") R> p_bspline <- plotfun(mod1_bspline, "B-splines with\ntwointernal knots") R> p_mspline1 <- plotfun(mod1_mspline1, "M-splines with\ntwointernal knots") R> p_mspline2 <- plotfun(mod1_mspline2, "M-splines with\nfiveinternal knots") and then combine the plots using the \codeplot_grid() function from the \pkgcowplot \proglangR package (Wilke:2019): {Schunk} {Sinput} R> p_combined <- plot_grid(p_exp, + p_weibull, + p_gompertz, + p_bspline, + p_mspline1, + p_mspline2, + ncol = 3)

Figure 1 shows the resulting plot with the estimated baseline hazard for each model and 95% posterior uncertainty limits. We can clearly see from the plot the additional flexibility the cubic spline models provide. They are able to capture at least two turning points in the hazard function, one around 1.5 years and another one around 4 years. It is also interesting to note that there appears to be very little change in the fit of the M-spline model when we increase the number of internal knots from two to five.

We can also compare the fit of these models using the \codeloo method for \codestansurv objects: {Schunk} {Sinput} R> compare_models(loo(mod1_exp), + loo(mod1_weibull), + loo(mod1_gompertz), + loo(mod1_bspline), + loo(mod1_mspline1), + loo(mod1_mspline2)) {Soutput} Model formulas: : NULL : NULL : NULL : NULL : NULL : NULL elpd_diff se_diff mod1_mspline1 0.0 0.0 mod1_bspline -0.4 1.7 mod1_mspline2 -1.6 1.5 mod1_weibull -18.0 5.3 mod1_gompertz -31.5 6.1 mod1_exp -36.3 6.0 where we see that models with a flexible parametric (spline-based) baseline hazard fit the data best followed by the standard parametric (Weibull, Gompertz, exponential) models. Roughly speaking, the B-spline and M-spline models seem to fit the data equally well since the differences in \codeelpd between the models are very small relative to their standard errors. Moreover, increasing the number of internal knots for the M-splines from two (the default) to five doesn’t seem to improve the fit (that is the default number of knots seems to provide sufficient flexibility for modelling the baseline hazard).

After fitting the survival model we often want to estimate the predicted survival function for individuals with different covariate patterns. Therefore, let us obtain the predicted survival function between 0 and 5 years for an individual in each of the prognostic groups. To do this we can use the \codeposterior_survfit() method for ‘\codestansurv’ objects and it’s associated \codeplot() method.

First let us construct the prediction (covariate) data: {Schunk} {Sinput} R> nd <- data.frame(group = c("Good", "Medium", "Poor")) R> head(nd) {Soutput} group 1 Good 2 Medium 3 Poor and then generate the posterior predictions: {Schunk} {Sinput} R> ps <- posterior_survfit(mod1, + newdata = nd, + times = 0, + extrapolate = TRUE, + control = list(edist = 5)) R> head(ps) {Soutput} stan_surv predictions num. individuals: 3 prediction type: event free probability standardised?: no conditional?: no

id cond_time time median ci_lb ci_ub 1 1 NA 0.0000 1.0000 1.0000 1.0000 2 1 NA 0.0505 0.9999 0.9993 1.0000 3 1 NA 0.1010 0.9996 0.9987 0.9999 4 1 NA 0.1515 0.9992 0.9979 0.9997 5 1 NA 0.2020 0.9987 0.9971 0.9995 6 1 NA 0.2525 0.9981 0.9960 0.9991

Here we note that the \codeid variable in the data frame of posterior predictions identifies which row of \codenewdata the predictions correspond to. For demonstration purposes we have also included other arguments in the \codeposterior_survfit() call, namely:

• •

  the \codetimes = 0 argument says that we want to predict at time $t=0$ (i.e. baseline) for each individual in the \codenewdata (this is the default anyway);

• •

  the \codeextrapolate = TRUE argument says that we want to extrapolate forward from time $t=0$ (this is also the default); and

• •

  the \codecontrol = list(edist = 5) identifies the control of the extrapolation; this is saying we wish to extrapolate the survival function forward from time $t=0$ for a distance of 5 time units (the default would have been to extrapolate as far as the largest event or censoring time in the estimation dataset, which is 7.28 years in the \codebcancer data).

Let us now plot the survival predictions. We will relabel the \codeid variable with meaningful labels identifying the covariate profile of each new individual in our prediction data: {Schunk}