Package 'flexsurv'

Description

flexsurv: Flexible parametric models for time-to-event data, including the generalized gamma, the generalized F and the Royston-Parmar spline model, and extensible to user-defined distributions.

Details

flexsurvreg fits parametric models for time-to-event (survival) data. Data may be right-censored, and/or left-censored, and/or left-truncated. Several built-in parametric distributions are available. Any user-defined parametric model can also be employed by supplying a list with basic information about the distribution, including the density or hazard and ideally also the cumulative distribution or hazard.

Covariates can be included using a linear model on any parameter of the distribution, log-transformed to the real line if necessary. This typically defines an accelerated failure time or proportional hazards model, depending on the distribution and parameter.

flexsurvspline fits the flexible survival model of Royston and Parmar (2002) in which the log cumulative hazard is modelled as a natural cubic spline function of log time. Covariates can be included on any of the spline parameters, giving either a proportional hazards model or an arbitrarily-flexible time-dependent effect. Alternative proportional odds or probit parameterisations are available.

Output from the models can be presented as survivor, cumulative hazard and hazard functions (summary.flexsurvreg). These can be plotted against nonparametric estimates (plot.flexsurvreg) to assess goodness-of-fit. Any other user-defined function of the parameters may be summarised in the same way.

Multi-state models for time-to-event data can also be fitted with the same functions. Predictions from those models can then be made using the functions pmatrix.fs, pmatrix.simfs, totlos.fs, totlos.simfs, or sim.fmsm, or alternatively by msfit.flexsurvreg followed by mssample or probtrans from the package mstate.

Distribution (“dpqr”) functions for the generalized gamma and F distributions are given in GenGamma, GenF (preferred parameterisations) and GenGamma.orig, GenF.orig (original parameterisations). flexsurv also includes the standard Gompertz distribution with unrestricted shape parameter, see Gompertz.

User guide

The flexsurv user guide vignette explains the methods in detail, and gives several worked examples. A further vignette flexsurv-examples gives a few more complicated examples, and users are encouraged to submit their own.

Author(s)

Christopher Jackson [email protected]

References

Jackson, C. (2016). flexsurv: A Platform for Parametric Survival Modeling in R. Journal of Statistical Software, 70(8), 1-33. doi:10.18637/jss.v070.i08

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.

Cox, C. (2008). The generalized $F$ distribution: An umbrella for parametric survival analysis. Statistics in Medicine 27:4301-4312.

Cox, C., Chu, H., Schneider, M. F. and Muñoz, A. (2007). Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Statistics in Medicine 26:4252-4374

See Also

Useful links:

• https://github.com/chjackson/flexsurv

• http://chjackson.github.io/flexsurv/

• Report bugs at https://github.com/chjackson/flexsurv/issues

Description

Defined as the sum of the AICs of the transition-specific models.

Usage

## S3 method for class 'fmsm'
AIC(object, ..., k = 2)

Arguments

Value

The sum of the AICs of the transition-specific models.

Description

Generic function for the second-order Akaike information criterion. The only current implementation of this in flexsurv is in AICc.flexsurvreg, see there for more details.

Usage

AICc(object, ...)

AICC(object, ...)

Arguments

Details

This can be spelt either as AICC or AICc.

Description

Second-order or "corrected" Akaike information criterion, often known as AICc (see, e.g. Burnham and Anderson 2002). This is defined as -2 log-likelihood + (2*p*n)/(n - p -1), whereas the standard AIC is defined as -2 log-likelihood + 2*p, where p is the number of parameters and n is the sample size. The correction is intended to adjust AIC for small-sample bias, hence it only makes a difference to the result for small n.

Usage

## S3 method for class 'flexsurvreg'
AICc(object, cens = TRUE, ...)

## S3 method for class 'flexsurvreg'
AICC(object, cens = TRUE, ...)

Arguments

Details

This can be spelt either as AICC or AICc.

Value

The second-order AIC of the fitted model.

References

Burnham, K. P., Anderson, D. R. (2002) Model Selection and Multimodel Inference: a practical information-theoretic approach. Second edition. Springer: New York.

See Also

BIC, AIC, BIC.flexsurvreg, nobs.flexsurvreg

Description

Given a fitted flexsurvmix model, return the Aalen-Johansen estimates of the probability of occupying each state at a series of times covering the observed data. State 1 represents not having experienced any of the competing events, while state 2 and any further states correspond to having experienced each of the competing events respectively. These estimates can be compared with the fitted probabilities returned by p_flexsurvmix to check the fit of a flexsurvmix model.

Usage

ajfit(x, newdata = NULL, tidy = TRUE)

Arguments

Details

This is only supported for models with no covariates or models containing only factor covariates.

For models with factor covariates, the Aalen-Johansen estimates are computed for the subsets of the data defined in newdata. If newdata is not supplied, then this function returns state occupancy probabilities for all possible combinations of the factor levels.

The Aalen-Johansen estimates are computed using survfit from the survival package (Therneau 2020).

References

Therneau T (2020). _A Package for Survival Analysis in R_. R package version 3.2-3, <URL: https://CRAN.R-project.org/package=survival>.

Description

This computes Aalen-Johansen estimates of the probability of occupying each state at a series of times, using ajfit. The equivalent estimates from the parametric model are then produced using p_flexsurvmix, and concatenated with the nonparametric estimates to form a tidy data frame. This data frame can then simply be plotted using ggplot.

Usage

ajfit_flexsurvmix(x, maxt = NULL, startname = "Start", B = NULL)

Arguments

Description

Computes both parametric and comparable Aalen-Johansen nonparametric estimates from a flexible parametric multi-state model, and returns them together in a tidy data frame. Only models with no covariates, or only factor covariates, are supported. If there are factor covariates, then the nonparametric estimates are computed for subgroups defined by combinations of the covariates. Another restriction of this function is that all transitions must have the same covariates on them.

Usage

ajfit_fmsm(x, maxt = NULL, newdata = NULL)

Arguments

Value

Tidy data frame containing both Aalen-Johansen and parametric estimates of state occupancy over time, and by subgroup if subgroups are included.

Description

Augment accepts a model object and a dataset and adds information about each observation in the dataset. Most commonly, this includes predicted values in the .fitted column, residuals in the .resid column, and standard errors for the fitted values in a .se.fit column. New columns always begin with a . prefix to avoid overwriting columns in the original dataset.

Usage

## S3 method for class 'flexsurvreg'
augment(
  x,
  data = NULL,
  newdata = NULL,
  type.predict = "response",
  type.residuals = "response",
  ...
)

Arguments

Details

If neither of data or newdata are specified, then model.frame(x) will be used. It is worth noting that model.frame(x) will include a Surv object and not the original time-to-event variables used when fitting the flexsurvreg object. If the original data is desired, specify data.

