Log hazard scale parameter.

Spline basis coefficients. These should sum to 1, otherwise they are normalised internally to sum to 1. Supplied either as a vector with one element per basis term, or a matrix with one column per basis term, and rows for alternative values of the coefficients (in vectorised usage of this function). If an array is supplied, it is collapsed into a matrix with number of columns equal to the final dimension of the array.

Locations of knots on the axis of time, supplied in increasing order. These include the two boundary knots.

In vectorised usage of these functions, the knots and degree must be the same for all alternative times and parameter values.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

logical; if TRUE (default), probabilities are P(X \le x), otherwise, P(X > x).

Probability of "cure", which defaults to zero. If this is non-zero, this defines a "mixture cure" version of the M-spline model.

Constant to be added to the cumulative hazard.

A data frame that defines the background hazard as a piecewise-constant function of time. This should have two columns, named "time" and "hazard". Each row gives the "background" hazard between the specified time and the next time. The first element of "time" should be 0, and the final row specifies the hazard at all times greater than the last element of "time".

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Vector of times.

Return log density or probability.

Constant to be added to the hazard, e.g. representing a "background" risk of death from causes other than the cause of interest.

Vector of probabilities.

Number of random numbers to simulate.

Discounting rate used to calculate the discounted mean or restricted mean survival time, using an exponential discounting function.

Method of numerical integration to obtain the restricted mean survival time from the survival function.

The default is method="gl", a Gauss-Legendre method with 100 nodes between zero and the maximum of t.

method="adaptive" uses the base R integrate function, which is much slower, but potentially more robust for badly-behaved survival functions.

Number of nodes for the Gauss-Legendre method.

Vector of times.

log hazard intercept, spline coefficients and cure probability before the start of the waning period ("treated")

log hazard intercept, spline coefficients and cure probability after the end of the waning period ("untreated")

Locations of knots on the axis of time, supplied in increasing order. These include the two boundary knots.

In vectorised usage of these functions, the knots and degree must be the same for all alternative times and parameter values.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

vector of two components giving start and stop of waning period

time resolution for piecewise constant hazard approximation in the waning period. If this is not specified, defaults to dividing the waning period into 10 pieces.

Constant to be added to the cumulative hazard.

A data frame that defines the background hazard as a piecewise-constant function of time. This should have two columns, named "time" and "hazard". Each row gives the "background" hazard between the specified time and the next time. The first element of "time" should be 0, and the final row specifies the hazard at all times greater than the last element of "time".

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Return log density or probability.

logical; if TRUE (default), probabilities are P(X \le x), otherwise, P(X > x).

Constant to be added to the hazard, e.g. representing a "background" risk of death from causes other than the cause of interest.

Vector of probabilities.

Number of random numbers to simulate.

Discounting rate used to calculate the discounted mean or restricted mean survival time, using an exponential discounting function.

Method of numerical integration to obtain the restricted mean survival time from the survival function.

The default is method="gl", a Gauss-Legendre method with 100 nodes between zero and the maximum of t.

method="adaptive" uses the base R integrate function, which is much slower, but potentially more robust for badly-behaved survival functions.

Number of nodes for the Gauss-Legendre method.

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then the MCMC samples are returned instead of being summarised as a median and 95% credible intervals.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

A time point, or vector of time points, to calculate the cumulative hazard at.

A data frame with two columns, named "time" and "hazard". Each row gives the "background" hazard between the specified time and the next time.

A vector of the same length as t indicating strata. There should then be

• a column in backhaz named stratum, which should be a factor with the same levels as strata (or analogous character vector).

• a row in backhaz indicating the background hazard at each time interval for each stratum in strata.

We assume that the background hazard is known at all times, and defined by a piecewise-constant function (step function) of time. The first element of "time" should be 0, and the final row specifies the hazard at all times greater than the last element of "time".

In additive hazards models, the background hazard is the hazard of death from causes other than the specific cause of interest. The overall hazard is defined as the sum of the background hazard and the cause-specific hazard. In survextrap, the cause-specific hazard is modelled with the M-spline model, and the background hazard is assumed to be known.

A fitted model object as returned by survextrap

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then the MCMC samples are returned instead of being summarised as a median and 95% credible intervals.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

A fitted model object as returned by survextrap

A data frame with two rows. The hazard ratio will be defined as hazard(second row) divided by hazard(first row). If the only covariate in the model is a factor with two levels, then newdata defaults to a data frame containing the levels of this factor, so that the hazard ratio for the second level versus the first level is computed. For any other models, newdata must be supplied explicitly.

Standardisation (with standardise_to) is not supported. This might be done by hand by using hazard(...,sample=TRUE) to obtain posterior samples for the two standardised hazards separately, then summarising by hand.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then the MCMC samples are returned instead of being summarised as a median and 95% credible intervals.

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

Quantiles which define the "low" and "high" values of a time-varying quantity (hazard in prior_haz_sd and the hazard ratio in prior_hr_sd). The ratio between the high and low values will be summarised, as a measure of time-dependence. By default, this is c(0.1, 0.9), so that the 10% and 90% quantiles are used respectively.

If TRUE then the MCMC samples are returned instead of being summarised as a median and 95% credible intervals.

