rstpm2 (version 1.3.2)

stpm2: Fully parametric generalised survival model

Description

This implements the generalised survival model g(S(t|x)) = eta, where g is a link function, S is survival, t is time, x are covariates and eta is a linear predictor. The main model assumption is that the time effects in the linear predictor are smooth. This extends the class of flexible parametric survival models developed by Royston and colleagues. The model has been extended to include relative survival, Gamma frailties and normal random effects.

Usage

stpm2(formula, data, 
      smooth.formula = NULL, smooth.args = NULL,
      df = 3, cure = FALSE, logH.args = NULL,
      logH.formula = NULL, tvc = NULL, tvc.formula =
      NULL, control = list(parscale = 1, maxit = 300),
      init = NULL, coxph.strata = NULL, weights = NULL,
      robust = FALSE, baseoff = FALSE, bhazard = NULL,
      timeVar = "", time0Var = "", use.gr = TRUE,
      reltol=1.0e-8, trace = 0,
      link.type=c("PH","PO","probit","AH"), 
      frailty = !is.null(cluster), cluster = NULL, logtheta=-6, nodes=9,
      RandDist=c("Gamma","LogN"), 
      contrasts = NULL,
      subset = NULL, ...)

Arguments

formula

a formula object, with the response on the left of a ~ operator, and the regression terms (excluding time) on the right. The response must be a survival object as returned by the Surv function. T

data

a data.frame in which to interpret the variables named in the formula argument. [at present: required]

smooth.formula

a formula for describing the time effects for the linear predictor, including the baseline and the time-dependent effects (default=NULL). Only one of df, smooth.formula, smooth.args, logH.args or l

smooth.args

a list describing the arguments for the nsx function for modelling the baseline time effect on the linear predictor scale (default=NULL). Use this or smooth.formula for changing the knot placement and specifying cure models.

df

an integer that describes the degrees of freedom for the ns function for modelling the baseline log-cumulative hazard (default=3).

logH.args

as per smooth.args. Deprecated.

logH.formula

as per smooth.formula. Deprecated.

tvc

a list with the names of the time-varying coefficients and the degrees of freedom (e.g. tvc=list(x=3) specifies x as a time-varying coefficient with 3 degrees of freedom).

tvc.formula

a formula for describing the time-varying coefficients. If a time-varying coefficient is being model, then only one of tvc and tvc.formula is required.

bhazard

a vector for the background hazard for relative survival estimation. At present, this does not use data and it is required for all individuals - although it is only used at the event times.

control

control argument passed to optim.

init

init should either be FALSE, such that initial values will be determined using Cox regression, or a numeric vector of initial values.

coxph.strata

variable in the data argument for stratification of the coxph model fit for estimating initial values.

weights

an optional vector of 'prior weights' to be used in the fitting process. Should be NULL or a numeric vector.

robust

Boolean used to determine whether to use a robust variance estimator.

baseoff

Boolean used to determine whether fully define the model using tvc.formula rather than combining logH.formula and tvc.formula

timeVar

variable defining the time variable. By default, this is determined from the survival object, however this may be ambiguous if two variables define the time

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

subset

an optional vector specifying a subset of observations to be used in the fitting process.

cure

logical for whether to estimate a cure model.

time0Var

string variable to determine the entry variable; useful for when more than one data variable is used in the entry time.

use.gr

logical indicating whether to use gradients in the calculation

link.type

type of link function.

reltol

relative tolerance for the model convergence

trace

logical for whether to provide trace information

frailty

logical for whether to fit a shared frailty model

cluster

string for the data variable that determines the cluster for the frailty

nodes

number of integration points for Gaussian quadrature

RandDist

type of distribution for the random effect or frailty

logtheta

initial value for log-theta used in the gamma shared frailty model

...

additional arguments to be passed to the mle2 .

Value

• An stpm2-class object that inherits from mle2-class.

Details

The implementation extends the mle2 object from the bbmle package. The model inherits all of the methods from the mle2 class.

The default linear predictor includes a time effect modelled using natural splines for log(time) with three degrees of freedom.

Examples

data(brcancer)
summary(fit <- stpm2(Surv(rectime,censrec==1)~hormon,data=brcancer,df=3))

## some predictions
head(predict(fit,se.fit=TRUE,type="surv"))
head(predict(fit,se.fit=TRUE,type="hazard"))

## some plots
plot(fit,newdata=data.frame(hormon=0),type="hazard")
plot(fit,newdata=data.frame(hormon=0),type="surv")

## the same model using logH.formula
summary(stpm2(Surv(rectime,censrec==1)~hormon,data=brcancer,logH.formula=~ns(log(rectime),df=3)))


## time-varying coefficient
summary(fit.tvc <- stpm2(Surv(rectime,censrec==1)~hormon,data=brcancer,df=3,
                     tvc=list(hormon=3)))
anova(fit,fit.tvc) # compare with and without tvc

## some more plots
plot(fit.tvc,newdata=data.frame(hormon=0),type="hr",var="hormon", ylim=c(0,2))
                                        # no lines method: use add=TRUE
plot(fit.tvc,newdata=data.frame(hormon=1),type="hr",var="hormon",
     add=TRUE,ci=FALSE,line.col=2)

plot(fit.tvc,newdata=data.frame(hormon=0),type="sdiff",var="hormon")

plot(fit.tvc,newdata=data.frame(hormon=0),type="hdiff",var="hormon")

plot(fit.tvc,newdata=data.frame(hormon=0),type="hazard")
plot(fit.tvc,newdata=data.frame(hormon=1),type="hazard",line.col=2,ci=FALSE,add=TRUE)

## compare number of knots
hormon0 <- data.frame(hormon=0)
plot(fit,type="hazard",newdata=hormon0)
AIC(fit)
for (df in 4:6) {
    fit.new <- stpm2(Surv(rectime,censrec==1)~hormon,data=brcancer,df=df)
    plot(fit.new,type="hazard",newdata=hormon0,add=TRUE,ci=FALSE,line.col=df)
    print(AIC(fit.new))
}