Value

A tibble containing data or newdata and possible additional columns:

• .fitted Fitted values of model

• .se.fit Standard errors of fitted values

• .resid Residuals (not present if newdata specified)

Examples

fit <- flexsurvreg(formula = Surv(futime, fustat) ~ age, data = ovarian, dist = "exp")
augment(fit, data = ovarian)

Description

Compute a basis for a natural cubic spline, by default using the parameterisation described by Royston and Parmar (2002). Used for flexible parametric survival models.

Usage

basis(knots, x, spline = "rp")

Arguments

Details

The exact formula for the basis is given in flexsurvspline.

Value

A matrix with one row for each ordinate and one column for each knot.

basis returns the basis, and dbasis returns its derivative with respect to x.

fss and dfss are the same, but with the order of the arguments swapped around for consistency with similar functions in other R packages.

Author(s)

Christopher Jackson <[email protected]>

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.

Wang W, Yan J (2021). Shape-Restricted Regression Splines with R Package splines2. Journal of Data Science, 19(3), 498-517.

See Also

flexsurvspline.

Description

Survival times of 686 patients with primary node positive breast cancer.

Usage

bc

Format

A data frame with 686 rows.

Source

German Breast Cancer Study Group, 1984-1989. Used as a reference dataset for the spline-based survival model of Royston and Parmar (2002), implemented here in flexsurvspline. Originally provided with the stpm (Royston 2001, 2004) and stpm2 (Lambert 2009, 2010) Stata modules.

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.

Sauerbrei, W. and Royston, P. (1999). Building multivariable prognostic and diagnostic models: transformation of the predictors using fractional polynomials. Journal of the Royal Statistical Society, Series A 162:71-94.

See Also

flexsurvspline

Description

Bayesian Information Criterion (BIC) for comparison of flexsurvreg models

Usage

## S3 method for class 'flexsurvreg'
BIC(object, cens = TRUE, ...)

Arguments

Details

There is no "official" definition of what the sample size should be for the use of BIC in censored survival analysis. BIC is based on an approximation to Bayesian model comparison using Bayes factors and an implicit vague prior. Informally, the sample size describes the number of "units" giving rise to a distinct piece of information (Kass and Raftery 1995). However censored observations provide less information than observed events, in principle. The default used here is the number of individuals, for consistency with more familiar kinds of statistical modelling. However if censoring is heavy, then the number of events may be a better represent the amount of information. Following these principles, the best approximation would be expected to be somewere in between.

AIC and BIC are intended to measure different things. Briefly, AIC measures predictive ability, whereas BIC is expected to choose the true model from a set of models where one of them is the truth. Therefore BIC chooses simpler models for all but the tiniest sample sizes ($log(n)>2$, $n>7$). AIC might be preferred in the typical situation where "all models are wrong but some are useful". AIC also gives similar results to cross-validation (Stone 1977).

Value

The BIC of the fitted model. This is minus twice the log likelihood plus p*log(n), where p is the number of parameters and n is the sample size of the data. If weights was supplied to flexsurv, the sample size is defined as the sum of the weights.

References

Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773-795.

Stone, M. (1977). An asymptotic equivalence of choice of model by cross‐validation and Akaike's criterion. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 44-47.

See Also

BIC, AIC, AICC.flexsurvreg, nobs.flexsurvreg

Description

Calculate a confidence interval for a model output by repeatedly replacing the parameters in a fitted model object with a draw from the multivariate normal distribution of the maximum likelihood estimates, then recalculating the output function.

Usage

bootci.fmsm(
  x,
  B,
  fn,
  cl = 0.95,
  attrs = NULL,
  cores = NULL,
  sample = FALSE,
  ...
)

Arguments

Value

A matrix with two rows, giving the upper and lower confidence limits respectively. Each row is a vector of the same length as the unlisted result of the function corresponding to fncall.

Examples

## How to use bootci.fmsm

## Write a function with one argument called x giving a fitted model,
## and returning some results of the model.  The results may be in any form.   

tmat <- rbind(c(NA,1,2),c(NA,NA,3),c(NA,NA,NA))
bexp <- flexsurvreg(Surv(Tstart, Tstop, status) ~ trans, data=bosms3, dist="exp")

summfn <- function(x, t){
 resp <-  flexsurv::pmatrix.fs(x, trans=tmat, t=t)
 rest <- flexsurv::totlos.fs(x, trans=tmat, t=t)
 list(resp, rest)
}

## Use bootci.fmsm to obtain the confidence interval
## The matrix columns are in the order of the unlisted results of the original
## summfn.  You will have to rearrange them into the format that you want.
## If summfn has any extra arguments, in this case \code{t}, make sure they are
## passed through via the ... argument to bootci.fmsm

bootci.fmsm(bexp, B=3, fn=summfn, t=10)
bootci.fmsm(bexp, B=3, fn=summfn, t=5)

Description

A dataset containing histories of bronchiolitis obliterans syndrome (BOS) from lung transplant recipients. BOS is a chronic decline in lung function, often observed after lung transplantation.

Format

A data frame containing a sequence of observed or censored transitions to the next stage of severity or death. It is grouped by patient and includes histories of 204 patients. All patients start in state 1 (no BOS) at six months after transplant, and may subsequently develop BOS or die.

bosms3 contains the data for a three-state model: no BOS, BOS or death. bosms4 uses a four-state representation: no BOS, mild BOS, moderate/severe BOS or death.

Details

The entry time of each patient into each stage of BOS was estimated by clinicians, based on their history of lung function measurements and acute rejection and infection episodes. BOS is only assumed to occur beyond six months after transplant. In the first six months the function of each patient's new lung stabilises. Subsequently BOS is diagnosed by comparing the lung function against the "baseline" value.

The same data are provided in the msm package, but in the native format of msm to allow Markov models to be fitted. In flexsurv, much more flexible models can be fitted.

Source

Papworth Hospital, U.K.

References

Heng. D. et al. (1998). Bronchiolitis Obliterans Syndrome: Incidence, Natural History, Prognosis, and Risk Factors. Journal of Heart and Lung Transplantation 17(12)1255–1263.

Description

Extract model coefficients from fitted flexible survival models. This presents all parameter estimates, transformed to the real line if necessary. For example, shape or scale parameters, which are constrained to be positive, are returned on the log scale.

Usage

## S3 method for class 'flexsurvreg'
coef(object, ...)

Arguments

Details

This matches the behaviour of coef.default for standard R model families such as glm, where intercepts in regression models are presented on the same scale as the covariate effects. Note that any parameter in a distribution fitted by flexsurvreg or flexsurvreg may be an intercept in a regression model.

Value

This returns the mod$res.t[,"est"] component of the fitted model object mod. See flexsurvreg, flexsurvspline for full documentation of all components.