A fitted model object as returned by survextrap

Vector of times. The restricted mean survival time up to each one of these times will be computed.

A data frame with two rows. The result will be the restricted mean for the covariates in the second row, minus the restricted mean for the covariates in the first row. If newdata is omitted for models where the only covariate is a factor with two levels, then this is taken from these levels. Otherwise newdata must be supplied explicitly.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then an MCMC sample is returned from the posterior of the output, rather than summary statistics.

Discounting rate used to calculate the discounted mean or restricted mean survival time, using an exponential discounting function.

Method of numerical integration to obtain the restricted mean survival time from the survival function.

The default is method="gl", a Gauss-Legendre method with 100 nodes between zero and the maximum of t.

method="adaptive" uses the base R integrate function, which is much slower, but potentially more robust for badly-behaved survival functions.

Number of nodes for the Gauss-Legendre method.

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

Discounting rate used to calculate the discounted mean or restricted mean survival time, using an exponential discounting function.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then an MCMC sample is returned from the posterior of the output, rather than summary statistics.

Other options (currently unused).

A numeric vector of times at which to evaluate the basis.

Spline knots

Spline degree

If TRUE, then the integrated M-spline (I-spline) basis is returned.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary, which improves smoothness at the boundary.

A list with components knots (vector of knots), degree (polynomial degree) and bsmooth (logical for smoothness constraint at boundary), defining an M-spline configuration.

If TRUE then the multinomial logit transform of the coefficients is returned. This is a vector of length one less than the number of coefficients, with the rth element defined by log(coefs[r+1] / coefs[1]).

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Location of the final spline knot.

Vector of observation times whose quantiles will be used to choose knot locations

A list with any or none of the following components: df, degree, bsmooth, knots, bknot, as documented in mspline_init.

Vector of observation times whose quantiles will be used to choose knot locations

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Location of the final spline knot.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Coefficients of the spline basis terms. These are normalised internally to sum to 1, if they do not already sum to 1.

Scale parameter. After computing the standard M-spline function as a weighted sum of the basis terms, the function is multiplied by scale. The log of the scale is the parameter called alpha in the results of a survextrap model, the intercept of the linear model on the log hazard.

Minimum plotting time. Defaults to zero.

Maximum plotting time. Defaults to the highest knot.

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Location of the final spline knot.

Minimum plotting time. Defaults to zero.

Maximum plotting time. Defaults to the highest knot.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

A survival formula in standard R formula syntax, with a call to Surv() on the left hand side.

Covariates included on the right hand side of the formula with be modelled with proportional hazards, or if nonprop is TRUE then a non-proportional hazards is used.

If data is omitted, so that the model is being fitted to external aggregate data alone, without individual data, then the formula should not include a Surv() call. The left-hand side of the formula will then be empty, and the right hand side specifies the covariates as usual. For example, formula = ~1 if there are no covariates.

Data frame containing variables in formula. Variables should be in a data frame, and not in the working environment.

This may be omitted, in which case external must be supplied. This allows a model to be fitted to external aggregate data alone, without any individual-level data.

If TRUE, a mixture cure model is used, where the "uncured" survival is defined by the M-spline model, and the cure probability is estimated.

Non-proportional hazards model specification. This is achieved by modelling the spline basis coefficients in terms of the covariates. See the methods vignette for more details.

If TRUE, then all covariates are modelled with non-proportional hazards, using the same model formula as formula.

If this is a formula, then this is assumed to define a model for the dependence of the basis coefficients on the covariates.

IF this is NULL or FALSE (the default) then any covariates are modelled with proportional hazards.

Background hazard, that is, for causes of death other than the cause of interest. This defines a "relative survival" or "additive hazards" model. The overall hazard that describes the all-cause survival data (given in the data and/or external argument) is then modelled as the sum of a cause-specific hazard and a background hazard.

The background hazard is assumed to be known, and the cause-specific hazard is modelled with the flexible parametric model.

The background hazard can be supplied in two forms. The meaning of predictions from the model depends on which of these is used.

(a) A data frame with columns "hazard" and "time", specifying the background hazard at all times as a piecewise-constant (step) function. Each row gives the background hazard between the specified time and the next time. The first element of "time" should be 0, and the final row specifies the hazard at all times greater than the last element of "time". Predictions from the model fitted by survextrap will then include this background hazard, because it is known at all times.

(b) The (quoted) name of a variable in the data giving the background hazard. For censored cases, the exact value does not matter. The predictions from survextrap will then describe the excess hazard or survival on top of this background. The overall hazard cannot be predicted in general, because the background hazard is only specified over a limited range of time.

If there is external data, and backhaz is supplied in form (b), then the user should also supply the background survival at the start and stop points in columns of the external data named "backsurv_start" and "backsurv_stop". That is, the probability of survival up to each of these times for someone alive at time 0. This should describe the same reference population as backhaz, though the package does not check for consistency between these.

If there are stratifying variables specified in backhaz_strata, then there should be multiple rows giving the background hazard value for each time period and stratifying variable.

If backhaz is NULL (the default) then no background hazard component is included in the model.

