Package description
The hdcuremodels and survival packages
should be loaded.
Data examples
The package includes two datasets: amltrain (Archer et
al, 2024) and amltest (Archer et al, 2024; Bamopoulos et al
2020). Both datasets include patients diagnosed with acute myeloid
leukemia (AML) who were cytogenetically normal at diagnosis along with
the same variables: cryr is the duration of complete
response (in years), relapse.death is a censoring variable
where 1 indicates the patient relapsed or died and 0 indicates the
patient was alive at last follow-up, and expression for 320 transcripts
measured using RNA-sequencing. The restriction to 320 transcripts was to
reduce run time. Therefore, results obtained with these data will not
precisely recapitulate those in the original publication (Archer et al,
2024). amltrain includes the 306 subjects that were used
for training the penalized MCM while amltest includes the
40 subjects that were used to test the penalized MCM.
We also included a function, generate_cure_data, that
allows the user to generate time-to-event data that includes a cured
fraction. Various parameters in this function will allow the user to
explore the impact of sample size (n), number of penalized
variables in the model (j), number of penalized variables
truly associated with the outcome (n_true), effect size
(referred to as the signal amplitude) of the penalized variables
(a), correlation among penalized variables
(rho), same_signs which if TRUE
indicates that the penalized incidence and latency coefficients should
have the same signs, number of unpenalized variables to coerce into the
model (nonp), itct_mean the incidence
intercept which controls the cure fraction, cens_ub the
upper bound for the Uniform[0, cens_ub] distribution for
generating censoring times, model which specifies the type
of regression model to use for generating event times,
alpha and lambda the numeric shape and rate
parameter values to specify when model = "weibull". When
a = 0 there is no relationship between the j
variables and the outcome. When rho = 0 the j
variables generated are uncorrelated while as rho
approaches 1 the variables generated are highly correlated. Note that
rho cannot be set equal to 1.
withr::local_seed(23)
data <- generate_cure_data(n = 200, j = 10, n_true = 5, a = 1.8, rho = 0.2)
training <- data$training
testing <- data$testing
head(training)
#> Time Censor U1 U2 X1 X2 X3
#> 1 0.3073157 1 -0.38092821 -1.1512299 0.2601014 -0.08664437 0.37085391
#> 2 9.9987949 0 0.27775133 0.3498232 -0.1313600 -0.14703388 0.11419749
#> 3 3.6234583 0 -0.09619218 -0.5486285 0.1070170 -0.55391727 -0.06514085
#> 4 1.7712017 0 1.00829091 0.4991512 0.4844801 0.78733289 -0.73085715
#> 5 1.9837930 0 1.35379114 -0.7851983 0.1338578 0.29613139 -0.26529766
#> 6 3.3668557 0 0.26424693 -0.2152782 -0.8450739 -0.65111161 -0.05986652
#> X4 X5 X6 X7 X8 X9
#> 1 0.80194354 0.10743195 -0.11153197 -0.43545676 -0.16495055 -0.6693095
#> 2 -0.59643655 -0.47107853 -0.69620473 -0.04103287 -0.21026165 0.4764264
#> 3 -0.70499475 0.47765142 0.60504607 1.01672281 0.71624649 0.2765438
#> 4 -0.15857135 0.41828145 0.17988730 -0.77016976 -0.05602750 0.2451666
#> 5 0.07002325 0.59285090 0.05405545 -0.18371201 0.69365239 -0.1808401
#> 6 -0.06927526 0.04030413 -0.01702558 0.41253565 0.04395145 -0.1051484
#> X10
#> 1 0.04322245
#> 2 -0.41536326
#> 3 0.30864241
#> 4 0.14679400
#> 5 -0.88188624
#> 6 -0.10867668The training and testing data.frames have the same structure which
includes Time, Censor, and covariates. Covariates prefixed with “U” are
the unpenalized covariates and is equal to the value passed to
nonp (default is 2). Variables prefixed with “X” are the
penalized covariates and is equal to the value passed to j.
To identify which penalized covariates are truly related to the
incidence and latency portions of the model, the user can inspect the
indices returned by nonzero_b and nonzero_beta
in the parameters component of the list.
data$parameters
#> $nonzero_b
#> [1] 2 3 6 8 10
#>
#> $nonzero_beta
#> [1] 1 3 6 8 10
#>
#> $b_u
#> [1] 0.1669904 -0.4123894
#>
#> $beta_u
#> [1] -0.3969023 0.3423121
#>
#> $b_p_nz
#> [1] 1.8 1.8 -1.8 -1.8 1.8
#>
#> $beta_p_nz
#> [1] 1.8 -1.8 1.8 1.8 1.8
#>
#> $itct
#> [1] 0.4489596In this example, the 2nd, 3rd, 6th, 8th, and 10th penalized covariates, corresponding to
names(training)[grep("^X", names(training))][data$parameters$nonzero_b]
#> [1] "X2" "X3" "X6" "X8" "X10"are truly associated with incidence while the 1st, 3rd, 6th, 8th, and 10th penalized covariates are truly associated with latency.
Assessing model assumptions for fitting a mixture cure model
The workflow for fitting a mixture cure model should include the
assessment of two assumptions: first, that a non-zero cure fraction is
present; second, that there is sufficient follow-up (Maller and Zhou,
1996). Inferential tests for assessing these two assumptions are
included in the hdcuremodels package. The functions
nonzerocure_test and sufficient_fu_test both
take a survfit object as their argument.
As can be seen from the Kaplan-Meier plot, there is a long-plateau
that does not drop down to zero. This may indicate the presence of a
cured fraction. We can test the null hypothesis that the cured fraction
is zero against the alternative hypothesis that the cured fraction is
not zero using nonzerocure_test (Maller and Zhou,
1996).
nonzerocure_test(km_train)
#> Warning in Surv(object$time, object$n.event): Invalid status value, converted
#> to NA
#> $proportion_susceptible
#> [1] 0.7146919
#>
#> $proportion_cured
#> [1] 0.2853081
#>
#> $p_value
#> [1] 0.001
#>
#> $time_95_percent_of_events
#> [1] 5.553847Given the small p-value we reject the null hypothesis and conclude
there is a non-zero cure fraction present. We can also extract the cured
fraction as the Kaplan-Meier estimate beyond the last observed event
(Goldman, 1991) using the cure_estimate function.
This estimate requires sufficiently long follow-up which can be
tested using the sufficient_fu_test function (Maller and
Zhou, 1996).
This function tests the null hypothesis of insufficient follow-up against the alternative that there is sufficient follow-up. Based on these results, we reject the null hypothesis and conclude there is sufficient follow-up. Having verified these two assumptions, we can now fit a mixture cure model.
Penalized mixture cure models
Fitting penalized mixture cure models using GMIFS
The primary function for fitting parametric models using the GMIFS
algorithm in the hdcuremodels package is
curegmifs. The function arguments are
args(curegmifs)
#> function (formula, data, subset, x_latency = NULL, model = c("weibull",
#> "exponential"), penalty_factor_inc = NULL, penalty_factor_lat = NULL,
#> epsilon = 0.001, thresh = 1e-05, scale = TRUE, maxit = 10000,
#> inits = NULL, verbose = TRUE, suppress_warning = FALSE, na.action = na.omit,
#> ...)
#> NULLThe curegmifs function accepts a model formula that
specifies the time-to-event outcome on the left-hand side of the
equation as a Surv object and any incidence predictor
variable(s) on the right-hand side of the equation. Note that at least
some incidence predictor variables must be included in order to fit a
penalized mixture cure model, otherwise, the survival
package functions should be used to fit time-to-event models that lack
an incidence component. The data parameter specifies the
name of the data.frame. Users should note that rows with missing data
are omitted (that is, only na.action = na.omit is
operational) therefore users may want to impute missing data prior to
calling this function. The optional subset parameter can be
used to limit model fitting to a subset of observations in the data. The
x_latency parameter specifies the variables to be included
in the latency portion of the model and can be either a matrix of
predictors, a model formula with the right hand side specifying the
latency variables, or the same data.frame passed to the
data parameter. Note that when using the model formula
syntax for x_latency it cannot handle
x_latency = ~ .. The curegmifs function can
fit either a either "weibull" or "exponential"
model, which is specified using the model parameter. Other
parameters include penalty_factor_inc which is an optional
numeric vector with length equal to the number of incidence variables,
where 1 indicates that variable should be penalized and 0 indicates that
variable is unpenalized (default is that all variables are penalized).
Likewise penalty_factor_lat is an optional numeric vector
with length equal to the number of latency variables, where 1 indicates
that variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Unpenalized
predictors are those that we want to coerce into the model (e.g., age)
so that no penalty is applied. By default the variables are centered and
scaled (scale = TRUE). The parameter epsilon
is the size of the coefficient update at each step (default = 0.001).
The GMIFS algorithm stops when either the difference between successive
log-likelihoods is less than thresh (default 1e-05) or the
algorithm has exceeded the maximum number of iterations
(maxit). Initialization is automatically provided by the
function though inits can be used to provide initial values
for the incidence intercept, unpenalized incidence and latency
coefficients, rate parameter, and shape parameter if fitting a Weibull
mixture cure model. By default verbose = TRUE so that
running information is echoed to the R console.
fitgmifs <- curegmifs(Surv(cryr, relapse.death) ~ .,
data = amltrain,
x_latency = amltrain, model = "weibull"
)Details of the GMIFS mixture cure model have been described in Fu et al, 2022.
Fitting penalized mixture cure models using E-M algorithm
The primary function for fitting penalized MCMs using the E-M
algorithm in the hdcuremodels package is
cureem. The function arguments are
args(cureem)
#> function (formula, data, subset, x_latency = NULL, model = c("cox",
#> "weibull", "exponential"), penalty = c("lasso", "MCP", "SCAD"),
#> penalty_factor_inc = NULL, penalty_factor_lat = NULL, thresh = 0.001,
#> scale = TRUE, maxit = NULL, inits = NULL, lambda_inc = 0.1,
#> lambda_lat = 0.1, gamma_inc = 3, gamma_lat = 3, na.action = na.omit,
#> ...)
#> NULLThe cureem function accepts a model formula that
specifies the time-to-event outcome on the left-hand side of the
equation as a Surv object and any incidence predictor
variable(s) on the right-hand side of the equation. Note that at least
some incidence predictor variables must be included in order to fit a
penalized mixture cure model, otherwise, the survival
package functions should be used to fit time-to-event models that lack
an incidence component. The data parameter specifies the
name of the data.frame and the optional subset parameter
can be used to limit model fitting to a subset of observations in the
data. The x_latency parameter specifies the variables to be
included in the latency portion of the model and can be either a matrix
of predictors, a model formula with the right hand side specifying the
latency variables, or the same data.frame passed to the
data parameter. Note that when using the model formula
syntax for x_latency it cannot handle
x_latency = ~ .. The cureem function can fit
one of three models which is specified using the model
parameter, which can be either "cox" (default),
"weibull", or "exponential". Other parameters
include penalty which can be "lasso",
"MCP", or "SCAD" when fitting a
"cox" model but must be "lasso" when fitting a
"weibull" or "exponential" model.
penalty_factor_inc is an optional numeric vector with
length equal to the number of incidence variables, where 1 indicates
that variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Likewise
penalty_factor_lat is an optional numeric vector with
length equal to the number of latency variables, where 1 indicates that
variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Unpenalized
predictors are those that we want to coerce into the model (e.g., age)
so that no penalty is applied. The iterative process stops when the
differences between successive expected penalized complete-data
log-likelihoods for both incidence and latency components are less than
thresh (default = 0.001). By default the variables are
centered and scaled (scale = TRUE). The user can specify
the maximum number of passes over the data for each lambda using
maxit, which defaults to 100 when
penalty = "lasso" and 1000 when either
penalty = "MCP" or penalty = "SCAD".
Initialization is automatically provided by the function though
inits can be used to provide initial values for the
incidence intercept, unpenalized indicidence and latency coefficients,
rate parameter (for Weibull and exponential MCM), and shape parameter
(for Weibull MCM). By default verbose = TRUE so that
running information is echoed to the R console. The user can also
specify the penalty parameter for the incidence
(lambda_inc) and latency (lambda_lat) portions
of the model and the \(\gamma\) penalty
when MCP or SCAD is used (gamma_inc and
gamma_lat).
Details of the E-M MCM have been described in the Supplementary Material of Archer et al, 2024.
Cross-validation
There is a function for performing cross-validation (CV)
corresponding to each of the two optimization methods. The primary
function for fitting cross-validated penalized MCMs using the E-M
algorithm in the hdcuremodels package is
cv_cureem. The function arguments are
args(cv_cureem)
#> function (formula, data, subset, x_latency = NULL, model = c("cox",
#> "weibull", "exponential"), penalty = c("lasso", "MCP", "SCAD"),
#> penalty_factor_inc = NULL, penalty_factor_lat = NULL, fdr_control = FALSE,
#> fdr = 0.2, grid_tuning = FALSE, thresh = 0.001, scale = TRUE,
#> maxit = NULL, inits = NULL, lambda_inc_list = NULL, lambda_lat_list = NULL,
#> nlambda_inc = NULL, nlambda_lat = NULL, gamma_inc = 3, gamma_lat = 3,
#> lambda_min_ratio_inc = 0.1, lambda_min_ratio_lat = 0.1, n_folds = 5,
#> measure_inc = c("c", "auc"), one_se = FALSE, cure_cutoff = 5,
#> parallel = FALSE, seed = NULL, verbose = TRUE, na.action = na.omit,
#> ...)
#> NULLThe cv_cureem function accepts a model formula that
specifies the time-to-event outcome on the left-hand side of the
equation as a Surv object and any incidence predictor
variable(s) on the right-hand side of the equation. Note that at least
some incidence predictor variables must be included in order to fit a
penalized mixture cure model, otherwise, the survival
package functions should be used to fit time-to-event models that lack
an incidence component. The data parameter specifies the
name of the data.frame and the optional subset parameter
can be used to limit model fitting to a subset of observations in the
data. The x_latency parameter specifies the variables to be
included in the latency portion of the model and can be either a matrix
of predictors, a model formula with the right hand side specifying the
latency variables, or the same data.frame passed to the
data parameter. Note that when using the model formula
syntax for x_latency it cannot handle
x_latency = ~ .. The cv_cureem function can
fit one of three models which is specified using the model
parameter, which can be either "cox" (default),
"weibull", or "exponential". Other parameters
include penalty which can be "lasso",
"MCP", or "SCAD" when fitting a
"cox" model but must be "lasso" when fitting a
"weibull" or "exponential" model.
penalty_factor-inc is an optional numeric vector with
length equal to the number of incidence variables, where 1 indicates
that variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Likewise
penalty_factor_lat is an optional numeric vector with
length equal to the number of latency variables, where 1 indicates that
variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Unpenalized
predictors are those that we want to coerce into the model (e.g., age)
so that no penalty is applied. The user can choose to use the model-X
knock-off procedure to control the false discovery rate (FDR) by
specifying fdr_control = TRUE and optionally changing the
FDR threshold (default fdr = 0.20) (Candes et al, 2018). To
identify the optimal \(\lambda\) for
the incidence and latency portions of the model, the user can set
grid_tuning = TRUE (default is that one value for \(\lambda\) is used in both portions of the
model). Other useful parameters for the cross-validation function
include n_folds, an integer specifying the number of folds
for the k-fold cross-validation procedure (default is 5);
measure_inc which specifies the evaluation criterion used
in selecting the optimal penalty which can be "c" for the
C-statistic using cure status weighting (Asano and Hirakawa, 2017) or
"auc" for cure prediction using mean score imputation
(Asano et al, 2014) (default is measure_inc = "c");
one_se is a logical variable that if TRUE then the one
standard error rule is used which selects the most parsimonious model
having evaluation criterion no more than one standard error worse than
that of the best evaluation criterion (default is FALSE); and
cure_cutoff which is a numeric value representing the
cutoff time used to represent subjects not experiencing the event by
this time are cured which is used to produce a proxy for the unobserved
cure status when calculating the C-statistic and AUC (default is 5). If
the logical parameter parallel is TRUE, cross-validation
will be performed using parallel processing which requires the
foreach and doMC R packages. To foster
reproducibility of cross-validation results, seed can be
set to an integer.
As with cureem, the iterative process stops when the
differences between successive expected penalized complete-data
log-likelihoods for both incidence and latency components are less than
thresh (default = 0.001). By default the variables are
centered and scaled (scale = TRUE). The user can specify
the maximum number of passes over the data for each lambda using
maxit, which defaults to 100 when
penalty = "lasso" and 1000 when either
penalty = "MCP" or penalty = "SCAD".
Initialization is automatically provided by the function though
inits can be used to provide initial values for the
incidence intercept, unpenalized indicidence and latency coefficients,
rate parameter (for Weibull and exponential MCM), and shape parameter
(for Weibull MCM). When model = "cox", inits
should also include a numeric vector for the latency survival
probabilities. Optionally, the user can supply a numeric vector to
search for the optimal penalty for the incidence portion
(lambda_inc_list) and a numeric vector to search for the
optimal penalty for the latency portion (lambda_lat_list)
of the model. By default the number of values to search for the optimal
incidence penalty is 10 which can be changed by specifying an integer
for nlambda_int and similarly for latency by specifying an
integer for nlambda_lat. If penalty is either
"MCP" or "SCAD", the user can optionally
specify the penalization parameter \(\gamma\) for the incidence
(gamma_inc) and latency (gamma_lat) portions
of the model. By default verbose = TRUE so that running
information is echoed to the R console. The user can also specify the
penalty parameter for the incidence (lambda_inc) and
latency (lambda_lat) portions of the model and the \(\gamma\) penalty when MCP or SCAD is used
(gamma_inc and gamma_lat).
fit_cv <- cv_cureem(Surv(Time, Censor) ~ .,
data = training,
x_latency = training, fdr_control = FALSE,
grid_tuning = FALSE, nlambda_inc = 10,
nlambda_lat = 10, n_folds = 2, seed = 23,
verbose = TRUE
)
#> Fold 1 out of 2 training...
#> Fold 2 out of 2 training...
#> Selected lambda for incidence: 0.052
#> Selected lambda for latency: 0.052
#> Maximum C-statistic: 0.717964448478618
#> Fitting a final model...Notice in the previous section describing cureem that
values were supplied for the \(\lambda\) penalty parameters for both the
incidence and latency portions of the model using
lambda_inc and lambda_lat. Those values were
determined from the following repeated 10-fold cross-validation where
the optimal \(\lambda\) for the
incidence portion was identified by fitting the models to maximize the
AUC while the optimal \(\lambda\) for
the latency portion was identified by fitting the models to maximize the
C-statistic. After the CV procedure the mode for each was taken. Because
the run time for the repeated 10-fold CV procedure was 5.65 hours, this
code chunk is not evaluated herein.
lambda_inc <- lambda_lat <- rep(0, 100)
for (k in 1:100) {
print(k)
coxem_auc_k <- cv_cureem(Surv(cryr, relapse.death) ~ .,
data = amltrain, x_latency = amltrain,
model = "cox", penalty = "lasso",
scale = TRUE, grid_tuning = TRUE,
nfolds = 10, nlambda_inc = 20,
nlambda_lat = 20, verbose = FALSE,
parallel = TRUE, measure_inc = "auc"
)
lambda_inc[k] <- coxem_auc_k$selected_lambda_inc
coxem_c_k <- cv_cureem(Surv(cryr, relapse.death) ~ .,
data = amltrain,
x_latency = amltrain, model = "cox",
penalty = "lasso", scale = TRUE,
grid_tuning = TRUE, nfolds = 10,
nlambda_inc = 20, nlambda_lat = 20,
verbose = FALSE, parallel = TRUE,
measure_inc = "c"
)
lambda_lat[k] <- coxem_c_k$selected_lambda_lat
}
table(lambda_inc)
table(lambda_lat)The primary function for fitting cross-validated penalized MCMs using
the GMIFS algorithm in the hdcuremodels package is
cv_curegmifs. The function arguments are
args(cv_curegmifs)
#> function (formula, data, subset, x_latency = NULL, model = c("weibull",
#> "exponential"), penalty_factor_inc = NULL, penalty_factor_lat = NULL,
#> fdr_control = FALSE, fdr = 0.2, epsilon = 0.001, thresh = 1e-05,
#> scale = TRUE, maxit = 10000, inits = NULL, n_folds = 5, measure_inc = c("c",
#> "auc"), one_se = FALSE, cure_cutoff = 5, parallel = FALSE,
#> seed = NULL, verbose = TRUE, na.action = na.omit, ...)
#> NULLThe cv_curegmifs function accepts a model formula that
specifies the time-to-event outcome on the left-hand side of the
equation as a Surv object and any incidence predictor
variable(s) on the right-hand side of the equation. Note that at least
some incidence predictor variables must be included in order to fit a
penalized mixture cure model, otherwise, the survival
package functions should be used to fit time-to-event models that lack
an incidence component. The data parameter specifies the
name of the data.frame and the optional subset parameter
can be used to limit model fitting to a subset of observations in the
data. The x_latency parameter specifies the variables to be
included in the latency portion of the model and can be either a matrix
of predictors, a model formula with the right hand side specifying the
latency variables, or the same data.frame passed to the
data parameter. Note that when using the model formula
syntax for x_latency it cannot handle
x_latency = ~ .. The cv_curegmifs function can
fit either a either "weibull" or "exponential"
model, which is specified using the model parameter. Other
parameters include penalty_factor_inc which is an optional
numeric vector with length equal to the number of incidence variables,
where 1 indicates that variable should be penalized and 0 indicates that
variable is unpenalized (default is that all variables are penalized).
Likewise penalty_factor_lat is an optional numeric vector
with length equal to the number of latency variables, where 1 indicates
that variable should be penalized and 0 indicates that variable is
unpenalized (default is that all variables are penalized). Unpenalized
predictors are those that we want to coerce into the model (e.g., age)
so that no penalty is applied. The user can choose to use the model-X
knock-off procedure to control the false discovery rate (FDR) by
specifying fdr_control = TRUE and optionally changing the
FDR threshold (default fdr = 0.20) (Candes et al, 2018). By
default the variables are centered and scaled
(scale = TRUE). The parameter epsilon is the
size of the coefficient update at each step (default = 0.001). The GMIFS
algorithm stops when either the difference between successive
log-likelihoods is less than thresh (default 1e-05) or the
algorithm has exceeded the maximum number of iterations
(maxit). Initialization is automatically provided by the
function though inits can be used to provide initial values
for the incidence intercept, unpenalized indicidence and latency
coefficients, rate parameter, and shape parameter if fitting a Weibull
mixture cure model. Other useful parameters for the cross-validation
function include n_folds, an integer specifying the number
of folds for the k-fold cross-validation procedure (default is 5);
measure_inc which specifies the evaluation criterion used
in selecting the optimal penalty which can be "c" for the
C-statistic using cure status weighting (Asano and Hirakawa, 2017) or
"auc" for cure prediction using mean score imputation
(Asano et al, 2014) (default is measure_inc = "c");
one_se is a logical variable that if TRUE then the one
standard error rule is used which selects the most parsimonious model
having evaluation criterion no more than one standard error worse than
that of the best evaluation criterion (default is FALSE); and
cure_cutoff which is a numeric value representing the
cutoff time used to represent subjects not experiencing the event by
this time are cured which is used to produce a proxy for the unobserved
cure status when calculating the C-statistic and AUC (default is 5). If
the logical parameter parallel is TRUE, cross-validation
will be performed using parallel processing which requires the
foreach and doMC R packages. To foster
reproducibility of cross-validation results, seed can be
set to an integer. By default verbose = TRUE so that
running information is echoed to the R console.
Other Package Functions
The four modeling functions cureem,
curegmifs, cv_cureem, and
cv_curegmifs all result in an object of class
mixturecure. Generic functions for resulting
mixturecure objects are available for extracting meaningful
results. The print function prints the first several
incidence and latency coefficients and the rate (exponential and
Weibull) and alpha (Weibull) when fitting a parametric MCM and returns
the fitted object invisible to the user.
The summary function prints the following output for a
model fit using either cureem or
curegmifs:
- the number of non-zero coefficients in the incidence portion of the fitted mixture cure model when using the minimum AIC
- the number of non-zero coefficients in the latency portions of the fitted mixture cure model when using the minimum AIC
- the step and value that maximizes the log-likelihood;
- the step and value that minimizes the AIC;
- the step and value that minimizes the modified AIC (mAIC);
- the step and value that minimizes the corrected AIC (cAIC);
- the step and value that minimizes the BIC;
- the step and value that minimizes the modified BIC (mBIC);
- the step and value that minimizes the extended BIC (EBIC).
summary(fitem)
#> Mixture cure model fit using the EM algorithm
#> Number of non-zero incidence covariates at minimum AIC: 112
#> Number of non-zero latency covariates at minimum AIC: 88
#> Optimal step for selected information criterion: EM algorithm
#> at step = 25 logLik = -1113.55183538453
#> at step = 12 AIC = 2634.47640955092
#> at step = 12 mAIC = 5510.63178303021
#> at step = 12 cAIC = 3415.28410185861
#> at step = 12 BIC = 3382.91701504335
#> at step = 12 mBIC = 5423.13688876735
#> at step = 12 EBIC = 3777.8145382339The summary function prints the following output for a
model fit using either cv_cureem or
cv_curegmifs when fdr_control = FALSE:
- the number of non-zero coefficients in the incidence portion of the fitted mixture cure model
- the number of non-zero coefficients in the latency portions of the fitted mixture cure model
summary(fit_cv)
#> Mixture cure model fit using the EM algorithm
#> using cross-validation
#> Number of non-zero incidence covariates: 5
#> Number of non-zero latency covariates: 6The summary function prints the following output for a
model fit using either cv_cureem or
cv_curegmifs when fdr_control = TRUE:
- Number of non-zero incidence covariates
- Number of non-zero latency covariates
For a cureem or curegmifs fitted
mixturecure object, the plot function provides
a trace of the coefficients’ paths by default. Because it is common for
the user to include at least some of the same variables in both the
incidence and latency portions of the MCM, when
label = TRUE a legend is provided with variables included
in the incidence portion prefixed with “I” and variables included in the
latency portion prefixed with “L.” Alternatively, the type
parameter can be used to specify any of the information criterion
(“logLik”, “AIC”, “cAIC”, “mAIC”, “BIC”, “mBIC”, “EBIC”). For a
cv_cureem or cv_curegmifs fitted
mixturecure object, a lollipop plot of the estimated
incidence and latency coefficients is produced. The height of each
lollipop represents the estimated coefficient for the optimal
cross-validated model.
Coefficient estimates can be extracted from the fitted model using
the coef for any of these model criteria (“logLik”, “AIC”,
“cAIC”, “mAIC”, “BIC”, “mBIC”, “EBIC”) or by specifying the step at
which the model is desired by specifying the model.select
parameter. For example,
is equivalent to
as demonstrated by comparing the results in each object:
names(coef_cAIC)
#> [1] "b0" "beta_inc" "beta_lat"
all.equal(coef_cAIC$rate, coef_12$rate)
#> [1] TRUE
all.equal(coef_cAIC$alpha, coef_12$alpha)
#> [1] TRUE
all.equal(coef_cAIC$b0, coef_12$b0)
#> [1] TRUE
all.equal(coef_cAIC$beta_inc, coef_12$beta_inc)
#> [1] TRUE
all.equal(coef_cAIC$beta_lat, coef_12$beta_lat)
#> [1] TRUEAgain, there are two sets of coefficients: those in the incidence
portion of the model (beta_inc) and those in the latency
portion of the model (beta_lat). Additionally,
b0 is the intercept in the incidence portion of the model.
Depending on the model fit, coef will return
rate (exponential and Weibull) and alpha
(Weibull).
Predictions can be extracted at a given step or information criterion
(“logLik”, “AIC”, “cAIC”, “mAIC”, “BIC”, “mBIC”, “EBIC”) using the
predict function with model_select
specified.
This returns three objects: p_uncured is the estimated
probability of being susceptible (\(\hat{p}(\mathbf{x})\)),
linear_latency is \(\hat{\boldsymbol{\beta}}\mathbf{w}\), while
latency_risk applies high risk and low risk labels using
zero as the cutpoint from the linear_latency vector.
Perhaps we want to apply the 0.5 threshold to p_uncured to
create Cured and Susceptible labels.
Then we can assess how well our MCM identified patients likely to be cured from those likely to be susceptible visually by examining the Kaplan-Meier curves.
We can assess how well our MCM identified higher versus lower risk patients among those predicted to be susceptible visually by examining the Kaplan-Meier curves.
km_suscept <- survfit(Surv(cryr, relapse.death) ~ train_predict$latency_risk, data = amltrain,
subset = (p_group == "Susceptible"))
Of course, we expect our model to perform well on our training data.
We can also assess how well our fitted MCM performs using the
independent test set amltest. In this case we use the
predict function with newdata specified.
Again we will apply the 0.5 threshold to p_uncured to
create Cured and Susceptible labels.
Then we can assess how well our MCM identified patients likely to be cured from those likely to be susceptible visually by examining the Kaplan-Meier curves.
km_suscept_test <- survfit(Surv(cryr, relapse.death) ~ test_predict$latency_risk,
data = amltest, subset = (test_p_group == "Susceptible"))
The hdcuremodels package also includes two functions for
assessing the performance of MCMs. The ability of the MCM to
discriminate between those cured (\(Y_i=0\)) versus those susceptible (\(Y_i=1\)) can be assessed by calculating the
mean score imputation area under the curve using the
auc_mcm function (Asano et al, 2014). In a MCM, when \(\delta_i=1\) we know that the subject
experienced the event. However, when \(\delta_i=0\) either the subject was cured
or the subject would have experienced the event if followed longer than
their censoring time. Therefore, for a cure_cutoff \(\tau\) (default is 5) the outcome \(Y_i\) is defined as \[
Y_i = \begin{cases}
0 \text{ if } T_i >\tau\\
1 \text{ if } T_i \le\tau \text{ and } \delta_i=1\\
\text{missing} \text{ if } T_i \le\tau \text{ and } \delta_i=0\\
\end{cases}.
\] The mean score imputation AUC lets \(Y_i = 1-\hat{p}(\mathbf{x}_i)\) for those
subjects with a missing outcome. The C-statistic for MCMs was adapted to
weight patients by their outcome (cured, susceptible, censored) and is
available in the concordance_mcm function (Asano &
Hirakawa, 2017). In both functions, if newdata is not
specified, the training data are used.
auc_mcm(fitem, model_select = "cAIC")
#> [1] 0.9690409
auc_mcm(fitem, newdata = amltest, model_select = "cAIC")
#> [1] 0.8049214
concordance_mcm(fitem, model_select = "cAIC")
#> [1] 0.8546535
concordance_mcm(fitem, newdata = amltest, model_select = "cAIC")
#> [1] 0.6987875Other generic functions return the number of observations for the
training data nobs, the log-likehood at the selected step
logLik, the family of the mixture cure model
family, and the model formula and data used for training
the incidence and latency components formula.