Parametric survival analysis using the flexsurvPlus package: understanding the theory

Introduction

Parametric survival models are often the preferred method of extrapolating survival data for use in economic models. The National Institute for Health and Care Excellence (NICE) Decision Support Unit (DSU) technical support document (TSD) 14 recommends that the Exponential, Weibull, Gompertz, log-logistic, log normal and Generalized Gamma parametric models should all be considered.[1] More recently, NICE also discusses more flexible models in NICE DSU TSD 21, however, more these models are not in the scope of this package.[2] The Canadian Agency for Drugs and Technologies in Health (CADTH) additionally specifies that the Gamma distribution must also be considered. This document therefore details the characteristics of each of these distributions and demonstrates how the parameters from each distribution, outputted using the flexsurvPlus package, can be implemented within an economic model.[3] The Generalized F distribution is not commonly used, however it has been included in this package in case more flexible modeling is required.

The flexsurvPlus package allows the inclusion of a treatment effect in the following three ways:

  • Separate models - Models fitted to each treatment arm separately

    • Similarly, a separate model can be fitted to data from a single treatment arm
  • Independent shape models - Models fitted to both/all treatment arms including a treatment covariate to model the effect of treatment on both the scale and shape parameter(s) of the distribution.

  • Common shape models - Models fitted to both/all treatment arms including a treatment covariate to model the effect of treatment on the scale parameter of the distribution. The shape parameter(s) of the distribution is common across treatments which reflects an assumption of proportional hazards or accelerated failure time between treatments depending on the distribution

This document details the parameters for each of these three model types. A separate vignette; “Fitting parametric survival models in R” details how to perform the analyses in R using the flexsurvPlus package.

Properties of the distributions

The properties of each parametric distribution are summarized in the Table below.

By default, covariates are placed on the “location” parameter of the distribution, typically the “scale” or “rate” parameter, through a linear model, or a log-linear model if this parameter must be positive. This gives an accelerated failure time model or a proportional hazards model (see Table below) depending on how the distribution is parametrized.The other parameters are ancillary parameters that determine the shape, variance, or higher moments of the distribution. These parameters impact the hazard function, which can take a variety of shapes depending on the distribution:[4,5]

  • the exponential distribution only supports a constant hazard;
  • the Weibull, Gompertz, and gamma distributions support monotonically increasing and decreasing hazards;
  • the log-logistic and log-normal distributions support arc-shaped and monotonically decreasing hazards;
  • and the generalized gamma distribution supports an arc-shaped, bathtub-shaped, monotonically increasing, and monotonically decreasing hazards.

The log-hazard function of the Weibull model is linear with respect to the log of time. The log-hazard function of the Gompertz model is linear with respect to time.

Special cases

  • The Weibull distribution reduces to an exponential distribution when the Weibull shape parameter equals 1 (a = 1) and the exponential rate parameter is 1 divided by the the Weibull scale parameter ($\frac{1}{b}$).
  • The Gompertz distribution reduces to an exponential distribution when the Gompertz shape parameter equals 0 (a = 0).
  • The Generalized Gamma distribution reduces to an exponential model (Q = σ = 1), Weibull model (Q = 1), gamma model (Q = σ) and log-normal model (Q = 0).
  • The Generalized F distribution reduces to the Generalized Gamma distribution (P = 1)

Initial diagnostics

Before fitting any models, the following initial diagnostic plots are helpful to understand which distribution might be the best fit to the data being modeled:

  • Kaplan–Meier plot
    • To display the survival function for each treatment
  • Log-cumulative hazard plot
    • Visual inspection of a log-cumulative hazard plot allows assessment of the proportional hazards assumption. In addition, they provide an alternative way to assess changes in hazards over time.
  • Smoothed hazard plot
    • Does the shape of the hazard plot indicate that a certain distribution is more appropriate than another? Possible shapes of the hazard function for each of the distributions that can be fitted with the flexsurvPlus package are described in Table X below. For example, if the hazard is constant, an exponential distribution might be the most appropriate fit to the data.