A character vector of names of variables that appear in backhaz that indicate strata, e.g. backhaz_strata = c("agegroup","sex"). This allows different background hazard values to be used for different subgroups. These variables must also appear in the datasets being modelled, that is, in data, external or both. Each row of those datasets should then have a corresponding row in backhaz which has the same values of the stratifying variables.

This is NULL by default, indicating no stratification of the background hazard.

If stratification is done, then backhaz must be supplied in form (a), as a data frame rather than a variable in the data.

External data as a data frame of aggregate survival counts with columns named:

start: Start time

stop: Follow-up time

n: Number of people alive at start

r: Number of those people who are still alive at stop

If there are covariates in formula, then the values they take in the external data must be supplied as additional columns in external. Therefore if there are external data, the covariates in formula and data should not be named start,stop,n or r.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Additional knots, other than those determined from the quantiles of the individual data. Typically used to add a maximum knot at the time that we want to extrapolate to.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Location of the final spline knot.

Vector of observation times whose quantiles will be used to choose knot locations

Basis coefficients

Hazard scale parameter

Best guess (prior median) for a typical hazard ratio

Upper limit of 95% credible interval for hazard ratio

Best guess (prior median) for a typical survival time

Upper limit of 95% credible interval for a survival time

A list with components knots (vector of knots), degree (polynomial degree) and bsmooth (logical for smoothness constraint at boundary), defining an M-spline configuration.

A fitted model object as returned by survextrap

"survival" for a plot of the survival function, "hazard" for the hazard function, against time.

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Additional arguments, passed on to plot_hazard and plot_survival.

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

If TRUE then credible intervals are drawn. Defaults to drawing the intervals if the plot shows the curve for only one covariate value.

X-axis label

Y-axis label

Passed to geom_line

Transparency for the credible interval ribbons

Show the locations of the spline knots as vertical lines

A fitted model object as returned by survextrap

A data frame with two rows. The hazard ratio will be defined as hazard(second row) divided by hazard(first row). If the only covariate in the model is a factor with two levels, then newdata defaults to a data frame containing the levels of this factor, so that the hazard ratio for the second level versus the first level is computed. For any other models, newdata must be supplied explicitly.

Standardisation (with standardise_to) is not supported. This might be done by hand by using hazard(...,sample=TRUE) to obtain posterior samples for the two standardised hazards separately, then summarising by hand.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

If TRUE then credible intervals are drawn. Defaults to drawing the intervals if the plot shows the curve for only one covariate value.

X-axis label

Y-axis label

Passed to geom_line

Transparency for the credible interval ribbons

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Location of the final spline knot.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Coefficients of the spline basis terms. These are normalised internally to sum to 1, if they do not already sum to 1.

Scale parameter. After computing the standard M-spline function as a weighted sum of the basis terms, the function is multiplied by scale. The log of the scale is the parameter called alpha in the results of a survextrap model, the intercept of the linear model on the log hazard.

Minimum plotting time. Defaults to zero.

Maximum plotting time. Defaults to the highest knot.

Show the positions of the knots, including the upper boundary

Show the "mean" around which each basis term is centred (defined as the mean of a random variable whose PDF is defined by the basis term).

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

If TRUE then a Kaplan-Meier curve of the observed data is plotted, using the results of survival::survfit() on the formula originally used for the survextrap fit. By default, this is only done when there are no covariates or one factor covariate.

The Kaplan-Meier estimates are returned in the km component of the fitted model object returned by survextrap, for use in hand-crafted plots like these.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

If TRUE then credible intervals are drawn. Defaults to drawing the intervals if the plot shows the curve for only one covariate value.

X-axis label

Y-axis label

Passed to geom_line

Transparency for the credible interval ribbons

Show the locations of the spline knots as vertical lines

A fitted model object as returned by survextrap

Other arguments (currently unused).

A fitted model object as returned by survextrap

A list of control parameters defining the spline model.

knots: Spline knots. If this is not supplied, then the number of knots is taken from df, and their location is taken from equally-spaced quantiles of the observed event times in the individual-level data.

add_knots: This is intended to be used when there are external data included in the model. External data are typically outside the time period covered by the individual data. add_knots would then be chosen to span the time period covered by the external data, so that the hazard trajectory can vary over that time.

If there are external data, and both knots and add_knots are omitted, then a default set of knots is chosen to span both the individual and external data, by taking the quantiles of a vector defined by concatenating the individual-level event times with the start and stop times in the external data.

df: Degrees of freedom, i.e. the number of parameters (or basis terms) intended to result from choosing knots based on quantiles of the data. The total number of parameters will then be df plus the number of additional knots specified in add_knots. df defaults to 10. This does not necessarily overfit, because the function is smoothed through the prior.

degree: Polynomial degree used for the basis function. The default is 3, giving a cubic. This can only be changed from 3 if bsmooth is FALSE.

bsmooth: If TRUE (on by default) the spline is smoother at the highest knot, by defining the derivative and second derivative at this point to be zero.

Spline basis coefficients that define the prior mean for the hazard function. By default, these are set to values that define a constant hazard function (see mspline_constant_coefs). They are normalised to sum to 1 internally (if they do not already).

Gamma prior for the standard deviation that controls the variability over time (or smoothness) of the hazard function. This should be a call to p_gamma(). The default is p_gamma(2,1). See prior_haz_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

