R/curesurv.R

In curesurv: Mixture and Non Mixture Parametric Cure Models to Estimate Cure Indicators

#' @title Fitting cure models using curesurv
#'
#' @description Fits the non-mixture cure model proposed by Boussari et \emph{al}. (2020),
#' or mixture cure model such as proposed by De Angelis et \emph{al}. (1999) with
#' the possibility to correct the background mortality as proposed by Phillips et \emph{al}. (2002)
#' in the net survival framework.
#'
#' ## Non-mixture cure model
#' ### The Boussari model
#'
#' This model allows for direct estimation of time-to-null-excess-hazard which can be interpreted as time-to-cure.
#' The parametrization offers various link functions for the covariates effects on the time-to-null-excess-hazard:
#' \ifelse{html}{\out{&tau;(z<sub>k</sub>) = g(&tau;<sub>0</sub> + z<sub>k</sub> &tau;<sub>k</sub>)}}{\eqn{\tau(z_k)=g(\tau_0+z_k\tau_k)}}.
#' If `link_tau=linear`, then g is the identity function. If `link_tau=loglinear` then g is the exponential function.
#' In this model, the cure proportion is expressed as:
#' \ifelse{html}{\out{&pi;(z;&theta;) = exp(&minus;g(&tau;<sub>0</sub> + z<sub>k</sub> &tau;<sub>k</sub>) \text{Beta}((\alpha_{0} + Z_{k} \alpha_{k}), \beta))}}{\eqn{\pi(z;\theta) = \exp(-g(\tau_0 + z_k \tau_k) \text{Beta}((\alpha_{0} + Z_{k} \alpha_{k}), \beta))}}.
#'
#' ## Mixture cure model
#'
#' The user can choose the survival function modeling the uncured patients net
#' survival among Weibull (default) and exponentiated Weibull. The parametrization
#' for weibull distribution is \ifelse{html}{\out{S<sub>u</sub>(t) = exp(&#8722;&#955;t<sup>&#947;</sup>)<sup>exp(&#948;Z)</sup>}}{\eqn{S_{u}(t) = (\exp\{-\lambda * (t) ^{\gamma}\}) ^ {\exp(\{\delta Z\})}}{ASCII}}.
#' The related hazard function is expressed as: \ifelse{html}{\out{&#955;<sub>u</sub>(t) = &#947;&#955;t<sup>&#947;&#8722; 1</sup> exp(&#948;z)}}{\eqn{\\lambda_{u}(t) = \\gamma\\lambda t^ {\\gamma - 1} \\exp(\\delta z)}{ASCII}}
#' The net survival and the excess hazard functions can be respectively expressed
#' as \ifelse{html}{\out{S<sub>E</sub>(t) = &#960;(z;&#946;) + (1 &#8722; &#960;(z;&#946;)) * S<sub>u</sub>(t)}}{\eqn{S_{E}(t) = \pi(z;\beta) + (1 - \pi(z;\beta)) S_{u}(t)}{ASCII}}.
#' and \ifelse{html}{\out{&#955;<sub>E</sub>(t)=[((1&#8722;&#960;(z;&#946;))fu(t))/(&#960;(z;&#946;) + (1 &#8722; &#960;(z;&#946;))S<sub>u</sub>(t))]}}{\eqn{\lambda_{E}(t)=\frac{(1-\pi(z;\beta))fu(t)}{\pi(z;\beta) + (1 - \pi(z;\beta))S_{u}(t)}}{ASCII}},
#' with \ifelse{html}{\out{&#960;(z;&#946;) = [1/((1 + exp(&#8722;[&#946;<sub>0</sub> + Z&#946;])))]}}{\eqn{\pi(z;\beta) =  \frac{1}{(1 + exp(-[\beta_0 + Z\beta]))}}{ASCII}}.
#'
#' ## Correction of background mortality
#'
#' Usually, in the net survival framework the expected hazard is directly obtained from life tables.
#' However some patients in cancer registries can have some factors impacting their expected mortality rates (such as comorbidities, deprivation)
#' that are not always accounted #'  for in the available life tables, and there is a need to account for this problem.
#' The correction proposed by Phillips et \emph{al} (2002) assumes that \ifelse{html}{\out{&#955;<sub>exp</sub>(t,z)=&#945;&#955;<sub>pop</sub>(t,z<sub>k</sub>)}}{\eqn{\lambda_{exp}(t,z)=\alpha\lambda_{pop}(t,z_k)}{ASCII}}
#' with \ifelse{html}{\out{&#955;<sub>exp</sub>(t,z)}}{\eqn{\lambda_{exp}(t,z)}{ASCII}} the patient expected hazard and \ifelse{html}{\out{&#955;<sub>pop</sub>(t,z<sub>k</sub>)}}{\eqn{\lambda_{pop}(t,z_k)}{ASCII}} the population hazard obtained from life table.
#'
#' @note Note that all these models can be fitted in the overall survival
#' setting.
#'
#'
#'
#' @param formula a formula object of the \code{\link[survival]{Surv}} function with the
#' response on the left of a \code{~} operator and the terms on the right. The
#' response must be a survival object as returned by the \code{\link[survival]{Surv}}
#' function (\code{time} in first and \code{status} in second).
#' @note \code{time} is OBLIGATORY in years
#'
#'
#'
#' @param data a data frame in which to interpret the variables named in the
#' formula
#'
#'
#'
#' @param pophaz corresponds to the name of the column in the data representing
#' the values of the population instantaneous mortality rates. If the
#' \code{pophaz} argument is not specified, overall survival is fitted.
#'
#'
#'
#' @param cumpophaz corresponds to the name of the column in the data
#'  representing the values of the instantaneous population cumulative mortality
#'  rates. If not specified, the model cannot be compared with model with
#'  \code{pophaz.alpha = TRUE} using AIC.
#'
#'
#'
#'
#' @param model To fit a mixture model, specify \code{model = "mixture"}.
#' To fit Time-To-Null Excess Hazard model the argument is \code{model = "tneh"}.
#'
#'
#'
#'
#' @param pophaz.alpha to be specified if user want an excess hazard model
#' with correction of mortality rates by a scale parameter
#'
#'
#'
#'
#'
#' @param dist For mixture model, it corresponds to the function used to fit
#' the uncured patients survival. By default, \code{("weib")} is used. Another
#' option is the exponentiated Weibull function \code{("eweib")}. For
#' non-mixture models, this argument corresponds to the name of the model. By
#' default, \code{("tneh")} is used to fit the time to null excess hazard model
#' proposed by Boussari et al..
#'
#'
#' @param init a list containing the vector of initial values \code{theta_init},
#'  the vector of upper bounds \code{theta_upper} and the vector of the
#'  lower bounds \code{theta_lower} for the parameters to estimate.
#'  For each elements of the list, give the name of the covariate followed by
#'  the vector of the fixed initials values
#'
#'
#' @param method_opt optimization method used in \code{optim} function. The default
#' algorithm is "L-BFGS-B".
#'
#'
#'
#' @param maxit_opt option for maximum of iteration in optimization function
#'
#'
#'
#'
#' @param gradient True if optimization process requires gradient to be provided
#'
#'
#'
#' @param hessian_varcov TRUE if user wants variance covariance matrix using
#' hessian function
#'
#'
#'
#' @param optim_func specify which function to be used for optimization purposes.
#'
#'
#'
#'
#' @param optimizer only use this argument when \code{optim_func="bbmle"}
#'
#'
#'
#'
#' @param optim_fixed to specify with parameter to not estimated in the
#' estimation process
#'
#'
#'
#' @param link_tau must be specified only for \code{model ="tneh"}.
#' Default is linear link \code{("linear")}. Another link is loglinear
#'  \code{("loglinear")}.
#'
#'
#'
#'
#' @param ncoor_des if null, the initial parameters are defaults. If else,
#' the initials parameters are obtained via coordinates descent algorithms
#'
#'
#'
#'
#' @param  trace Non-negative integer corresponding to the \code{trace} argument
#' as in \code{optim}
#'
#'
#'
#'
#'
#' @param iter_eps this parameter only works when \code{ncoor_des = "iter"}; It
#'    allows to run coordinates descent algorithm until the stooping criteria
#'    equal at least to the specified value.
#'
#'
#'
#'
#' @param clustertype related to cluster type in \code{marqLevAlg} package
#'
#'
#'
#' @param nproc number of processors for parallel computing as in \code{marqLevAlg}
#'
#'
#'
#' @param nvalues number of set of initial values when using multiple
#'  initials values
#'
#' @param subset an expression indicating which subset of the data should be used
#'  in the modeling. All observations are included by default
#'
#'
#' @param na.action as in the coxph function, a missing-data filter function.
#'
#'
#' @param sign_delta only used for mixture cure rate models to specify if the
#' effects or minus the effects of covariates acting on uncured survival to be
#' considered. Default will be sign_delta = "1". The alternative is
#' sign_delta = "-1".
#'
#'
#' @param  ... additional parameters such \code{z_alpha},
#'  and \code{z_tau}. For more details, use the help function.
#'
#'
#'
#'
#' @keywords curesurv
#'
#' @return An object of class \code{curesurv}.
#' This object is a list containing the following components:
#'
#'   \item{iter_coords}{number of iterations performed to obtain initial values
#'    of the parameters in tneh model only}
#'   \item{coefficients}{estimates found for the model}
#'   \item{estimates}{estimates in the appropriate scale for the model}
#'   \item{loglik}{corresponds to the log-likelihood computed; if only the
#'   pophaz is provided, the log-likelihood doesn't correspond to the total
#'   log-likelihood. The part of the cumulative population hazard is a constant
#'   and is dropped for the computation as presented in Esteve et \emph{al.} (1990);
#'   The total log-likekihood is calculed if the user specifies a column name equal
#'   expected cumulative mortality (cumpophaz)}
#'   \item{iteractions}{the number iterations attained to estimate the parameters of the related model}
#'   \item{evaluations}{the number of times the log-likelihood function was evaluated until to reach the convergence}
#'   \item{convergence}{an integer code as in optim when \code{L-BFGS-B} method
#' is used in optim.}
#'
#'   \item{message}{a character string returned by the optimizer}
#'   \item{varcov}{the variance covariance matrix of the parameters estimated}
#'   \item{varcov_star}{the variance covariance matrix of the coefficients of the model of interest}
#'   \item{std_err}{the standard errors of the estimated parameters}
#'   \item{std_err_star}{the standard errors of the coefficients of the model of interest}
#'   \item{AIC}{the Akaike information criteria from the model of interest}
#'   \item{n.events}{the number of events in the dataset. Events are considered}
#'   \item{n.obs}{the number of observations in the dataset.}
#'   \item{model}{if fitted model is a mixture model, it returns "mixture".
#'    If fitted model is Time-To-Null Excess Hazard model, it returns "nmixture".}
#'   \item{Terms}{the representation of the terms in the model}
#'   \item{pophaz.alpha}{logical value to indicate if fitted cure model requires
#'    correction of mortality rates by a scale parameter}
#'   \item{pophaz}{corresponds to the the population instantaneous mortality rates.}
#'   \item{cumpophaz}{corresponds to the population cumulative mortality rates.}
#'   \item{frailtyhp}{a booleen to be specified if a frailty correction is needed for the population hazard.}
#'   \item{dist}{For mixture model, it corresponds to the function used to fit
#'   the uncured patients survival. By default, ("weib") is used. Another option
#'   is the exponentiated Weibull function ("eweib"). For non-mixture models,
#'   this argument corresponds to the name of the model. By default, ("tneh")
#'   is used to fit the time to null excess hazard model proposed by Boussari
#'   et al.}
#'   \item{xmax}{maximum follow-up time to evaluate the TTC}
#'   \item{z_tau}{Covariates acting on parameter tau in non mixture cure model tneh}
#'   \item{link_tau}{returned  only for model ="tneh"; returned by default is
#'   "linear" or "loglinear" for linear or loglinear link function of covariates
#'    acting on tau parameter.}
#'   \item{z_alpha}{Covariates acting on parameter alpha in non mixture cure model tneh}
#'   \item{z_c}{Covariates acting on cure fraction in mixture cure model}
#'   \item{z_ucured}{covariates acting on survival of uncured in mixture cure model}
#'   \item{z_pcured}{Covariates acting on cure fraction in mixture cure model}
#'   \item{z_ucured}{covariates acting on survival of uncured in mixture cure model}
#'   \item{data}{the dataset used to run the model}
#'   \item{call}{the function call based on model}
#'   \item{formula}{the formula as a formula object}
#'
#' @author Juste Goungounga, Judith Breaud, Olayide Boussari, Laura Botta, Valerie Jooste
#' @references Boussari O, Bordes L, Romain G, Colonna M, Bossard N, Remontet L,
#'  Jooste V. Modeling excess hazard with time-to-cure as a parameter.
#'   Biometrics. 2021 Dec;77(4):1289-1302. doi: 10.1111/biom.13361.
#'    Epub 2020 Sep 12. PMID: 32869288.
#' (\href{https://pubmed.ncbi.nlm.nih.gov/32869288/}{pubmed})
#'
#'
#'  Boussari O, Romain G, Remontet L, Bossard N, Mounier M, Bouvier AM,
#'  Binquet C, Colonna M, Jooste V. A new approach to estimate time-to-cure from
#'  cancer registries data. Cancer Epidemiol. 2018 Apr;53:72-80.
#'  doi: 10.1016/j.canep.2018.01.013. Epub 2018 Feb 4. PMID: 29414635.
#' (\href{https://pubmed.ncbi.nlm.nih.gov/29414635/}{pubmed})
#'
#'
#'  Phillips N, Coldman A, McBride ML. Estimating cancer prevalence using
#'  mixture models for cancer survival. Stat Med. 2002 May 15;21(9):1257-70.
#'  doi: 10.1002/sim.1101. PMID: 12111877.
#'  (\href{https://pubmed.ncbi.nlm.nih.gov/12111877/}{pubmed})
#'
#'
#'   De Angelis R, Capocaccia R, Hakulinen T, Soderman B, Verdecchia A. Mixture
#'   models for cancer survival analysis: application to population-based data
#'   with covariates. Stat Med. 1999 Feb 28;18(4):441-54.
#'   doi: 10.1002/(sici)1097-0258(19990228)18:4<441::aid-sim23>3.0.co;2-m.
#'   PMID: 10070685.
#'   (\href{https://pubmed.ncbi.nlm.nih.gov/10070685/}{pubmed})
#'
#'
#'   Botta L, Caffo O, Dreassi E, Pizzoli S, Quaglio F, Rugge M, Valsecchi MG.
#'   A new cure model that corrects for increased risk of non-cancer death:
#'   analysis of reliability and robustness, and application to real-life data.
#'   BMC Med Res Methodol. 2023 Mar 25;23(1):70.
#'   doi: 10.1186/s12874-023-01876-x. PMID: N/A.
#'   (\href{https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-023-01876-x}{pubmed})
#'
#'
#' @examples
#'
#' library("curesurv")
#' library("survival")
#'
#'
#'
#' # Net survival setting
#' # Mixture cure model with Weibull function for the uncured patients survival:
#' # no covariate
#'
#' theta_init2 <- rep(0, 3)
#' theta_lower2 <- c(-Inf,-Inf,-Inf)
#' theta_upper2 <- c(Inf, Inf, Inf)
#'
#'
#' fit_m0_ml <- curesurv(Surv(time_obs, event) ~ 1 | 1,
#'              pophaz = "ehazard",
#'              cumpophaz = "cumehazard",
#'              model = "mixture", dist = "weib",
#'              data = testiscancer,
#'              init = list(theta_init = theta_init2,
#'              theta_lower = theta_lower2,
#'              theta_upper = theta_upper2),
#'              method_opt = "L-BFGS-B")
#' fit_m0_ml
#'
#'
#'
#' # Mixture cure model with Weibull function for the uncured patients survival:
#' #standardized age as covariate
#'
#'
#' fit_m2_ml <- curesurv(Surv(time_obs, event) ~ age_cr | age_cr,
#'                    pophaz = "ehazard",
#'                    cumpophaz = "cumehazard",
#'                    model = "mixture", dist = "weib",
#'                    data = testiscancer,
#'                    method_opt = "L-BFGS-B")
#'
#'  fit_m2_ml
#'
#'
#'
#' ## Non mixture cure model
#' ### TNEH Null model
#' #### loglinear effect of covariates on time-to-null excess hazard
#'
#' theta_init2 <- rep(0, 3)
#' theta_lower2 <- c(-Inf,-Inf,-Inf)
#' theta_upper2 <- c(Inf, Inf, Inf)
#'
#' fit_m0_mult_tneh <- curesurv(Surv(time_obs, event) ~ 1,
#'                           pophaz = "ehazard",
#'                           cumpophaz = "cumehazard",
#'                           model = "nmixture",
#'                           dist = "tneh", link_tau = "loglinear",
#'                           data = testiscancer,
#'                           init = list(theta_init = theta_init2,
#'                                       theta_lower = theta_lower2,
#'                                       theta_upper = theta_upper2),
#'                           method_opt = "L-BFGS-B")
#'
#'
#' fit_m0_mult_tneh
#'
#' #### Additive parametrization
#' theta_init2 <- c(1, 6, 6)
#' theta_lower2 <- c(0,1,0)
#' theta_upper2 <- c(Inf, Inf, Inf)
#'
#' fit_m0_ad_tneh <- curesurv(Surv(time_obs, event) ~ 1,
#'                           pophaz = "ehazard",
#'                           cumpophaz = "cumehazard",
#'                           model = "nmixture",
#'                           dist = "tneh", link_tau = "linear",
#'                           data = testiscancer,
#'                           init = list(theta_init = theta_init2,
#'                                       theta_lower = theta_lower2,
#'                                       theta_upper = theta_upper2),
#'                           method_opt = "L-BFGS-B")
#'
#'
#'
#'  fit_m0_ad_tneh
#'
#' #### Additive parametrization, with covariates
#' fit_m1_ad_tneh <- curesurv(Surv(time_obs, event) ~ z_alpha(age_cr) +
#'                           z_tau(age_cr),
#'                           pophaz = "ehazard",
#'                           cumpophaz = "cumehazard",
#'                           model = "nmixture",
#'                           dist = "tneh", link_tau = "linear",
#'                           data = testiscancer,
#'                           method_opt = "L-BFGS-B")
#'
#'
#'
#'  fit_m1_ad_tneh
#'
#'
#'
#'
#' @seealso [curesurv::predict.curesurv()], [curesurv::print.curesurv()], `browseVignettes("curesurv")`
#'
#'
#' @import survival stats Formula
#'
#' @export
curesurv <- function(formula,
                     data,
                     pophaz = NULL,
                     cumpophaz = NULL,
                     pophaz.alpha = FALSE,
                     model = "nmixture",
                     dist = "weib",
                     link_tau = "linear",
                     ncoor_des = NULL,
                     init = NULL,
                     maxit_opt = 10000,
                     gradient = FALSE,
                     hessian_varcov = TRUE,
                     optim_func = "optim",
                     optimizer = "optim",
                     method_opt = "L-BFGS-B",
                     trace = 0,
                     nvalues = 10,
                     iter_eps = 1e-8,
                     optim_fixed = NULL,
                     clustertype = NULL,
                     nproc = 1,
                     subset,
                     na.action,
                     sign_delta,
                     ...) {

  Call <- match.call()
  other.Args <- list(...)
  m <- match.call(expand.dots = FALSE)
  indx <- match(c("formula", "data", "subset", "na.action"),
                names(Call),
                nomatch = 0)

  temp <- Call[c(1, indx)]
  temp[[1]] <- as.name("model.frame")
  f <- Formula::Formula(formula)
  special <- c("strata")
  Terms <- if (missing(data)) {
    terms(f, special)
  } else {
    terms(f, special, data = data)
  }
  temp$formula <- Terms
  rate_data <- ratedata(pophaz = pophaz, cumpophaz = cumpophaz, data = data)
  modelout <- cureformula(f = f,
                          temp = temp,
                          Terms = Terms,
                          model = model,
                          link_tau = link_tau,
                          dist = dist,
                          pophaz = rate_data$pophaz,
                          cumpophaz = rate_data$cumpophaz,
                          method_opt = method_opt,
                          pophaz.alpha = pophaz.alpha,
                          maxit_opt = maxit_opt,
                          optim_func = optim_func,
                          gradient = gradient,
                          hessian_varcov = hessian_varcov,
                          ncoor_des = ncoor_des,
                          trace = trace,
                          optim_fixed = optim_fixed,
                          optimizer = optimizer,
                          init = init,
                          nvalues = nvalues,
                          iter_eps = iter_eps,
                          clustertype = clustertype,
                          nproc = nproc,
                          subset = subset,
                          na.action = na.action,
                          sign_delta = sign_delta)
  modelout$data <- data
  modelout$Call <- Call
  modelout$formula<-formula
  class(modelout) <- c("curesurv","list")
  return(modelout)
}