Properties of common parametric survival distributions
Distribution Number of parameters Location parameter Ancillary parameter Survival function: S(t) Hazard function: h(t) Possible shapes of the hazard function Treatment effect model assumption
Exponential 1: b rate (b) NA ebt b Constant PH and AFT
Weibull 2: a, b scale (b) shape (a) $e^{(\frac{-t}{b})^{a}}$ $a\frac{t}{b}^{a -1}$ Constant, monotonically increasing, monotonically decreasing AFT
Gompertz 2: a, b rate (b) shape (a) $e^{\frac{-b}{a} (e^{a t}-1)}$ beat Constant, monotonically increasing, monotonically decreasing PH
Log-logistic 2: a, b scale (b) shape (a) $\frac{1}{1 + (\frac{t}{b})^a}$ $\frac{\frac{a}{b}(\frac{t}{b})^{a - 1}}{1 + (\frac{t}{b})^a}$ Monotonically decreasing, initial increase the decrease AFT
Log-normal 2: μ, σ meanlog (μ) sdlog (σ) $1 - \Phi \frac{log t - \mu}{\sigma}$ $\frac{\frac{1}{\sigma \sqrt{2 \pi}}t^{-1}e[-\frac{(log t - \mu)^2}{2 \sigma^2}]}{S(t)}$ Initially increases then decreases AFT
Gamma 2: a, b rate (b) shape (a) $1 - \frac{\gamma(a, bt)}{\Gamma(a)}$ $\frac{b^at^{a -1}e^{-bt}}{\Gamma(a) S(t)}$ Constant, monotonically increasing, monotonically decreasing AFT
Generalized gamma 3: μ, σ, Q mu (μ) sigma (σ), Q $1 - \frac{\gamma(Q^{-2}, Q^{-2}e^{Q \frac{log(t) – \mu}{\sigma}})}{\Gamma(Q^{-2})}; \ if \ Q \neq 0$  $1 - \Phi(\frac{logt – \mu}{\sigma}); \ if \ Q = 0$ $\frac{|Q|(Q^{-2})^{Q^{-2}}}{\sigma t \Gamma(Q^{-2}) S(t)} \exp\left[Q^{-2}\left(Q\frac{logt – \mu}{\sigma}-e^{Q\frac{logt – \mu}{\sigma}}\right)\right]$ Constant, monotonically increasing, monotonically decreasing, bathtub, arc-shaped AFT
Generalized F 4: μ, σ, Q, P mu (μ) sigma (σ), Q, P $$ \int_{0}^{s_{1}(s_{2}+s_{1}e^{w})^{-1}} \frac{x^{s_{2}-1}(1-x)^{s_{1}-1}}{B(s_{2}, s_{1})} \ dx $$ $$ \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})S(t)}$$ Decreasing, arc-shaped AFT

Exponential

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “exp”

uses the flexsurv package to fit two exponential models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "exp") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "exp")

The estimated rate parameter from each of these models can be then used within the following Excel formula to estimate survival at different time points:

EXP(-rate*time)

This is equivalent to running the following code in R:

summary( flexsurvreg exponential model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “exp”

uses the flexsurv package to fit one exponential model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "exp")

This model outputs a rate parameter and treatment coefficient. Rate parameters for each treatment are estimated in run_PSM as follows:

Rate parameter for intervention and reference treatments from a common shape model
Rate (intervention treatment) Rate (reference treatment)
Rate parameter + treatment coefficient Rate parameter

Survival can then be estimated in Excel using these two rate parameters, as for the separate models, with formula:

EXP(-rate*time)

This is equivalent to running the following code in R:

summary( flexsurvreg exponential model , type="survival", B=0, t=time)

Weibull

The Weibull model may be parametrized as either a proportional hazards model or as an accelerated failure time model; however, the flexsurv package parametrizes the Weibull as an accelerated failure time model.

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “weibull”

uses the flexsurv package to fit two weibull models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "weibull") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "weibull")

Estimated shape and scale parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

1-WEIBULL(time,shape,scale,TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg weibull model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “weibull”

uses the flexsurv package to fit one weibull model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "weibull")

This model outputs shape and scale parameters and a treatment coefficient on the scale parameter. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Shape and scale parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
Shape Shape parameter Shape parameter
Scale Scale parameter + treatment coefficient Scale parameter

Survival can then be estimated in Excel using the shape and scale parameters specific to each treatment with formula:

1-WEIBULL(time,shape,scale,TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg weibull model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “weibull”

uses the flexsurv package to fit one weibull model with a covariate for treatment on both the scale and shape parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + shape(ARM), data = data with a treatment variable named “ARM” , dist = "weibull")

This model outputs shape and scale parameters and a treatment coefficient for both shape and scale. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Shape and scale parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
Shape Shape parameter + treatment coefficient for shape Shape parameter
Scale Scale parameter + treatment coefficient for scale Scale parameter

