Package {rpsurv}

rpsurv: Fast Royston-Parmar Flexible Parametric Survival Models

Description

Fits Royston-Parmar flexible parametric survival models with the likelihood and gradient evaluated in parallel C++ ('RcppParallel').

Author(s)

Maintainer: Imad El Badisy elbadisyimad@gmail.com

See Also

Useful links:

• https://github.com/ielbadisy/rpsurv

• Report bugs at https://github.com/ielbadisy/rpsurv/issues

German breast cancer data

Description

Recurrence-free survival for 686 women in a randomized trial of hormonal therapy for breast cancer (Sauerbrei and Royston), used throughout the package's documentation and vignette as the baseline rpsurv() fitting example. Identical to rstpm2::brcancer, bundled here so examples do not require rstpm2 to be installed.

Usage

brcancer

Format

A data frame with 686 rows and 14 variables:

id

subject identifier

hormon

hormonal therapy (0 = no, 1 = yes)

x1

age, years

x2

menopausal status

x3

tumour size, mm

x4

tumour grade

x5

number of positive lymph nodes

x6

progesterone receptor, fmol

x7

estrogen receptor, fmol

rectime

recurrence-free survival time, days

censrec

event indicator (1 = recurrence/death, 0 = censored)

x4a

tumour grade >= 2

x4b

tumour grade == 3

x5e

exp(-0.12 * x5)

Source

https://www.stata-press.com/data/r11/brcancer.dta, as redistributed by the rstpm2 package.

Cox-Snell residual diagnostic plot

Description

Plots the Nelson-Aalen cumulative hazard of the Cox-Snell residuals against the residuals themselves. Under a correctly specified model the residuals are a censored sample from a unit exponential distribution, so the points should lie close to the y = x reference line. plot(fit, type = "coxsnell") draws the same plot.

Usage

coxsnell_plot(object, ...)

Arguments

Details

For counting-process data (Surv(start, stop, status)), each row's residual is H(stop) - H(start). Given survival to start, it is unit exponential and censored at its observed value when the row has no event, so the same check applies row by row. Tied residuals are handled by the Nelson-Aalen estimator.

Value

Called for its side effect of drawing the diagnostic plot. Invisibly returns a data.frame with columns time (the sorted distinct Cox-Snell residuals) and cumhaz (their Nelson-Aalen cumulative hazard), the coordinates of the plotted step function.

Default knot placement for a Royston-Parmar spline

Description

Places boundary knots at the min/max of x and interior knots at equally spaced centiles of x restricted to event times, matching the default behaviour of rstpm2::stpm2() / flexsurv::flexsurvspline().

Usage

default_knots(x, df)

Arguments

Details

When event times are heavily concentrated near one end of the range (e.g. an accelerating hazard piling events up before administrative censoring), a centile can coincide with a boundary knot. An interior knot placed exactly at a boundary knot makes that spline basis column identically zero (see rcs_basis()), which is a silent rank deficiency: such knots are dropped, reducing the realised spline df below df (with a warning) rather than fitting a structurally unidentified model.

Stanford heart transplant data

Description

Survival of patients on the waiting list for the Stanford heart transplant program (Crowley and Hu, 1977), already in counting-process (start, stop, event) format with a time-varying transplant covariate: a patient's row is split at the moment of transplantation, if it occurred. This makes it a natural, real (rather than constructed) illustration of rpsurv()'s support for genuine time-varying covariates via Surv(start, stop, status). Identical to survival::heart, bundled here under the same name for convenience.

Usage

heart

Format

A data frame with 172 rows and 8 variables:

start

interval start time, days

stop

interval stop time, days

event

event indicator at stop (1 = death, 0 = censored)

age

age minus 48 years

year

year of acceptance into the program, in years since 1967

surgery

prior bypass surgery (0 = no, 1 = yes)

transplant

transplant status during this interval (0 = no, 1 = yes)

id

patient identifier

Source

Crowley, J. and Hu, M. (1977) Covariance analysis of heart transplant survival data. Journal of the American Statistical Association, 72, 27-36, as redistributed by the survival package.

Compare the fitted survival curve against the Kaplan-Meier estimate

Description

A calibration / goodness-of-fit diagnostic: overlays the model-based survival curve on the nonparametric Kaplan-Meier estimate, overall or by groups of a covariate. Close agreement supports the chosen spline df and scale; systematic divergence suggests more baseline df, a different scale, or a tve term is needed. plot(fit, type = "km") draws the same plot.

Usage

km_compare_plot(
  object,
  by = NULL,
  col = NULL,
  n_grid = 200,
  legend = TRUE,
  ...
)

Arguments

Details