Author(s)

C. H. Jackson [email protected]

See Also

flexsurvreg, flexsurvspline.

Description

Cox-Snell residuals from a parametric survival model

Usage

coxsnell_flexsurvreg(x)

Arguments

Value

A data frame with a column called est giving the Cox-Snell residual, defined as the fitted cumulative hazard at each data point. fitted cumulative hazard at the given observed data point, and other columns indicating the observation time, observed event status, and covariate values defining the data at this point.

The cumulative hazards est should form a censored sample from an Exponential(1). Therefore to check the fit of the model, plot a nonparametric estimate of the cumulative hazard curve against a diagonal line through the origin, which is the theoretical cumulative hazard trajectory of the Exponential(1).

Examples

fitg <- flexsurvreg(formula = Surv(futime, fustat) ~ age, data = ovarian, dist = "gengamma")
  cs <- coxsnell_flexsurvreg(fitg)
  
  ## Model appears to fit well, with some small sample noise 
  surv <- survfit(Surv(cs$est, ovarian$fustat) ~ 1)
  plot(surv, fun="cumhaz")
  abline(0, 1, col="red")

Description

In a mixture model for competing events, an individual can experience one of a set of different events. We specify a model for the probability that they will experience each event before the others, and a model for the time to the event conditionally on that event occurring first.

Usage

flexsurvmix(
  formula,
  data,
  event,
  dists,
  pformula = NULL,
  anc = NULL,
  partial_events = NULL,
  initp = NULL,
  inits = NULL,
  fixedpars = NULL,
  dfns = NULL,
  method = "direct",
  em.control = NULL,
  optim.control = NULL,
  aux = NULL,
  sr.control = survreg.control(),
  integ.opts,
  hess.control = NULL,
  ...
)

Arguments

Details

This differs from the more usual "competing risks" models, where we specify "cause-specific hazards" describing the time to each competing event. This time will not be observed for an individual if one of the competing events happens first. The event that happens first is defined by the minimum of the times to the alternative events.

The flexsurvmix function fits a mixture model to data consisting of a single time to an event for each individual, and an indicator for what type of event occurs for that individual. The time to event may be observed or censored, just as in flexsurvreg, and the type of event may be known or unknown. In a typical application, where we follow up a set of individuals until they experience an event or a maximum follow-up time is reached, the event type is known if the time is observed, and the event type is unknown when follow-up ends and the time is right-censored.

The model is fitted by maximum likelihood, either directly or by using an expectation-maximisation (EM) algorithm, by wrapping flexsurvreg to compute the likelihood or to implement the E and M steps.

Some worked examples are given in the package vignette about multi-state modelling, which can be viewed by running vignette("multistate", package="flexsurv").

Value

List of objects containing information about the fitted model. The important one is res, a data frame containing the parameter estimates and associated information.

References

Jackson, C. H. and Tom, B. D. M. and Kirwan, P. D. and Mandal, S. and Seaman, S. R. and Kunzmann, K. and Presanis, A. M. and De Angelis, D. (2022) A comparison of two frameworks for multi-state modelling, applied to outcomes after hospital admissions with COVID-19. Statistical Methods in Medical Research 31(9) 1656-1674.

Larson, M. G., & Dinse, G. E. (1985). A mixture model for the regression analysis of competing risks data. Journal of the Royal Statistical Society: Series C (Applied Statistics), 34(3), 201-211.

Lau, B., Cole, S. R., & Gange, S. J. (2009). Competing risk regression models for epidemiologic data. American Journal of Epidemiology, 170(2), 244-256.

Description

Parametric modelling or regression for time-to-event data. Several built-in distributions are available, and users may supply their own.

Usage

flexsurvreg(
  formula,
  anc = NULL,
  data,
  weights,
  bhazard,
  rtrunc,
  subset,
  na.action,
  dist,
  inits,
  fixedpars = NULL,
  dfns = NULL,
  aux = NULL,
  cl = 0.95,
  integ.opts = NULL,
  sr.control = survreg.control(),
  hessian = TRUE,
  hess.control = NULL,
  ...
)

Arguments

Details

Parameters 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.

The usage of flexsurvreg is intended to be similar to survreg in the survival package.

Value

A list of class "flexsurvreg" containing information about the fitted model. Components of interest to users may include:

Custom distributions

flexsurvreg is intended to be easy to extend to handle new distributions. To define a new distribution for use in flexsurvreg, construct a list with the following elements:

"name"

A string naming the distribution. If this is called "dist", for example, then there must be visible in the working environment, at least, either

a) a function called ddist which defines the probability density,

or

b) a function called hdist which defines the hazard.

Ideally, in case a) there should also be a function called pdist which defines the probability distribution or cumulative density, and in case b) there should be a function called Hdist defining the cumulative hazard. If these additional functions are not provided, flexsurv attempts to automatically create them by numerically integrating the density or hazard function. However, model fitting will be much slower, or may not even work at all, if the analytic versions of these functions are not available.

The functions must accept vector arguments (representing different times, or alternative values for each parameter) and return the results as a vector. The function Vectorize may be helpful for doing this: see the example below. These functions may be in an add-on package (see below for an example) or may be user-written. If they are user-written they must be defined in the global environment, or supplied explicitly through the dfns argument to flexsurvreg. The latter may be useful if the functions are created dynamically (as in the source of flexsurvspline) and thus not visible through R's scoping rules.

Arguments other than parameters must be named in the conventional way – for example x for the first argument of the density function or hazard, as in dnorm(x, ...) and q for the first argument of the probability function. Density functions should also have an argument log, after the parameters, which when TRUE, computes the log density, using a numerically stable additive formula if possible.

Additional functions with names beginning with "DLd" and "DLS" may be defined to calculate the derivatives of the log density and log survival probability, with respect to the parameters of the distribution. The parameters are expressed on the real line, for example after log transformation if they are defined as positive. The first argument must be named t, representing the time, and the remaining arguments must be named as the parameters of the density function. The function must return a matrix with rows corresponding to times, and columns corresponding to the parameters of the distribution. The derivatives are used, if available, to speed up the model fitting with optim.

"pars"

Vector of strings naming the parameters of the distribution. These must be the same names as the arguments of the density and probability functions.

"location"

Name of the main parameter governing the mean of the distribution. This is the default parameter on which covariates are placed in the formula supplied to flexsurvreg.

"transforms"

List of R functions which transform the range of values taken by each parameter onto the real line. For example, c(log, log) for a distribution with two positive parameters.

"inv.transforms"

List of R functions defining the corresponding inverse transformations. Note these must be lists, even for single parameter distributions they should be supplied as, e.g. c(exp) or list(exp).

"inits"

