cox_cure {intsurv}

Cox Cure Rate Model

Description

For right-censored data, fit a regular Cox cure rate model (Kuk and Chen, 1992; Sy and Taylor, 2000) via an EM algorithm. For right-censored data with uncertain event status, fit the Cox cure model proposed by Wang et al. (2020).

Usage

cox_cure(
  surv_formula,
  cure_formula,
  time,
  event,
  data,
  subset,
  contrasts = NULL,
  bootstrap = 0,
  firth = FALSE,
  surv_start = NULL,
  cure_start = NULL,
  surv_offset = NULL,
  cure_offset = NULL,
  em_max_iter = 200,
  em_rel_tol = 1e-05,
  surv_max_iter = 30,
  surv_rel_tol = 1e-05,
  cure_max_iter = 30,
  cure_rel_tol = 1e-05,
  tail_completion = c("zero", "exp", "zero-tau"),
  tail_tau = NULL,
  pmin = 1e-05,
  early_stop = TRUE,
  verbose = FALSE,
  ...
)

cox_cure.fit(
  surv_x,
  cure_x,
  time,
  event,
  cure_intercept = TRUE,
  bootstrap = 0,
  firth = FALSE,
  surv_start = NULL,
  cure_start = NULL,
  surv_offset = NULL,
  cure_offset = NULL,
  surv_standardize = TRUE,
  cure_standardize = TRUE,
  em_max_iter = 200,
  em_rel_tol = 1e-05,
  surv_max_iter = 30,
  surv_rel_tol = 1e-05,
  cure_max_iter = 30,
  cure_rel_tol = 1e-05,
  tail_completion = c("zero", "exp", "zero-tau"),
  tail_tau = NULL,
  pmin = 1e-05,
  early_stop = TRUE,
  verbose = FALSE,
  ...
)

Arguments

surv_formula

A formula object starting with ~ for the model formula in survival model part. For Cox model, no intercept term is included even if an intercept is specified or implied in the model formula. A model formula with an intercept term only is not allowed.

cure_formula

A formula object starting with ~ for the model formula in incidence model part. For logistic model, an intercept term is included by default and can be excluded by adding + 0 or - 1 to the model formula.

time

A numeric vector for the observed survival times.

event

A numeric vector for the event indicators. NA's are allowed and represent uncertain event status.

data

An optional data frame, list, or environment that contains the covariates and response variables (time and event), included in the model. If they are not found in data, the variables are taken from the environment of the specified formula, usually the environment from which this function is called.

subset

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

contrasts

An optional list, whose entries are values (numeric matrices or character strings naming functions) to be used as replacement values for the contrasts replacement function and whose names are the names of columns of data containing factors. See contrasts.arg of model.matrix.default for details.

bootstrap

An integer representing the number of bootstrap samples for estimating standard errors of the coefficient estimates. The bootstrap procedure will not run with bootstrap = 0 by default. If bootstrap > 0, the specified number of bootstrap samples will be used in estimating standard errors.

firth

A logical value indicating whether to use Firth's bias-reduction method (Firth, 1993) in the logistic model component. The default value is FALSE for fitting a regular logistic model. Notice that this argument is experimental and only available for regular Cox cure rate model currently.

surv_start, cure_start

An optional numeric vector representing the starting values for the survival model component or the incidence model component. If surv_start = NULL is specified, the starting values will be obtained from fitting a regular Cox to events only. Similarly, if cure_start = NULL is specified, the starting values will be obtained from fitting a regular logistic model to the non-missing event indicators.

surv_offset, cure_offset

An optional numeric vector representing the offset term in the survival model compoent or the incidence model component. The function will internally try to find values of the specified variable in the data first. Alternatively, one or more offset terms can be specified in the formula (by stats::offset()). If more than one offset terms are specified, their sum will be used.

em_max_iter

A positive integer specifying the maximum iteration number of the EM algorithm. The default value is 200.

em_rel_tol

A positive number specifying the tolerance that determines the convergence of the EM algorithm in terms of the convergence of the covariate coefficient estimates. The tolerance is compared with the relative change between estimates from two consecutive iterations, which is measured by ratio of the L1-norm of their difference to the sum of their L1-norm. The default value is 1e-5.

surv_max_iter, cure_max_iter

A positive integer specifying the maximum iteration number of the M-step routine related to the survival model component or the incidence model component. The default value is 200.

surv_rel_tol, cure_rel_tol

A positive number specifying the tolerance that determines the convergence of the M-step related to the survival model component or the incidence model component in terms of the convergence of the covariate coefficient estimates. The tolerance is compared with the relative change between estimates from two consecutive iterations, which is measured by ratio of the L1-norm of their difference to the sum of their L1-norm. The default value is 1e-5.

tail_completion

A character string specifying the tail completion method for conditional survival function. The available methods are "zero" for zero-tail completion after the largest event times (Sy and Taylor, 2000), "exp" for exponential-tail completion (Peng, 2003), and "zero-tau" for zero-tail completion after a specified tail_tau. The default method is the zero-tail completion proposed by Sy and Taylor (2000).

tail_tau

A numeric number specifying the time of zero-tail completion. It will be used only if tail_completion = "zero-tau". A reasonable choice must be a time point between the largest event time and the largest survival time.