The Kaplan-Meier curve estimates the marginal survival of the group, so the model curve is the matching population-averaged (direct-adjusted) estimand, not the survival curve at the mean covariate profile (which is a different, generally higher or lower, curve when covariates vary):

• Right-censored data: the mean over the group's subjects of their predicted survival curves, ⁠mean_i S(t | x_i)⁠. Kaplan-Meier estimates the same curve when censoring does not depend on the covariates; when it does, the two can differ slightly even for a correct model.

• Counting-process data (Surv(start, stop, status), left truncation or time-varying covariates): the Kaplan-Meier curve is then a product of conditional survival over the rows at risk, so the model curve is exp(-A(t)), where A cumulates the average model hazard of the rows at risk at each time (the model counterpart of the Nelson-Aalen estimator, computed on a fine time grid). A mean of S(t | x_i) from time 0 would not be comparable, because subjects enter the sample only after surviving to their entry time.

Value

Called for its side effect of drawing the comparison plot. Invisibly returns a data.frame with columns group, curve ("km" or "model"), time and surv, the plotted coordinates.

Examples

fit <- rpsurv(survival::Surv(time, status) ~ trt + celltype, data = veteran)
km_compare_plot(fit, by = "celltype")

Non-linear effect component of a covariate

Description

Returns the fitted non-linear effect g(x) - g(ref) of a covariate named in nle, on the linear-predictor scale, with a delta-method 95% confidence interval. Because the spline enters the linear predictor additively and does not depend on time, this is an exact linear combination of the coefficients. On the hazard scale, and when var is not also in tve, it is the log hazard ratio of x versus ref at every time. When var is also in tve, it is the time-constant component of the effect, separate from the time-varying component returned by tvecurve().

Usage

nlecurve(object, var, at = NULL, ref = NULL, se.fit = TRUE)

Arguments

Value

a data frame of class "nlecurve" with columns x and est, and optionally se, lower, and upper; see plot.nlecurve().

Examples

set.seed(1)
n <- 500
x <- runif(n, 0, 3)
dat <- data.frame(time = rexp(n, 0.2 * exp((x - 1.5)^2 / 2)),
                  status = 1, x = x)
fit <- rpsurv(survival::Surv(time, status) ~ x, data = dat, nle = "x")
g <- nlecurve(fit, "x", ref = 1.5)
head(g)
plot(g)

Plot a non-linear effect

Description

Draws the output of nlecurve(): g(x) - g(ref) against the covariate, with its 95% confidence band and a reference line at no effect.

Usage

## S3 method for class 'nlecurve'
plot(
  x,
  exp = !identical(attr(x, "scale"), "normal"),
  ci = TRUE,
  log = if (exp) "y" else "",
  col = 1,
  ref.line = TRUE,
  add = FALSE,
  xlab = NULL,
  ylab = NULL,
  rug = TRUE,
  ...
)

Arguments

Value

Called for its side effect. Invisibly returns x.

Examples

fit <- rpsurv(survival::Surv(time, status) ~ trt + karno, data = veteran,
              nle = "karno")
plot(nlecurve(fit, "karno", ref = 50))

Plot a Royston-Parmar model

Description

Draws predicted survival, hazard or cumulative hazard curves for covariate profiles, or one of the two goodness-of-fit diagnostics: type = "coxsnell" calls coxsnell_plot() and type = "km" calls km_compare_plot().

Usage

## S3 method for class 'rpsurv'
plot(
  x,
  newdata = NULL,
  type = c("survival", "hazard", "cumhaz", "coxsnell", "km"),
  times = NULL,
  ci = TRUE,
  col = NULL,
  legend = NULL,
  by = NULL,
  ...
)

Arguments

Value

Called for its side effect of drawing the plot. For the diagnostics, invisibly returns the value of coxsnell_plot() or km_compare_plot(). For predicted curves, invisibly returns the data.frame of predictions produced by predict.rpsurv(), with columns id (row of newdata), time, est, and, when ci = TRUE, lower and upper for the 95% confidence band. These are exactly the plotted values.

Examples

fit <- rpsurv(survival::Surv(time, status) ~ trt + celltype, data = veteran)
plot(fit, newdata = data.frame(trt = 1, celltype = c("squamous", "adeno")))
plot(fit, type = "coxsnell")
plot(fit, type = "km", by = "celltype")

Plot a time-varying effect

Description

Draws the output of tvecurve(): the effect against time, with its 95% confidence band and a reference line at no effect.

Usage

## S3 method for class 'tvecurve'
plot(
  x,
  exp = attr(x, "type") != "link",
  ci = TRUE,
  log = if (exp) "y" else "",
  col = 1,
  ref.line = TRUE,
  add = FALSE,
  xlab = "Time",
  ylab = NULL,
  ...
)