The default "random_walk", specifies a random walk prior for the multinomial-logit spline coefficients, based on logistic distributions. See the methods vignette for full details.

The alternative "exchangeable" uses independent logistic priors on the multinomial-logit spline coefficients, conditionally on a common smoothing variance parameter. Note this is the method explained in the original survextrap paper (Jackson, BMC Med Res 2023). The random walk model is shown to perform better in Timmins et al (2025).

In non-proportional hazards models, setting smooth_model also determines whether an exchangeable or random walk model is used for the non-proportionality parameters (\delta).

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

A survival formula in standard R formula syntax, with a call to Surv() on the left hand side.

Covariates included on the right hand side of the formula with be modelled with proportional hazards, or if nonprop is TRUE then a non-proportional hazards is used.

If data is omitted, so that the model is being fitted to external aggregate data alone, without individual data, then the formula should not include a Surv() call. The left-hand side of the formula will then be empty, and the right hand side specifies the covariates as usual. For example, formula = ~1 if there are no covariates.

If TRUE, a mixture cure model is used, where the "uncured" survival is defined by the M-spline model, and the cure probability is estimated.

Non-proportional hazards model specification. This is achieved by modelling the spline basis coefficients in terms of the covariates. See the methods vignette for more details.

If TRUE, then all covariates are modelled with non-proportional hazards, using the same model formula as formula.

If this is a formula, then this is assumed to define a model for the dependence of the basis coefficients on the covariates.

IF this is NULL or FALSE (the default) then any covariates are modelled with proportional hazards.

A data frame with one row, containing variables in the model formulae. Samples will then be drawn, for any covariate-dependent parameters, with covariates set to the values given here.

Prior for the standard deviation parameters that smooth the non-proportionality effects over time in non-proportional hazards models. This should be a call to p_gamma() or a list of calls to p_gamma() with one component per covariate, as in prior_loghr. See prior_hr_sd for a way to calibrate this to represent a meaningful belief.

Minimum plotting time. Defaults to zero.

Maximum plotting time. Defaults to the highest knot.

Number of simulations to draw

Quantiles which define the "low" and "high" values of a time-varying quantity (hazard in prior_haz_sd and the hazard ratio in prior_hr_sd). The ratio between the high and low values will be summarised, as a measure of time-dependence. By default, this is c(0.1, 0.9), so that the 10% and 90% quantiles are used respectively.

Quantiles used to summarise the implied prior distributions of the simulated quantities.

A data frame with one row, containing "reference" values of variables in the model formulae. The hazard ratio between the hazards at newdata and newdata0 will be returned.

A list of control parameters defining the spline model.

knots: Spline knots. If this is not supplied, then the number of knots is taken from df, and their location is taken from equally-spaced quantiles of the observed event times in the individual-level data.

add_knots: This is intended to be used when there are external data included in the model. External data are typically outside the time period covered by the individual data. add_knots would then be chosen to span the time period covered by the external data, so that the hazard trajectory can vary over that time.

If there are external data, and both knots and add_knots are omitted, then a default set of knots is chosen to span both the individual and external data, by taking the quantiles of a vector defined by concatenating the individual-level event times with the start and stop times in the external data.

df: Degrees of freedom, i.e. the number of parameters (or basis terms) intended to result from choosing knots based on quantiles of the data. The total number of parameters will then be df plus the number of additional knots specified in add_knots. df defaults to 10. This does not necessarily overfit, because the function is smoothed through the prior.

degree: Polynomial degree used for the basis function. The default is 3, giving a cubic. This can only be changed from 3 if bsmooth is FALSE.

bsmooth: If TRUE (on by default) the spline is smoother at the highest knot, by defining the derivative and second derivative at this point to be zero.

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

Quantiles used to summarise the implied prior distributions of the simulated quantities.

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

Quantiles used to summarise the implied prior distributions of the simulated quantities.

Sample size of the simulated dataset. Each observation in the dataset is generated from a model with the same parameters. These parameters are generated from a single simulation from the prior distribution.

If TRUE, then one value of the parameter vector is drawn from the prior, followed by n individual-level times given this common prior value. If FALSE, then to produce each sampled individual time, a different sample from the prior is used.

A list of control parameters defining the spline model.

knots: Spline knots. If this is not supplied, then the number of knots is taken from df, and their location is taken from equally-spaced quantiles of the observed event times in the individual-level data.

add_knots: This is intended to be used when there are external data included in the model. External data are typically outside the time period covered by the individual data. add_knots would then be chosen to span the time period covered by the external data, so that the hazard trajectory can vary over that time.

If there are external data, and both knots and add_knots are omitted, then a default set of knots is chosen to span both the individual and external data, by taking the quantiles of a vector defined by concatenating the individual-level event times with the start and stop times in the external data.

df: Degrees of freedom, i.e. the number of parameters (or basis terms) intended to result from choosing knots based on quantiles of the data. The total number of parameters will then be df plus the number of additional knots specified in add_knots. df defaults to 10. This does not necessarily overfit, because the function is smoothed through the prior.

degree: Polynomial degree used for the basis function. The default is 3, giving a cubic. This can only be changed from 3 if bsmooth is FALSE.