Survival can then be estimated in Excel using the shape and scale parameters specific to each treatment with formula:

1-WEIBULL(time,shape,scale,TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg weibull model , type="survival", B=0, t=time)

Gompertz

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “gompertz”

uses the flexsurv package to fit two gompertz models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "gompertz") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "gompertz")

Estimation of the rate and shape parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

EXP((-1/shape) * rate * (EXP(shape*time)-1))

This is equivalent to running the following code in R:

summary( flexsurvreg gompertz model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “gompertz”

uses the flexsurv package to fit one gompertz model with a covariate for treatment on the rate parameter. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "gompertz")

This model outputs rate and shape parameters and a treatment coefficient. Rate and shape parameters for each treatment are estimated in run_PSM as follows:

Rate and shape parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
Rate Rate parameter + treatment coefficient Rate parameter
Shape Shape parameter Shape parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

EXP((-1/shape) * rate * (EXP(shape*time)-1))

This is equivalent to running the following code in R:

summary( flexsurvreg gompertz model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “gompertz”

uses the flexsurv package to fit one gompertz model with a covariate for treatment on both the rate and shape parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + shape(ARM), data = data with a treatment variable named “ARM” , dist = "gompertz")

This model outputs rate and shape parameters and a treatment coefficient for both rate and shape. Rate and shape parameters for each treatment are estimated in run_PSM as follows:

Rate and shape parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
Rate Rate parameter + treatment coefficient on rate Rate parameter
Shape Shape parameter + treatment coeffecient on shape Shape parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

EXP((-1/shape) * rate * (EXP(shape*time)-1))

This is equivalent to running the following code in R:

summary( flexsurvreg gompertz model , type="survival", B=0, t=time)

Log-logistic

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “llogis”

uses the flexsurv package to fit two log-logistic models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "llogis") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "llogis")

Estimation of the shape and scale parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

1 /(1+((time / scale)^shape))

This is equivalent to running the following code in R:

summary( flexsurvreg log-logistic model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “llogis”

uses the flexsurv package to fit one log-logistic model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "llogis")

This model outputs shape and scale parameters and a treatment coefficient. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Shape and scale parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
Shape Shape parameter Shape parameter
Scale Scale parameter + treatment coefficient Scale parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

1 /(1+((time / scale)^shape))

This is equivalent to running the following code in R:

summary( flexsurvreg log-logistic model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “llogis”

uses the flexsurv package to fit one log-logistic model with a covariate for treatment on both the shape and scale parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + shape(ARM), data = data with a treatment variable named “ARM” , dist = "llogis")

This model outputs shape and scale parameters and a treatment coefficient for both shape and scale. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Shape and scale parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
Shape Shape parameter + treatment coefficient on shape Shape parameter
Scale Scale parameter + treatment coefficient on scale Scale parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

1 /(1+((time / scale)^shape))

This is equivalent to running the following code in R:

summary( flexsurvreg log-logistic model , type="survival", B=0, t=time)

Log-normal

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “lnorm”

uses the flexsurv package to fit two log-normal models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "lnorm") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "lnorm")

Estimation of the meanlog and sdlog parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

1-LOGNORMDIST(time, meanlog, sdlog))

This is equivalent to running the following code in R:

summary( flexsurvreg log-normal model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “lnorm”

uses the flexsurv package to fit one log-normal model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "lnorm")

This model outputs meanlog and sdlog parameters and a treatment coefficient. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Meanlog and sdlog parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
Meanlog Meanlog parameter + treatment coefficient Meanlog parameter
Sdlog Sdlog parameter Sdlog parameter

Survival can then be estimated in Excel using the meanlog and sdlog parameters specific to each treatment with formula:

1-LOGNORMDIST(time, meanlog, sdlog))

This is equivalent to running the following code in R:

summary( flexsurvreg log-normal model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “lnorm”

uses the flexsurv package to fit one log-normal model with a covariate for treatment on both the meanlog and sdlog parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event==1) ~ ARM + sdlog(ARM), data = data with a treatment variable named “ARM” , dist = "lnorm")

This model outputs meanlog and sdlog parameters and a treatment coefficient for both meanlog and sdlog. Shape and scale parameters for each treatment are estimated in run_PSM as follows:

Meanlog and sdlog parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
Meanlog Meanlog parameter + treatment coefficient on meanlog Meanlog parameter
Sdlog Sdlog parameter + treatment coefficient on sdlog Sdlog parameter