Arguments

Details

The default y label states what is drawn. For type = "hr" it is the hazard ratio of var = ref + 1 versus var = ref (other covariates at 0). It is not a ratio "per unit" of var: with a time-varying effect the derivative of the linear predictor depends on var, so the hazard ratio changes with ref, and off the hazard scale it also depends on the other covariates. For type = "link" it is beta(t), the effect per unit of var on the model scale: a log cumulative hazard ratio (scale = "hazard"), a log cumulative odds ratio ("odds"), or a shift of the probit of the failure probability ("normal"); exponentiated, the first two are a cumulative hazard ratio and a cumulative odds ratio.

Value

Called for its side effect. Invisibly returns x.

Examples

fit <- rpsurv(survival::Surv(time, status) ~ trt + karno, data = veteran,
              tve = "karno")
plot(tvecurve(fit, "karno", ref = 50))

Predict from a fitted Royston-Parmar model

Description

Predict from a fitted Royston-Parmar model

Usage

## S3 method for class 'rpsurv'
predict(
  object,
  newdata = NULL,
  newdata0 = NULL,
  times = NULL,
  type = c("survival", "hazard", "cumhaz", "link", "hr", "sdiff"),
  se.fit = FALSE,
  ...
)

Arguments

Value

a long-format data frame with columns id, time, est, and optionally lower/upper.

Examples

dat <- data.frame(
  time = rexp(200, 0.2), status = rbinom(200, 1, 0.7), x1 = rbinom(200, 1, 0.5)
)
fit <- rpsurv(survival::Surv(time, status) ~ x1, data = dat, tve = "x1")
tt <- seq(0.5, 4, length.out = 20)

hr <- predict(fit,
  newdata = data.frame(x1 = 1), newdata0 = data.frame(x1 = 0),
  times = tt, type = "hr", se.fit = TRUE
)
head(hr)

sd <- predict(fit,
  newdata = data.frame(x1 = 1), newdata0 = data.frame(x1 = 0),
  times = tt, type = "sdiff", se.fit = TRUE
)
head(sd)

Restricted cubic spline basis (Durrleman-Simon / Royston-Parmar parameterisation)

Description

Builds the basis used by Royston & Parmar (2002) flexible parametric models: a linear term plus one term per interior knot, each term being a natural (restricted) cubic spline component that is linear beyond the boundary knots.

Usage

rcs_basis(x, knots, derivative = FALSE)

Arguments

Value

a matrix with length(knots) - 1 columns.

Restricted cubic spline basis, computed in C++

Description

Same basis as rcs_basis() (Durrleman-Simon / Royston-Parmar parameterisation): a linear term plus one term per interior knot.

Usage

rcs_basis_cpp(x, knots, derivative)

Arguments

Value

an ⁠n x (length(knots) - 1)⁠ matrix.

Residuals for a Royston-Parmar model

Description

Residuals for a Royston-Parmar model

Usage

## S3 method for class 'rpsurv'
residuals(object, type = c("coxsnell", "martingale", "deviance"), ...)

Arguments

Value

numeric vector of residuals, one per observation used to fit the model.

Restricted mean survival time difference between two covariate profiles

Description

Integrates the survival difference sdiff = S(t | newdata) - S(t | newdata0) from 0 to tau by the trapezoidal rule, with an optional delta-method standard error and 95% confidence interval.

Usage

rmst_diff(object, newdata, newdata0, tau, se.fit = FALSE, n_grid = 200)

Arguments

Value

a one-row data frame with columns est and, if se.fit = TRUE, se, lower, upper.

Examples

dat <- data.frame(
  time = rexp(200, 0.2), status = rbinom(200, 1, 0.7), x1 = rbinom(200, 1, 0.5)
)
fit <- rpsurv(survival::Surv(time, status) ~ x1, data = dat, tve = "x1")
rmst_diff(fit, newdata = data.frame(x1 = 1), newdata0 = data.frame(x1 = 0),
          tau = 4, se.fit = TRUE)

Build the design matrices for a Royston-Parmar model

Description

Constructs eta's design matrix X and the design matrix for ⁠d(eta)/d(log t)⁠, dX, which share the same column layout: ⁠[intercept, baseline spline(log t), covariates, tve spline(log t):covariate]⁠. A covariate named in nle enters the covariate block as its restricted cubic spline basis in x (first column linear) instead of a single linear column; it does not depend on time, so its dX columns are zero.

Usage

rp_design(
  log_time,
  cov_data,
  knots,
  tve = NULL,
  tve_knots = NULL,
  nle = NULL,
  nle_knots = NULL
)