A function of the observed survival times t (including right-censoring times, and using the halfway point for interval-censored times) which returns a vector of reasonable initial values for maximum likelihood estimation of each parameter. For example, function(t){ c(1, mean(t)) } will always initialize the first of two parameters at 1, and the second (a scale parameter, for instance) at the mean of t.

For example, suppose we want to use an extreme value survival distribution. This is available in the CRAN package eha, which provides conventionally-defined density and probability functions called eha::dEV and eha::pEV. See the Examples below for the custom list in this case, and the subsequent command to fit the model.

Author(s)

Christopher Jackson <[email protected]>

References

Jackson, C. (2016). flexsurv: A Platform for Parametric Survival Modeling in R. Journal of Statistical Software, 70(8), 1-33. doi:10.18637/jss.v070.i08

Cox, C. (2008) The generalized $F$ distribution: An umbrella for parametric survival analysis. Statistics in Medicine 27:4301-4312.

Cox, C., Chu, H., Schneider, M. F. and Muñoz, A. (2007) Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Statistics in Medicine 26:4252-4374

Jackson, C. H. and Sharples, L. D. and Thompson, S. G. (2010) Survival models in health economic evaluations: balancing fit and parsimony to improve prediction. International Journal of Biostatistics 6(1):Article 34.

Nelson, C. P., Lambert, P. C., Squire, I. B., & Jones, D. R. (2007). Flexible parametric models for relative survival, with application in coronary heart disease. Statistics in medicine, 26(30), 5486-5498.

See Also

flexsurvspline for flexible survival modelling using the spline model of Royston and Parmar.

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

Examples

## Compare generalized gamma fit with Weibull
fitg <- flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, dist="gengamma")
fitg
fitw <- flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian, dist="weibull")
fitw
plot(fitg)
lines(fitw, col="blue", lwd.ci=1, lty.ci=1)
## Identical AIC, probably not enough data in this simple example for a
## very flexible model to be worthwhile.

## Custom distribution
## make "dEV" and "pEV" from eha package (if installed)
##   available to the working environment
if (require("eha")) {
custom.ev <- list(name="EV",
                      pars=c("shape","scale"),
                      location="scale",
                      transforms=c(log, log),
                      inv.transforms=c(exp, exp),
                      inits=function(t){ c(1, median(t)) })
fitev <- flexsurvreg(formula = Surv(futime, fustat) ~ 1, data = ovarian,
                    dist=custom.ev)
fitev
lines(fitev, col="purple", col.ci="purple")
}


## Custom distribution: supply the hazard function only
hexp2 <- function(x, rate=1){ rate } # exponential distribution
hexp2 <- Vectorize(hexp2)
custom.exp2 <- list(name="exp2", pars=c("rate"), location="rate",
                    transforms=c(log), inv.transforms=c(exp),
                    inits=function(t)1/mean(t))
flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist=custom.exp2)
flexsurvreg(Surv(futime, fustat) ~ 1, data = ovarian, dist="exp")
## should give same answer

Description

This function estimates the distribution of the time between an initial and final event, in situations where individuals are only observed if they have experienced both events before a certain time, thus they are right-truncated at this time. The time of the initial event provides information about the time from initial to final event, given the truncated observation scheme, and initial events are assumed to occur with an exponential growth rate.

Usage

flexsurvrtrunc(
  t,
  tinit,
  rtrunc,
  tmax,
  data = NULL,
  method = "joint",
  dist,
  theta = NULL,
  fixed.theta = TRUE,
  inits = NULL,
  fixedpars = NULL,
  dfns = NULL,
  integ.opts = NULL,
  cl = 0.95,
  optim_control = list()
)

Arguments

Details

Covariates are not currently supported.

Note that flexsurvreg, with an rtrunc argument, can fit models for a similar kind of data, but those models ignore the information provided by the time of the initial event.

A nonparametric estimator of survival under right-truncation is also provided in survrtrunc. See Seaman et al. (2020) for a full comparison of the alternative models.

References

Seaman, S., Presanis, A. and Jackson, C. (2020) Estimating a Time-to-Event Distribution from Right-Truncated Data in an Epidemic: a Review of Methods

See Also

flexsurvreg, survrtrunc.

Examples

set.seed(1) 
## simulate time to initial event
X <- rexp(1000, 0.2)
## simulate time between initial and final event
T <- rgamma(1000, 2, 10) 

## right-truncate to keep only those with final event time
## before tmax
tmax <- 40
obs <- X + T < tmax 
rtrunc <- tmax - X
dat <- data.frame(X, T, rtrunc)[obs,]

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gamma", theta=0.2)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gamma", theta=0.2, fixed.theta=FALSE)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gamma", theta=0.2, inits=c(1, 8))

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gamma", theta=0.2, method="final")

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gamma", fixed.theta=TRUE)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="weibull", fixed.theta=TRUE)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="lnorm", fixed.theta=TRUE)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gengamma", fixed.theta=TRUE)

flexsurvrtrunc(t=T, rtrunc=rtrunc, tinit=X, tmax=40, data=dat,
                dist="gompertz", fixed.theta=TRUE)

Description

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

Usage

flexsurvspline(
  formula,
  data,
  weights,
  bhazard,
  rtrunc,
  subset,
  k = 0,
  knots = NULL,
  bknots = NULL,
  scale = "hazard",
  timescale = "log",
  spline = "rp",
  ...
)

Arguments

Details

This function works as a wrapper around flexsurvreg by dynamically constructing a custom distribution using dsurvspline, psurvspline and unroll.function.

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.

The spline coefficients $\gamma_j: j=1, 2 \ldots$, which are called the "ancillary parameters" above, may also be modelled as linear functions of covariates $\mathbf{z}$, as

$\gamma_j(\mathbf{z}) = \gamma_{j0} + \gamma_{j1}z_1 + \gamma_{j2}z_2 + ...$

giving a model where the effects of covariates are arbitrarily flexible functions of time: a non-proportional hazards or odds model.

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,\boldsymbol{\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)$.

Value

A list of class "flexsurvreg" with the same elements as described in flexsurvreg, and including extra components describing the spline model. See in particular:

Author(s)

Christopher Jackson <[email protected]>

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.

Wang W, Yan J (2021). Shape-Restricted Regression Splines with R Package splines2. Journal of Data Science, 19(3), 498-517.

Jackson, C. (2016). flexsurv: A Platform for Parametric Survival Modeling in R. Journal of Statistical Software, 70(8), 1-33. doi:10.18637/jss.v070.i08

See Also

flexsurvreg for flexible survival modelling using general parametric distributions.

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

Examples

## 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, lwd=3, ci=FALSE)

## Simple Weibull model fits much less well

splw <- flexsurvspline(Surv(recyrs, censrec) ~ group, data=bc, k=0, scale="hazard")
lines(splw, col="blue", ci=FALSE)