bsmooth: If TRUE (on by default) the spline is smoother at the highest knot, by defining the derivative and second derivative at this point to be zero.

Right-censoring time to impose on the simulated event times.

Spline basis coefficients that define the prior mean for the hazard function. By default, these are set to values that define a constant hazard function (see mspline_constant_coefs). They are normalised to sum to 1 internally (if they do not already).

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

Gamma prior for the standard deviation that controls the variability over time (or smoothness) of the hazard function. This should be a call to p_gamma(). The default is p_gamma(2,1). See prior_haz_sd for a way to calibrate this to represent a meaningful belief.

A data frame with one row, containing variables in the model formulae. Samples will then be drawn, for any covariate-dependent parameters, with covariates set to the values given here.

A model formula with no response, defining the covariates on the hazard scale.

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

Prior for the standard deviation parameters that smooth the non-proportionality effects over time in non-proportional hazards models. This should be a call to p_gamma() or a list of calls to p_gamma() with one component per covariate, as in prior_loghr. See prior_hr_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline cure probability. This should be a call to p_beta(). The default is a uniform prior, p_beta(1,1). Baseline is defined by the mean of continuous covariates and the reference level of factor covariates.

A list of control parameters defining the spline model.

knots: Spline knots. If this is not supplied, then the number of knots is taken from df, and their location is taken from equally-spaced quantiles of the observed event times in the individual-level data.

add_knots: This is intended to be used when there are external data included in the model. External data are typically outside the time period covered by the individual data. add_knots would then be chosen to span the time period covered by the external data, so that the hazard trajectory can vary over that time.

If there are external data, and both knots and add_knots are omitted, then a default set of knots is chosen to span both the individual and external data, by taking the quantiles of a vector defined by concatenating the individual-level event times with the start and stop times in the external data.

df: Degrees of freedom, i.e. the number of parameters (or basis terms) intended to result from choosing knots based on quantiles of the data. The total number of parameters will then be df plus the number of additional knots specified in add_knots. df defaults to 10. This does not necessarily overfit, because the function is smoothed through the prior.

degree: Polynomial degree used for the basis function. The default is 3, giving a cubic. This can only be changed from 3 if bsmooth is FALSE.

bsmooth: If TRUE (on by default) the spline is smoother at the highest knot, by defining the derivative and second derivative at this point to be zero.

Spline basis coefficients that define the prior mean for the hazard function. By default, these are set to values that define a constant hazard function (see mspline_constant_coefs). They are normalised to sum to 1 internally (if they do not already).

Gamma prior for the standard deviation that controls the variability over time (or smoothness) of the hazard function. This should be a call to p_gamma(). The default is p_gamma(2,1). See prior_haz_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

The default "random_walk", specifies a random walk prior for the multinomial-logit spline coefficients, based on logistic distributions. See the methods vignette for full details.

The alternative "exchangeable" uses independent logistic priors on the multinomial-logit spline coefficients, conditionally on a common smoothing variance parameter. Note this is the method explained in the original survextrap paper (Jackson, BMC Med Res 2023). The random walk model is shown to perform better in Timmins et al (2025).

In non-proportional hazards models, setting smooth_model also determines whether an exchangeable or random walk model is used for the non-proportionality parameters (\delta).

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

A model formula with no response, defining the covariates on the hazard scale.

A model formula with no response, giving any covariates on the cure proportion.

A model formula with no response, defining any covariates affecting the spline basis coefficients, which gives a nonproportional hazards model.

A data frame with one row, containing variables in the model formulae. Samples will then be drawn, for any covariate-dependent parameters, with covariates set to the values given here.

A data frame with one row, containing "reference" values of variables in the model formulae. The hazard ratio between the hazards at newdata and newdata0 will be returned.

Prior for the standard deviation parameters that smooth the non-proportionality effects over time in non-proportional hazards models. This should be a call to p_gamma() or a list of calls to p_gamma() with one component per covariate, as in prior_loghr. See prior_hr_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline cure probability. This should be a call to p_beta(). The default is a uniform prior, p_beta(1,1). Baseline is defined by the mean of continuous covariates and the reference level of factor covariates.

Priors for log odds ratios on cure probabilities. This should be a call to p_normal() or p_t(). The default is p_normal(0,2.5).

Number of simulations to draw

Vector of knot locations. If not supplied, df has to be specified. One of two rules is then used to choose the knot locations. If bknot is specified, a set of equally spaced knots between zero and bknot is used. Otherwise if obstimes is supplied, the knots are chosen as equally spaced quantiles of obstimes.

The number of knots (excluding zero) is df - degree + 1 if bsmooth is TRUE, or df - degree - 1 otherwise.

Desired number of basis terms, or "degrees of freedom" in the spline. If knots is not supplied, the number of knots is then chosen to satisfy this.

Spline polynomial degree. Can only be changed from the default of 3 if bsmooth is FALSE.

If TRUE then the function is constrained to also have zero derivative and second derivative at the boundary.

Spline basis coefficients that define the prior mean for the hazard function. By default, these are set to values that define a constant hazard function (see mspline_constant_coefs). They are normalised to sum to 1 internally (if they do not already).