pmin

A numeric number specifying the minimum value of probabilities for sake of numerical stability. The default value is 1e-5.

early_stop

A logical value specifying whether to stop the iteration once the negative log-likelihood unexpectedly increases, which may suggest convergence on likelihood, or indicate numerical issues or implementation bugs. The default value is TRUE.

verbose

A logical value. If TRUE, a verbose information will be given along iterations for tracing the convergence. The default value is FALSE.

...

Other arguments for future usage. A warning will be thrown if any invalid argument is specified.

surv_x

A numeric matrix for the design matrix of the survival model component.

cure_x

A numeric matrix for the design matrix of the cure rate model component. The design matrix should exclude an intercept term unless we want to fit a model only including the intercept term. In that case, we need further set cure_intercept = FALSE to not standardize the intercept term.

cure_intercept

A logical value specifying whether to add an intercept term to the cure rate model component. If TRUE by default, an intercept term is included.

surv_standardize

A logical value specifying whether to standardize the covariates for the survival model part. If FALSE, the covariates will be standardized internally to have mean zero and standard deviation one.

cure_standardize

A logical value specifying whether to standardize the covariates for the cure rate model part. If TRUE (by default), the covariates will be standardized internally to have mean zero and standard deviation one.

Value

cox_cure object for regular Cox cure rate model or cox_cure_uncer object for Cox cure rate model with uncertain events.

References

Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38.

Kuk, A. Y. C., & Chen, C. (1992). A mixture model combining logistic regression with proportional hazards regression. Biometrika, 79(3), 531–541.

Peng, Y. (2003). Estimating baseline distribution in proportional hazards cure models. Computational Statistics & Data Analysis, 42(1-2), 187–201.

Sy, J. P., & Taylor, J. M. (2000). Estimation in a Cox proportional hazards cure model. Biometrics, 56(1), 227–236.

Wang, W., Luo, C., Aseltine, R. H., Wang, F., Yan, J., & Chen, K. (2020). Suicide Risk Modeling with Uncertain Diagnostic Records. arXiv preprint arXiv:2009.02597.

See Also

cox_cure_net for regularized Cox cure rate model with elastic-net penalty.

Examples

library(intsurv)

### regular Cox cure rate model ======================================
## 1. simulate right-censored data with a cure fraction
set.seed(123)
n_obs <- 2e2
p <- 5
x_mat <- matrix(rnorm(n_obs * p), nrow = n_obs, ncol = p)
colnames(x_mat) <- paste0("x", seq_len(p))
cure_beta <- rep(0.5, p)
b0 <- - 1
expit <- binomial()$linkinv
ncure_prob <- expit(as.numeric(b0 + x_mat %*% cure_beta))
is_cure <- 1 - rbinom(n_obs, size = 1, prob = ncure_prob)
surv_beta <- rep(0.5, p)
risk_score <- as.numeric(x_mat %*% surv_beta)
event_time <- rexp(n_obs, exp(as.numeric(x_mat %*% surv_beta)))
censor_time <- 10
event <- ifelse(event_time < censor_time & ! is_cure, 1, 0)
obs_time <- ifelse(event > 0, event_time, censor_time)

## model-fitting from given design matrices
fit1 <- cox_cure.fit(x_mat, x_mat, obs_time, event, bootstrap = 30)
summary(fit1)

## coefficient estimates from both model parts
coef(fit1)

## or a particular part
coef(fit1, "surv")
coef(fit1, "cure")

## weighted concordance index (C-index)
fit1$model$c_index
## which also can be computed as follows
cIndex(time = obs_time, event = event,
       risk_score = fit1$fitted$surv_xBeta,
       weight = ifelse(event > 0, 1, fit1$fitted$susceptible_prob))

## 2. create a toy dataset
toy_dat <- data.frame(time = obs_time, status = event)
toy_dat$group <- cut(abs(x_mat[, 1L]), breaks = c(0, 0.5, 1, 3, Inf),
                     labels = LETTERS[1:4])
toy_dat <- cbind(toy_dat, as.data.frame(x_mat[, - 1L, drop = FALSE]))

## model-fitting from given model formula
fit2 <- cox_cure(~ x3 + x4 + group, ~ group + x3 + offset(x2),
                 time = time, event = status, surv_offset = x2,
                 data = toy_dat, subset = group != "D", bootstrap = 30)
summary(fit2)

## get BIC's
BIC(fit1)
BIC(fit2)
BIC(fit1, fit2)


### Cox cure rate model with uncertain event status ==================
## simulate sample data
sim_dat <- simData4cure(nSubject = 200, max_censor = 10, lambda_censor = 0.1,
                        survMat = x_mat, cureMat = x_mat, b0 = 1)
table(sim_dat$case)
table(sim_dat$obs_event, useNA = "ifany")

## use formula
fit3 <- cox_cure(~ x1 + x2 + x3, ~ z1 + z2 + z3,
                 time = obs_time, event = obs_event, data = sim_dat)
summary(fit3)

## use design matrix
fit4 <- cox_cure.fit(x_mat, x_mat, time = sim_dat$obs_time,
                     event = sim_dat$obs_event)
summary(fit4)

## get BIC's
BIC(fit3, fit4)