## Alternative way of fitting the Weibull

## Not run: 
splw2 <- flexsurvreg(Surv(recyrs, censrec) ~ group, data=bc, dist="weibull")

## End(Not run)

Description

Constructor for a mixture multi-state model based on flexsurvmix

Usage

fmixmsm(...)

Arguments

Value

A list of flexsurvmix objects, with the following attribute(s):

pathways A list of all potential pathways until absorption, for models without cycles. For models with cycles this will have an element has_cycle=TRUE, plus the pathways discovered before the function found the cycle and gave up.

Description

Construct a multi-state model from a set of parametric survival models

Usage

fmsm(..., trans)

Arguments

Value

A list containing the objects given in ..., and with attributes "trans" and "statenames" defining the transition structure matrix and state names, and with list components named to describe the transitions they correspond to. If any of the arguments in ... are named, then these are used to define the transition names, otherwise default names are chosen based on the state names.

Description

Density, distribution function, hazards, quantile function and random generation for the generalized F distribution, using the reparameterisation by Prentice (1975).

Usage

dgenf(x, mu = 0, sigma = 1, Q, P, log = FALSE)

pgenf(q, mu = 0, sigma = 1, Q, P, lower.tail = TRUE, log.p = FALSE)

Hgenf(x, mu = 0, sigma = 1, Q, P)

hgenf(x, mu = 0, sigma = 1, Q, P)

qgenf(p, mu = 0, sigma = 1, Q, P, lower.tail = TRUE, log.p = FALSE)

rgenf(n, mu = 0, sigma = 1, Q, P)

Arguments

Details

If $y \sim F(2s_1, 2s_2)$, and $w =$$\log(y)$ then $x = \exp(w\sigma + \mu)$ has the original generalized F distribution with location parameter $\mu$, scale parameter $\sigma>0$ and shape parameters $s_1,s_2$.

In this more stable version described by Prentice (1975), $s_1,s_2$ are replaced by shape parameters $Q,P$, with $P>0$, and

$s_1 = 2(Q^2 + 2P + Q\delta)^{-1}, \quad s_2 = 2(Q^2 + 2P - Q\delta)^{-1}$

equivalently

$Q = \left(\frac{1}{s_1} - \frac{1}{s_2}\right)\left(\frac{1}{s_1} + \frac{1}{s_2}\right)^{-1/2}, \quad P = \frac{2}{s_1 + s_2}$

Define $\delta = (Q^2 + 2P)^{1/2}$, and $w = (\log(x) - \mu)\delta /\sigma$, then the probability density function of $x$ is

$f(x) = \frac{\delta (s_1/s_2)^{s_1} e^{s_1 w}}{\sigma x (1 + s_1 e^w/s_2) ^ {(s_1 + s_2)} B(s_1, s_2)}$

The original parameterisation is available in this package as dgenf.orig, for the sake of completion / compatibility. With the above definitions,

dgenf(x, mu=mu, sigma=sigma, Q=Q, P=P) = dgenf.orig(x, mu=mu, sigma=sigma/delta, s1=s1, s2=s2)

The generalized F distribution with P=0 is equivalent to the generalized gamma distribution dgengamma, so that dgenf(x, mu, sigma, Q, P=0) equals dgengamma(x, mu, sigma, Q). The generalized gamma reduces further to several common distributions, as described in the GenGamma help page.

The generalized F distribution includes the log-logistic distribution (see eha::dllogis) as a further special case:

dgenf(x, mu=mu, sigma=sigma, Q=0, P=1) = eha::dllogis(x, shape=sqrt(2)/sigma, scale=exp(mu))

The range of hazard trajectories available under this distribution are discussed in detail by Cox (2008). Jackson et al. (2010) give an application to modelling oral cancer survival for use in a health economic evaluation of screening.

Value

dgenf gives the density, pgenf gives the distribution function, qgenf gives the quantile function, rgenf generates random deviates, Hgenf retuns the cumulative hazard and hgenf the hazard.

Note

The parameters Q and P are usually called $q$ and $p$ in the literature - they were made upper-case in these R functions to avoid clashing with the conventional arguments q in the probability function and p in the quantile function.

Author(s)

Christopher Jackson <[email protected]>

References

R. L. Prentice (1975). Discrimination among some parametric models. Biometrika 62(3):607-614.

Cox, C. (2008). The generalized $F$ distribution: An umbrella for parametric survival analysis. Statistics in Medicine 27:4301-4312.

Jackson, C. H. and Sharples, L. D. and Thompson, S. G. (2010). Survival models in health economic evaluations: balancing fit and parsimony to improve prediction. International Journal of Biostatistics 6(1):Article 34.

See Also

GenF.orig, GenGamma

Description

Density, distribution function, quantile function and random generation for the generalized F distribution, using the less flexible original parameterisation described by Prentice (1975).

Usage

dgenf.orig(x, mu = 0, sigma = 1, s1, s2, log = FALSE)

pgenf.orig(q, mu = 0, sigma = 1, s1, s2, lower.tail = TRUE, log.p = FALSE)

Hgenf.orig(x, mu = 0, sigma = 1, s1, s2)

hgenf.orig(x, mu = 0, sigma = 1, s1, s2)

qgenf.orig(p, mu = 0, sigma = 1, s1, s2, lower.tail = TRUE, log.p = FALSE)

rgenf.orig(n, mu = 0, sigma = 1, s1, s2)

Arguments

Details

If $y \sim F(2s_1, 2s_2)$, and $w = \log(y)$ then $x = \exp(w\sigma + \mu)$ has the original generalized F distribution with location parameter $\mu$, scale parameter $\sigma>0$ and shape parameters $s_1>0,s_2>0$. The probability density function of $x$ is

$f(x | \mu, \sigma, s_1, s_2) = \frac{(s_1/s_2)^{s_1} e^{s_1 w}}{\sigma x (1 + s_1 e^w/s_2) ^ {(s_1 + s_2)} B(s_1, s_2)}$

where $w = (\log(x) - \mu)/\sigma$, and $B(s_1,s_2) = \Gamma(s_1)\Gamma(s_2)/\Gamma(s_1+s_2)$ is the beta function.

As $s_2 \rightarrow \infty$, the distribution of $x$ tends towards an original generalized gamma distribution with the following parameters:

dgengamma.orig(x, shape=1/sigma, scale=exp(mu) / s1^sigma, k=s1)

See GenGamma.orig for how this includes several other common distributions as special cases.

The alternative parameterisation of the generalized F distribution, originating from Prentice (1975) and given in this package as GenF, is preferred for statistical modelling, since it is more stable as $s_1$ tends to infinity, and includes a further new class of distributions with negative first shape parameter. The original is provided here for the sake of completion and compatibility.