Gamma prior for the standard deviation that controls the variability over time (or smoothness) of the hazard function. This should be a call to p_gamma(). The default is p_gamma(2,1). See prior_haz_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

The default "random_walk", specifies a random walk prior for the multinomial-logit spline coefficients, based on logistic distributions. See the methods vignette for full details.

The alternative "exchangeable" uses independent logistic priors on the multinomial-logit spline coefficients, conditionally on a common smoothing variance parameter. Note this is the method explained in the original survextrap paper (Jackson, BMC Med Res 2023). The random walk model is shown to perform better in Timmins et al (2025).

In non-proportional hazards models, setting smooth_model also determines whether an exchangeable or random walk model is used for the non-proportionality parameters (\delta).

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

A survival formula in standard R formula syntax, with a call to Surv() on the left hand side.

Covariates included on the right hand side of the formula with be modelled with proportional hazards, or if nonprop is TRUE then a non-proportional hazards is used.

If data is omitted, so that the model is being fitted to external aggregate data alone, without individual data, then the formula should not include a Surv() call. The left-hand side of the formula will then be empty, and the right hand side specifies the covariates as usual. For example, formula = ~1 if there are no covariates.

If TRUE, a mixture cure model is used, where the "uncured" survival is defined by the M-spline model, and the cure probability is estimated.

Non-proportional hazards model specification. This is achieved by modelling the spline basis coefficients in terms of the covariates. See the methods vignette for more details.

If TRUE, then all covariates are modelled with non-proportional hazards, using the same model formula as formula.

If this is a formula, then this is assumed to define a model for the dependence of the basis coefficients on the covariates.

IF this is NULL or FALSE (the default) then any covariates are modelled with proportional hazards.

A data frame with one row, containing variables in the model formulae. Samples will then be drawn, for any covariate-dependent parameters, with covariates set to the values given here.

A data frame with one row, containing "reference" values of variables in the model formulae. The hazard ratio between the hazards at newdata and newdata0 will be returned.

Prior for the standard deviation parameters that smooth the non-proportionality effects over time in non-proportional hazards models. This should be a call to p_gamma() or a list of calls to p_gamma() with one component per covariate, as in prior_loghr. See prior_hr_sd for a way to calibrate this to represent a meaningful belief.

Minimum plotting time. Defaults to zero.

Maximum plotting time. Defaults to the highest knot.

Number of simulations to draw

Prior location. For the normal distribution, this is the mean. Defaults to 0

Prior scale. For the normal distribution, this is the standard deviation. Defaults to 2.5.

Prior degrees of freedom (only for Student t distribution).

First shape parameter (for Beta distribution, defaults to 1).

Second shape parameter (for Beta distribution, defaults to 1).

Shape parameter (for Gamma distribution, defaults to 2).

Rate parameter (for Gamma distribution, defaults to 1).

A fitted model object as returned by survextrap

Vector of times. The restricted mean survival time up to each one of these times will be computed.

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.

Discounting rate used to calculate the discounted mean or restricted mean survival time, using an exponential discounting function.

Method of numerical integration to obtain the restricted mean survival time from the survival function.

The default is method="gl", a Gauss-Legendre method with 100 nodes between zero and the maximum of t.

method="adaptive" uses the base R integrate function, which is much slower, but potentially more robust for badly-behaved survival functions.

Number of nodes for the Gauss-Legendre method.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then an MCMC sample is returned from the posterior of the output, rather than summary statistics.

Data frame describing a population.

Number of draws from the population distribution used per MCMC sample from the parameters when random=TRUE. With the default of 1, the value of the covariate vector X is essentially treated as if it were an additional parameter in the Bayesian model, drawn by Monte Carlo independently of the remaining parameters.

By default this is FALSE, indicating that standardised samples should be obtained by concatenating the posterior samples for each covariate value in the standard population. The sample from the standardised posterior of parameters then has size niter times the number of rows in newdata, where niter is the number of MCMC iterations used in the original survextrap fit. Computing the resulting output function (e.g. RMST which uses numerical integration) can then be computationally intensive if this sample size is large.

A quicker alternative is to sample a random row of the standard population for each MCMC iteration. The standardised sample from the posterior then has size niter. This is specified by using random=TRUE. If this is used, then the result depends on the random number seed, and it should be checked that the results are stable to within the required number of significant figures. If not, run survextrap with more MCMC iterations or increase nstd here.

A fitted model object as returned by survextrap

A list of functions to calculate different posterior summaries from the MCMC sample. This is passed to posterior::summarise_draws. If the list is named, then the names will be used for the columns of the output.

See the examples below for different ways this can be used.

Defaults to list(median = median, ~quantile(.x, probs=c(0.025, 0.975)), sd = sd, rhat = posterior::rhat, ess_bulk = posterior::ess_bulk)

Many useful such functions are provided with the posterior package.

Summary functions can also be supplied in separate arguments here. They will then be added to those supplied in summ_fns.

A survival formula in standard R formula syntax, with a call to Surv() on the left hand side.

Covariates included on the right hand side of the formula with be modelled with proportional hazards, or if nonprop is TRUE then a non-proportional hazards is used.

If data is omitted, so that the model is being fitted to external aggregate data alone, without individual data, then the formula should not include a Surv() call. The left-hand side of the formula will then be empty, and the right hand side specifies the covariates as usual. For example, formula = ~1 if there are no covariates.