Negative log-likelihood and gradient for a Royston-Parmar model

Description

Supports left truncation / counting-process data (genuine time-varying covariates): each row contributes logS(stop) - logS(entry) to the log-likelihood, Xentry and hasEntry giving the design row and flag for the entry (left-truncation) time. Pass hasEntry all-zero (and Xentry a matching-shape dummy matrix) for standard right-censored data.

Usage

rp_negloglik_grad_cpp(beta, X, dX, Xentry, hasEntry, logtime, status, scale)

Arguments

Value

a list with value (negative log-likelihood) and gradient.

Fit a Royston-Parmar flexible parametric survival model

Description

Fit a Royston-Parmar flexible parametric survival model

Usage

rpsurv(
  formula,
  data,
  df = 4,
  knots = NULL,
  tve = NULL,
  tve.df = 3,
  nle = NULL,
  nle.df = 3,
  scale = c("hazard", "odds", "normal"),
  control = list()
)

Arguments

Value

an object of class "rpsurv".

Time-varying effect vs. time-varying covariate

These are two distinct extensions and rpsurv() supports both, separately or combined:

• Time-varying effect (non-proportional hazards): a covariate's value is fixed for a subject but its association with the outcome changes over time, e.g. a treatment effect that fades. Request this with tve.

• Time-varying covariate: a covariate's value itself changes during follow-up, e.g. a lab measurement updated at clinic visits. Represent this in the data as counting-process (start, stop] intervals, one row per interval with the covariate value held constant within it, and fit with Surv(start, stop, status) ~ .... rpsurv() handles the corresponding left-truncated likelihood automatically.

Non-linear effect vs. time-varying effect

tve lets a covariate's effect change with time; nle lets it change with the covariate's own value. Each has its own argument and spline degrees of freedom, so the two components are fitted, reported, and extracted separately: tvecurve() returns beta(t) and nlecurve() returns g(x) - g(ref).

Estimation

The log-likelihood and its analytic gradient are maximised with BFGS over a linear reparameterisation of the coefficients in which the design matrix has orthogonal, unit-scale columns (from its QR decomposition). This keeps the problem well conditioned when covariates on large scales multiply the time spline basis in tve terms. The Hessian is obtained by central differences of the analytic gradient in that parameterisation, and the coefficients and their variance-covariance matrix are mapped back exactly to the scale of the original design.

Time-varying effect component of a covariate

Description

Returns the fitted time-varying effect of a covariate named in tve, with a delta-method 95% confidence interval.

Usage

tvecurve(
  object,
  var,
  times = NULL,
  type = c("hr", "link"),
  ref = 0,
  se.fit = TRUE
)

Arguments

Details

• type = "hr" (default): the log of the instantaneous hazard ratio at each time for var = ref + 1 versus var = ref, as in predict.rpsurv() with type = "hr". Other covariates are set to 0; on the hazard scale they cancel unless they also have a time-varying effect. If var is also in nle, its non-linear component is held fixed, so the result is the time-varying component only (see nlecurve() for the other one).

• type = "link": the coefficient function beta(t) on the linear-predictor (log cumulative hazard, log cumulative odds, or probit) scale, the main effect plus its spline in log time; when var is also in nle the main effect belongs to g, so beta(t) is the spline part only.

Value

a data frame of class "tvecurve" with columns time and est, and optionally se, lower, and upper; see plot.tvecurve().

Examples

fit <- rpsurv(survival::Surv(time, status) ~ trt + karno, data = veteran,
              tve = "karno", nle = "karno")
tvecurve(fit, "karno", times = c(30, 90, 180), ref = 50)
plot(tvecurve(fit, "karno", ref = 50))

Veterans' Administration lung cancer trial

Description

A randomized trial comparing two treatment regimens for lung cancer (Kalbfleisch and Prentice). celltype is a well-known example of a non-proportional-hazards effect: its hazard ratio changes materially over follow-up, which makes this dataset a natural illustration of rpsurv()'s tve (time-varying effect) argument. Identical to survival::veteran, bundled here under the same name for convenience.

Usage

veteran

Format

A data frame with 137 rows and 8 variables:

trt

treatment (1 = standard, 2 = test)

celltype

tumour cell type: squamous, smallcell, adeno, large

time

survival time, days

status

event indicator (1 = dead, 0 = censored)

karno

Karnofsky performance score (0-100)

diagtime

months from diagnosis to randomization

age

age, years

prior

prior therapy (0 = no, 10 = yes)

Source

Kalbfleisch, J. and Prentice, R. (1980) The Statistical Analysis of Failure Time Data. Wiley, New York, as redistributed by the survival package.