Large diffs are not rendered by default.

@@ -12,17 +12,21 @@
#' @param event Status indicator, normally 0 = alive, 1 = dead.
#' Must be of the same length with the number of rows as \code{x}.
#' @param ddist Data frame version of x, made by \code{\link[rms]{datadist}}.
#' @param lambda Value of the penalty parameter lambda in
#' \code{\link[glmnet]{glmnet}} or \code{\link[ncvreg]{ncvsurv}}.
#' Required except when \code{model.type == "flasso"}.
#' We will use the selected variables at the provided \code{lambda} to
#' build the nomogram, and make predictions.
#' See the example for choosing a proper lambda value extracted
#' from cross-validation results.
#' @param pred.at Time point at which to plot nomogram prediction axis.
#' @param fun.at Function values to label on axis.
#' @param funlabel Label for \code{fun} axis.
#'
#' @note We will try to use the value of the selected penalty
#' parameter (e.g. lambda, alpha) in the model object fitted by
#' \code{\link[glmnet]{glmnet}}, \code{\link[ncvreg]{ncvsurv}}, or
#' \code{\link[penalized]{penalized}}. The selected variables under
#' the penalty parameter will be used to build the nomogram
#' and make predictions. This means, if you use routines other
#' than \code{hdcox.*} functions to fit the model, please
#' make sure that there is only one set of necessary parameters
#' (e.g. only a single lambda should be in glmnet objects intead
#' of a vector of multiple lambdas) in the fitted model object.
#'
#' @export hdnom.nomogram
#'
#' @importFrom rms ols nomogram
@@ -39,87 +43,27 @@
#' y = Surv(time, event)
#'
#' # Fit penalized Cox model with lasso penalty
#' lassofit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#' fit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' # Prepare data needed for plotting nomogram
#' x.df = as.data.frame(x)
#' dd = datadist(x.df)
#' options(datadist = "dd")
#'
#' # Generate hdnom.nomogram objects and plot nomogram
#' nom = hdnom.nomogram(lassofit$lasso_model, model.type = "lasso",
#' x, time, event, x.df,
#' lambda = lassofit$lasso_best_lambda, pred.at = 365 * 2,
#' nom = hdnom.nomogram(fit$lasso_model, model.type = "lasso",
#' x, time, event, x.df, pred.at = 365 * 2,
#' funlabel = "2-Year Overall Survival Probability")
#'
#' print(nom)
#' plot(nom)
# ### Testing fused lasso models ###
# library("penalized")
# library("survival")
# library("rms")
#
# # Load imputed SMART data
# data(smart)
# x = as.matrix(smart[, -c(1, 2)])
# x = x[1:500, ]
# time = smart$TEVENT[1:500]
# event = smart$EVENT[1:500]
# x.df = as.data.frame(x)
# dd = datadist(x.df)
# options(datadist = "dd")
#
# # Fit penalized Cox model (fused lasso penalty) with penalized package
# set.seed(1010)
# cvfit = optL1(response = Surv(time, event), penalized = x, fusedl = TRUE,
# standardize = TRUE, model = 'cox', fold = 3, trace = TRUE)
#
# # Generate hdnom.nomogram objects and plot nomogram
# nom = hdnom.nomogram(cvfit$fullfit, model.type = "flasso",
# x, time, event, x.df, pred.at = 365 * 2,
# funlabel = "2-Year Overall Survival Probability")
#
# print(nom)
# plot(nom)
#
# ### Testing SCAD models ###
# library("ncvreg")
#
# cvfit = cv.ncvsurv(x, Surv(time, event),
# model = 'cox', penalty = 'SCAD',
# alpha = 1, nfolds = 5,
# seed = 1010, trace = TRUE)
# fit = ncvsurv(x, Surv(time, event), model = 'cox', penalty = 'SCAD', alpha = 1)
#
# nom = hdnom.nomogram(fit, model.type = "scad", x, time, event, x.df,
# lambda = cvfit$lambda.min, pred.at = 365 * 2,
# funlabel = "2-Year Overall Survival Probability")
#
# print(nom)
# plot(nom)
#
# ### Testing Mnet models ###
# cvfit = cv.ncvsurv(x, Surv(time, event),
# model = 'cox', penalty = 'MCP',
# alpha = 0.5, nfolds = 5,
# seed = 1010, trace = TRUE)
# fit = ncvsurv(x, Surv(time, event), model = 'cox', penalty = 'MCP', alpha = 0.5)
#
# # Generate hdnom.nomogram objects and plot nomogram
# nom = hdnom.nomogram(fit, model.type = "mnet", x, time, event, x.df,
# lambda = cvfit$lambda.min, pred.at = 365 * 2,
# funlabel = "2-Year Overall Survival Probability")
#
# print(nom)
# plot(nom)
hdnom.nomogram = function(object,
model.type = c('lasso', 'alasso', 'flasso',
'enet', 'aenet',
'mcp', 'mnet',
'scad', 'snet'),
x, time, event, ddist,
lambda = NULL, pred.at = NULL,
fun.at = NULL, funlabel = NULL) {
pred.at = NULL, fun.at = NULL, funlabel = NULL) {

model.type = match.arg(model.type)

@@ -133,13 +77,16 @@ hdnom.nomogram = function(object,
if (!all(c('coxnet', 'glmnet') %in% class(object)))
stop('object class must be "glmnet" and "coxnet"')

if (is.null(lambda)) stop('Missing argument lambda')
if (length(object$'lambda') != 1L)
stop('There should be one and only one lambda in the model object')

glmnet_pred_lp = as.vector(predict(object, newx = x, s = lambda, type = 'link'))
glmnet_pred_lp = as.vector(
predict(object, newx = x, s = object$'lambda', type = 'link'))

all_vars = rownames(object$beta)
all_vars = rownames(object$'beta')
selected_vars =
all_vars[which(as.vector(abs(coef(object, s = lambda)) > .Machine$double.eps))]
all_vars[which(as.vector(
abs(coef(object, s = object$'lambda')) > .Machine$double.eps))]
ols_formula = paste('glmnet_pred_lp ~',
paste(paste('`', selected_vars, '`', sep = ''),
collapse = ' + '))
@@ -154,7 +101,7 @@ hdnom.nomogram = function(object,
survtime_at_idx = names(survtime_at)

survcurve = glmnet.survcurve(object = object, time = time, event = event,
x = x, survtime = survtime_ones, lambda = lambda)
x = x, survtime = survtime_ones)

}

@@ -163,11 +110,9 @@ hdnom.nomogram = function(object,
if (!all(c('ncvsurv', 'ncvreg') %in% class(object)))
stop('object class must be "ncvreg" and "ncvsurv"')

if (is.null(lambda)) stop('Missing argument lambda')

ncvreg_pred_lp = predict(object, X = x, type = 'link')

all_vars = rownames(object$beta)
all_vars = rownames(object$'beta')
selected_vars =
all_vars[which(as.vector(abs(coef(object)) > .Machine$double.eps))]
ols_formula = paste('ncvreg_pred_lp ~',
@@ -184,7 +129,7 @@ hdnom.nomogram = function(object,
survtime_at_idx = names(survtime_at)

survcurve = ncvreg.survcurve(object = object, time = time, event = event,
x = x, survtime = survtime_ones, lambda = lambda)
x = x, survtime = survtime_ones)

}

@@ -193,7 +138,6 @@ hdnom.nomogram = function(object,
if (!('penfit' %in% class(object)))
stop('object class must be "penfit"')

lambda = object@lambda1
penalized_pred_lp = as.numeric(object@lin.pred)

all_vars = colnames(x)
@@ -228,9 +172,12 @@ hdnom.nomogram = function(object,
if (is.null(funlabel))
funlabel = paste('Overall Survival Probability at Time', pred.at)

nom = list('ols_fit' = ols_fit, 'lambda' = lambda, 'survcurve' = survcurve,
'bhfun' = bhfun, 'pred.at' = pred.at,
'fun.at' = fun.at, 'funlabel' = funlabel)
nom = list('ols_fit' = ols_fit,
'survcurve' = survcurve,
'bhfun' = bhfun,
'pred.at' = pred.at,
'fun.at' = fun.at,
'funlabel' = funlabel)

class(nom) = 'hdnom.nomogram'

@@ -246,10 +193,10 @@ hdnom.nomogram = function(object,
#' linear predictors for all samples
#'
#' @keywords internal
glmnet.survcurve = function(object, time, event, x, survtime, lambda) {
glmnet.survcurve = function(object, time, event, x, survtime) {

lp = as.numeric(predict(object, newx = data.matrix(x),
s = lambda, type = 'link'))
s = object$'lambda', type = 'link'))
basesurv = glmnet.basesurv(time, event, lp, sort(survtime))
p = exp(exp(lp) %*% (-t(basesurv$cumulative_base_hazard)))
colnames(p) = names(sort(survtime))
@@ -302,7 +249,7 @@ glmnet.basesurv = function(time, event, lp,
#' linear predictors for all samples
#'
#' @keywords internal
ncvreg.survcurve = function(object, time, event, x, survtime, lambda) {
ncvreg.survcurve = function(object, time, event, x, survtime) {

lp = as.numeric(predict(object, X = data.matrix(x), type = 'link'))
basesurv = ncvreg.basesurv(time, event, lp, sort(survtime))
@@ -422,8 +369,8 @@ print.hdnom.nomogram = function(x, ...) {
if (class(x) != 'hdnom.nomogram')
stop('object class must be "hdnom.nomogram"')

nom = nomogram(fit = x$ols_fit, fun = x$bhfun,
fun.at = x$fun.at, funlabel = x$funlabel,
nom = nomogram(fit = x$'ols_fit', fun = x$'bhfun',
fun.at = x$'fun.at', funlabel = x$'funlabel',
lp = TRUE, vnames = 'labels', ...)

print(nom)
@@ -450,8 +397,8 @@ plot.hdnom.nomogram = function(x, ...) {
if (class(x) != 'hdnom.nomogram')
stop('object class must be "hdnom.nomogram"')

nom = nomogram(fit = x$ols_fit, fun = x$bhfun,
fun.at = x$fun.at, funlabel = x$funlabel,
nom = nomogram(fit = x$'ols_fit', fun = x$'bhfun',
fun.at = x$'fun.at', funlabel = x$'funlabel',
lp = TRUE, vnames = 'labels', ...)

plot(nom)
@@ -24,6 +24,8 @@
#' From the built \emph{adaptive lasso} or \emph{adaptive elastic-net} model.
#' @param gamma Value of the model parameter gamma for
#' MCP/SCAD/Mnet/Snet models.
#' @param lambda1 Value of the penalty parameter lambda1 for fused lasso model.
#' @param lambda2 Value of the penalty parameter lambda2 for fused lasso model.
#' @param method Validation method.
#' Could be \code{"bootstrap"}, \code{"cv"}, or \code{"repeated.cv"}.
#' @param boot.times Number of repetitions for bootstrap.
@@ -69,27 +71,27 @@
#' y = Surv(time, event)
#'
#' # Fit penalized Cox model with lasso penalty
#' lassofit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#' fit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' # Model validation by bootstrap with time-dependent AUC
#' # Normally boot.times should be set to 200 or more,
#' # we set it to 3 here only to save example running time.
#' val.boot = hdnom.validate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "bootstrap", boot.times = 3,
#' tauc.type = "UNO", tauc.time = seq(0.25, 2, 0.25) * 365,
#' seed = 1010)
#'
#' # Model validation by 5-fold cross-validation with time-dependent AUC
#' val.cv = hdnom.validate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "cv", nfolds = 5,
#' tauc.type = "UNO", tauc.time = seq(0.25, 2, 0.25) * 365,
#' seed = 1010)
#'
#' # Model validation by repeated cross-validation with time-dependent AUC
#' val.repcv = hdnom.validate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "repeated.cv", nfolds = 5, rep.times = 3,
#' tauc.type = "UNO", tauc.time = seq(0.25, 2, 0.25) * 365,
#' seed = 1010)
@@ -122,7 +124,7 @@
#
# set.seed(1010)
# val.boot = hdnom.validate(x, time, event, model.type = "flasso",
# lambda = 60,
# lambda1 = 5, lambda2 = 2,
# method = "bootstrap", boot.times = 10,
# tauc.type = "UNO", tauc.time = seq(0.25, 2, 0.25) * 365,
# seed = 1010)
@@ -156,6 +158,7 @@ hdnom.validate = function(x, time, event,
'mcp', 'mnet',
'scad', 'snet'),
alpha, lambda, pen.factor = NULL, gamma,
lambda1, lambda2,
method = c('bootstrap', 'cv', 'repeated.cv'),
boot.times = NULL, nfolds = NULL, rep.times = NULL,
tauc.type = c("CD", "SZ", "UNO"), tauc.time,
@@ -215,7 +218,7 @@ hdnom.validate = function(x, time, event,
tauc[[i]] =
penalized.validate.internal(
x_tr = x_tr, x_te = x_te, y_tr = y_tr, y_te = y_te,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
tauc.type = tauc.type, tauc.time = tauc.time
)
}
@@ -271,7 +274,7 @@ hdnom.validate = function(x, time, event,
tauc[[i]] =
penalized.validate.internal(
x_tr = x_tr, x_te = x_te, y_tr = y_tr, y_te = y_te,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
tauc.type = tauc.type, tauc.time = tauc.time
)
}
@@ -332,7 +335,7 @@ hdnom.validate = function(x, time, event,
tauc[[j]][[i]] =
penalized.validate.internal(
x_tr = x_tr, x_te = x_te, y_tr = y_tr, y_te = y_te,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
tauc.type = tauc.type, tauc.time = tauc.time
)
}
@@ -377,7 +380,8 @@ hdnom.validate = function(x, time, event,
class(tauc) = c('hdnom.validate',
'penalized.validate.bootstrap')
attr(tauc, 'model.type') = model.type
attr(tauc, 'lambda') = lambda
attr(tauc, 'lambda1') = lambda1
attr(tauc, 'lambda2') = lambda2
attr(tauc, 'boot.times') = boot.times
attr(tauc, 'tauc.type') = tauc.type
attr(tauc, 'tauc.time') = tauc.time
@@ -418,7 +422,8 @@ hdnom.validate = function(x, time, event,
class(tauc) = c('hdnom.validate',
'penalized.validate.cv')
attr(tauc, 'model.type') = model.type
attr(tauc, 'lambda') = lambda
attr(tauc, 'lambda1') = lambda1
attr(tauc, 'lambda2') = lambda2
attr(tauc, 'nfolds') = nfolds
attr(tauc, 'tauc.type') = tauc.type
attr(tauc, 'tauc.time') = tauc.time
@@ -461,7 +466,8 @@ hdnom.validate = function(x, time, event,
class(tauc) = c('hdnom.validate',
'penalized.validate.repeated.cv')
attr(tauc, 'model.type') = model.type
attr(tauc, 'lambda') = lambda
attr(tauc, 'lambda1') = lambda1
attr(tauc, 'lambda2') = lambda2
attr(tauc, 'nfolds') = nfolds
attr(tauc, 'rep.times') = rep.times
attr(tauc, 'tauc.type') = tauc.type
@@ -596,11 +602,12 @@ ncvreg.validate.internal = function(x_tr, x_te, y_tr, y_te, model.type,
#'
#' @keywords internal
penalized.validate.internal = function(x_tr, x_te, y_tr, y_te,
lambda,
lambda1, lambda2,
tauc.type, tauc.time) {

samp_fit = penalized(response = y_tr, penalized = x_tr,
lambda1 = lambda, lambda2 = 0,
lambda1 = lambda1, lambda2 = lambda2,
maxiter = 25, epsilon = 1e-3, # for faster convergence, consistent with `hdcox.flasso()`
fusedl = TRUE, standardize = TRUE, model = 'cox')

lp_tr = as.vector(samp_fit@lin.pred)
@@ -748,7 +755,8 @@ print.hdnom.validate = function(x, ...) {
cat('Validation method: bootstrap\n')
cat('Bootstrap samples:', attr(x, 'boot.times'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Time-dependent AUC type:', attr(x, 'tauc.type'), '\n')
cat('Evaluation time points for tAUC:', attr(x, 'tauc.time'))
},
@@ -759,7 +767,8 @@ print.hdnom.validate = function(x, ...) {
cat('Validation method: k-fold cross-validation\n')
cat('Cross-validation folds:', attr(x, 'nfolds'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Time-dependent AUC type:', attr(x, 'tauc.type'), '\n')
cat('Evaluation time points for tAUC:', attr(x, 'tauc.time'))
},
@@ -771,7 +780,8 @@ print.hdnom.validate = function(x, ...) {
cat('Cross-validation folds:', attr(x, 'nfolds'), '\n')
cat('Cross-validation repeated times:', attr(x, 'rep.times'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Time-dependent AUC type:', attr(x, 'tauc.type'), '\n')
cat('Evaluation time points for tAUC:', attr(x, 'tauc.time'))
}
@@ -22,7 +22,10 @@
#' model fits on the resampled data. From the Cox model you have built.
#' @param pen.factor Penalty factors to apply to each coefficient.
#' From the built \emph{adaptive lasso} or \emph{adaptive elastic-net} model.
#' @param gamma Value of the model parameter gamma for MCP/SCAD/Mnet/Snet models.
#' @param gamma Value of the model parameter gamma for
#' MCP/SCAD/Mnet/Snet models.
#' @param lambda1 Value of the penalty parameter lambda1 for fused lasso model.
#' @param lambda2 Value of the penalty parameter lambda2 for fused lasso model.
#' @param method Calibration method.
#' Options including \code{"fitting"}, \code{"bootstrap"}, \code{"cv"},
#' and \code{"repeated.cv"}.
@@ -53,11 +56,11 @@
#' y = Surv(time, event)
#'
#' # Fit Cox model with lasso penalty
#' lassofit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#' fit = hdcox.lasso(x, y, nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' # Model calibration by fitting the original data directly
#' cal.fitting = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "fitting",
#' pred.at = 365 * 9, ngroup = 5,
#' seed = 1010)
@@ -66,21 +69,21 @@
#' # Normally boot.times should be set to 200 or more,
#' # we set it to 3 here only to save example running time.
#' cal.boot = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "bootstrap", boot.times = 3,
#' pred.at = 365 * 9, ngroup = 5,
#' seed = 1010)
#'
#' # Model calibration by 5-fold cross-validation
#' cal.cv = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "cv", nfolds = 5,
#' pred.at = 365 * 9, ngroup = 5,
#' seed = 1010)
#'
#' # Model calibration by repeated cross-validation
#' cal.repcv = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$lasso_best_lambda,
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "repeated.cv", nfolds = 3, rep.times = 3,
#' pred.at = 365 * 9, ngroup = 5,
#' seed = 1010)
@@ -113,7 +116,7 @@
#
# set.seed(1010)
# cal.fitting = hdnom.calibrate(x, time, event, model.type = "flasso",
# lambda = 60,
# lambda1 = 5, lambda2 = 2,
# method = "fitting",
# pred.at = 365 * 9, ngroup = 5,
# seed = 1010)
@@ -131,7 +134,7 @@
# seed = 1010)
#
# cal.repcv = hdnom.calibrate(x, time, event, model.type = "flasso",
# lambda = 60,
# lambda1 = 5, lambda2 = 2,
# method = "repeated.cv", nfolds = 5, rep.times = 3,
# pred.at = 365 * 9, ngroup = 5,
# seed = 1010)
@@ -157,6 +160,7 @@ hdnom.calibrate = function(x, time, event,
'mcp', 'mnet',
'scad', 'snet'),
alpha, lambda, pen.factor = NULL, gamma,
lambda1, lambda2,
method = c('fitting', 'bootstrap', 'cv', 'repeated.cv'),
boot.times = NULL, nfolds = NULL, rep.times = NULL,
pred.at, ngroup = 5,
@@ -195,7 +199,7 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
pred_list = penalized.calibrate.internal.pred(
x_tr = x, x_te = x, y_tr = Surv(time, event),
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
pred.at = pred.at
)
}
@@ -258,7 +262,7 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
pred_list_list[[i]] = penalized.calibrate.internal.pred(
x_tr = x_tr, x_te = x_te, y_tr = y_tr,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
pred.at = pred.at
)
}
@@ -331,7 +335,7 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
pred_list_list[[i]] = penalized.calibrate.internal.pred(
x_tr = x_tr, x_te = x_te, y_tr = y_tr,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
pred.at = pred.at
)
}
@@ -412,7 +416,7 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
pred_list_list[[j]][[i]] = penalized.calibrate.internal.pred(
x_tr = x_tr, x_te = x_te, y_tr = y_tr,
lambda = lambda,
lambda1 = lambda1, lambda2 = lambda2,
pred.at = pred.at
)
}
@@ -466,23 +470,24 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('lasso', 'alasso', 'enet', 'aenet')) {
class(prob) = c('hdnom.calibrate',
'glmnet.calibrate.fitting')
attr(prob, 'alpha') = alpha
attr(prob, 'lambda') = lambda
attr(prob, 'alpha') = alpha
attr(prob, 'lambda') = lambda
attr(prob, 'pen.factor') = pen.factor
}

if (model.type %in% c('mcp', 'mnet', 'scad', 'snet')) {
class(prob) = c('hdnom.calibrate',
'ncvreg.calibrate.fitting')
attr(prob, 'alpha') = alpha
attr(prob, 'alpha') = alpha
attr(prob, 'lambda') = lambda
attr(prob, 'gamma') = gamma
attr(prob, 'gamma') = gamma
}

if (model.type %in% c('flasso')) {
class(prob) = c('hdnom.calibrate',
'penalized.calibrate.fitting')
attr(prob, 'lambda') = lambda
attr(prob, 'lambda1') = lambda1
attr(prob, 'lambda2') = lambda2
}

},
@@ -510,7 +515,8 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
class(prob) = c('hdnom.calibrate',
'penalized.calibrate.bootstrap')
attr(prob, 'lambda') = lambda
attr(prob, 'lambda1') = lambda1
attr(prob, 'lambda2') = lambda2
attr(prob, 'boot.times') = boot.times
}

@@ -539,7 +545,8 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
class(prob) = c('hdnom.calibrate',
'penalized.calibrate.cv')
attr(prob, 'lambda') = lambda
attr(prob, 'lambda1') = lambda1
attr(prob, 'lambda2') = lambda2
attr(prob, 'nfolds') = nfolds
}

@@ -570,7 +577,8 @@ hdnom.calibrate = function(x, time, event,
if (model.type %in% c('flasso')) {
class(prob) = c('hdnom.calibrate',
'penalized.calibrate.repeated.cv')
attr(prob, 'lambda') = lambda
attr(prob, 'lambda1') = lambda1
attr(prob, 'lambda2') = lambda2
attr(prob, 'nfolds') = nfolds
attr(prob, 'rep.times') = rep.times
}
@@ -705,11 +713,12 @@ ncvreg.calibrate.internal.pred = function(x_tr, x_te, y_tr,
#'
#' @keywords internal
penalized.calibrate.internal.pred = function(x_tr, x_te, y_tr,
lambda,
lambda1, lambda2,
pred.at) {

object = penalized(response = y_tr, penalized = x_tr,
lambda1 = lambda, lambda2 = 0,
lambda1 = lambda1, lambda2 = lambda2,
maxiter = 25, epsilon = 1e-3, # for faster convergence, consistent with `hdcox.flasso()`
fusedl = TRUE, standardize = TRUE, model = 'cox')

lp = as.vector(data.matrix(x_tr) %*% as.matrix(object@penalized))
@@ -913,7 +922,8 @@ print.hdnom.calibrate = function(x, ...) {
cat('Random seed:', attr(x, 'seed'), '\n')
cat('Calibration method: fitting\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Calibration time point:', attr(x, 'pred.at'), '\n')
cat('Number of groups formed for calibration:', attr(x, 'ngroup'), '\n')
},
@@ -924,7 +934,8 @@ print.hdnom.calibrate = function(x, ...) {
cat('Calibration method: bootstrap\n')
cat('Bootstrap samples:', attr(x, 'boot.times'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Calibration time point:', attr(x, 'pred.at'), '\n')
cat('Number of groups formed for calibration:', attr(x, 'ngroup'), '\n')
},
@@ -935,7 +946,8 @@ print.hdnom.calibrate = function(x, ...) {
cat('Calibration method: k-fold cross-validation\n')
cat('Cross-validation folds:', attr(x, 'nfolds'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Calibration time point:', attr(x, 'pred.at'), '\n')
cat('Number of groups formed for calibration:', attr(x, 'ngroup'), '\n')
},
@@ -947,7 +959,8 @@ print.hdnom.calibrate = function(x, ...) {
cat('Cross-validation folds:', attr(x, 'nfolds'), '\n')
cat('Cross-validation repeated times:', attr(x, 'rep.times'), '\n')
cat('Model type:', attr(x, 'model.type'), '\n')
cat('Fused lasso model lambda:', attr(x, 'lambda'), '\n')
cat('Fused lasso model lambda1:', attr(x, 'lambda1'), '\n')
cat('Fused lasso model lambda2:', attr(x, 'lambda2'), '\n')
cat('Calibration time point:', attr(x, 'pred.at'), '\n')
cat('Number of groups formed for calibration:', attr(x, 'ngroup'), '\n')
}
@@ -1036,23 +1049,27 @@ summary.hdnom.calibrate = function(object, ...) {
},

penalized.calibrate.fitting = {
attr(object, 'lambda') = NULL
attr(object, 'lambda1') = NULL
attr(object, 'lambda2') = NULL
},

penalized.calibrate.bootstrap = {
attr(object, 'lambda') = NULL
attr(object, 'boot.times') = NULL
attr(object, 'lambda1') = NULL
attr(object, 'lambda2') = NULL
attr(object, 'boot.times') = NULL
},

penalized.calibrate.cv = {
attr(object, 'lambda') = NULL
attr(object, 'nfolds') = NULL
penalized.calibrate.cv = {
attr(object, 'lambda1') = NULL
attr(object, 'lambda2') = NULL
attr(object, 'nfolds') = NULL
},

penalized.calibrate.repeated.cv = {
attr(object, 'lambda') = NULL
attr(object, 'nfolds') = NULL
attr(object, 'rep.times') = NULL
attr(object, 'lambda1') = NULL
attr(object, 'lambda2') = NULL
attr(object, 'nfolds') = NULL
attr(object, 'rep.times') = NULL
}

)
@@ -57,12 +57,12 @@
#' time_new = smart$TEVENT[1001:2000]
#' event_new = smart$EVENT[1001:2000]
#'
#' # Fit Cox model by lasso penalization
#' lassofit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#' # Fit Cox model with lasso penalty
#' fit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' # External validation with time-dependent AUC
#' val.ext =
#' hdnom.external.validate(lassofit, x, time, event,
#' hdnom.external.validate(fit, x, time, event,
#' x_new, time_new, event_new,
#' tauc.type = "UNO",
#' tauc.time = seq(0.25, 2, 0.25) * 365)
@@ -41,12 +41,12 @@
#' time_new = smart$TEVENT[1001:2000]
#' event_new = smart$EVENT[1001:2000]
#'
#' # Fit Cox model by lasso penalization
#' lassofit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#' # Fit Cox model with lasso penalty
#' fit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' # External calibration
#' cal.ext =
#' hdnom.external.calibrate(lassofit, x, time, event,
#' hdnom.external.calibrate(fit, x, time, event,
#' x_new, time_new, event_new,
#' pred.at = 365 * 5, ngroup = 5)
#'
@@ -37,21 +37,21 @@
#' time_new = smart$TEVENT[1001:2000]
#' event_new = smart$EVENT[1001:2000]
#'
#' # Fit Cox model by lasso penalization
#' lassofit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#' # Fit Cox model with lasso penalty
#' fit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' ### Internal calibration
#' # Internal calibration
#' cal.int = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$'lasso_best_lambda',
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "cv", nfolds = 5,
#' pred.at = 365 * 9, ngroup = 3)
#'
#' hdnom.kmplot(cal.int, group.name = c('High risk', 'Medium risk', 'Low risk'),
#' time.at = 1:6 * 365)
#'
#' ### External calibration
#' # External calibration
#' cal.ext =
#' hdnom.external.calibrate(lassofit, x, time, event,
#' hdnom.external.calibrate(fit, x, time, event,
#' x_new, time_new, event_new,
#' pred.at = 365 * 5, ngroup = 3)
#'
@@ -29,20 +29,20 @@
#' time_new = smart$TEVENT[1001:2000]
#' event_new = smart$EVENT[1001:2000]
#'
#' # Fit Cox model by lasso penalization
#' lassofit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#' # Fit Cox model with lasso penalty
#' fit = hdcox.lasso(x, Surv(time, event), nfolds = 5, rule = "lambda.1se", seed = 11)
#'
#' ### Internal calibration
#' # Internal calibration
#' cal.int = hdnom.calibrate(x, time, event, model.type = "lasso",
#' alpha = 1, lambda = lassofit$'lasso_best_lambda',
#' alpha = 1, lambda = fit$lasso_best_lambda,
#' method = "cv", nfolds = 5,
#' pred.at = 365 * 9, ngroup = 3)
#'
#' hdnom.logrank(cal.int)
#'
#' ### External calibration
#' # External calibration
#' cal.ext =
#' hdnom.external.calibrate(lassofit, x, time, event,
#' hdnom.external.calibrate(fit, x, time, event,
#' x_new, time_new, event_new,
#' pred.at = 365 * 5, ngroup = 3)
#'
@@ -52,12 +52,12 @@ hdnom.logrank = function(object) {
if (!(any(c('hdnom.calibrate', 'hdnom.external.calibrate') %in% class(object))))
stop('object class must be "hdnom.calibrate" or "hdnom.external.calibrate"')

time = attr(object, 'surv.time')
time = attr(object, 'surv.time')
event = attr(object, 'surv.event')
grp = attr(object, 'risk.group')
grp = attr(object, 'risk.group')

sdiff = survdiff(formula('Surv(time, event) ~ grp'))
pval = pchisq(sdiff$'chisq', length(sdiff$'n') - 1L, lower.tail = FALSE)
pval = pchisq(sdiff$'chisq', length(sdiff$'n') - 1L, lower.tail = FALSE)
sdiff$'pval' = pval

return(sdiff)
@@ -151,7 +151,8 @@ hdnom.compare.validate =
tauclist[[i]] =
hdnom.validate(
x, time, event, model.type = 'flasso',
alpha = 1, lambda = cvfit$'flasso_best_lambda',
lambda1 = cvfit$'flasso_best_lambda1',
lambda2 = cvfit$'flasso_best_lambda2',
method = method, boot.times = boot.times, nfolds = nfolds, rep.times = rep.times,
tauc.type = tauc.type, tauc.time = tauc.time,
seed = seed, trace = trace)
@@ -382,7 +383,7 @@ plot.hdnom.compare.validate =
geom_point(data = df, aes_string(x = 'Time', y = 'Median',
colour = 'Model', fill = 'Model')) +
geom_line(data = df, aes_string(x = 'Time', y = 'Median',
colour = 'Model', fill = 'Model'),
colour = 'Model'),
linetype = 'dashed') +
scale_x_continuous(breaks = df$'Time') +
scale_colour_manual(values = col_pal) +
@@ -399,7 +400,7 @@ plot.hdnom.compare.validate =
geom_point(data = df, aes_string(x = 'Time', y = 'Median',
colour = 'Model', fill = 'Model')) +
geom_line(data = df, aes_string(x = 'Time', y = 'Median',
colour = 'Model', fill = 'Model'),
colour = 'Model'),
linetype = 'dashed') +
geom_ribbon(data = df, aes_string(ymin = 'Qt25', ymax = 'Qt75',
colour = 'Model', fill = 'Model'),
@@ -133,7 +133,8 @@ hdnom.compare.calibrate =
problist[[i]] =
hdnom.calibrate(
x, time, event, model.type = 'flasso',
alpha = 1, lambda = cvfit$'flasso_best_lambda',
lambda1 = cvfit$'flasso_best_lambda1',
lambda2 = cvfit$'flasso_best_lambda2',
method = method, boot.times = boot.times, nfolds = nfolds, rep.times = rep.times,
pred.at = pred.at, ngroup = ngroup,
seed = seed, trace = trace)
@@ -1,8 +1,8 @@
#' Nomograms for High-Dimensional Data with Penalized Cox Models
#' Benchmarking and Visualization Toolkit for Penalized Cox Models
#'
#' Nomograms for High-Dimensional Data with Penalized Cox Models
#' Benchmarking and Visualization Toolkit for Penalized Cox Models
#'
#' The vignette can be opened with \code{vignette('hdnom')}.
#' Browse vignettes with \code{browseVignettes(package = "hdnom")}.
#'
#' \tabular{ll}{ Package: \tab hdnom\cr
#' Type: \tab Package\cr
BIN -1.12 KB vignettes/aenetfit.rds
Binary file not shown.
BIN +1.13 KB vignettes/fit.rds
Binary file not shown.
@@ -89,13 +89,13 @@ Fit a penalized Cox model by adaptive elastic-net regularization,
with `hdcox.aenet`:

```{r, eval = FALSE}
# Enable parallel parameter tuning
# enable parallel parameter tuning
suppressMessages(library("doParallel"))
registerDoParallel(detectCores())
aenetfit = hdcox.aenet(x, y, nfolds = 10, rule = "lambda.1se",
seed = c(5, 7), parallel = TRUE)
names(aenetfit)
fit = hdcox.aenet(x, y, nfolds = 10, rule = "lambda.1se",
seed = c(5, 7), parallel = TRUE)
names(fit)
```

```
@@ -105,14 +105,14 @@ names(aenetfit)
```

```{r, echo = FALSE}
aenetfit = readRDS("aenetfit.rds")
fit = readRDS("fit.rds")
```

Adaptive elastic-net includes two estimation steps. The random seed used
for parameter tuning, the selected best $\alpha$, the selected best $\lambda$,
the model fitted for each estimation step, and the penalty factor for the model
coefficients in the second estimation step are all stored in the model object
`aenetfit`.
`fit`.

# Nomogram Plotting

@@ -121,10 +121,10 @@ about the model, namely, the model object and parameters,
from the result of the last step:

```{r}
fit = aenetfit$aenet_model
alpha = aenetfit$aenet_best_alpha
lambda = aenetfit$aenet_best_lambda
adapen = aenetfit$pen_factor
model = fit$aenet_model
alpha = fit$aenet_best_alpha
lambda = fit$aenet_best_lambda
adapen = fit$pen_factor
```

To plot the nomogram, first we make `x` available as a `datadist` object
@@ -137,8 +137,8 @@ x.df = as.data.frame(x)
dd = datadist(x.df)
options(datadist = "dd")
nom = hdnom.nomogram(fit, model.type = "aenet", x, time, event, x.df,
lambda = lambda, pred.at = 365 * 2,
nom = hdnom.nomogram(model, model.type = "aenet",
x, time, event, x.df, pred.at = 365 * 2,
funlabel = "2-Year Overall Survival Probability")
plot(nom)
```
@@ -231,7 +231,7 @@ event_new = smart$EVENT[1001:2000]
# External validation with time-dependent AUC
val.ext =
hdnom.external.validate(aenetfit, x, time, event,
hdnom.external.validate(fit, x, time, event,
x_new, time_new, event_new,
tauc.type = "UNO",
tauc.time = seq(0.25, 2, 0.25) * 365)
@@ -293,7 +293,7 @@ use `hdnom.external.calibrate()`:

```{r, fig.width = 8, fig.height = 8, out.width = 600, out.height = 600}
cal.ext =
hdnom.external.calibrate(aenetfit, x, time, event,
hdnom.external.calibrate(fit, x, time, event,
x_new, time_new, event_new,
pred.at = 365 * 5, ngroup = 3)
@@ -423,7 +423,7 @@ the `smart` dataset as the new samples, and predict their overall
survival probability from one year to ten years:

```{r}
predict(aenetfit, x, y, newx = x[101:105, ], pred.at = 1:10 * 365)
predict(fit, x, y, newx = x[101:105, ], pred.at = 1:10 * 365)
```

# Customize Color Palette