Data frame containing variables in formula. Variables should be in a data frame, and not in the working environment.

This may be omitted, in which case external must be supplied. This allows a model to be fitted to external aggregate data alone, without any individual-level data.

External data as a data frame of aggregate survival counts with columns named:

start: Start time

stop: Follow-up time

n: Number of people alive at start

r: Number of those people who are still alive at stop

If there are covariates in formula, then the values they take in the external data must be supplied as additional columns in external. Therefore if there are external data, the covariates in formula and data should not be named start,stop,n or r.

If TRUE, a mixture cure model is used, where the "uncured" survival is defined by the M-spline model, and the cure probability is estimated.

Non-proportional hazards model specification. This is achieved by modelling the spline basis coefficients in terms of the covariates. See the methods vignette for more details.

If TRUE, then all covariates are modelled with non-proportional hazards, using the same model formula as formula.

If this is a formula, then this is assumed to define a model for the dependence of the basis coefficients on the covariates.

IF this is NULL or FALSE (the default) then any covariates are modelled with proportional hazards.

Prior for the baseline log hazard scale parameter (alpha or log(eta)). This should be a call to a prior constructor function, such as p_normal(0,1) or p_t(0,2,2). Supported prior distribution families are normal (parameters mean and SD) and t distributions (parameters location, scale and degrees of freedom). The default is a normal distribution with mean 0 and standard deviation 20.

Note that eta is not in itself a hazard, but it is proportional to the hazard (see the vignette for the full model specification).

"Baseline" is defined by the continuous covariates taking a value of zero and factor covariates taking their reference level. To use a different baseline, the data should be transformed appropriately beforehand, so that a value of zero has a different meaning. For continuous covariates, it helps for both computation and interpretation to define the value of zero to denote a typical value in the data, e.g. the mean.

Priors for log hazard ratios. This should be a call to p_normal() or p_t(). A list of calls can also be provided, to give different priors to different coefficients, where the name of each list component matches the name of the coefficient, e.g. list("age45-59" = p_normal(0,1), "age60+" = p_t(0,2,3))

The default is p_normal(0,2.5) for all coefficients.

Gamma prior for the standard deviation that controls the variability over time (or smoothness) of the hazard function. This should be a call to p_gamma(). The default is p_gamma(2,1). See prior_haz_sd for a way to calibrate this to represent a meaningful belief.

Prior for the baseline cure probability. This should be a call to p_beta(). The default is a uniform prior, p_beta(1,1). Baseline is defined by the mean of continuous covariates and the reference level of factor covariates.

Priors for log odds ratios on cure probabilities. This should be a call to p_normal() or p_t(). The default is p_normal(0,2.5).

Prior for the standard deviation parameters that smooth the non-proportionality effects over time in non-proportional hazards models. This should be a call to p_gamma() or a list of calls to p_gamma() with one component per covariate, as in prior_loghr. See prior_hr_sd for a way to calibrate this to represent a meaningful belief.

Background hazard, that is, for causes of death other than the cause of interest. This defines a "relative survival" or "additive hazards" model. The overall hazard that describes the all-cause survival data (given in the data and/or external argument) is then modelled as the sum of a cause-specific hazard and a background hazard.

The background hazard is assumed to be known, and the cause-specific hazard is modelled with the flexible parametric model.

The background hazard can be supplied in two forms. The meaning of predictions from the model depends on which of these is used.

(a) A data frame with columns "hazard" and "time", specifying the background hazard at all times as a piecewise-constant (step) function. Each row gives the background hazard between the specified time and the next time. The first element of "time" should be 0, and the final row specifies the hazard at all times greater than the last element of "time". Predictions from the model fitted by survextrap will then include this background hazard, because it is known at all times.

(b) The (quoted) name of a variable in the data giving the background hazard. For censored cases, the exact value does not matter. The predictions from survextrap will then describe the excess hazard or survival on top of this background. The overall hazard cannot be predicted in general, because the background hazard is only specified over a limited range of time.

If there is external data, and backhaz is supplied in form (b), then the user should also supply the background survival at the start and stop points in columns of the external data named "backsurv_start" and "backsurv_stop". That is, the probability of survival up to each of these times for someone alive at time 0. This should describe the same reference population as backhaz, though the package does not check for consistency between these.

If there are stratifying variables specified in backhaz_strata, then there should be multiple rows giving the background hazard value for each time period and stratifying variable.

If backhaz is NULL (the default) then no background hazard component is included in the model.

A character vector of names of variables that appear in backhaz that indicate strata, e.g. backhaz_strata = c("agegroup","sex"). This allows different background hazard values to be used for different subgroups. These variables must also appear in the datasets being modelled, that is, in data, external or both. Each row of those datasets should then have a corresponding row in backhaz which has the same values of the stratifying variables.

This is NULL by default, indicating no stratification of the background hazard.

If stratification is done, then backhaz must be supplied in form (a), as a data frame rather than a variable in the data.

A list of control parameters defining the spline model.

knots: Spline knots. If this is not supplied, then the number of knots is taken from df, and their location is taken from equally-spaced quantiles of the observed event times in the individual-level data.