Value

dgenf.orig gives the density, pgenf.orig gives the distribution function, qgenf.orig gives the quantile function, rgenf.orig generates random deviates, Hgenf.orig retuns the cumulative hazard and hgenf.orig the hazard.

Author(s)

Christopher Jackson <[email protected]>

References

R. L. Prentice (1975). Discrimination among some parametric models. Biometrika 62(3):607-614.

See Also

GenF, GenGamma.orig, GenGamma

Description

Density, distribution function, hazards, quantile function and random generation for the generalized gamma distribution, using the parameterisation originating from Prentice (1974). Also known as the (generalized) log-gamma distribution.

Usage

dgengamma(x, mu = 0, sigma = 1, Q, log = FALSE)

pgengamma(q, mu = 0, sigma = 1, Q, lower.tail = TRUE, log.p = FALSE)

Hgengamma(x, mu = 0, sigma = 1, Q)

hgengamma(x, mu = 0, sigma = 1, Q)

qgengamma(p, mu = 0, sigma = 1, Q, lower.tail = TRUE, log.p = FALSE)

rgengamma(n, mu = 0, sigma = 1, Q)

Arguments

Details

If $\gamma \sim Gamma(Q^{-2}, 1)$ , and $w = log(Q^2 \gamma) / Q$, then $x = \exp(\mu + \sigma w)$ follows the generalized gamma distribution with probability density function

$f(x | \mu, \sigma, Q) = \frac{|Q|(Q^{-2})^{Q^{-2}}}{\sigma x \Gamma(Q^{-2})} \exp(Q^{-2}(Qw - \exp(Qw)))$

This parameterisation is preferred to the original parameterisation of the generalized gamma by Stacy (1962) since it is more numerically stable near to $Q=0$ (the log-normal distribution), and allows $Q<=0$. The original is available in this package as dgengamma.orig, for the sake of completion and compatibility with other software - this is implicitly restricted to Q>0 (or k>0 in the original notation). The parameters of dgengamma and dgengamma.orig are related as follows.

dgengamma.orig(x, shape=shape, scale=scale, k=k) =

dgengamma(x, mu=log(scale) + log(k)/shape, sigma=1/(shape*sqrt(k)), Q=1/sqrt(k))

The generalized gamma distribution simplifies to the gamma, log-normal and Weibull distributions with the following parameterisations:

The properties of the generalized gamma and its applications to survival analysis are discussed in detail by Cox (2007).

The generalized F distribution GenF extends the generalized gamma to four parameters.

Value

dgengamma gives the density, pgengamma gives the distribution function, qgengamma gives the quantile function, rgengamma generates random deviates, Hgengamma retuns the cumulative hazard and hgengamma the hazard.

Author(s)

Christopher Jackson <[email protected]>

References

Prentice, R. L. (1974). A log gamma model and its maximum likelihood estimation. Biometrika 61(3):539-544.

Farewell, V. T. and Prentice, R. L. (1977). A study of distributional shape in life testing. Technometrics 19(1):69-75.

Lawless, J. F. (1980). Inference in the generalized gamma and log gamma distributions. Technometrics 22(3):409-419.

Cox, C., Chu, H., Schneider, M. F. and Muñoz, A. (2007). Parametric survival analysis and taxonomy of hazard functions for the generalized gamma distribution. Statistics in Medicine 26:4252-4374

Stacy, E. W. (1962). A generalization of the gamma distribution. Annals of Mathematical Statistics 33:1187-92

See Also

GenGamma.orig, GenF, Lognormal, GammaDist, Weibull.

Description

Density, distribution function, hazards, quantile function and random generation for the generalized gamma distribution, using the original parameterisation from Stacy (1962).

Usage

dgengamma.orig(x, shape, scale = 1, k, log = FALSE)

pgengamma.orig(q, shape, scale = 1, k, lower.tail = TRUE, log.p = FALSE)

Hgengamma.orig(x, shape, scale = 1, k)

hgengamma.orig(x, shape, scale = 1, k)

qgengamma.orig(p, shape, scale = 1, k, lower.tail = TRUE, log.p = FALSE)

rgengamma.orig(n, shape, scale = 1, k)

Arguments

Details

If $w \sim Gamma(k,1)$, then $x = \exp(w/shape + \log(scale))$ follows the original generalised gamma distribution with the parameterisation given here (Stacy 1962). Defining shape$=b>0$, scale$=a>0$, $x$ has probability density

$f(x | a, b, k) = \frac{b}{\Gamma(k)} \frac{x^{bk - 1}}{a^{bk}}$

$\exp(-(x/a)^b)$

The original generalized gamma distribution simplifies to the gamma, exponential and Weibull distributions with the following parameterisations:

Also as k tends to infinity, it tends to the log normal (as in dlnorm) with the following parameters (Lawless, 1980):

dlnorm(x, meanlog=log(scale) + log(k)/shape, sdlog=1/(shape*sqrt(k)))

For more stable behaviour as the distribution tends to the log-normal, an alternative parameterisation was developed by Prentice (1974). This is given in dgengamma, and is now preferred for statistical modelling. It is also more flexible, including a further new class of distributions with negative shape k.

The generalized F distribution GenF.orig, and its similar alternative parameterisation GenF, extend the generalized gamma to four parameters.

Value

dgengamma.orig gives the density, pgengamma.orig gives the distribution function, qgengamma.orig gives the quantile function, rgengamma.orig generates random deviates, Hgengamma.orig retuns the cumulative hazard and hgengamma.orig the hazard.

Author(s)

Christopher Jackson <[email protected]>

References

Stacy, E. W. (1962). A generalization of the gamma distribution. Annals of Mathematical Statistics 33:1187-92.

Prentice, R. L. (1974). A log gamma model and its maximum likelihood estimation. Biometrika 61(3):539-544.

Lawless, J. F. (1980). Inference in the generalized gamma and log gamma distributions. Technometrics 22(3):409-419.

See Also

GenGamma, GenF.orig, GenF, Lognormal, GammaDist, Weibull.

Description

Evaluate baseline time-to-event distribution parameters given covariate values in a flexsurvmix model

Usage

get_basepars(x, newdata, event)

Arguments

Description

Glance accepts a model object and returns a tibble with exactly one row of model summaries.

Usage

## S3 method for class 'flexsurvreg'
glance(x, ...)

Arguments

Value

A one-row tibble containing columns:

• N Number of observations used in fitting

• events Number of events

• censored Number of censored events

• trisk Total length of time-at-risk (i.e. follow-up)

• df Degrees of freedom (i.e. number of estimated parameters)

• logLik Log-likelihood

• AIC Akaike's "An Information Criteria"

• BIC Bayesian Information Criteria

Examples