Survival can then be estimated in Excel using the meanlog and sdlog parameters specific to each treatment with formula:

1-LOGNORMDIST(time, meanlog, sdlog))

This is equivalent to running the following code in R:

summary( flexsurvreg log-normal model , type="survival", B=0, t=time)

Gamma

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “gamma”

uses the flexsurv package to fit two gamma models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "gamma") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "gamma")

Estimation of the rate and shape parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

1-GAMMA.DIST(time, shape, 1/(rate),TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg gamma model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “gamma”

uses the flexsurv package to fit one gamma model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "gamma")

This model outputs rate and shape parameters and a treatment coefficient. Rate and shape parameters for each treatment are estimated in run_PSM as follows:

Rate and shape parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
Rate Rate parameter + treatment coefficient Rate parameter
Shape Shape parameter Shape parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

1-GAMMA.DIST(time, shape, 1/(rate),TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg gamma model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “gamma”

uses the flexsurv package to fit one gamma model with a covariate for treatment on both the rate and shape parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + shape(ARM), data = data with a treatment variable named “ARM” , dist = "gamma")

This model outputs rate and shape parameters and a treatment coefficient for both rate and shape. Rate and shape parameters for each treatment are estimated in run_PSM as follows:

Rate and shape parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
Rate Rate parameter + treatment coefficient on rate Rate parameter
Shape Shape parameter + treatment coeffecient on shape Shape parameter

Survival can then be estimated in Excel using the rate and shape parameters specific to each treatment with formula:

1-GAMMA.DIST(time, shape, 1/(rate),TRUE)

This is equivalent to running the following code in R:

summary( flexsurvreg gamma model , type="survival", B=0, t=time)

Generalized gamma

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “gengamma”

uses the flexsurv package to fit two generalized gamma models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "gengamma") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "gengamma")

Estimation of the mu, sigma and Q parameters from each of these models can be then used within the following Excel formula to estimate survival at different time points:

IF(Q<0,GAMMADIST((-Q^-2) * EXP(-Q* -((LN(time)-(mu))/sigma)),-Q^-2,1,1),1-GAMMADIST((-Q^-2) * EXP(-Q * -((LN(time)-(mu))/sigma)),-Q^-2,1,1))

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized Gamma model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “gengamma”

uses the flexsurv package to fit one Generalized Gamma model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "gengamma")

This model outputs mu, sigma and Q parameters and a treatment coefficient. Mu, sigma and Q parameters for each treatment are estimated in run_PSM as follows:

Mu, sigma and Q parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
mu mu coefficient + treatment coefficient mu coefficient
Scale sigma coefficient sigma coefficient
Q Q coefficient Q coefficient

Survival can then be estimated in Excel using the mu, sigma and Q parameters specific to each treatment with formula:

IF(Q<0,GAMMADIST((-Q^-2) * EXP(-Q* -((LN(time)-(mu))/sigma)),-Q^-2,1,1),1-GAMMADIST((-Q^-2) * EXP(-Q * -((LN(time)-(mu))/sigma)),-Q^-2,1,1))

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized Gamma model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “gengamma”

uses the flexsurv package to fit one Generalized Gamma model with a covariate for treatment on both the mu, sigma and Q parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + sigma(ARM) + Q(ARM), data = data with a treatment variable named “ARM” , dist = "gengamma")

This model outputs mu, sigma and Q parameters and a treatment coefficient for mu, sigma and Q. Parameters for each treatment are estimated in run_PSM as follows:

Mu, sigma and Q parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
mu mu coefficient + treatment coefficient for mu mu coefficient
Scale sigma coefficient + treatment coefficient for sigma sigma coefficient
Q Q coefficient + treatment coefficient for Q Q coefficient

Survival can then be estimated in Excel using the mu, sigma and Q parameters specific to each treatment with formula:

IF(Q<0,GAMMADIST((-Q^-2) * EXP(-Q* -((LN(time)-(mu))/sigma)),-Q^-2,1,1),1-GAMMADIST((-Q^-2) * EXP(-Q * -((LN(time)-(mu))/sigma)),-Q^-2,1,1))

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized Gamma model , type="survival", B=0, t=time)

Generalized F

Separate models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Separate”
  • distr = “genf”

uses the flexsurv package to fit two generalized F models (one to data for the intervention treatment and one to data for the reference treatment). The R code for the fitted models is as follows:

flexsurvreg(formula = Surv(Time, Event)~1, data = data for intervention treatment , dist = "genf") flexsurvreg(formula = Surv(Time, Event)~1, data = data for reference treatment , dist = "genf")

Estimation of the mu, sigma, Q and P parameters from each of these models can be then used within the following Excel formulas to estimate survival at different time points:

BETA.DIST(h_m2/(h_m2 + h_m1*time^(h_d/sigma)/EXP(h_d/sigma*mu)),h_m2,h_m1,TRUE)

  • h_m1 = 2/(Q*Q + 2*P + Q*SQRT(Q*Q + 2*P))
  • h_m2 = 2/(Q*Q + 2*P - Q*SQRT(Q*Q + 2*P))
  • h_d = SQRT(Q*Q + 2*P)

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized F model , type="survival", B=0, t=time)

Common shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Common shape”
  • distr = “genf”

uses the flexsurv package to fit one Generalized F model with a covariate for treatment. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event)~ARM, data = data with a treatment variable named “ARM” , dist = "genf")

This model outputs mu, sigma Q and P parameters and a treatment coefficient. Mu, sigma Q and P parameters for each treatment are estimated in run_PSM as follows:

Mu, sigma, Q and P parameters for intervention and reference treatments from a common shape model
Parameter Intervention Reference
mu mu coefficient + treatment coefficient mu coefficient
Sigma sigma coefficient sigma coefficient
Q Q coefficient Q coefficient
P P coefficient P coefficient

Survival can then be estimated in Excel using the mu, sigma, Q and P parameters specific to each treatment with formula:

BETA.DIST(h_m2/(h_m2 + h_m1*time^(h_d/sigma)/EXP(h_d/sigma*mu)),h_m2,h_m1,TRUE)

  • h_m1 = 2/(Q*Q + 2*P + Q*SQRT(Q*Q + 2*P))
  • h_m2 = 2/(Q*Q + 2*P - Q*SQRT(Q*Q + 2*P))
  • h_d = SQRT(Q*Q + 2*P)

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized f model , type="survival", B=0, t=time)

Independent shape models

The run_PSM function in the flexsurvPlus package with the following arguments:

  • model.type = “Independent shape”
  • distr = “genf”

uses the flexsurv package to fit one Generalized F model with a covariate for treatment on both the mu, sigma, Q and P parameters. The R code for the fitted model is as follows:

flexsurvreg(formula = Surv(Time, Event) ~ ARM + sigma(ARM) + Q(ARM) + P(ARM), data = data with a treatment variable named “ARM” , dist = "genf")

This model outputs mu, sigma, Q and P parameters and a treatment coefficient for mu, sigma, Q and P. Parameters for each treatment are estimated in run_PSM as follows:

Mu, sigma, Q and P parameters for intervention and reference treatments from an independent shape model
Parameter Intervention Reference
mu mu coefficient + treatment coefficient for mu mu coefficient
Scale sigma coefficient + treatment coefficient for sigma sigma coefficient
Q Q coefficient + treatment coefficient for Q Q coefficient
P P coefficient + treatment coefficient for P P coefficient

Survival can then be estimated in Excel using the mu, sigma, Q and P parameters specific to each treatment with formula:

BETA.DIST(h_m2/(h_m2 + h_m1*time^(h_d/sigma)/EXP(h_d/sigma*mu)),h_m2,h_m1,TRUE)

  • h_m1 = 2/(Q*Q + 2*P + Q*SQRT(Q*Q + 2*P))
  • h_m2 = 2/(Q*Q + 2*P - Q*SQRT(Q*Q + 2*P))
  • h_d = SQRT(Q*Q + 2*P)

This is equivalent to running the following code in R:

summary( flexsurvreg Generalized F model , type="survival", B=0, t=time)

References

[1]
Latimer N. NICE DSU technical support document 14: Survival analysis for economic evaluations alongside clinical trials-extrapolation with patient-level data. Sheffield: Report by the Decision Support Unit 2011;2013.
[2]
Rutherford MJ, Lambert PC, Sweeting MJ, Pennington B, Crowther MJ, Abrams KR, et al. NICE DSU TECHNICAL SUPPORT DOCUMENT 21: Flexible methods for survival analysis 2020.
[3]
CADTH. Procedures for CADTH drug reimbursement reviews 2020.
[4]
Incerti D. Parametric survival modeling n.d.
[5]
Cox C. The generalized f distribution: An umbrella for parametric survival analysis. Statistics in Medicine 2008;27:4301–12.