add_knots: This is intended to be used when there are external data included in the model. External data are typically outside the time period covered by the individual data. add_knots would then be chosen to span the time period covered by the external data, so that the hazard trajectory can vary over that time.

If there are external data, and both knots and add_knots are omitted, then a default set of knots is chosen to span both the individual and external data, by taking the quantiles of a vector defined by concatenating the individual-level event times with the start and stop times in the external data.

df: Degrees of freedom, i.e. the number of parameters (or basis terms) intended to result from choosing knots based on quantiles of the data. The total number of parameters will then be df plus the number of additional knots specified in add_knots. df defaults to 10. This does not necessarily overfit, because the function is smoothed through the prior.

degree: Polynomial degree used for the basis function. The default is 3, giving a cubic. This can only be changed from 3 if bsmooth is FALSE.

bsmooth: If TRUE (on by default) the spline is smoother at the highest knot, by defining the derivative and second derivative at this point to be zero.

Any extra knots beyond those chosen from the individual-level data (or supplied in knots). All other spline specifications are set to their defaults, as described in mspline. For example, add_knots = 10 is a shorthand for mspline = list(add_knots = 10).

The default "random_walk", specifies a random walk prior for the multinomial-logit spline coefficients, based on logistic distributions. See the methods vignette for full details.

The alternative "exchangeable" uses independent logistic priors on the multinomial-logit spline coefficients, conditionally on a common smoothing variance parameter. Note this is the method explained in the original survextrap paper (Jackson, BMC Med Res 2023). The random walk model is shown to perform better in Timmins et al (2025).

In non-proportional hazards models, setting smooth_model also determines whether an exchangeable or random walk model is used for the non-proportionality parameters (\delta).

Smoothing variance parameter estimation.

"bayes": the smoothing parameter is estimated by full Bayes (the default).

"eb": empirical Bayes is used.

Alternatively, if a number is supplied here, then the smoothing parameter is fixed to this number.

Spline basis coefficients that define the prior mean for the hazard function. By default, these are set to values that define a constant hazard function (see mspline_constant_coefs). They are normalised to sum to 1 internally (if they do not already).

Method from rstan used to fit the model.

If fit_method="mcmc" then a sample from the posterior is drawn using Markov Chain Monte Carlo sampling, via rstan's rstan::sampling() function. This is the default. It is the most accurate but the slowest method.

If fit_method="opt", then instead of an MCMC sample from the posterior, survextrap returns the posterior mode calculated using optimisation, via rstan's rstan::optimizing() function. A sample from a normal approximation to the (real-line-transformed) posterior distribution is drawn in order to obtain credible intervals. This is useful for model development, while using MCMC for the "final answer".

If fit_method="vb", then variational Bayes methods are used, via rstan's rstan::vb() function. This is labelled as "experimental" by rstan. It might give a better approximation to the posterior than fit_method="opt", and is faster than MCMC, in particular for large datasets, but has not been investigated in depth for these models.

Compute leave-one-out cross-validation statistics. This is done by default. Set to FALSE to not compute them. If these statistics are computed, then they are returned in the loo and loo_external components of the object returned by survextrap. loo describes the fit of the model to the individual-level data, and loo_external describes fit to the external data.

See the "examples" vignette for more explanation of these.

Additional arguments to supply to control the Stan fit, passed to the appropriate rstan function, depending on which is chosen through the fit_method argument.

A fitted model object as returned by survextrap

Data frame of covariate values to compute the output for. If there are covariates in the model and this is not supplied, the following default is used:

(a) if the only covariate is one factor variable, then the output is computed for each level of this factor.

(b) if there are multiple covariates, or any numeric covariates, then the output is computed at the mean of each numeric covariate in the original data, and at the baseline level of each factor covariate.

Note: caution is required about how treatment groups (for example) are stored in your data. If these are coded as numeric (0/1), then if newdata is not specified, only one output will be shown. This relates to the average value of this numeric variable over the data, which doesn't correspond to either of the treatment groups. To avoid this, a treatment group should be stored as a factor.

Vector of times at which to compute the estimates.

Maximum time at which to compute the estimates. If t is supplied, then this is ignored. If t is not supplied, then t is set to a set of 100 equally spaced time points from 0 to tmax. If both tmax and t are not supplied, then tmax is set to the maximum follow up time in the data.

Number of MCMC iterations to use to compute credible intervals. Set to a low value to make this function quicker, at the cost of some approximation error (which may not be important for plotting or model development).

A list of functions to use to summarise the posterior sample. This is passed to posterior::summarise_draws. By default this is list(median=median, ~quantile(.x, probs=c(0.025, 0.975))). If the list is named, then the names will be used for the columns of the output.

If TRUE then the MCMC samples are returned instead of being summarised as a median and 95% credible intervals.

Data frame of covariate values defining the "untreated" group for use in treatment waning models. See Survmspline_wane.

Vector of two numbers, defining the time period over which the hazard is interpolated between the hazard of the "treated" group (taken from newdata) and the hazard of the "untreated" group (taken from newdata0). Optional - if this is not supplied, then no waning is assumed.

Number of intervals defining the piecewise constant approximation to the hazard during the waning period.