fitg <- flexsurvreg(formula = Surv(futime, fustat) ~ age, data = ovarian, dist = "gengamma")
glance(fitg)

Description

Density, distribution function, hazards, quantile function and random generation for the Gompertz distribution with unrestricted shape.

Usage

dgompertz(x, shape, rate = 1, log = FALSE)

pgompertz(q, shape, rate = 1, lower.tail = TRUE, log.p = FALSE)

qgompertz(p, shape, rate = 1, lower.tail = TRUE, log.p = FALSE)

rgompertz(n, shape = 1, rate = 1)

hgompertz(x, shape, rate = 1, log = FALSE)

Hgompertz(x, shape, rate = 1, log = FALSE)

Arguments

Details

The Gompertz distribution with shape parameter $a$ and rate parameter $b$ has probability density function

$f(x | a, b) = be^{ax}\exp(-b/a (e^{ax} - 1))$

and hazard

$h(x | a, b) = b e^{ax}$

The hazard is increasing for shape $a>0$ and decreasing for $a<0$. For $a=0$ the Gompertz is equivalent to the exponential distribution with constant hazard and rate $b$.

The probability distribution function is

$F(x | a, b) = 1 - \exp(-b/a (e^{ax} - 1))$

Thus if $a$ is negative, letting $x$ tend to infinity shows that there is a non-zero probability $\exp(b/a)$ of living forever. On these occasions qgompertz and rgompertz will return Inf.

Value

dgompertz gives the density, pgompertz gives the distribution function, qgompertz gives the quantile function, hgompertz gives the hazard function, Hgompertz gives the cumulative hazard function, and rgompertz generates random deviates.

Note

Some implementations of the Gompertz restrict $a$ to be strictly positive, which ensures that the probability of survival decreases to zero as $x$ increases to infinity. The more flexible implementation given here is consistent with streg in Stata.

The functions eha::dgompertz and similar available in the package eha label the parameters the other way round, so that what is called the shape there is called the rate here, and what is called 1 / scale there is called the shape here. The terminology here is consistent with the exponential dexp and Weibull dweibull distributions in R.

Author(s)

Christopher Jackson <[email protected]>

References

Stata Press (2007) Stata release 10 manual: Survival analysis and epidemiological tables.

See Also

dexp

Description

Hazard and cumulative hazard functions for distributions which are built into flexsurv, and whose distribution functions are in base R.

Usage

hexp(x, rate = 1, log = FALSE)

Hexp(x, rate = 1, log = FALSE)

hgamma(x, shape, rate = 1, log = FALSE)

Hgamma(x, shape, rate = 1, log = FALSE)

hlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)

Hlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)

hweibull(x, shape, scale = 1, log = FALSE)

Hweibull(x, shape, scale = 1, log = FALSE)

Arguments

Details

For the exponential and the Weibull these are available analytically, and so are programmed here in numerically stable and efficient forms.

For the gamma and log-normal, these are simply computed as minus the log of the survivor function (cumulative hazard) or the ratio of the density and survivor function (hazard), so are not expected to be robust to extreme values or quick to compute.

Value

Hazard (functions beginning 'h') or cumulative hazard (functions beginning 'H').

Author(s)

Christopher Jackson <[email protected]>

See Also

dexp,dweibull,dgamma,dlnorm,dgompertz,dgengamma,dgenf

Description

Hazard ratio as a function of time from a parametric survival model

Usage

hr_flexsurvreg(
  x,
  newdata = NULL,
  t = NULL,
  start = 0,
  ci = TRUE,
  B = 1000,
  cl = 0.95,
  na.action = na.pass
)

Arguments

Value

A data frame with estimate and confidence limits for the hazard ratio, and one row for each of the times requested in t.

Description

Add fitted survival (or hazard or cumulative hazard) curves from a flexsurvreg model fit to an existing plot.

Usage

## S3 method for class 'flexsurvreg'
lines(
  x,
  newdata = NULL,
  X = NULL,
  type = "survival",
  t = NULL,
  est = TRUE,
  ci = NULL,
  B = 1000,
  cl = 0.95,
  col = "red",
  lty = 1,
  lwd = 2,
  col.ci = NULL,
  lty.ci = 2,
  lwd.ci = 1,
  ...
)

Arguments

Details

Equivalent to plot.flexsurvreg(...,add=TRUE).

Author(s)

C. H. Jackson [email protected]

See Also

flexsurvreg

Description

Density, distribution function, hazards, quantile function and random generation for the log-logistic distribution.

Usage

dllogis(x, shape = 1, scale = 1, log = FALSE)

pllogis(q, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE)

qllogis(p, shape = 1, scale = 1, lower.tail = TRUE, log.p = FALSE)

rllogis(n, shape = 1, scale = 1)

hllogis(x, shape = 1, scale = 1, log = FALSE)

Hllogis(x, shape = 1, scale = 1, log = FALSE)

Arguments

Details

The log-logistic distribution with shape parameter $a>0$ and scale parameter $b>0$ has probability density function

$f(x | a, b) = (a/b) (x/b)^{a-1} / (1 + (x/b)^a)^2$

and hazard

$h(x | a, b) = (a/b) (x/b)^{a-1} / (1 + (x/b)^a)$

for $x>0$. The hazard is decreasing for shape $a\leq 1$, and unimodal for $a > 1$.

The probability distribution function is

$F(x | a, b) = 1 - 1 / (1 + (x/b)^a)$

If $a > 1$, the mean is $b c / sin(c)$, and if $a > 2$ the variance is $b^2 * (2*c/sin(2*c) - c^2/sin(c)^2)$, where $c = \pi/a$, otherwise these are undefined.

Value

dllogis gives the density, pllogis gives the distribution function, qllogis gives the quantile function, hllogis gives the hazard function, Hllogis gives the cumulative hazard function, and rllogis generates random deviates.

Note

Various different parameterisations of this distribution are used. In the one used here, the interpretation of the parameters is the same as in the standard Weibull distribution (dweibull). Like the Weibull, the survivor function is a transformation of $(x/b)^a$ from the non-negative real line to [0,1], but with a different link function. Covariates on $b$ represent time acceleration factors, or ratios of expected survival.

The same parameterisation is also used in eha::dllogis in the eha package.

Author(s)

Christopher Jackson <[email protected]>

References

Stata Press (2007) Stata release 10 manual: Survival analysis and epidemiological tables.

See Also

dweibull

Description

Log likelihood from a flexsurvreg model

Usage

## S3 method for class 'flexsurvreg'
logLik(object, ...)

Arguments

Value

Log-likelihood (numeric) with additional attributes df (degrees of freedom, or number of parameters that were estimated), and number of observations nobs (including observed events and censored observations).

Description

Mean and restricted mean survival time functions for distributions which are built into flexsurv.

Usage

