Learn R Programming

flexsurv (version 0.1.5)

flexsurvspline: Flexible survival regression using the Royston/Parmar spline model.

Description

Flexible parametric modelling of time-to-event data using the spline model of Royston and Parmar (2002).

Usage

flexsurvspline(formula, data, k=0, knots=NULL, scale="hazard", inits=NULL,
               fixedpars=NULL, cl=0.95,...)

Arguments

formula
A formula expression in conventional R linear modelling syntax. The response must be a survival object as returned by the Surv function, and any covariates are given on the right-hand side. For ex
data
A data frame in which to find variables supplied in formula. If not given, the variables should be in the working environment.
k
Number of knots in the spline. k=0 will give a Weibull, log-logistic or lognormal model, if "scale" is "hazard", "odds" or "normal" respectively. k is equivalent to
knots
Locations of knots on the axis of log time. If not specified, these are chosen as equally-spaced quantiles, for example, at the median with one knot, or at the 33% and 67% quantiles of log time with two knots. Either k or
scale
If "hazard", the log cumulative hazard is modelled as a spline function of log time.

If "odds", the log cumulative odds is modelled as a spline function of log time.

If "normal", $-\Phi^{-1}(S(t)

inits
A numeric vector giving initial values for each unknown parameter. If not specified, default initial values are chosen by estimating the baseline survival at each observed death time from the equivalent Cox model, transforming to the log c
fixedpars
Vector of indices of parameters whose values will be fixed at their initial values during optimisation. The indices are ordered with the intercept "gamma0" first, then the remaining spline coefficients "gamma1","gamma2"...<
cl
Width of symmetric confidence intervals for maximum likelihood estimates, by default 0.95.
...
Optional arguments to the general-purpose Roptimisation routine optim. See flexsurvreg for examples.

Value

  • A list of class "flexsurvreg" with the following elements.
  • callA copy of the function call, for use in post-processing.
  • kNumber of knots.
  • knotsLocation of knots on the log time axis.
  • resMatrix of maximum likelihood estimates and confidence limits. Spline coefficients are labelled "gamma...", and covariate effects are labelled with the names of the covariates.

    Coefficients gamma1,gamma2,... here are the equivalent of s0,s1,... in Stata streg, and gamma0 is the equivalent of the xb constant term. To reproduce results, use the noorthog option in Stata, since no orthogonalisation is performed on the spline basis here. In the Weibull model, for example, gamma0,gamma1 are -shape log(scale), shape respectively in dweibull or flexsurvreg notation, or (-Intercept/scale, 1/scale) in survreg notation.

    In the log-logistic model with shape a and scale b (as in dllogis from the eha package), 1/b^a is equivalent to exp(gamma0), and a is equivalent to gamma1.

    In the log-normal model with log-scale mean mu and standard deviation sigma, -mu/sigma is equivalent to gamma0 and 1/sigma is equivalent to gamma1.

  • loglikThe maximised log-likelihood. This will differ from Stata, where the sum of the log uncensored survival times is added to the log-likelihood in survival models, to remove dependency on the time scale.
  • AICAkaike's information criterion (-2*log likelihood + 2*number of estimated parameters)

Details

In the spline-based survival model of Royston and Parmar (2002), a transformation $g(S(t,z))$ of the survival function is modelled as a natural cubic spline function of log time $x = \log(t)$ plus linear effects of covariates $z$.

$$g(S(t,z)) = s(x, \bm{\gamma}) + \bm{\beta}^T \mathbf{z}$$

The proportional hazards model (scale="hazard") defines $g(S(t,\mathbf{z})) = \log(-\log(S(t,\mathbf{z}))) = \log(H(t,\mathbf{z}))$, the log cumulative hazard.

The proportional odds model (scale="odds") defines $g(S(t,\mathbf{z})) = \log(S(t,\mathbf{z})^{-1} - 1)$, the log cumulative odds. The probit model (scale="normal") defines $g(S(t,\mathbf{z})) = -\Phi^{-1}(S(t,\mathbf{z}))$, where $\Phi^{-1}()$ is the inverse normal distribution function qnorm.

With no knots, the spline reduces to a linear function, and these models are equivalent to Weibull, log-logistic and lognormal models respectively. Natural cubic splines are cubic splines constrained to be linear beyond boundary knots $k_{min},k_{max}$. The spline function is defined as

$$s(x,\bm{\gamma}) = \gamma_0 + \gamma_1 x + \gamma_2 v_1(x) + \ldots + \gamma_{m+1} v_m(x)$$

where $v_j(x)$ is the $j$th basis function

$$v_j(x) = (x - k_j)^3_+ - \lambda_j(x - k_{min})^3_+ - (1 - \lambda_j) (x - k_{max})^3_+$$

$$\lambda_j = \frac{k_{max} - k_j}{k_{max} - k_{min}}$$

and $(x - a)_+ = max(0, x - a)$. Parameters $\bm{\gamma},\bm{\beta}$ are estimated by maximum likelihood using the algorithms available in the standard R optim function. Parameters defined to be positive are estimated on the log scale. Confidence intervals are estimated from the Hessian at the maximum, and transformed back to the original scale of the parameters.

References

Royston, P. and Parmar, M. (2002). Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Statistics in Medicine 21(1):2175-2197.

See Also

flexsurvreg for flexible survival modelling using fully parametric distributions including the generalized F and gamma.

plot.flexsurvreg and lines.flexsurvreg to plot fitted survival, hazards and cumulative hazards from models fitted by flexsurvspline and flexsurvreg.

Examples

Run this code
data(bc)
bc$recyrs <- bc$rectime/365

## Best-fitting model to breast cancer data from Royston and Parmar (2002) 
## One internal knot (2 df) and cumulative odds scale 
spl <- flexsurvspline(Surv(recyrs, censrec) ~ group, data=bc, k=1, scale="odds")

## Fitted survival 
plot(spl)

## Simple Weibull model fits much less well 
splw <- flexsurvspline(Surv(recyrs, censrec) ~ group, data=bc, k=0, scale="hazard")
lines(splw, col="blue")

## Alternative way of fitting the Weibull 
splw2 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data=bc, dist="weibull")

Run the code above in your browser using DataLab