Introduction

Parametric cure models are a type of parametric survival model model in which it is assumed that there are a proportion of subjects who will not experience the event. In a mixture cure model, these ‘cured’ and ‘uncured’ subjects are modeled separately, with the cured individuals subject to no excess risk and the uncured individuals subject to excess risk modeled using a parametric survival distribution. In a non-mixture model, a parametric survival distribution is scaled such that survival asymptotically approaches the cure fraction.

Mixture Cure Model

The following code fits a mixture cure model to the bc dataset from flexsurv using a Weibull distribution and a logistic link function for the cure fraction:

library(flexsurvcure)
cure_model <- flexsurvcure(Surv(rectime, censrec)~group, data=bc, link="logistic", dist="weibullPH", mixture=T)
print(cure_model)

## Call:
## flexsurvcure(formula = Surv(rectime, censrec) ~ group, data = bc, 
##     dist = "weibullPH", link = "logistic", mixture = T)
## 
## Estimates: 
##              data mean  est        L95%       U95%       se         exp(est) 
## theta               NA   6.73e-01   5.84e-01   7.52e-01         NA         NA
## shape               NA   1.55e+00   1.38e+00   1.74e+00   9.07e-02         NA
## scale               NA   1.61e-05   5.10e-06   5.11e-05   9.50e-06         NA
## groupMedium   3.34e-01  -1.23e+00  -1.74e+00  -7.09e-01   2.64e-01   2.93e-01
## groupPoor     3.32e-01  -3.48e+00  -5.56e+00  -1.40e+00   1.06e+00   3.08e-02
##              L95%       U95%     
## theta               NA         NA
## shape               NA         NA
## scale               NA         NA
## groupMedium   1.75e-01   4.92e-01
## groupPoor     3.85e-03   2.47e-01
## 
## N = 686,  Events: 299,  Censored: 387
## Total time at risk: 771400
## Log-likelihood = -2580.012, df = 5
## AIC = 5170.025

Model results can be displayed graphically using the plot S3 method:

plot(cure_model)

Predicted survival probabilities can also be generated using the summary S3 method:

summary(cure_model, t=seq(from=0,to=3000,by=1000), type="survival", tidy=T)

##    time        est        lcl       ucl  group
## 1     0 1.00000000 1.00000000 1.0000000   Good
## 2  1000 0.83251996 0.78865790 0.8685350   Good
## 3  2000 0.71313032 0.64476832 0.7770906   Good
## 4  3000 0.67959584 0.60453748 0.7562965   Good
## 5     0 1.00000000 1.00000000 1.0000000 Medium
## 6  1000 0.68057306 0.62450014 0.7301375 Medium
## 7  2000 0.45286672 0.38259264 0.5312609 Medium
## 8  3000 0.38890796 0.30637734 0.4867416 Medium
## 9     0 1.00000000 1.00000000 1.0000000   Poor
## 10 1000 0.51803653 0.46954409 0.6370861   Poor
## 11 2000 0.17446456 0.12979119 0.4138290   Poor
## 12 3000 0.07796116 0.03881832 0.3724180   Poor

More complex models may be fitted by adding covariates to the parametric distribution used to model the uncured individuals. This is done by passing a list of formula, named according to the parameters affected, through the anc argument:

cure_model_complex <- flexsurvcure(Surv(rectime, censrec)~group, data=bc, link="logistic", dist="weibullPH", mixture=T, anc=list(scale=~group))
print(cure_model_complex)

## Call:
## flexsurvcure(formula = Surv(rectime, censrec) ~ group, data = bc, 
##     dist = "weibullPH", link = "logistic", mixture = T, anc = list(scale = ~group))
## 
## Estimates: 
##                     data mean  est        L95%       U95%       se       
## theta                      NA   2.62e-02         NA         NA         NA
## shape                      NA   1.49e+00         NA         NA         NA
## scale                      NA   4.65e-06         NA         NA         NA
## groupMedium          3.34e-01  -1.78e+00         NA         NA         NA
## groupPoor            3.32e-01   2.07e+00         NA         NA         NA
## scale(groupMedium)   3.34e-01   8.29e-01         NA         NA         NA
## scale(groupPoor)     3.32e-01   2.11e+00         NA         NA         NA
##                     exp(est)   L95%       U95%     
## theta                      NA         NA         NA
## shape                      NA         NA         NA
## scale                      NA         NA         NA
## groupMedium          1.69e-01         NA         NA
## groupPoor            7.96e+00         NA         NA
## scale(groupMedium)   2.29e+00         NA         NA
## scale(groupPoor)     8.22e+00         NA         NA
## 
## N = 686,  Events: 299,  Censored: 387
## Total time at risk: 771400
## Log-likelihood = -2571.857, df = 7
## AIC = 5157.713

plot(cure_model_complex)

Non-Mixture Cure Model

Non-mixture cure models can be fit by passing mixture=FALSE to flexsurvcure:

library(flexsurvcure)
cure_model_nmix <- flexsurvcure(Surv(rectime, censrec)~group, data=bc, link="loglog", dist="weibullPH", mixture=F)
print(cure_model_nmix)

## Call:
## flexsurvcure(formula = Surv(rectime, censrec) ~ group, data = bc, 
##     dist = "weibullPH", link = "loglog", mixture = F)
## 
## Estimates: 
##              data mean  est       L95%      U95%      se        exp(est)
## theta              NA   6.35e-01  7.31e-01  5.17e-01        NA        NA
## shape              NA   1.72e+00  1.53e+00  1.92e+00  1.01e-01        NA
## scale              NA   3.07e-06  9.19e-07  1.03e-05  1.89e-06        NA
## groupMedium  3.34e-01   8.35e-01  4.99e-01  1.17e+00  1.71e-01  2.31e+00
## groupPoor    3.32e-01   1.63e+00  1.31e+00  1.95e+00  1.64e-01  5.09e+00
##              L95%      U95%    
## theta              NA        NA
## shape              NA        NA
## scale              NA        NA
## groupMedium  1.65e+00  3.22e+00
## groupPoor    3.69e+00  7.02e+00
## 
## N = 686,  Events: 299,  Censored: 387
## Total time at risk: 771400
## Log-likelihood = -2567.8, df = 5
## AIC = 5145.6