mean_exp(rate = 1)

rmst_exp(t, rate = 1, start = 0)

mean_gamma(shape, rate = 1)

rmst_gamma(t, shape, rate = 1, start = 0)

rmst_genf(t, mu, sigma, Q, P, start = 0)

mean_genf(mu, sigma, Q, P)

rmst_genf.orig(t, mu, sigma, s1, s2, start = 0)

mean_genf.orig(mu, sigma, s1, s2)

rmst_gengamma(t, mu = 0, sigma = 1, Q, start = 0)

mean_gengamma(mu = 0, sigma = 1, Q)

rmst_gengamma.orig(t, shape, scale = 1, k, start = 0)

mean_gengamma.orig(shape, scale = 1, k)

rmst_gompertz(t, shape, rate = 1, start = 0)

mean_gompertz(shape, rate = 1)

mean_llogis(shape = 1, scale = 1)

rmst_llogis(t, shape = 1, scale = 1, start = 0)

mean_lnorm(meanlog = 0, sdlog = 1)

rmst_lnorm(t, meanlog = 0, sdlog = 1, start = 0)

mean_weibull(shape, scale = 1)

rmst_weibull(t, shape, scale = 1, start = 0)

rmst_weibullPH(t, shape, scale = 1, start = 0)

mean_weibullPH(shape, scale = 1)

Arguments

Details

For the exponential, Weibull, log-logistic, lognormal, and gamma, mean survival is provided analytically. Restricted mean survival for the exponential distribution is also provided analytically. Mean and restricted means for other distributions are calculated via numeric integration.

Value

mean survival (functions beginning 'mean') or restricted mean survival (functions beginning 'rmst_').

Author(s)

Christopher Jackson <[email protected]>

See Also

dexp,dweibull,dgamma,dlnorm,dgompertz,dgengamma,dgenf

Description

This returns the mean of each event-specific parametric time-to-event distribution in the mixture model, which is the mean time to event conditionally on that event being the one that happens.

Usage

mean_flexsurvmix(x, newdata = NULL, B = NULL)

Arguments

Value

Mean times to next event conditionally on each alternative event, given the specified covariate values.

Description

Calculate the mean time from the start of the process to a final (or "absorbing") state in a mixture multi-state model. Models with cycles are not supported.

Usage

meanfinal_fmixmsm(x, newdata = NULL, final = FALSE, B = NULL)

Arguments

Value

A data frame of mean times to absorption, by covariate values and pathway (or by final state)

Description

Returns a list of data frames, with each component containing the data that were used for the model fit for that mixture component.

Usage

## S3 method for class 'flexsurvmix'
model.frame(formula, ...)

Arguments

Value

A list of data frames

Description

Extract the data from a model fitted with flexsurvreg.

Usage

## S3 method for class 'flexsurvreg'
model.frame(formula, ...)

## S3 method for class 'flexsurvreg'
model.matrix(object, par = NULL, ...)

Arguments

Value

model.frame returns a data frame with all the original variables used for the model fit.

model.matrix returns a design matrix for a part of the model that includes covariates. The required part is indicated by the "par" argument (see above).

Author(s)

C. H. Jackson [email protected]

See Also

flexsurvreg, model.frame, model.matrix.

Description

Cumulative transition-specific intensity/hazard functions for fully-parametric multi-state or competing risks models, using a piecewise-constant approximation that will allow prediction using the functions in the mstate package.

Usage

msfit.flexsurvreg(
  object,
  t,
  newdata = NULL,
  variance = TRUE,
  tvar = "trans",
  trans,
  B = 1000
)

Arguments

Value

An object of class "msfit", in the same form as the objects used in the mstate package. The msfit method from mstate returns the equivalent cumulative intensities for Cox regression models fitted with coxph.

Author(s)

C. H. Jackson [email protected]

References

Liesbeth C. de Wreede, Marta Fiocco, Hein Putter (2011). mstate: An R Package for the Analysis of Competing Risks and Multi-State Models. Journal of Statistical Software, 38(7), 1-30. doi:10.18637/jss.v038.i07

Mandel, M. (2013). "Simulation based confidence intervals for functions with complicated derivatives." The American Statistician 67(2):76-81

See Also

flexsurv provides alternative functions designed specifically for predicting from parametric multi-state models without calling mstate. These include pmatrix.fs and pmatrix.simfs for the transition probability matrix, and totlos.fs and totlos.simfs for expected total lengths of stay in states. These are generally more efficient than going via mstate.

Examples

## 3 state illness-death model for bronchiolitis obliterans
## Compare clock-reset / semi-Markov multi-state models

## Simple exponential model (reduces to Markov)

bexp <- flexsurvreg(Surv(years, status) ~ trans,
                    data=bosms3, dist="exp")
tmat <- rbind(c(NA,1,2),c(NA,NA,3),c(NA,NA,NA))
mexp <- msfit.flexsurvreg(bexp, t=seq(0,12,by=0.1),
                          trans=tmat, tvar="trans", variance=FALSE)

## Cox semi-parametric model within each transition

bcox <- coxph(Surv(years, status) ~ strata(trans), data=bosms3)

if (require("mstate")){

mcox <- mstate::msfit(bcox, trans=tmat)

## Flexible parametric spline-based model 

bspl <- flexsurvspline(Surv(years, status) ~ trans + gamma1(trans),
                       data=bosms3, k=3)
mspl <- msfit.flexsurvreg(bspl, t=seq(0,12,by=0.1),
                         trans=tmat, tvar="trans", variance=FALSE)

## Compare fit: exponential model is OK but the spline is better

plot(mcox, lwd=1, xlim=c(0, 12), ylim=c(0,4))
cols <- c("black","red","green")
for (i in 1:3){
    lines(mexp$Haz$time[mexp$Haz$trans==i], mexp$Haz$Haz[mexp$Haz$trans==i],
             col=cols[i], lwd=2, lty=2)
    lines(mspl$Haz$time[mspl$Haz$trans==i], mspl$Haz$Haz[mspl$Haz$trans==i],
             col=cols[i], lwd=3)
}
legend("topright", lwd=c(1,2,3), lty=c(1,2,1),
   c("Cox", "Exponential", "Flexible parametric"), bty="n")

}

## Fit a list of models, one for each transition
## More computationally efficient, but only valid if parameters
## are different between transitions.

## Not run: 
bexp.list <- vector(3, mode="list")
for (i in 1:3) { 
  bexp.list[[i]] <- flexsurvreg(Surv(years, status) ~ 1, subset=(trans==i),
                                data=bosms3, dist="exp")
}

## The list of models can be passed to this and other functions,
## as if it were a single multi-state model. 

msfit.flexsurvreg(bexp.list, t=seq(0,12,by=0.1), trans=tmat)

## End(Not run)