R/Score.R

Defines functions AUC.competing.risks AUC.survival AUC.binary auRoc.factor auRoc.numeric delongtest AireTrap Brier.competing.risks Brier.survival Brier.binary Score Score.list

Documented in Score Score.list

# {{{ roxy header
##' Methods to score the predictive performance of risk markers and risk prediction models
##'
##' The function implements a toolbox for the risk prediction modeller:
##' all tools work for the three outcomes: (1) binary (uncensored),
##' (2) right censored time to event without competing risks,
##' (3) right censored time to event with competing risks
##'
##' Computed are the (time-dependent) Brier score and the (time-dependent)
##' area under the ROC curve for a list of risk prediction models either in
##' external validation data or in the learning data using bootstrap
##' cross-validation. The function optionally provides results for plotting (time-point specific)
##' ROC curves, for (time-point specific) calibration curves and for (time-point specific) retrospective boxplots.
##'
##' For uncensored binary outcome the Delong-Delong test is used to contrast AUC of rival models.
##' In right censored survival data (with and without competing risks)
##' the p-values correspond to Wald tests based on standard errors obtained with an estimate of the influence function
##' as described in detail in the appendix of Blanche et al. (2015).
##' @title Score risk predictions
##' @aliases Score Score.list
##' @param object List of risk predictions (see details and examples).
##' @param formula A formula which identifies the outcome (left hand
##'     side). E.g., \code{Y ~ 1} for binary and \code{Hist(time,status) ~ 1} for time-to-event outcome.
##' In right censored data, the right hand side of the
##'     formula is used to estimate the inverse probability of censoring weights (IPCW) model.
##' @param data \code{data.frame} or \code{data.table} in which the formula can be
##'     interpreted.
##' @param metrics Character vector specifying which metrics to
##'     apply. Case does not matter. Choices are \code{"AUC"} and \code{"Brier"}.
##' @param summary Character vector specifying which summary
##'     statistics to apply to the predicted risks. Choices are \code{"risks"}, \code{"IPA"},
##'     \code{"riskQuantile"} and \code{"ibs"}. Can be all \code{c("risks","IPA","riskQuantile","ibs")} or a subset thereof.
##'     \itemize{
##'     \item \code{"risks"} adds the predicted risks to the output.
##'     \item \code{"ipa"} computes the index of prediction accuracy (AKA R-squared) based on Brier scores for model vs null model
##'     \item \code{"riskQuantile"} calculates
##'     time-point specific boxplots for the
##'     predicted risks (or biomarker values) conditional on the outcome at the time-point.
##'     \item \code{"ibs"} calculates integrated Brier scores across the time points at which the Brier score is computed. This works only with
##'     time-to-event outcome and the results depend on the argument \code{times}.
##'     }
##'     Set to \code{NULL} to avoid estimation of summary statistics.
##' @param plots Character vector specifying for which plots to put data into the result.
##'     Currently implemented are \code{"ROC"}, \code{"Calibration"} and \code{"boxplot"}.
##'     In addition, one can plot AUC and Brier score as function of time as soon as
##'     \code{times} has at least two different values.
##' @param cause Event of interest. Used for binary outcome \code{Y}
##'     to specify that risks are risks of the event \code{Y=event}
##'     and for competing risks outcome to specify the cause of
##'     interest.
##' @param times For survival and competing risks outcome: list of
##'     prediction horizons. All times which are greater than the
##'     maximal observed time in the data set are automatically removed.
##'     Note that the object returned by the function may become huge when
##'     the prediction performance is estimated at many prediction horizons.
##' @param landmarks Not yet implemented.
##' @param use.event.times If \code{TRUE} merge all unique event times with
##'     the vector given by argument \code{times}.
##' @param null.model If \code{TRUE} fit a risk prediction model which ignores
##'     the covariates and predicts the same value for all subjects. The model is fitted using \code{data}
##' and the left hand side of \code{formula}. For binary outcome this is just the empirical prevalence. For (right censored) time to event outcome, the null models are
##' equal to the Kaplan-Meier estimator (no competing risks) and the Aalen-Johansen estimator (with competing risks).
##' @param se.fit Logical or \code{0} or \code{1}. If \code{FALSE} or \code{0} do not calculate standard errors.
##' @param conservative Logical, only relevant in right censored data. If \code{TRUE} ignore
##'        variability of the estimate of the inverse probability of censoring weights when calculating standard
##'        errors for prediction performance parameters. This can potentially reduce computation time and memory usage
##'        at a usually very small expense of a slightly higher standard error.
##' @param multi.split.test Logical or \code{0} or \code{1}. If \code{FALSE} or \code{0} do not calculate multi-split tests. This argument is ignored when \code{split.method} is "none".
##' @param conf.int Either logical or a numeric value between 0 and 1. In right censored data,
##'     confidence intervals are based on Blanche et al (see references). Setting \code{FALSE} prevents the
##'     computation confidence intervals. \code{TRUE} means compute 95 percent confidence
##'     intervals and corresponding p-values for AUC and Brier score. If set to 0.87, the
##'     level of significance is 13 percent. So, do not set it to 0.87.
##' @param contrasts Either logical or a list of contrasts. A list of contrasts defines which risk prediction models (markers)
##'    should be contrasted with respect to their prediction performance.
##'   If \code{TRUE} do all possible comparisons. For
##'     example, when \code{object} is a list with two risk prediction models and
##'     \code{null.model=TRUE} setting \code{TRUE} is equivalent to
##'     \code{list(c(0,1,2),c(1,2))} where \code{c(0,1,2)} codes for the
##'     two comparisons: 1 vs 0 and 2 vs 0 (positive integers refer to
##'     elements of \code{object}, 0 refers to the benchmark null
##'     model which ignores the covariates).  This again is equivalent
##'     to explicitly setting \code{list(c(0,1),c(0,2),c(1,2))}. A
##'     more complex example: Suppose \code{object} has 7 elements and you
##'     want to do the following 3 comparisons: 6 vs 3, 2 vs 5 and 2
##'     vs 3, you should set \code{contrasts=c(6,3),c(2,5,3)}.
##' @param probs Quantiles for retrospective summary statistics of the
##'     predicted risks. This affects the result of the function \code{boxplot.Score}.
##' @param cens.method Method for dealing with right censored
##'     data. Either \code{"ipcw"} or \code{"pseudo"}.
##'     Here IPCW refers to inverse probability of censoring weights and \code{pseudo} for jackknife pseudo values.
##'     Right now pseudo values  are only used for calibration curves.
##' @param cens.model Model for estimating inverse probability of
##'     censored weights. Implemented are the Kaplan-Meier method (\code{"km"}) and
##' Cox regression (\code{"cox"}) both applied to the censored times. If the right hand side of \code{formula} does not specify covariates,
##' the Kaplan-Meier method is used even if this argument is set to \code{"cox"}.
##' @param split.method Method for cross-validation. Right now the only choice is \code{bootcv} in which case bootstrap learning sets
##' are drawn with our without replacement (argument \code{M}) from \code{data}. The data not included in the current bootstrap learning
##' set are used as validation set to compute the prediction performance.
##' @param B Number of bootstrap sets for cross-validation.
##' @param M Size of subsamples for bootstrap cross-validation. If specified it
##'     has to be an integer smaller than the size of \code{data}.
##' @param seed Super seed for setting training data seeds when
##'     randomly splitting (bootstrapping) the data during cross-validation.
##' @param trainseeds Seeds for training models during cross-validation.
##' @param parallel The type of parallel operation to be used (if any). If missing, the default is \code{"no"}.
##' @param ncpus integer: number of processes to be used in parallel operation.
##' @param cl An optional \code{parallel} or \code{snow} cluster for use if \code{parallel = "snow"}. If not supplied, a cluster on the local machine is created for the duration of the \code{Score} call.
##' @param progress.bar Style for \code{txtProgressBar}. Can be 1,2,3 see \code{help(txtProgressBar)} or NULL to avoid the progress bar.
##' @param keep list of characters (not case sensitive) which determines additional output.
##' \code{"residuals"} provides Brier score residuals and
##' \code{"splitindex"} provides sampling index used to split the data into training and validation sets.
##' \code{"vcov"} provides the variance-covariance matrix of the estimated parameters.
##' @param predictRisk.args
##'  A list of argument-lists to control how risks are predicted.
##'  The names of the lists should be the S3-classes of the \code{object}.
##'  The argument-lists are then passed on to the S3-class specific predictRisk method.
##'  For example, if your object contains one or several random forest model fitted with the function randomForestSRC::rfsrc then you can
##'  specify additional arguments for the function riskRegression::predictRisk.rfsrc which will pass
##'  these on to the function randomForestSRC::predict.rfsrc. A specific example in this case would be
##'  \code{list(rfsrc=list(na.action="na.impute"))}.
##'
##'  A more flexible approach is to write a new predictRisk S3-method. See Details.
##' @param debug Logical. If \code{TRUE} indicate landmark in progress of the program.
##' @param useEventTimes obsolete.
##' @param nullModel obsolete.
##' @param censMethod obsolete.
##' @param censModel obsolete.
##' @param splitMethod obsolete.
##' @param ... Named list containging additional arguments that are passed on to the \code{predictRisk} methods corresponding to object. See examples.
##' @return List with scores and assessments of contrasts, i.e.,
##'     tests and confidence limits for performance and difference in performance (AUC and Brier),
##'     summaries and plots. Most elements are in\code{data.table} format.
##' @details
##' This function works with one or multiple models that predict the risk of an event R(t|X) for a subject
##' characterized by predictors X at time t.
##' With binary endpoints (outcome 0/1 without time component) the risk is simply R(X).
##' In case of a survival object
##' without competing risks the function still works with predicted event probabilities, i.e., R(t|X)=1-S(t|X) where
##' S(t|X) is the predicted survival chance for subject X at time t.
##'
##' The already existing predictRisk methods (see methods(predictRisk)) may not cover all models and methods
##' for predicting risks. But users can quickly extend the package as explained in detail in Mogensen et al. (2012) for
##' the predecessors \code{pec::predictSurvProb} and \code{pec::predictEventProb} which have been unified as
##' \code{riskRegression::predictRisk}.
##'
##' Bootstrap Crossvalidation (see also Gerds & Schumacher 2007 and Mogensen et al. 2012)
##'
##' B=10, M (not specified or M=NROW(data))
##' Training of each of the models in each of 10 bootstrap data sets (learning data sets).
##' Learning data sets are obtained by sampling \code{NROW(data)} subjects of the data set
##' with replacement. There are roughly \code{.632*NROW(data)} subjects in the learning data (inbag)
##' and \code{.368*NROW(data)} subjects not in the validation data sets (out-of-bag).
##'
##' These are used to estimate the scores: AUC, Brier, etc. Reported are averages across the 10 splits.
##'
##' ## Bootstrap with replacement
##' \code{
##' set.seed(13)
##' N=17
##' data = data.frame(id=1:N, y=rbinom(N,1,.3),x=rnorm(N))
##' boot.index = sample(1:N,size=N,replace=TRUE)
##' boot.index
##' inbag = 1:N %in% boot.index
##' outofbag = !inbag
##' learn.data = data[inbag]
##' val.data = data[outofbag]
##' riskRegression:::getSplitMethod("bootcv",B=10,N=17)$index
##' }
##' NOTE: the number .632 is the expected probability to draw one subject (for example subject 1) with
##' replacement from the data, which does not depend on the sample size:
##' \code{B=10000}
##' \code{N=137}
##' \code{mean(sapply(1:B, function(b){match(1,sample(1:N,size=N,replace=TRUE),nomatch=0)}))}
##' \code{N=30}
##' \code{mean(sapply(1:B, function(b){match(1,sample(1:N,size=N,replace=TRUE),nomatch=0)}))}
##' \code{N=300}
##' \code{mean(sapply(1:B, function(b){match(1,sample(1:N,size=N,replace=TRUE),nomatch=0)}))}
##'
##'
##' ## Bootstrap without replacement (training size set to be 70 percent of data)
##' B=10, M=.7
##'
##' Training of each of the models in each of 10 bootstrap data sets (learning data sets).
##' Learning data sets are obtained by sampling \code{round(.8*NROW(data))} subjects of the data set
##' without replacement. There are \code{NROW(data)-round(.8*NROW(data))} subjects not in the learning data sets.
##' These are used to estimate the scores: AUC, Brier, etc. Reported are averages across the 10 splits.
##' \code{
##' set.seed(13)
##' N=17
##' data = data.frame(id=1:N, y=rbinom(N,1,.3),x=rnorm(N))
##' boot.index = sample(1:N,size=M,replace=FALSE)
##' boot.index
##' inbag = 1:N %in% boot.index
##' outofbag = !inbag
##' learn.data = data[inbag]
##' val.data = data[outofbag]
##' riskRegression:::getSplitMethod("bootcv",B=10,N=17,M=.7)$index
##' }
##'
##' @examples
##' # binary outcome
##' library(lava)
##' set.seed(18)
##' learndat <- sampleData(48,outcome="binary")
##' testdat <- sampleData(40,outcome="binary")
##'
##' ## score logistic regression models
##' lr1 = glm(Y~X1+X2+X7+X9,data=learndat,family=binomial)
##' lr2 = glm(Y~X3+X5,data=learndat,family=binomial)
##' Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5)"=lr2),formula=Y~1,data=testdat)
##'
##' ## ROC curve and calibration plot
##' xb=Score(list("LR(X1+X2+X7+X9)"=lr1,"LR(X3+X5+X6)"=lr2),formula=Y~1,
##'          data=testdat,plots=c("calibration","ROC"))
##' \dontrun{plotROC(xb)
##' plotCalibration(xb)
##' }
##'
##' ## compute AUC for a list of continuous markers
##' markers = as.list(testdat[,.(X6,X7,X8,X9,X10)])
##' Score(markers,formula=Y~1,data=testdat,metrics=c("auc"))
##'
##' # cross-validation
##' \dontrun{
##'     learndat=sampleData(400,outcome="binary")
##'     lr1a = glm(Y~X6,data=learndat,family=binomial)
##'     lr2a = glm(Y~X7+X8+X9,data=learndat,family=binomial)
##'     ## bootstrap cross-validation
##'     x1=Score(list("LR1"=lr1a,"LR2"=lr2a),formula=Y~1,data=learndat,split.method="bootcv",B=100)
##'     x1
##'     ## leave-one-out and leave-pair-out bootstrap
##'     x2=Score(list("LR1"=lr1a,"LR2"=lr2a),formula=Y~1,data=learndat,
##'              split.method="loob",
##'              B=100,plots="calibration")
##'     x2
##' }
##' # survival outcome
##'
##' # Score Cox regression models
##' \dontrun{library(survival)
##' library(rms)
##' library(prodlim)
##' set.seed(18)
##' trainSurv <- sampleData(100,outcome="survival")
##' testSurv <- sampleData(40,outcome="survival")
##' cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=trainSurv, y=TRUE, x = TRUE)
##' cox2 = coxph(Surv(time,event)~X3+X5+X6,data=trainSurv, y=TRUE, x = TRUE)
##' xs=Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'          formula=Surv(time,event)~1,data=testSurv,conf.int=FALSE,times=c(5,8))
##' xs
##' }
##'
##' # Integrated Brier score
##' \dontrun{
##' xs=Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'          formula=Surv(time,event)~1,data=testSurv,conf.int=FALSE,
##'          summary="ibs",
##'          times=sort(unique(testSurv$time)))
##' }
##'
##' # time-dependent AUC for list of markers
##' \dontrun{survmarkers = as.list(testSurv[,.(X6,X7,X8,X9,X10)])
##' Score(survmarkers,
##'       formula=Surv(time,event)~1,metrics="auc",data=testSurv,
##'       conf.int=TRUE,times=c(5,8))
##'
##' # compare models on test data
##' Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'       formula=Surv(time,event)~1,data=testSurv,conf.int=TRUE,times=c(5,8))
##' }
##' # crossvalidation models in traindata
##' \dontrun{
##'     library(survival)
##'     set.seed(18)
##'     trainSurv <- sampleData(400,outcome="survival")
##'     cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=trainSurv, y=TRUE, x = TRUE)
##'     cox2 = coxph(Surv(time,event)~X3+X5+X6,data=trainSurv, y=TRUE, x = TRUE)
##'     x1 = Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'                formula=Surv(time,event)~1,data=trainSurv,conf.int=TRUE,times=c(5,8),
##'                split.method="loob",B=100,plots="calibration")
##'
##'     x2= Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'               formula=Surv(time,event)~1,data=trainSurv,conf.int=TRUE,times=c(5,8),
##'               split.method="bootcv",B=100)
##' }
##'
##' # restrict number of comparisons
##' \dontrun{
##'     Score(list("Cox(X1+X2+X7+X9)"=cox1,"Cox(X3+X5+X6)"=cox2),
##'           formula=Surv(time,event)~1,data=trainSurv,contrasts=TRUE,
##'           null.model=FALSE,conf.int=TRUE,times=c(5,8),split.method="bootcv",B=3)
##'
##'     # competing risks outcome
##'     set.seed(18)
##'     trainCR <- sampleData(40,outcome="competing.risks")
##'     testCR <- sampleData(40,outcome="competing.risks")
##'     library(riskRegression)
##'     library(cmprsk)
##'     # Cause-specific Cox regression
##'     csc1 = CSC(Hist(time,event)~X1+X2+X7+X9,data=trainCR)
##'     csc2 = CSC(Hist(time,event)~X3+X5+X6,data=trainCR)
##'     # Fine-Gray regression
##'     fgr1 = FGR(Hist(time,event)~X1+X2+X7+X9,data=trainCR,cause=1)
##'     fgr2 = FGR(Hist(time,event)~X3+X5+X6,data=trainCR,cause=1)
##'     Score(list("CSC(X1+X2+X7+X9)"=csc1,"CSC(X3+X5+X6)"=csc2,
##'                "FGR(X1+X2+X7+X9)"=fgr1,"FGR(X3+X5+X6)"=fgr2),
##'           formula=Hist(time,event)~1,data=testCR,se.fit=1L,times=c(5,8))
##' }
##'
##'
##'
##' \dontrun{
##'     # reproduce some results of Table IV of Blanche et al. Stat Med 2013
##'     data(Paquid)
##'     ResPaquid <- Score(list("DSST"=-Paquid$DSST,"MMSE"=-Paquid$MMSE),
##'                        formula=Hist(time,status)~1,
##'                        data=Paquid,
##'                        null.model = FALSE,
##'                        conf.int=TRUE,
##'                        metrics=c("auc"),
##'                        times=c(3,5,10),
##'                        plots="ROC")
##'     ResPaquid
##'     plotROC(ResPaquid,time=5)
##' }
##' \dontrun{
##' # parallel options
##' # by erikvona: Here is a generic example of using future
##' # and doFuture, works great with the current version:
##' library(riskRegression)
##' library(future)
##' library(foreach)
##' library(doFuture)
##' library(survival)
##' # Register all available cores for parallel operation
##' plan(multiprocess, workers = availableCores())
##' registerDoFuture()
##' trainSurv <- sampleData(400,outcome="survival")
##' cox1 = coxph(Surv(time,event)~X1+X2+X7+X9,data=trainSurv,
##'              y=TRUE, x = TRUE)
##' # Bootstrapping on multiple cores
##' x1 = Score(list("Cox(X1+X2+X7+X9)"=cox1),
##'      formula=Surv(time,event)~1,data=trainSurv, times=c(5,8), 
##'      parallel = "as.registered", split.method="bootcv",B=100)
##' }
##'
##'
##' 
##' @author Thomas A Gerds \email{tag@@biostat.ku.dk} and Paul Blanche \email{paul.blanche@@univ-ubs.fr}
##' @references
##'
##' Ulla B. Mogensen, Hemant Ishwaran, Thomas A. Gerds (2012).
##' Evaluating Random Forests for Survival Analysis Using Prediction Error
##' Curves. Journal of Statistical Software, 50(11), 1-23. URL
##' http://www.jstatsoft.org/v50/i11/.
##'
##' Paul Blanche, Cecile Proust-Lima, Lucie Loubere, Claudine Berr, Jean- Francois Dartigues, and
##' Helene Jacqmin-Gadda. Quantifying and comparing dynamic predictive accuracy of joint models
##' for longitudinal marker and time-to-event in presence of censoring and competing risks.
##' Biometrics, 71 (1):102--113, 2015.
##'
##' P. Blanche, J-F Dartigues, and H. Jacqmin-Gadda. Estimating and comparing
##' time-dependent areas under receiver operating characteristic curves for
##' censored event times with competing risks. Statistics in Medicine,
##' 32(30):5381--5397, 2013.
##'
#' E. Graf et al.  (1999), Assessment and comparison of prognostic
#' classification schemes for survival data. Statistics in Medicine, vol 18,
#' pp= 2529--2545.
#'
#' Efron, Tibshirani (1997) Journal of the American Statistical Association 92,
#' 548--560 Improvement On Cross-Validation: The .632+ Bootstrap Method.
#'
#' Gerds, Schumacher (2006), Consistent estimation of the expected Brier score
#' in general survival models with right-censored event times. Biometrical
#' Journal, vol 48, 1029--1040.
#'
#' Thomas A. Gerds, Martin Schumacher (2007) Efron-Type Measures of Prediction
#' Error for Survival Analysis Biometrics, 63(4), 1283--1287
#' doi:10.1111/j.1541-0420.2007.00832.x
#'
#' Martin Schumacher, Harald Binder, and Thomas Gerds. Assessment of survival
#' prediction models based on microarray data. Bioinformatics, 23(14):1768-74,
#' 2007.
#'
#' Mark A. van de Wiel, Johannes Berkhof, and Wessel N. van Wieringen Testing
#' the prediction error difference between 2 predictors Biostatistics (2009)
#' 10(3): 550-560 doi:10.1093/biostatistics/kxp011
#'
#' Michael W Kattan and Thomas A Gerds. The index of prediction accuracy: an
#' intuitive measure useful for evaluating risk prediction models. Diagnostic
#' and Prognostic Research, 2(1):7, 2018.
##'
##' @export
##'
# }}}
# {{{ Score.list
Score.list <- function(object,
                       formula,
                       data,
                       metrics=c("auc","brier"),
                       summary=NULL, # riskQuantile, risk
                       plots= NULL, # c("roc","calibration","boxplot","p-values"),
                       cause,
                       times,
                       landmarks,
                       use.event.times=FALSE,
                       null.model=TRUE,
                       se.fit=TRUE,
                       conservative=FALSE,
                       multi.split.test=FALSE,
                       conf.int=.95,
                       contrasts=TRUE,
                       probs=c(0,0.25,0.5,0.75,1),
                       cens.method="ipcw",
                       cens.model="cox",
                       split.method,
                       B,
                       M,
                       seed,
                       trainseeds,
                       parallel=c("no","multicore","snow","as.registered"),
                       ncpus=1,
                       cl=NULL,
                       progress.bar=3,
                       keep,
                       predictRisk.args,
                       debug=0L,
                       useEventTimes,
                       nullModel,
                       censMethod,
                       censModel,
                       splitMethod,
                       ...){
    if (!missing(useEventTimes)) {
        warning("(Spelling of) argument 'useEventTimes' is obsolete and will be disabled in future versions of this function.\nUse 'use.event.times' instead.")
        use.event.times <- useEventTimes
    }
    if (!missing(nullModel)) {
        warning("(Spelling of) argument 'nullModel' is obsolete and will be disabled in future versions of this function.\nUse 'null.model' instead.")
        null.model <- nullModel
    }
    if (!missing(splitMethod)) {
        warning("(Spelling of) argument 'splitMethod' is obsolete and will be disabled in future versions of this function.\nUse 'split.method' instead.")
        split.method <- splitMethod
    }
    if (!missing(censMethod)) {
        warning("(Spelling of) argument 'censMethod' is obsolete and will be disabled in future versions of this function.\nUse 'cens.method' instead.")
        cens.method <- censMethod
    }
    if (!missing(censModel)) {
        warning("(Spelling of) argument 'censModel' is obsolete and will be disabled in future versions of this function.\nUse 'cens.model' instead.")
        cens.model <- censModel
    }

    se.conservative=IPCW=IF.AUC.conservative=IF.AUC0=IF.AUC=IC0=Brier=AUC=casecontrol=se=nth.times=time=status=ID=WTi=risk=IF.Brier=lower=upper=crossval=b=time=status=model=reference=p=model=pseudovalue=ReSpOnSe=residuals=event=j=NULL

    # }}}
    theCall <- match.call()
    # ----------------------------find metrics and plots ----------------------
    # {{{

    ## Metrics <- lapply(metrics,grep,c("AUC","Brier"),ignore.case=TRUE,value=TRUE)
    plots[grep("^box",plots,ignore.case=TRUE)] <- "boxplot"
    summary[grep("^riskQuantile",summary,ignore.case=TRUE)] <- "riskQuantile"
    # force riskQuantile to make boxplots
    if (("boxplot"%in% plots) && !("riskQuantile" %in% summary))
        summary <- c(summary,"riskQuantile")
    summary[grep("^risk$|^risks$",summary,ignore.case=TRUE)] <- "risks"
    metrics[grep("^auc$",metrics,ignore.case=TRUE)] <- "AUC"
    metrics[grep("^brier$",metrics,ignore.case=TRUE)] <- "Brier"
    plots[grep("^roc$",plots,ignore.case=TRUE)] <- "ROC"
    plots[grep("^cal",plots,ignore.case=TRUE)] <- "Calibration"
    if (length(posIBS <- grep("^ibs$|^crps$",summary,ignore.case=TRUE))>0){
        summary <- summary[-posIBS]
        if (!("Brier" %in% metrics)) metrics <- c(metrics,"Brier")
        ibs <- TRUE
    }else{
        ibs <- FALSE
    }
    ## IPA
    if (length(posRR <- grep("^ipa$|^rr$|^r2|rsquared$",summary,ignore.case=TRUE))>0){
        if (!null.model) stop("Need the null model to compute IPA/R^2 but argument 'null.model' is FALSE.")
        summary <- summary[-posRR]
        if (!("Brier" %in% metrics)) metrics <- c(metrics,"Brier")
        if (!null.model) {
            null.model <- TRUE
            warning("Value of argument 'null.model' ignored as the null model is needed to compute IPA/R^2.")
        }
        ipa <- TRUE
    }else{
        ipa <- FALSE
    }
    ## Plots <- lapply(plots,grep,c("Roc","Cal"),ignore.case=TRUE,value=TRUE)
    if ("ROC" %in% plots) {
        ## add AUC if needed
        if (!("AUC" %in% metrics)) metrics <- c(metrics,"AUC")
    }
    if ("Calibration" %in% plots) {
        ## add pseudo if needed
        if (!("pseudo" %in% cens.method)) cens.method <- c(cens.method,"pseudo")
    }

    # }}}
    # -----------------parse other arguments and prepare data---------
    # {{{ censoring model arguments
    if (length(grep("^km|^kaplan|^marg",cens.model,ignore.case=TRUE))>0){
        cens.model <- "KaplanMeier"
    } else{
        if (length(attr(terms(formula),"factors"))>0){
            cens.model <- "cox"
        } else{
            cens.model <- "KaplanMeier"
        }
    }
    # }}}
    # {{{ Response
    if (missing(formula)){stop("Argument formula is missing.")}
    formula.names <- try(all.names(formula),silent=TRUE)
    if (!(formula.names[1]=="~")
        ||
        (match("$",formula.names,nomatch=0)+match("[",formula.names,nomatch=0)>0)){
        stop("Invalid specification of formula.\n Could be that you forgot the right hand side:\n ~covariate1 + covariate2 + ...?\nNote that any subsetting, ie data$var or data[,\"var\"], is not supported by this function.")
    }
    if (missing(data)){stop("Argument data is missing.")}
    data <- data.table(data)
    responseFormula <- stats::update(formula,~1)
    ## if (missing(event)) event <- 1
    responsevars <- all.vars(responseFormula)
    if (!all(responsevars%in%names(data)))stop("Response variable(s) ",paste(responsevars,collapse=", ")," not found in data.")
    response <- getResponse(formula=responseFormula,
                            data=data,
                            cause=cause,
                            vars=responsevars)
    response.dim <- NCOL(response)
    response.type <- attr(response,"model")
    if (response.type %in% c("survival","competing.risks")){
        byvars <- c("model","times")
    } else{
        byvars <- c("model")
    }
    states <- attr(response,"states")
    if (missing(cause)||is.null(cause)){
        if (response.type=="binary"){
            cause <- attr(response,"event")
        }
        else{
            cause <- states[[1]]
        }
    }
    position.cause <- match(cause,states,nomatch=0)
    if (position.cause==0) stop(paste0("Requested cause: ",cause,". Available causes: ", paste(states,collapse=",")))
    # add null model and find names for the object
    if (null.model==TRUE){
        nullobject <- getNullModel(formula=formula,data=data,response.type=response.type)
    } else{
        nullobject <- NULL
    }
    ## put ReSpOnSe for binary and (time, event, status) in the first column(s)
    ## but first rename to avoid problems with pre-existing variables with same name(s)
    ## data[,eval(responsevars):=NULL]
    ## data have to be ordered when ipcw is called
    if (response.type=="binary"){
        if ("event" %in% names(data)){
            setnames(data,"event","protectedName.event")
        }
        data[,event:=factor(ReSpOnSe,
                            levels=1:(length(states)+1),
                            labels=c(states,"event-free"))]
    }
    if (response.type %in% c("survival","competing.risks")){
        test.names <- match(colnames(response),names(data),nomatch=0)
        if (sum(test.names)>0){
            protected.names <- names(data)[test.names]
            setnames(data,protected.names,paste0("protectedName.",protected.names))
        }
        data <- cbind(response,data)
        N <- as.numeric(NROW(data))
        neworder <- data[,order(time,-status)]
        data.table::setorder(data,time,-status)
    }else{
        data <- cbind(response,data)
        N <- as.numeric(NROW(data))
        neworder <- 1:N
    }
    ## add ID variable for merging purposes and because output has long format
    ## data[["ID"]]=1:N
    data[,ID:=1:N]
    if (response.type=="survival")
        formula <- stats::update(formula,"prodlim::Hist(time,status)~.")
    if (response.type=="competing.risks")
        formula <- stats::update(formula,"prodlim::Hist(time,event)~.")
    ## stop("Dont know how to predict response of type ",response.type))
    cens.type <- attr(response,"cens.type")
    if (is.null(cens.type)) cens.type <- "uncensored"
    rm(response)
    # }}}
    # {{{ SplitMethod & parallel stuff
    if (!missing(seed)) {
        ## message("Random seed set to control split of data: seed=",seed)
        set.seed(seed)
    }
    split.method <- getSplitMethod(split.method=split.method,B=B,N=N,M=M)
    B <- split.method$B
    splitIndex <- split.method$index
    do.resample <- !(is.null(splitIndex))
    if (split.method$internal.name!="noplan"){
        if (missing(parallel)) parallel <- "no"
        parallel <- match.arg(parallel)
        ## Additionally, I'd like the process to be adjusted when a cluster is
        ## passed or `as.registered` is specified to ignore the `ncpus` argument
        ## (very counter-intuitive that when I'm setting up the cluster I still
        ## need to specify ncpus > 1):
        if (ncpus<=1 && is.null(cl) && parallel != "as.registered") {            
            parallel <- "no"
        }
        ## if (ncpus <- pmin(ncpus,parallel::detectCores()))
        switch(parallel,"no"={
            foreach::registerDoSEQ()
        },"multicore"={
            if(.Platform$OS.type == "windows"){stop("multicore does not work on windows. set argument parallel to 'snow'.")}
            loadNamespace("parallel")
            if (is.null(cl)) cl <- parallel::makeForkCluster(nnodes=ncpus)
            doParallel::registerDoParallel(cl=cl,cores=ncpus)
            on.exit(parallel::stopCluster(cl))
        },"snow"={
            loadNamespace("parallel")
            if (is.null(cl)) cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
            on.exit(parallel::stopCluster(cl))
            doParallel::registerDoParallel(cl=cl,cores=ncpus)
        })
        if (split.method$name[1]=="BootCv" && multi.split.test[1]==TRUE){
            if ("AUC" %in% metrics) {
                warning("Cannot do multi-split test with AUC yet. Forced multi.split.test=FALSE")
                multi.split.test=FALSE
            }else{
                if (se.fit==TRUE & conservative==FALSE){
                    warning("Cannot do deal with conservative=FALSE when also multi.split.test=TRUE. Forced conservative=TRUE.")
                }
                conservative=TRUE
            }
        }
        if (se.fit==TRUE){
            if (response.type=="competing.risks")
                warning("Under construction. Check devtools::install_github('tagteam/riskRegression') for progress.")
        }
        if ("ROC" %in% plots){
            warning("Cannot (not yet) do ROC analysis in combination with internal validation\n. Check devtools::install_github('tagteam/riskRegression') for progress.")
        }
    }
    # }}}
    # {{{ Checking the ability of the elements of object to predict risks
    # {{{ number of models and their labels
    NF <- length(object)
    # }}}
    allmethods <- utils::methods(predictRisk)
    if (is.null(names(object))){
        names(object) <- sapply(object,function(o)class(o)[1])}
    else {
            names(object)[(names(object)=="")] <- sapply(object[(names(object)=="")],function(o)class(o)[1])
        }
    names.object <- names(object) <- make.unique(names(object))
    ## sanity checks
    object.classes <- lapply(object,function(x)class(x))
    lapply(1:NF,function(f){
        name <- names.object[[f]]
        if (any(c("integer","factor","numeric","matrix") %in% object.classes[[f]])){
            if (split.method$internal.name!="noplan")
                stop(paste0("Cannot crossvalidate performance of deterministic risk predictions:",name))
        }
        if (c("factor") %in% object.classes[[f]]){
            if (length(grep("^brier$",metrics,ignore.case=TRUE)>0) || length(grep("^cal",plots,ignore.case=TRUE)>0)){
                stop(paste0("Cannot compute Brier score or calibration plots for predictions that are coded as factor: ",name))
            }
        }
        candidateMethods <- paste("predictRisk",unlist(object.classes[[f]]),sep=".")
        if (all(match(candidateMethods,allmethods,nomatch=0)==0)){
            stop(paste("Cannot find function (S3-method) called ",
                       candidateMethods,
                       sep=""))
        }
    })
    # }}}
    # {{{ additional arguments for predictRisk methods
    if (!missing(predictRisk.args)){
        if (!(all(names(predictRisk.args) %in% unlist(object.classes))))
            stop(paste0("Argument predictRisk.args should be a list whose names match the S3-classes of the argument object.
  For example, if your object contains a random forest model fitted with the function randomForestSRC::rfsrc then you can
  specify additional arguments for the function riskRegression::predictRisk.rfsrc which will pass
  these on to the function randomForestSRC::predict.rfsrc. A specific example in this case would be
  predictRisk.args=list(\"rfsrc\"=list(na.action = \"na.impute\"). \n\nThe classes of your current object are: ",
  paste(unlist(object.classes),collapse=", ")))
    }else{
        predictRisk.args <- NULL
    }
    # }}}
    # {{{ add null model and check resampling ability
    if (!is.null(nullobject)) {
        mlevs <- 0:NF
        mlabels <- c(names(nullobject),names(object))
    } else{
        mlevs <- 1:NF
        mlabels <- names(object)
    }
    if (do.resample){
        nix <- lapply(1:length(object),function(f){
            fit <- object[[f]]
            if(is.null(fit$call))
                stop(paste("model",names(object)[f],"does not have a call argument."))
        })
    }
    # }}}
    # {{{ resolve keep statements
    if (!missing(keep) && is.character(keep)){
        if("residuals" %in% tolower(keep)) keep.residuals=TRUE else keep.residuals = FALSE
        if("vcov" %in% tolower(keep)) keep.vcov=TRUE else keep.vcov = FALSE
        if ("splitindex" %in% tolower(keep)) keep.splitindex=TRUE else keep.splitindex = FALSE
        if ("cv" %in% tolower(keep)) keep.cv=TRUE else keep.cv = FALSE
    }else{
        keep.residuals=FALSE
        keep.vcov=FALSE
        keep.cv=FALSE
        keep.splitindex=FALSE
    }
    # }}}
    # {{{ resolve se.fit and contrasts
    if (missing(se.fit)){
        if (is.logical(conf.int)[[1]] && conf.int[[1]]==FALSE
            || conf.int[[1]]<=0
            || conf.int[[1]]>1)
            se.fit <- 0L
        else
            se.fit <- 1L
    }else{
        stopifnot(is.logical(se.fit)||se.fit[[1]]==0||se.fit[[1]]==1)
    }
    if (split.method$internal.name=="noplan") multi.split.test <- FALSE
    if (se.fit[1]==1L) {
        if (is.numeric(conf.int) && conf.int[1]<1 && conf.int[1]>0)
            alpha <- 1-conf.int
        else
            alpha <- .05
    }else{
        alpha <- NA
    }
    # allow user to write contrasts=0 instead of contrasts=FALSE
    if (length(contrasts)==1 && length(contrasts[[1]])==1 && contrasts==0) contrasts = FALSE
    if ((NF+length(nullobject))<=1) dolist <- NULL
    else{
        if ((is.logical(contrasts) && contrasts[1]==FALSE)){
            dolist <- NULL
        } else{
            if (is.logical(contrasts) && contrasts[1]==TRUE){
                if (is.null(nullobject)){
                    dolist <- lapply(1:(NF-1),function(i){c(i,(i+1):NF)})
                } else{
                    dolist <- lapply(0:(NF-1),function(i){c(i,(i+1):NF)})
                }
            }else{
                dolist <- contrasts
                if (!is.list(contrasts)) contrasts <- list(contrasts)
                if (!(all(sapply(dolist,function(x){
                    if (!(length(x)>1)) stop("All elements of the list contrasts must contain at least two elements, i.e., the two models to be compared.")
                    all(x<=NF) && all(x>=0)
                }))))
                    stop(paste("Argument contrasts should be a list of positive integers possibly mixed with 0 that refer to elements of object.\nThe object has ",NF,"elements but "))
            }
        }
    }
    # }}}
    # {{{ Evaluation landmarks and horizons (times)
    if (response.type %in% c("survival","competing.risks")){
        ## in case of a tie, events are earlier than right censored
        eventTimes <- unique(data[,time])
        maxtime <- eventTimes[length(eventTimes)]
        if (debug==TRUE){
            cat("\nThe maxtime is set at:",maxtime,"\n")
        }
        include.times <- NULL
        if (missing(landmarks)){
            start <- 0
            if (missing(times)){
                if (use.event.times==TRUE)
                    times <- unique(c(start,eventTimes))
                else{
                    ## times <- seq(start,maxtime,(maxtime - start)/100)
                    times <- median(eventTimes)
                }
            } else{
                if (use.event.times==TRUE)
                    times <- sort(c(start,unique(times),eventTimes))
                else
                    ## times <- sort(unique(c(start,times)))
                    times <- sort(unique(times))
            }
            (if (any(times>maxtime))
                 message(paste0("Upper limit of followup is ",
                                maxtime,"\nResults at evaluation time(s) beyond this time point are not computed.")))
            ## need to save indices to modify matrix input
            include.times <- times<=maxtime
            times <- times[include.times]
            NT <-  length(times)
            if (NT==0)
                stop("No evaluation time before end of followup.")
        }
        else{
            stop("Landmark updating not yet implemented.")
        }
    } else{
        if (!missing(times) && (!is.null(times)) && (times[1]!=FALSE)) warning("Function 'Score': Response type is not time-to-event: argument 'times' will be ignored.",call.=FALSE)
        times <- NULL
        NT <- 1
    }
    # }}}
    # -----------------Dealing with censored data outside the loop -----------------------
    # {{{
    if (response.type %in% c("survival","competing.risks")){
        if (cens.type=="rightCensored"){
            if (se.fit[1]>0L && ("AUC" %in% metrics) && (conservative[1]==TRUE)) {
                ## FIXME: need conservative formula for AUC
                warning("Cannot do conservative==TRUE with AUC yet.")
            }
            if ((se.fit[[1]]>0L) && ("AUC" %in% metrics) && (cens.model[[1]]=="cox")){
                if (!(split.method$name %in% c("LeaveOneOutBoot","BootCv"))){
                    warning("Cannot (not yet) estimate standard errors for AUC with Cox IPCW.\nTherefore, force cens.model to be marginal.")
                    cens.model <- "KaplanMeier"
                }
            }
            Weights <- getCensoringWeights(formula=formula,
                                           data=data,
                                           times=times,
                                           cens.model=cens.model,
                                           response.type=response.type,
                                           ## FIXME: need conservative formula for AUC
                                           influence.curve=(se.fit[[1]]==TRUE && (conservative[[1]]==0L || ("AUC" %in% metrics))))
            ##
            ## if cens.model is marginal then IC is a matrix (ntimes,newdata)
            ## if cens.model is Cox then IC is an array (nlearn, ntimes, newdata)
            ## IC is an array with dimension (nlearn, times, newdata)
            ## IC_G(t,z;x_k)
        }else {
            if (cens.type=="uncensored"){
                cens.method <- "none"
                cens.model <- "none"
                if ("auc" %in% metrics){
                    if (se.fit==TRUE) {
                        warning("Standard error for AUC with uncensored time-to-event outcome not yet implemented.")
                        se.fit <- FALSE
                    }
                }
                Weights <- NULL
            }
            else{
                stop("Cannot handle this type of censoring.")
            }
        }
        if ("pseudo" %in% cens.method){
            if (cens.type[1]=="rightCensored"
                && (response.type %in% c("survival","competing.risks"))){
                margForm <- update(formula,paste(".~1"))
                margFit <- prodlim::prodlim(margForm,data=data)
                ## position.cause is the result of match(cause, states)
                jack <- data.table(ID=data[["ID"]],
                                   times=rep(times,rep(N,NT)),
                                   pseudovalue=c(prodlim::jackknife(margFit,cause=position.cause,times=times)))
                if (response.type=="survival") jack[,pseudovalue:=1-pseudovalue]
            }
        }
    } else{
        Weights <- NULL
    }
    # }}}
    # -----------------performance program----------------------
    # {{{ define getPerformanceData, setup data in long-format, response, pred, weights, model, times, loop
    # {{{ header
    getPerformanceData <- function(testdata,
                                   testweights,
                                   traindata=NULL,
                                   trainseed=NULL){
        ID=NULL
        # inherit everything else from parent frame: object, nullobject, NF, NT, times, cause, response.type, etc.
        Brier=IPA=IBS=NULL
        looping <- !is.null(traindata)
        ## if (!looping) b=0
        N <- as.numeric(NROW(testdata))
        # split data vertically into response and predictors X
        response <- testdata[,1:response.dim,with=FALSE]
        response[,ID:=testdata[["ID"]]]
        setkey(response,ID)
        X <- testdata[,-c(1:response.dim),with=FALSE]
        ## restore sanity
        setnames(X,sub("^protectedName.","",names(X)))
        if (debug) if (looping) message(paste0("Loop round: ",b))
        if (debug) message("extracted test set and prepared output object")
        # }}}
        # {{{ collect pred as long format data.table
        args <- switch(response.type,"binary"={list(newdata=X)},
                       "survival"={list(newdata=X,times=times)},
                       "competing.risks"={list(newdata=X,times=times,cause=cause)},
                       stop("Unknown response.type."))
        ## remove our cbinded response (see above) from traindata to avoid clash when model uses Hist(time,status)
        ## where status has 0,1,2 but now event history response has status=0,1
        ## the original response is still there
        if(!is.null(traindata)){
            trainX <- traindata[,-c(1:response.dim),with=FALSE]
            ## restore sanity
            setnames(trainX,sub("^protectedName.","",names(trainX)))
            ## trainX <- copy(traindata)
            trainX[,ID:=NULL]
        }
        pred <- data.table::rbindlist(lapply(mlevs, function(f){
            ## add object specific arguments to predictRisk methods
            if (f[1]>0 && (length(extra.args <- unlist(lapply(object.classes[[f]],function(cc){predictRisk.args[[cc]]})))>0)){
                args <- c(args,extra.args)
            }
            ## predictions given as numeric values
            if (f[1]!=0 && any(c("integer","factor","numeric","matrix") %in% object.classes[[f]])){
                ## sort predictions by ID
                if (!is.null(dim(object[[f]]))) {## input matrix
                    if (response.type=="binary"){
                        p <- do.call("predictRisk", c(list(object=c(object[[f]])),args))[neworder]
                    }else{
                        if(!is.null(include.times)){ ## remove columns at times beyond max time
                            p <- c(do.call("predictRisk",c(list(object=object[[f]][,include.times,drop=FALSE]),args))[neworder,])
                        } else{
                            p <- c(do.call("predictRisk",c(list(object=object[[f]]),args))[neworder,])
                        }
                    }
                }
                else{ ## either binary or only one time point
                    p <- do.call("predictRisk", c(list(object=object[[f]]),args))[neworder]
                }
            }else{
                # predictions given as model which needs training in crossvalidation loops
                if (looping){
                    set.seed(trainseed)
                    if (f==0) model.f=nullobject[[1]] else model.f=object[[f]]
                    model.f$call$data <- trainX
                    # browser(skipCalls = 1)
                    trained.model <- try(eval(model.f$call),silent=TRUE)
                    if ("try-error" %in% class(trained.model)){
                        message(paste0("Failed to fit model ",f,ifelse(looping,paste0(" in cross-validation step ",b)),":"))
                        try(eval(model.f$call),silent=FALSE)
                        stop()
                    }
                } else{
                    if (f==0)
                        trained.model <- nullobject[[1]]
                    else
                        trained.model <- object[[f]]
                }
                p <- c(do.call("predictRisk", c(list(object=trained.model),args)))
                if (f[1]==0 && (response.type[1]!="binary")) {## glm predicts the same value for all subjects
                    p <- rep(p,rep(N,NT))
                }
            }
            if (response.type%in%c("survival","competing.risks")){
                out <- data.table(ID=testdata[["ID"]],model=f,risk=p,times=rep(times,rep(N,NT)))
                setkey(out,model,times,ID)
                out
            } else {
                out <- data.table(ID=testdata[["ID"]],model=f,risk=p)
                setkey(out,model,ID)
                out
            }
        }))
        if (debug) message("trained the model(s) and extracted the predictions")
        # }}}
        # {{{ merge with Weights (IPCW inner loop)
        if (response.type %in% c("survival","competing.risks")){
            if (cens.type=="rightCensored"){
                Weights <- testweights
                ## add subject specific weights
                response[,WTi:=Weights$IPCW.subject.times]
            } else {
                if (cens.type=="uncensored"){
                    Weights <- list(IPCW.times=rep(1,NT),IPCW.subject.times=matrix(1,ncol=NT,nrow=N))
                    Weights$method <- "marginal"
                    response[,WTi:=1]
                } else{
                    stop("Cannot handle this type of censoring.")
                }
            }
            ## add time point specific weights
            if (Weights$method=="marginal"){
                Wt <- data.table(times=times,Wt=Weights$IPCW.times)
                ## OBS: many digits in times may cause merge problems
                pred <- merge(pred,Wt,by=c("times"))
            }else{
                Wt <- data.table(times=rep(times,rep(N,NT)),
                                 Wt=c(Weights$IPCW.times),
                                 ID=testdata$ID)
                pred <- merge(pred,Wt,by=c("ID","times"))
            }
            if (debug) message("merged the weights with input for performance metrics")
        } else {
            ## if (response.type=="binary")
            Weights <- NULL
        }
        if (debug) message("added weights to predictions")
        # }}}
        # {{{ merge with response
        DT=merge(response,pred,by="ID")
        DT
    }
    # }}}
    # }}}
    # {{{ define function to test performance
    computePerformance <- function(DT,
                                   N,
                                   se.fit,
                                   conservative,
                                   cens.model,
                                   multi.split.test,
                                   keep.residuals,
                                   keep.vcov){
        IPA=IBS=Brier=NULL
        # inherit everything else from parent frame: summary, metrics, plots, alpha, probs, dolist, et
        out <- vector(mode="list",
                      length=length(c(summary,metrics,plots)))
        names(out) <- c(summary,metrics,plots)
        input <- list(DT=DT,
                      N=N,
                      NT=NT,
                      NF=NF,
                      alpha=alpha,
                      se.fit=se.fit,
                      conservative=conservative,
                      cens.model=cens.model,
                      multi.split.test=multi.split.test,
                      keep.residuals=keep.residuals,
                      keep.vcov=keep.vcov,
                      ## DT.residuals=DT.residuals,
                      dolist=dolist,Q=probs,ROC=FALSE,MC=Weights$IC)
        if (response.type=="competing.risks") {
            input <- c(input,list(cause=cause,states=states))
        }
        # {{{ collect data for summary statistics
        for (s in summary){
            if (s=="risks") {
                out[[s]] <- list(score=copy(input$DT)[,model:=factor(model,levels=mlevs,mlabels)],
                                 contrasts=NULL)
            } else{
                out[[s]] <- do.call(paste(s,response.type,sep="."),input)
                if (NROW(out[[s]]$score)>0){
                    out[[s]]$score[,model:=factor(model,levels=mlevs,mlabels)]
                }
                if (NROW(out[[s]]$contrasts)>0){
                    out[[s]]$contrasts[,model:=factor(model,levels=mlevs,mlabels)]
                    out[[s]]$contrasts[,reference:=factor(reference,levels=mlevs,mlabels)]
                }
            }
        }
        # }}}
        # {{{ collect data for calibration plots
        if ("Calibration" %in% plots){
            if (response.type[[1]]=="binary" || cens.type[[1]]=="uncensored")
                out[["Calibration"]]$plotframe <- DT[model!=0]
            else{
                out[["Calibration"]]$plotframe <- merge(jack,DT[model!=0],by=c("ID","times"))
            }
            out[["Calibration"]]$plotframe[,model:=factor(model,levels=mlevs,mlabels)]
        }
        # }}}
        ## make sure that Brier score comes first, so that we can remove the null.model afterwards
        for (m in sort(metrics,decreasing=TRUE)){
            if (m=="AUC" && ("ROC" %in% plots)){
                input <- replace(input, "ROC",TRUE)
                ## call AUC method
                out[[m]] <- do.call(paste(m,response.type,sep="."),input)
                out[["ROC"]]$plotframe <- out[[m]]$ROC
                out[["ROC"]]$plotframe[,model:=factor(model,levels=mlevs,mlabels)]
                out[[m]]$ROC <- NULL
            }else{
                input <- replace(input, "ROC",FALSE)
                ## call Brier or AUC method
                out[[m]] <- do.call(paste(m,response.type,sep="."),input)
            }
            if (!is.null(out[[m]]$score)){
                out[[m]]$score[,model:=factor(model,levels=mlevs,mlabels)]
            }
            ## set model and reference in model comparison results
            if (!is.null(out[[m]]$contrasts)>0){
                out[[m]]$contrasts[,model:=factor(model,levels=mlevs,mlabels)]
                out[[m]]$contrasts[,reference:=factor(reference,levels=mlevs,mlabels)]
            }
        }
        ## summary should be after metrics because IBS and IPA/R^2 depends on Brier score
        if (ibs){
            Dint <- function(x,y,range,na.omit=FALSE){
                if (is.null(range)) range=c(x[1],x[length(x)])
                ##   integrate a step function f with
                ##   values y=f(x) between range[1] and range[2]
                start <- max(range[1],min(x))
                Stop <- min(range[2],max(x))
                if ((Stop-start)<=0)
                    return(0)
                else{
                    Y=y[x>=start & x<Stop]
                    X=x[x>=start & x<Stop]
                    if (na.omit){
                        X=X[!is.na(Y)]
                        Y=Y[!is.na(Y)]
                    } else if (any(is.na(Y))|| any(is.na(X))){
                        return(NA)
                    }
                    return(1/(Stop-start) * sum(Y*diff(c(X,Stop))))
                }
            }
            if (response.type!="binary"){
                out[["Brier"]][["score"]][,IBS:=sapply(times,function(t){
                    Dint(x=c(0,times),y=c(0,Brier),range=c(0,t))
                }),by=c("model")]
            }
        }
        if (ipa == TRUE){
            if (response.type=="binary")
                out[["Brier"]][["score"]][,IPA:=1-Brier/Brier[model=="Null model"]]
            else
                out[["Brier"]][["score"]][,IPA:=1-Brier/Brier[model=="Null model"],by=times]
        }
        out
    }
    # }}}
    # -----------------apparent nosplit performance---------------------
    # {{{
    missing.predictions <- "Don't know yet"
    if (split.method$internal.name %in% c("noplan",".632+")){
        DT <- getPerformanceData(testdata=data,
                                 testweights=Weights,
                                 traindata=NULL,
                                 trainseed=NULL)
        if (any(is.na(DT[["risk"]]))){
            missing.predictions <- DT[,list("Missing.values"=sum(is.na(risk))),by=byvars]
            warning("Missing values in the predicted risk. See `missing.predictions' in output list.")
        }else{
            missing.predictions <- "None"
        }
        noSplit <- computePerformance(DT,
                                      N=N,
                                      se.fit=se.fit,
                                      conservative=conservative,
                                      cens.model,
                                      multi.split.test=multi.split.test,
                                      keep.residuals=keep.residuals,
                                      keep.vcov=keep.vcov)
        if (debug) message("computed apparent performance")
    }
    # }}}
    # -----------------crossvalidation performance---------------------
    # {{{
    # {{{ bootstrap re-fitting
    if (split.method$internal.name%in%c("BootCv","LeaveOneOutBoot","crossval")){
        if (missing(trainseeds)||is.null(trainseeds)){
            if (!missing(seed)) set.seed(seed)
            if (split.method$internal.name == "crossval"){
                trainseeds <- sample(1:1000000,size=B*split.method$k,replace=FALSE)
            }else{
                trainseeds <- sample(1:1000000,size=B,replace=FALSE)
            }
        }
        if (parallel=="snow") exports <- c("data","split.method","Weights","N","trainseeds") else exports <- NULL
        if (!is.null(progress.bar)){
            message("Running crossvalidation ...")
            if (!(progress.bar %in% c(1,2,3))) progress.bar <- 3
            if (B==1 && split.method$internal.name == "crossval")
                pb <- txtProgressBar(max = split.method$k, style = progress.bar,width=20)
            else
                pb <- txtProgressBar(max = B, style = progress.bar,width=20)
        }
        `%dopar%` <- foreach::`%dopar%`
        ## k-fold-CV
        if (split.method$internal.name == "crossval"){
            DT.B <- rbindlist(foreach::foreach (b=1:B,.export=exports,.packages="data.table",.errorhandling="pass") %dopar%{ ## repetitions of k-fold to avoid Monte-Carlo error
                index.b <- split.method$index[,b] ## contains a sample of the numbers 1:k with replacement
                if((B>1) && !is.null(progress.bar)){setTxtProgressBar(pb, b)}
                DT.b <- rbindlist(lapply(1:split.method$k,function(fold){
                    traindata=data[index.b!=fold]
                    testids <- index.b==fold # (1:N)[index.b!=fold]
                    if((B==1) && !is.null(progress.bar)){setTxtProgressBar(pb, fold)}
                    ## NOTE: subset.data.table preserves order ## So we need to use subset??
                    testdata <- subset(data,testids)
                    if (cens.type=="rightCensored"){
                        testweights <- Weights
                        # Need to check what's expected that testids is here and below:
                        testweights$IPCW.subject.times <- subset(testweights$IPCW.subject.times,testids)
                        if (Weights$dim>0){
                            testweights$IPCW.times <- subset(testweights$IPCW.times,testids)
                        }
                    } else {
                        testweights <- NULL
                    }
                    ## predicted risks of model trained without this fold
                    ## evaluated and added to this fold
                    DT.fold <- getPerformanceData(testdata=testdata,
                                                  testweights=testweights,
                                                  traindata=traindata,
                                                  trainseed=trainseeds[[(b-1)*split.method$k+fold]])
                    return(DT.fold)
                }))
                DT.b[,b:=b]
                DT.b
            })
        }else{# either LeaveOneOutBoot or BootCv
            DT.B <- rbindlist(foreach::foreach (b=1:B,.export=exports,.packages="data.table",.errorhandling="pass") %dopar%{
                if(!is.null(progress.bar)){setTxtProgressBar(pb, b)}
                ## DT.B <- rbindlist(lapply(1:B,function(b){
                traindata=data[split.method$index[,b]]
                ## setkey(traindata,ID)
                testids <- (match(1:N,unique(split.method$index[,b]),nomatch=0)==0)
                ## NOTE: subset.data.table preserves order
                testdata <- subset(data,testids)
                if (cens.type=="rightCensored"){
                    testweights <- Weights
                    testweights$IPCW.subject.times <- subset(testweights$IPCW.subject.times,testids)
                    if (Weights$dim>0){
                        testweights$IPCW.times <- subset(testweights$IPCW.times,testids)
                    }
                } else {
                    testweights <- NULL
                }
                DT.b <- getPerformanceData(testdata=testdata,
                                           testweights=testweights,
                                           traindata=traindata,
                                           trainseed=trainseeds[[b]])
                DT.b[,b:=b]
                DT.b
            })
        }
        if (any(is.na(DT.B[["risk"]]))){
            missing.predictions <- DT.B[,list("Missing.values"=sum(is.na(risk))),by=byvars]
            warning("Missing values in the predicted risk. See `missing.predictions' in output list.")
        }
        ## FIXME: subset influence curves
        ## case se.fit=1 we need only p-values for multi-split tests
        ## se.fit.cv <- se.fit*2
        ## cb <- computePerformance(DT.b,
        ## se.fit=se.fit.cv,
        ## multi.split.test=multi.split.test,
        ## keep.residuals=keep.residuals)
        ## cb
        if (debug) message("setup data for cross-validation performance")
        Response <- data[,c(1:response.dim),with=FALSE]
        Response[,ID:=data[["ID"]]]
        Response.names <- names(Response)
        Response.names <- Response.names[Response.names!="ID"]
        setkey(Response,ID)
        ## ## Show format for the data in DT.B
        ## cat(paste("\nDT.B for method:", split.method$name, "\n"))
        ## print(DT.B)
        # }}}
        # {{{ Leave-one-out bootstrap
        ## start clause split.method$name=="LeaveOneOutBoot"
        if (split.method$name=="LeaveOneOutBoot" | split.method$internal.name =="crossval"){  ## Testing if the crossval works in this loop
            ## if (split.method$name=="LeaveOneOutBoot" ){
            crossvalPerf <- lapply(metrics,function(m){
                # {{{ AUC LOOB
                if (m=="AUC"){
                    if (response.type=="binary")
                        auc.loob <- data.table(model=mlevs)
                    else
                        auc.loob <- data.table(expand.grid(times=times,model=mlevs))
                    auc.loob[,AUC:=as.numeric(NA)]
                    ## for each pair of individuals sum the concordance of the bootstraps where *both* individuals are out-of-bag
                    ## divide by number of times the pair is out-of-bag later
                    if (se.fit==TRUE){
                        aucDT <- NULL
                    }
                    if (response.type=="binary"){
                        NT <- 1
                    }
                    for (s in 1:NT){
                        t <- times[s]
                        if (response.type=="binary"){
                            ## the following indices have to be logical!!!
                            cases.index <- Response[,ReSpOnSe==1]
                            controls.index <- !cases.index
                            cc.status <- factor(cases.index,levels=c(TRUE,FALSE),labels=c("case","control"))
                        }else{
                            if (response.type=="survival"){
                                ## event of interest before times
                                ## the following indices have to be logical!!!
                                cases.index <- Response[,time<=t & status==1]
                                controls.index <- Response[,time>t]
                                cc.status <- factor(rep("censored",N),levels=c("censored","case","control"))
                                cc.status[cases.index] <- "case"
                                cc.status[controls.index] <- "control"
                            }
                            else{ ## competing risks
                                ## event of interest before times
                                ## the following indices have to be logical!!!
                                cases.index <- Response[,time<=t & status==1]
                                controls.index <- Response[,time>t | status==2]
                                cc.status <- factor(rep("censored",N),levels=c("censored","case","control"))
                                cc.status[cases.index] <- "case"
                                cc.status[controls.index] <- "control"
                            }
                        }
                        if (cens.type=="rightCensored"){
                            ## IPCW
                            weights.cases <- cases.index/Weights$IPCW.subject.times
                            if (Weights$method=="marginal"){
                                weights.controls <- controls.index/Weights$IPCW.times[s]
                            }else{
                                weights.controls <- controls.index/Weights$IPCW.times[,s]
                            }
                            weightMatrix <- outer(weights.cases[cases.index], weights.controls[controls.index], "*")
                        }else{ ## uncensored
                            weights.cases <- cases.index/1
                            weights.controls <- controls.index/1
                            weightMatrix <- outer(weights.cases[cases.index], weights.controls[controls.index], "*")
                        }
                        Phi <- (1/N^2)*sum(weights.cases[cases.index])*sum(weights.controls[controls.index])
                        which.cases <- (1:N)[cases.index]
                        which.controls <- (1:N)[controls.index]
                        for (mod in mlevs){
                            Ib <- matrix(0, sum(cases.index), sum(controls.index))
                            auc <- matrix(0, sum(cases.index), sum(controls.index))
                            for (u in 1:B){## cannot use b as running index because b==b does not work in data.table
                                ## test <- DT.B[model==mod&times==t&b==u]
                                oob <- match(1:N,unique(split.method$index[,u]),nomatch=0)==0
                                ## to use the cpp function AUCijFun we
                                ## need a vector of length equal to the number of cases (which.cases) for the current time point
                                ## which has arbitrary values in places where subjects are inbag and the predicted risks
                                ## for those out-of-bag. need another vector for controls.
                                riskset <- data.table::data.table(ID=1:N,casecontrol=cc.status,oob=oob)
                                data.table::setkey(riskset,ID)
                                if (is.null(t)){
                                    oob.risk <- DT.B[model==mod&b==u,data.table::data.table(ID,risk)]
                                }else{
                                    oob.risk <- DT.B[model==mod&times==t&b==u,data.table::data.table(ID,risk)]
                                }
                                data.table::setkey(oob.risk,ID)
                                riskset <- oob.risk[riskset]
                                riskset[is.na(risk),risk:=-9]
                                Ib.ij <- outer((cases.index*oob)[which.cases],(controls.index*oob)[which.controls],"*")
                                auc.ij <- AUCijFun(riskset[casecontrol=="case",risk],riskset[casecontrol=="control",risk])*Ib.ij
                                ## Ib.ij is 1 when the pair out of bag
                                ## print(head(oob))
                                ## print(auc.ij[1:5,1:5])
                                auc <- auc+auc.ij
                                Ib <- Ib + Ib.ij
                            }
                            auc <- (auc*weightMatrix)/Ib
                            # FIXME: why are there NA's?
                            auc[is.na(auc)] <- 0
                            ## Leave-one-pair-out bootstrap estimate of AUC
                            aucLPO <- (1/N^2)*sum(colSums(auc))*(1/Phi)
                            if (is.null(t)){
                                auc.loob[model==mod,AUC:=aucLPO]
                            }else{
                                auc.loob[times==t&model==mod,AUC:=aucLPO]
                            }
                            if (se.fit==1L){
                                ## ## First part of influence function
                                ic0Case <- rowSums(auc)
                                ic0Control <- colSums(auc)
                                ic0 <- (1/(Phi*N))*c(ic0Case, ic0Control)-2*aucLPO
                                id.cases <- data[["ID"]][cc.status=="case"]
                                id.controls <- data[["ID"]][cc.status=="control"]
                                id.censored <- data[["ID"]][cc.status=="censored"]
                                if (is.null(t)){
                                    this.aucDT <- data.table(model=mod,ID = c(id.cases,id.controls), IF.AUC0 = ic0)
                                }else{
                                    this.aucDT <- data.table(model=mod,times=t,ID = c(id.cases,id.controls,id.censored), IF.AUC0 = c(ic0, rep(-2*aucLPO,length(id.censored))))
                                }
                                aucDT <- rbindlist(list(aucDT,this.aucDT),use.names=TRUE,fill=TRUE)
                                if (response.type[[1]]=="binary" || cens.type[[1]]=="uncensored"){
                                    icPhi <- (aucLPO/Phi)*((weights.cases-(1/N)*sum(weights.cases))*(1/N)*sum(weights.controls)+(weights.controls-(1/N)*sum(weights.controls))*(1/N)*sum(weights.controls))-2*aucLPO
                                    if (is.null(t)){
                                        data.table::setkey(aucDT,model,ID)
                                        aucDT[model==mod,IF.AUC:=IF.AUC0-icPhi]
                                        auc.loob[model==mod,se:=sd(aucDT[["IF.AUC"]])/sqrt(N)]
                                    }else{
                                        data.table::setkey(aucDT,model,times,ID)
                                        aucDT[times==t&model==mod,IF.AUC:=IF.AUC0-icPhi]
                                        auc.loob[times==t&model==mod,se:=sd(aucDT[["IF.AUC"]])/sqrt(N)]
                                    }
                                }else{
                                    ic.weights <- matrix(0,N,N)
                                    if (cens.type[1]=="rightCensored" && (conservative[1]==FALSE)) {
                                        ## ## Influence function for G - i.e. censoring survival distribution
                                        if  (cens.model=="cox"){
                                            k=0 ## counts subject-times with event before t
                                            for (i in 1:N){
                                                if (i %in% id.cases){
                                                    k=k+1
                                                    ic.weights[i,] <- Weights$IC$IC.subject[i,k,]/(Weights$IPCW.subject.times[i])
                                                }else{
                                                    if (i %in% id.controls){ ## min(T,C)>t
                                                        ic.weights[i,] <- Weights$IC$IC.times[i,s,]/(Weights$IPCW.times[i,s])
                                                    }
                                                }
                                            }
                                        }else{
                                            k=0 ## counts subject-times with event before t
                                            for (i in 1:N){
                                                if (i %in% id.cases){
                                                    ## FIXME: need to check IC
                                                    pos.i <- sindex(jump.times=unique(data[["time"]]),eval.times=data[["time"]][i])
                                                    ic.weights[i,] <- Weights$IC[pos.i,]/(Weights$IPCW.subject.times[i])
                                                }else{
                                                    if (i %in% id.controls){ ## min(T,C)>t
                                                        ic.weights[i,] <- Weights$IC[s,]/(Weights$IPCW.times[s])
                                                    }
                                                }
                                            }
                                        }
                                        ## ## Part of influence function related to Weights
                                        ic.weightsCase <- as.numeric(rowSumsCrossprod(as.matrix(rowSums(auc)), ic.weights[which.cases,], 0))
                                        ic.weightsControl <- as.numeric(rowSumsCrossprod(as.matrix(colSums(auc)), ic.weights[which.controls,], 0))
                                        ic.weightsCC <- (1/(Phi*N^2))*(ic.weightsCase+ic.weightsControl)
                                    }
                                    ## ## Part of influence function related to Phi
                                    ## icPhiCase <- colMeans(ic.weights[which.cases,])
                                    icPhiCase <- as.numeric(rowSumsCrossprod(as.matrix(weights.cases[which.cases]),ic.weights[which.cases,],0))
                                    icPhiControl <- as.numeric(rowSumsCrossprod(as.matrix(weights.controls[which.controls]),ic.weights[which.controls,],0))
                                    icPhi <- (aucLPO/Phi)*((weights.cases-(1/N)*icPhiCase)*(1/N)*sum(weights.controls)+(weights.controls-(1/N)*icPhiControl)*(1/N)*sum(weights.cases)) - 2*aucLPO
                                    ## ## Combine all parts of influence function
                                    ## ic1 <- data.table(ID=data[["ID"]], "ic.weightsCC" = ic.weightsCC, "icPhi" = icPhi)
                                    data.table::setkey(aucDT,model,times,ID)
                                    if(conservative==TRUE){
                                        aucDT[model==mod&times==t, IF.AUC:=IF.AUC0-icPhi]
                                    }else{
                                        aucDT[model==mod&times==t, IF.AUC:=IF.AUC0-ic.weightsCC-icPhi]
                                        aucDT[model==mod&times==t, IF.AUC.conservative:=IF.AUC0-icPhi]
                                    }
                                    aucDT[,IF.AUC0:=NULL]
                                    auc.loob[model==mod&times==t,se:= sd(aucDT[model==mod&times==t,IF.AUC])/sqrt(N)]
                                    auc.loob[model==mod&times==t,se.conservative:=sd(aucDT[model==mod&times==t,IF.AUC.conservative])/sqrt(N)]
                                    ## testSE <- sqrt(sum(ic[["ic"]]^2))/N
                                }
                            }
                        }
                    }
                    if (se.fit==1L){
                        auc.loob[,lower:=pmax(0,AUC-qnorm(1-alpha/2)*se)]
                        auc.loob[,upper:=pmin(1,AUC + qnorm(1-alpha/2)*se)]
                    }
                    output <- list(score=auc.loob)
                    if (length(dolist)>0L){
                        if (se.fit==FALSE){
                            if (match("times",byvars,nomatch=0)){
                                contrasts.AUC <- auc.loob[,getComparisons(data.table(x=AUC,model=model),
                                                                          NF=NF,
                                                                          N=N,
                                                                          alpha=alpha,
                                                                          dolist=dolist,
                                                                          se.fit=FALSE),by=list(times)]
                            } else{
                                contrasts.AUC <- auc.loob[,getComparisons(data.table(x=AUC,model=model),
                                                                          NF=NF,
                                                                          N=N,
                                                                          alpha=alpha,
                                                                          dolist=dolist,
                                                                          se.fit=FALSE)]
                            }
                        }else{
                            aucDT <- merge(aucDT,auc.loob,by=byvars)
                            if (match("times",byvars,nomatch=0)){
                                contrasts.AUC <- aucDT[,getComparisons(data.table(x=AUC,IF=IF.AUC,model=model),
                                                                       NF=NF,
                                                                       N=N,
                                                                       alpha=alpha,
                                                                       dolist=dolist,
                                                                       se.fit=TRUE),by=list(times)]
                            } else{
                                contrasts.AUC <- aucDT[,getComparisons(data.table(x=AUC,IF=IF.AUC,model=model),
                                                                       NF=NF,
                                                                       N=N,
                                                                       alpha=alpha,
                                                                       dolist=dolist,
                                                                       se.fit=TRUE)]
                            }
                        }
                        setnames(contrasts.AUC,"delta","delta.AUC")
                        output <- c(output,list(contrasts=contrasts.AUC))
                    }
                    if (!is.null(output$score)){
                        output$score[,model:=factor(model,levels=mlevs,mlabels)]
                    }
                    ## set model and reference in model comparison results
                    if (!is.null(output$contrasts)>0){
                        output$contrasts[,model:=factor(model,levels=mlevs,mlabels)]
                        output$contrasts[,reference:=factor(reference,levels=mlevs,mlabels)]
                    }
                    return(output)
                }
                # }}}
                # {{{ Brier LOOB
                if (m=="Brier"){
                    ## sum across bootstrap samples where subject i is out of bag
                    if (cens.type=="rightCensored"){
                        if (response.type=="survival"){
                            ## event of interest before times
                            DT.B[time<=times & status==1,residuals:=(1-risk)^2/WTi]
                        }
                        else{ ## competing risks
                            ## event of interest before times
                            DT.B[time<=times & status==1 & event==cause,residuals:=(1-risk)^2/WTi]
                            ## competing event before times
                            DT.B[time<=times & status==1 &event!=cause,residuals:=(0-risk)^2/WTi]
                        }
                        ## right censored before times
                        DT.B[time<=times & status==0,residuals:=0]
                        ## no event at times
                        DT.B[time>times,residuals:=(risk)^2/Wt]
                    }else{
                        DT.B[,residuals:=(ReSpOnSe-risk)^2]
                    }
                    ## for each individual sum the residuals of the bootstraps where this individual is out-of-bag
                    ## divide by number of times out-off-bag later
                    ## DT.B <- DT.B[,data.table::data.table(residuals=sum(residuals)),by=c(byvars,"ID")]
                    DT.B <- DT.B[,data.table::data.table(risk=mean(risk),residuals=sum(residuals)),by=c(byvars,"ID")]
                    ## get denominator
                    if (split.method$name=="LeaveOneOutBoot"){
                        Ib <- split.method$B-tabulate(unlist(apply(split.method$index,2,unique)))
                        ## REMOVE ME
                        ## Ib <- Ib[order(data$ID)]
                        if (any(Ib==0)) {
                            warning("Some subjects are never out of bag.\nResults are unreliable until You increase the number of bootstrap replications (argument 'B').")
                            Ib.include <- Ib!=0
                            Ib <- Ib[Ib.include]
                            ## don't subset residuals, they are only
                            ## available for those with Ib.include==1
                            ## DT.B <- DT.B[Ib.include]
                        } else Ib.include <- NULL
                    }else{## cv-k crossvalidation
                        Ib <- rep(split.method$B,N)
                    }
                    ## Ib is the count of how many times subject i is out of bag
                    ## the order of Ib matches the order of ID
                    ## within groups defined by times and model (byvars)
                    data.table::setkeyv(DT.B,c(byvars,"ID"))
                    DT.B[,residuals:=residuals/Ib]
                    ## leave-one-out bootstrap estimate
                    DT.B[,Brier:=mean(residuals),by=byvars]
                    ## standard error via influence function
                    if (se.fit==1L){
                        ## influence function when censoring model is known or data are uncensored
                        DT.B[,IC0:=residuals-Brier]
                        ## se.brier <- DT.B[,list(se=sd(IC0, na.rm=TRUE)/sqrt(N)),by=byvars]
                        ## DT.B[,Brier:=NULL]
                        if (cens.type[1]=="rightCensored" && conservative[1]!=TRUE){
                            ## this is a new DT.B
                            DT.B <- cbind(data[,1:response.dim,with=FALSE],DT.B)
                            DT.B[,nth.times:=as.numeric(factor(times))]
                            WW <- data.table(ID=1:N,WTi=Weights$IPCW.subject.times,key="ID")
                            DT.B <- merge(WW,DT.B,by=ID)
                            ## DT.B[,WTi:=rep(Weights$IPCW.subject.times,NF+length(nullobject))]
                            if (Weights$method=="marginal"){
                                Wt <- data.table(times=times,Wt=Weights$IPCW.times)
                                ## OBS: many digits in times may cause merge problems
                                DT.B <- merge(DT.B,Wt,by=c("times"))
                            }else{
                                Wt <- data.table(times=rep(times,rep(N,NT)),
                                                 Wt=c(Weights$IPCW.times),
                                                 ID=data$ID)
                                DT.B <- merge(DT.B,Wt,by=c("ID","times"))
                            }
                            if (cens.type=="uncensored"){
                                DT.B[,IF.Brier:= residuals]
                                score.loob <- DT.B[,data.table(Brier=sum(residuals)/N,
                                                               se=sd(residuals)/sqrt(N),
                                                               se.conservative=sd(residuals)/sqrt(N)),by=byvars]
                            }else{
                                DT.B[,IF.Brier:=getInfluenceCurve.Brier(t=times[1],
                                                                        time=time,
                                                                        IC0,
                                                                        residuals=residuals,
                                                                        WTi=WTi,
                                                                        Wt=Wt,
                                                                        IC.G=Weights$IC,
                                                                        cens.model=cens.model,
                                                                        nth.times=nth.times[1]),by=byvars]
                                score.loob <- DT.B[,data.table(Brier=sum(residuals)/N,
                                                               se=sd(IF.Brier)/sqrt(N),
                                                               se.conservative=sd(IC0)/sqrt(N)),by=byvars]
                            }
                        }else{
                            ## either conservative == TRUE or binary or uncensored
                            score.loob <- DT.B[,data.table(Brier=sum(residuals)/N,se=sd(IC0)/sqrt(N)),
                                               by=byvars]
                            setnames(DT.B,"IC0","IF.Brier")
                        }
                        score.loob[,lower:=pmax(0,Brier-qnorm(1-alpha/2)*se)]
                        score.loob[,upper:=pmin(1,Brier + qnorm(1-alpha/2)*se)]
                    } else{
                        score.loob <- DT.B[,list(Brier=sum(residuals)/N), by=byvars]
                    }
                    if (ipa==TRUE){
                        if (response.type=="binary")
                            score.loob[,IPA:=1-Brier/Brier[model==0]]
                        else
                            score.loob[,IPA:=1-Brier/Brier[model==0],by=times]
                    }
                    data.table::setkeyv(score.loob,byvars)
                    ## data.table::setkey(DT.B,model,times)
                    ## DT.B <- DT.B[score.loob]
                    if (length(dolist)>0L){
                        if (se.fit==FALSE){
                            if (match("times",byvars,nomatch=0))
                                contrasts.Brier <- DT.B[,getComparisons(data.table(x=Brier,model=model),
                                                                        NF=NF,
                                                                        N=N,
                                                                        alpha=alpha,
                                                                        dolist=dolist,
                                                                        se.fit=FALSE),by=list(times)]
                            else
                                contrasts.Brier <- DT.B[,getComparisons(data.table(x=Brier,model=model),
                                                                        NF=NF,
                                                                        N=N,
                                                                        alpha=alpha,
                                                                        dolist=dolist,
                                                                        se.fit=FALSE)]
                        } else{
                            if (match("times",byvars,nomatch=0))
                                contrasts.Brier <- DT.B[,getComparisons(data.table(x=Brier,IF=IF.Brier,model=model),
                                                                        NF=NF,
                                                                        N=N,
                                                                        alpha=alpha,
                                                                        dolist=dolist,
                                                                        se.fit=TRUE),by=list(times)]
                            else
                                contrasts.Brier <- DT.B[,getComparisons(data.table(x=Brier,IF=IF.Brier,model=model),
                                                                        NF=NF,
                                                                        N=N,
                                                                        alpha=alpha,
                                                                        dolist=dolist,
                                                                        se.fit=TRUE)]
                        }
                        setnames(contrasts.Brier,"delta","delta.Brier")
                        output <- list(score=score.loob,contrasts=contrasts.Brier)
                    } else{
                        output <- list(score=score.loob)
                    }
                    if (keep.residuals) {
                        DT.B[,model:=factor(model,levels=mlevs,mlabels)]
                        if (all(c("Wt","WTi")%in%names(DT.B))){
                            DT.B[,IPCW:=1/WTi]
                            DT.B[time>=times,IPCW:=1/Wt]
                            DT.B[time<times & status==0,IPCW:=0]
                            output <- c(output,list(residuals=DT.B[,c("ID",names(response),"model","times","risk","residuals","IPCW"),with=FALSE]))
                        }else{
                            output <- c(output,list(residuals=DT.B[,c("ID",names(response),"model","times","risk","residuals"),with=FALSE]))
                        }
                    }
                    if (!is.null(output$score)){
                        output$score[,model:=factor(model,levels=mlevs,mlabels)]
                        if (response.type%in%c("survival","competing.risks"))
                            setkey(output$score,model,times)
                        else
                            setkey(output$score,model)
                    }
                    ## set model and reference in model comparison results
                    if (!is.null(output$contrasts)>0){
                        output$contrasts[,model:=factor(model,levels=mlevs,mlabels)]
                        output$contrasts[,reference:=factor(reference,levels=mlevs,mlabels)]
                        if (response.type%in%c("survival","competing.risks"))
                            setkey(output$score,model,times)
                        else
                            setkey(output$score,model)
                    }
                    return(output)
                }
                # }}}
            })
            names(crossvalPerf) <- metrics
            # }}}
        } ## end clause split.method$name=="LeaveOneOutBoot"
        if (FALSE) {  ## Dropping this for now, trying BootCv method below
            if (split.method$internal.name=="crossval"){
                N <- splitMethod$N # sample size
                k <- splitMethod$k # k-fold cross validation
                B <- splitMethod$B # number of times we repeat k-fold
                ## here we need a "copy" of the split.method$name=="LeaveOneOutBoot" clause (see above)
                ## and then fix possible bugs/problems which could be related to ordering of data and similar
                ## if it turns out that there are almost no bugs, then it may make sense to use the LeaveOneOutBoot clause
                ## also for k-fold ...

            }}
        if (split.method$name=="BootCv"){
            ## | split.method$internal.name=="crossval"){
                                        # {{{ bootcv
            if (parallel=="snow") exports <- c("DT.B","N.b","cens.model","multi.split.test") else exports <- NULL
            if (!is.null(progress.bar)){
                if (!(progress.bar %in% c(1,2,3))) progress.bar <- 3
                pb <- txtProgressBar(max = B, style = progress.bar,width=20)
            }
            crossval <- foreach(j=1:B,.export=exports,.packages="data.table",.errorhandling="pass") %dopar%{
                ## crossval <- lapply(1:B,function(j){
                DT.b <- DT.B[b==j]
                N.b <- length(unique(DT.b[["ID"]]))
                if(!is.null(progress.bar)){
                    setTxtProgressBar(pb, j)
                }
                computePerformance(DT.b,
                                   N=N.b,
                                   se.fit=FALSE,
                                   conservative=TRUE, ## cannot subset IC yet
                                   cens.model=cens.model,
                                   multi.split.test=multi.split.test,
                                   keep.residuals=FALSE,
                                   keep.vcov=FALSE)
            }
            crossvalPerf <- lapply(metrics,function(m){
                ## score
                if (length(crossval[[1]][[m]]$score)>0){
                    cv.score <- data.table::rbindlist(lapply(crossval,function(x){x[[m]]$score}))
                    if (se.fit==TRUE){
                        bootcv.score <- cv.score[,data.table::data.table(mean(.SD[[m]],na.rm=TRUE),
                                                                         lower=quantile(.SD[[m]],alpha/2,na.rm=TRUE),
                                                                         upper=quantile(.SD[[m]],(1-alpha/2),na.rm=TRUE)),by=byvars,.SDcols=m]
                        data.table::setnames(bootcv.score,c(byvars,m,"lower","upper"))
                    }else{
                        bootcv.score <- cv.score[,data.table::data.table(mean(.SD[[m]],na.rm=TRUE)),by=byvars,.SDcols=m]
                        data.table::setnames(bootcv.score,c(byvars,m))
                    }
                }else{
                    cv.score <- NULL
                    bootcv.score <- NULL
                }
                ## contrasts and multi-split test
                if (length(crossval[[1]][[m]]$contrasts)>0){
                    cv.contrasts <- data.table::rbindlist(lapply(crossval,function(x){x[[m]]$contrasts}))
                    delta.m <- paste0("delta.",m)
                    bootcv.contrasts <- switch(as.character(se.fit+3*multi.split.test),
                                               "4"={
                                                   cv.contrasts[,data.table::data.table(mean(.SD[[delta.m]],na.rm=TRUE),
                                                                                        lower=quantile(.SD[[delta.m]],alpha/2,na.rm=TRUE),
                                                                                        upper=quantile(.SD[[delta.m]],(1-alpha/2),na.rm=TRUE),
                                                                                        p=median(.SD[["p"]],na.rm=TRUE)),by=c(byvars,"reference"),.SDcols=c(delta.m,"p")]
                                               },
                                               "1"={cv.contrasts[,data.table::data.table(mean(.SD[[delta.m]],na.rm=TRUE),
                                                                                         lower=quantile(.SD[[delta.m]],alpha/2,na.rm=TRUE),
                                                                                         upper=quantile(.SD[[delta.m]],(1-alpha/2),na.rm=TRUE)),by=c(byvars,"reference"),.SDcols=c(delta.m)]
                                               },
                                               "3"={cv.contrasts[,data.table::data.table(mean(.SD[[delta.m]],na.rm=TRUE),
                                                                                         p=median(.SD[["p"]],na.rm=TRUE)),by=c(byvars,"reference"),.SDcols=c(delta.m,"p")]
                                               },
                                               "0"={
                                                   bootcv.contrasts <- cv.contrasts[,data.table::data.table(mean(.SD[[delta.m]],na.rm=TRUE)),by=c(byvars,"reference"),.SDcols=delta.m]
                                               })
                    data.table::setnames(bootcv.contrasts,"V1",delta.m)
                }else{
                    cv.contrasts <- NULL
                    bootcv.contrasts <- NULL
                }
                out <- list(score=bootcv.score,contrasts=bootcv.contrasts)
                if (keep.cv)
                    out <- c(out,list(cv.score=cv.score,cv.contrasts=cv.contrasts))
                out
            })
            names(crossvalPerf) <- metrics
            if (ipa==TRUE){
                if (response.type=="binary")
                    crossvalPerf[["Brier"]][["score"]][,IPA:=1-Brier/Brier[model=="Null model"]]
                else
                    crossvalPerf[["Brier"]][["score"]][,IPA:=1-Brier/Brier[model=="Null model"],by=times]
            }
        }
        if (!is.null(progress.bar)){
            cat("\n")
        }
                                        # }}}
                                        # {{{ collect data for plotRisk
        if ("risks"%in% summary){
            crossvalPerf[["risks"]]$score <- DT.B[,{
                c(.SD[1,Response.names,with=FALSE],
                  list(risk=mean(risk),
                       sd.risk=sd(risk),
                       oob=.N))},.SDcols=c(Response.names,"risk"),by=c(byvars,"ID")]
            crossvalPerf[["risks"]]$score[,model:=factor(model,levels=mlevs,mlabels)]
            setcolorder(crossvalPerf[["risks"]]$score,c("ID",byvars,Response.names,"risk","sd.risk"))
        }
                                        # }}}
                                        # {{{ collect data for calibration plots
        if ("Calibration" %in% plots){
            if (keep.residuals[[1]]==TRUE && split.method$name[[1]]=="LeaveOneOutBoot"){
                crossvalPerf[["Calibration"]]$plotframe <- crossvalPerf$Brier$Residuals[model!=0,]
            } else{
                ## there are no residuals in this case. residuals are only available for LOOB!
                crossvalPerf[["Calibration"]]$plotframe <- DT.B[model!=0,{
                    c(.SD[1,Response.names,with=FALSE],
                      list(risk=mean(risk),
                           oob=.N))},.SDcols=c(Response.names,"risk"),by=c(byvars,"ID")]
                setcolorder(crossvalPerf[["Calibration"]]$plotframe,c("ID",byvars,Response.names,"risk"))
            }
            crossvalPerf[["Calibration"]]$plotframe[,model:=factor(model,levels=mlevs,mlabels)]
            if (keep.residuals[[1]]==FALSE && split.method$name[[1]]=="LeaveOneOutBoot"){
                crossvalPerf$Brier$Residuals <- NULL
            }
            if (cens.type=="rightCensored")
                crossvalPerf[["Calibration"]]$plotframe <- merge(jack,crossvalPerf[["Calibration"]]$plotframe,by=c("ID","times"))
        }
                                        # }}}
    }
                                        # }}}
                                        #------------------output-----------------------------------
                                        # {{{ enrich the output object
    if (split.method$internal.name=="noplan"){
        if (keep.residuals==TRUE){
            noSplit$Brier$residuals[,model:=factor(model,levels=mlevs,mlabels)]
        }
        output <- noSplit
    } else{
        output <- crossvalPerf
    }
    models <- mlevs
    names(models) <- mlabels
    if (keep.vcov){
        lab.models <- mlabels
        if (null.model){
            names(lab.models) <- paste0("model=",0:(length(mlabels)-1))
        }
        else{
            names(lab.models) <- paste0("model=",1:(length(mlabels)))
        }
        if (!is.null(output$Brier$vcov))
            attr(output$Brier$vcov,"models") <- lab.models
        if (!is.null(output$AUC$vcov))
            if (response.type=="binary"){
                if (NCOL(output$AUC$vcov)==length(mlabels)){
                    colnames(output$AUC$vcov) <- mlabels
                    rownames(output$AUC$vcov) <- mlabels
                }else{
                    colnames(output$AUC$vcov) <- mlabels[-1]
                    rownames(output$AUC$vcov) <- mlabels[-1]
                }
            }else{
                for (ml in 1:length(mlabels)){
                    colnames(output$AUC$vcov) <- gsub(paste0("model=",ml),paste0("model=",mlabels[ml]),colnames(output$AUC$vcov))
                    rownames(output$AUC$vcov) <- gsub(paste0("model=",ml),paste0("model=",mlabels[ml]),rownames(output$AUC$vcov))}
            }
        attr(output$AUC$vcov,"models") <- lab.models
    }
    if (null.model==TRUE) nm <- names(models)[1] else nm <- NULL
    if (!keep.splitindex) split.method$index <- "split index was not saved"
    output <- c(output,list(response.type=response.type,
                            dolist=dolist,
                            cause=cause,
                            states=states,
                            null.model=nm,
                            models=models,
                            cens.type=cens.type,
                            censoringHandling=cens.method,
                            split.method=split.method,
                            metrics=metrics,
                            times=times,
                            alpha=alpha,
                            plots=plots,
                            summary=summary,
                            missing.predictions=missing.predictions,
                            call=theCall))
    for (p in c(plots)){
        output[[p]]$plotmethod <- p
        class(output[[p]]) <- paste0("score",p)
    }
    for (m in c(metrics,summary)){
        output[[m]]$metrics <- m
        class(output[[m]]) <- paste0("score",m)
    }
    class(output) <- "Score"
    output
                                        # }}}
}

##' @export
Score <- function(object,...){
    UseMethod("Score",object=object)
}

# {{{ Brier score

Brier.binary <- function(DT,
                         se.fit,
                         conservative=FALSE,
                         cens.model="none",
                         keep.vcov=FALSE,
                         multi.split.test,
                         alpha,
                         N,
                         NT,
                         NF,
                         dolist,
                         keep.residuals=FALSE,
                         ...){
    residuals=Brier=risk=model=ReSpOnSe=lower=upper=se=ID=IF.Brier=NULL
    DT[,residuals:=(ReSpOnSe-risk)^2]
    if (se.fit==1L){
        data.table::setkey(DT,model,ID)
        score <- DT[,data.table::data.table(Brier=sum(residuals)/N,se=sd(residuals)/sqrt(N)),by=list(model)]
        score[,lower:=pmax(0,Brier-qnorm(1-alpha/2)*se)]
        score[,upper:=pmin(1,Brier + qnorm(1-alpha/2)*se)]
        ## Influence function IPA
        ## Brier.null <- score[model==0][["Brier"]]
        ## Brier.model <- score[model!=0][["Brier"]]
        ## IC.ipa <- -(1/Brier.null)* DT[model!=0][["residuals"]] + DT[model==0][["residuals"]]*Brier.model/(Brier.null)^2
        if (keep.vcov){
            DT[,Brier:=sum(residuals)/N,by=list(model)]
            DT[,IF.Brier:=residuals-Brier]
        }
    }else{
        ## no se.fit
        score <- DT[,data.table(Brier=sum(residuals)/N),by=list(model)]
    }
    data.table::setkey(score,model)
    if (length(dolist)>0){
        ## merge with Brier score
        data.table::setkey(DT,model)
        data.table::setkey(score,model)
        DT <- DT[score]
        if (se.fit[[1]]==TRUE || multi.split.test[[1]]==TRUE){
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,
                                                             IF=residuals,
                                                             model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=multi.split.test,
                                                  se.fit=se.fit)]
        }else{
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=FALSE,
                                                  se.fit=FALSE)]
        }
        setnames(contrasts.Brier,"delta","delta.Brier")
        output <- list(score=score,contrasts=contrasts.Brier)
    }else{
        output <- list(score=score)
    }
    if (keep.vcov[1] && se.fit[1]==TRUE){
        output <- c(output,list(vcov=getVcov(DT,"IF.Brier")))
    }
    if (keep.residuals) {
        output <- c(output,list(residuals=DT[,data.table::data.table(ID,ReSpOnSe,model,risk,residuals)]))
    }
    output
}

Brier.survival <- function(DT,MC,se.fit,conservative,cens.model,keep.vcov=FALSE,multi.split.test,alpha,N,NT,NF,dolist,keep.residuals=FALSE,...){
    IC0=IPCW=nth.times=ID=time=times=raw.Residuals=risk=Brier=residuals=WTi=Wt=status=setorder=model=IF.Brier=data.table=sd=lower=qnorm=se=upper=NULL
    ## compute 0/1 outcome:
    DT[time<=times & status==1,residuals:=(1-risk)^2/WTi]
    DT[time<=times & status==0,residuals:=0]
    DT[time>times,residuals:=(risk)^2/Wt]
    if (se.fit[[1]]==1L || multi.split.test[[1]]==TRUE){
        ## data.table::setorder(DT,model,times,time,-status)
        data.table::setorder(DT,model,times,ID)
        DT[,nth.times:=as.numeric(factor(times))]
        DT[,IC0:=residuals-mean(residuals),by=list(model,times)]
        if (conservative==TRUE){
            score <- DT[,data.table(Brier=sum(residuals)/N,
                                    se=sd(IC0)/sqrt(N)),by=list(model,times)]
            setnames(DT,"IC0","IF.Brier")
        }else{
            if (cens.model=="none"){
                DT[,IF.Brier:=residuals]
                score <- DT[,data.table(Brier=sum(residuals)/N,
                                        se=sd(residuals)/sqrt(N),
                                        se.conservative=sd(residuals)),by=list(model,times)]
            }else{
                DT[,IF.Brier:=getInfluenceCurve.Brier(t=times[1],
                                                      time=time,
                                                      IC0,
                                                      residuals=residuals,
                                                      WTi=WTi,
                                                      Wt=Wt,
                                                      IC.G=MC,
                                                      cens.model=cens.model,
                                                      nth.times=nth.times[1]),by=list(model,times)]
                score <- DT[,data.table(Brier=sum(residuals)/N,
                                        se=sd(IF.Brier)/sqrt(N),
                                        se.conservative=sd(IC0)/sqrt(N)),by=list(model,times)]
            }
        }
        if (se.fit==TRUE){
            score[,lower:=pmax(0,Brier-qnorm(1-alpha/2)*se)]
            score[,upper:=pmin(1,Brier + qnorm(1-alpha/2)*se)]
        }
    }else{
        ## no se.fit
        score <- DT[,data.table(Brier=sum(residuals)/N),by=list(model,times)]
    }
    data.table::setkey(score,model,times)
    if (length(dolist)>0L){
        ## merge with Brier score
        data.table::setkey(DT,model,times)
        ## data.table::setkey(score,model,times)
        DT <- DT[score]
        if (se.fit[[1]]==TRUE || multi.split.test[[1]]==TRUE){
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,IF=IF.Brier,model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=multi.split.test,
                                                  se.fit=se.fit),by=list(times)]
        }else{
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=FALSE,
                                                  se.fit=FALSE),by=list(times)]
        }
        setnames(contrasts.Brier,"delta","delta.Brier")
        output <- list(score=score,contrasts=contrasts.Brier)
    } else{
        output <- list(score=score)
    }
    if (keep.vcov[1] && se.fit[1]==TRUE){
        output <- c(output,list(vcov=getVcov(DT,"IF.Brier",times=TRUE)))
    }
    if (keep.residuals) {
        if (all(c("Wt","WTi")%in%names(DT))){
            DT[,IPCW:=1/WTi]
            DT[time>=times,IPCW:=1/Wt]
            DT[time<times & status==0,IPCW:=0]
            output <- c(output,list(residuals=DT[,c("ID","time","status","model","times","risk","residuals","IPCW"),with=FALSE]))
        }else{
            output <- c(output,list(residuals=DT[,c("ID","time","status","model","times","risk","residuals"),with=FALSE]))
        }
    }
    output
}

Brier.competing.risks <- function(DT,MC,se.fit,conservative,cens.model,keep.vcov=FALSE,multi.split.test,alpha,N,NT,NF,dolist,keep.residuals=FALSE,cause,states,...){
    IC0=nth.times=ID=time=times=event=Brier=raw.Residuals=risk=residuals=WTi=Wt=status=setorder=model=IF.Brier=data.table=sd=lower=qnorm=se=upper=NULL
    ## compute 0/1 outcome:
    thecause <- match(cause,states,nomatch=0)
    if (length(thecause)==0) stop("Cannot identify cause of interest")
    DT[time<=times & status==1 & event==thecause,residuals:=(1-risk)^2/WTi]
    DT[time<=times & status==1 & event!=thecause,residuals:=(risk)^2/WTi]
    DT[time<=times & status==0,residuals:=0]
    DT[time>times,residuals:=(risk)^2/Wt]
    ## deal with censored observations before t
    DT[time<=times & status==0,residuals:=0]
    if (se.fit[[1]]==1L || multi.split.test[[1]]==TRUE){
        ## data.table::setorder(DT,model,times,time,-status)
        data.table::setorder(DT,model,times,ID)
        DT[,nth.times:=as.numeric(factor(times))]
        DT[,IC0:=residuals-mean(residuals),by=list(model,times)]
        if (conservative){
            score <- DT[,data.table(Brier=sum(residuals)/N,
                                    se=sd(IC0)/sqrt(N)),by=list(model,times)]
        }else{
            if (cens.model=="none"){
                DT[,IF.Brier:=residuals]
                score <- DT[,data.table(Brier=sum(residuals)/N,
                                        se=sd(residuals)/sqrt(N),
                                        se.conservative=sd(residuals)),by=list(model,times)]
            }else{
                DT[,IF.Brier:=getInfluenceCurve.Brier(t=times[1],
                                                      time=time,
                                                      IC0,
                                                      residuals=residuals,
                                                      WTi=WTi,
                                                      Wt=Wt,
                                                      IC.G=MC,
                                                      cens.model=cens.model,
                                                      nth.times=nth.times[1]),by=list(model,times)]
                score <- DT[,data.table(Brier=sum(residuals)/N,
                                        se=sd(IF.Brier)/sqrt(N),
                                        se.conservative=sd(IC0)/sqrt(N)),by=list(model,times)]
            }
        }
        if (se.fit==TRUE){
            score[,lower:=pmax(0,Brier-qnorm(1-alpha/2)*se)]
            score[,upper:=pmin(1,Brier + qnorm(1-alpha/2)*se)]
        }
    }else{
        ## no se.fit
        score <- DT[,data.table(Brier=sum(residuals)/N),by=list(model,times)]
    }
    data.table::setkey(score,model,times)
    if (length(dolist)>0){
        data.table::setkey(DT,model,times)
        ## merge with Brier score
        DT <- DT[score]
        data.table::setkey(score,model,times)
        if (se.fit[[1]]==TRUE || multi.split.test[[1]]==TRUE){
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,IF=IF.Brier,model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=multi.split.test,
                                                  se.fit=se.fit),by=list(times)]
        }else{
            contrasts.Brier <- DT[,getComparisons(data.table(x=Brier,model=model),
                                                  NF=NF,
                                                  N=N,
                                                  alpha=alpha,
                                                  dolist=dolist,
                                                  multi.split.test=FALSE,
                                                  se.fit=FALSE),by=list(times)]
        }
        setnames(contrasts.Brier,"delta","delta.Brier")
        output <- list(score=score,contrasts=contrasts.Brier)
    } else{
        output <- list(score=score)
    }
    if (keep.residuals) {
        output <- c(output,list(residuals=DT[,c("ID","time","status","model","times","risk","residuals"),with=FALSE]))
    }
    if (keep.vcov[1] && se.fit[1]==TRUE){
        output <- c(output,list(vcov=getVcov(DT,"IF.Brier",times=TRUE)))
    }
    output
}

# }}}

# {{{ AUC

## helper functions
AireTrap <- function(FP,TP,N){
    N <- length(FP)
    sum((FP-c(0,FP[-N]))*((c(0,TP[-N])+TP)/2))
}

## do not want to depend on Daim as they turn marker to ensure auc > 0.5
delongtest <-  function(risk,
                        score,
                        dolist,
                        response,
                        cause,
                        alpha,
                        multi.split.test,
                        se.fit,
                        keep.vcov) {
    cov=lower=upper=p=AUC=se=lower=upper=NULL
    Cases <- response == cause
    Controls <- response != cause
    nControls <- sum(!Cases)
    nCases <- sum(Cases)
    ## if (nbCases==0 || nbControls ==0 || length(unique(risk))==1) return(rep(0,n))
    nauc <- ncol(risk)
    auc <- score[["AUC"]]
    modelnames <- score[["model"]]
    riskcontrols <- as.matrix(risk[Controls,])
    riskcases <- as.matrix(risk[Cases,])
    V10 <- matrix(0, nrow = nCases, ncol = nauc)
    V01 <- matrix(0, nrow = nControls, ncol = nauc)
    tmn <- t(riskcontrols)
    tmp <- t(riskcases)
    for (i in 1:nCases) {
        V10[i, ] <- rowSums(tmn < tmp[, i]) + 0.5 * rowSums(tmn == tmp[, i])
    }
    for (i in 1:nControls) {
        V01[i, ] <- rowSums(tmp > tmn[, i]) + 0.5 * rowSums(tmp == tmn[, i])
    }
    V10 <- V10/nControls
    V01 <- V01/nCases
    W10 <- cov(V10)
    W01 <- cov(V01)
    S <- W10/nCases + W01/nControls
    se.auc <- sqrt(diag(S))
    score <- data.table(model=colnames(risk),AUC=auc)
    if (se.fit==1L){
        score[,se:=se.auc]
        score[,lower:=pmax(0,AUC-qnorm(1-alpha/2)*se)]
        score[,upper:=pmin(1,AUC+qnorm(1-alpha/2)*se)]
        setcolorder(score,c("model","AUC","se","lower","upper"))
    }else{
        setcolorder(score,c("model","AUC"))
    }
    names(auc) <- 1:nauc
    ## q1 <- auc/(2 - auc)
    ## q2 <- 2 * auc^2/(1 + auc)
    ## aucvar <- (auc * (1 - auc) + (nCases - 1) * (q1 - auc^2) + (nControls - 1) * (q2 - auc^2))/(nCases * nControls)
    if (length(dolist)>0){
        ncomp <- nauc * (nauc - 1)/2
        delta.AUC <- numeric(ncomp)
        se <- numeric(ncomp)
        model <- numeric(ncomp)
        reference <- numeric(ncomp)
        ctr <- 1
        Qnorm <- qnorm(1 - alpha/2)
        for (i in 1:(nauc - 1)) {
            for (j in (i + 1):nauc) {
                delta.AUC[ctr] <- auc[j]-auc[i]
                if (se.fit==TRUE){
                    ## cor.auc[ctr] <- S[i, j]/sqrt(S[i, i] * S[j, j])
                    LSL <- t(c(1, -1)) %*% S[c(j, i), c(j, i)] %*% c(1, -1)
                    ## print(c(1/LSL,rms::matinv(LSL)))
                    se[ctr] <- sqrt(LSL)
                }
                ## tmpz <- (delta.AUC[ctr]) %*% rms::matinv(LSL) %*% delta.AUC[ctr]
                ## tmpz <- (delta.AUC[ctr]) %*% (1/LSL) %*% delta.AUC[ctr]
                model[ctr] <- modelnames[j]
                reference[ctr] <- modelnames[i]
                ctr <- ctr + 1
            }
        }
        if (se.fit[[1]]==TRUE||multi.split.test[[1]]==TRUE){
            deltaAUC <- data.table(model,reference,delta.AUC=as.vector(delta.AUC),se)
        }else{
            deltaAUC <- data.table(model,reference,delta.AUC=as.vector(delta.AUC))
        }
        if (se.fit==TRUE){
            deltaAUC[,lower:=delta.AUC-Qnorm*se]
            deltaAUC[,upper:=delta.AUC+Qnorm*se]
        }
        if (se.fit[[1]]==TRUE||multi.split.test[[1]]==TRUE){
            deltaAUC[,p:=2*pnorm(abs(delta.AUC)/se,lower.tail=FALSE)]
        }
        out <- list(score = score, contrasts = deltaAUC)
    }else{
        out <- list(score = score, contrasts = NULL)
    }
    if (keep.vcov) {
        out <- c(out,list(vcov=S))
    }
    out
}

auRoc.numeric <- function(X,D,breaks,ROC){
    if (is.null(breaks)) breaks <- rev(sort(unique(X))) ## need to reverse when high X is concordant with {response=1}
    TPR <- c(prodlim::sindex(jump.times=X[D==1],eval.times=breaks,comp="greater",strict=FALSE)/sum(D==1))
    FPR <- c(prodlim::sindex(jump.times=X[D==0],eval.times=breaks,comp="greater",strict=FALSE)/sum(D==0))
    if (ROC==TRUE)
        data.table(risk=breaks,TPR,FPR)
    else
        0.5 * sum(diff(c(0,FPR,0,1)) * (c(TPR,0,1) + c(0,TPR,0)))
}
auRoc.factor <- function(X,D,ROC){
    TPR <- (sum(D==1)-table(X[D==1]))/sum(D==1)
    FPR <- table(X[D==0])/sum(D==0)
    if (ROC==TRUE)
        data.table(cbind(risk=c(sort(unique(X))),TPR,FPR))
    else
        0.5 * sum(diff(c(0,FPR,0,1)) * (c(TPR,0,1) + c(0,TPR,0)))
}

AUC.binary <- function(DT,breaks=NULL,se.fit,conservative=FALSE,cens.model="none",keep.vcov=FALSE,multi.split.test,alpha,N,NT,NF,dolist,ROC,...){
    model=risk=ReSpOnSe=FPR=TPR=ID=NULL
    aucDT <- DT[model>0]
    dolist <- dolist[sapply(dolist,function(do){match("0",do,nomatch=0L)})==0]
    data.table::setkey(aucDT,model,ID)
    if (is.factor(DT[["risk"]])){
        score <- aucDT[,auRoc.factor(risk,ReSpOnSe,ROC=ROC),by=list(model)]
    }
    else{
        score <- aucDT[,auRoc.numeric(risk,ReSpOnSe,breaks=NULL,ROC=ROC),by=list(model)]
    }
    if (ROC==FALSE){
        setnames(score,"V1","AUC")
        output <- list(score=score)
    } else{
        AUC <- score[,list(AUC=0.5 * sum(diff(c(0,FPR,0,1)) * (c(TPR,0,1) + c(0,TPR,0)))),by=list(model)]
        ROC <- score
        output <- list(score=AUC,ROC=ROC)
    }
    if (length(dolist)>0 || (se.fit[[1]]==1L)){
        xRisk <- data.table::dcast(aucDT[],ID~model,value.var="risk")[,-1,with=FALSE]
        delong.res <- delongtest(risk=xRisk,
                                 score=output$score,
                                 dolist=dolist,
                                 response=aucDT[model==model[1],ReSpOnSe],
                                 cause="1",
                                 alpha=alpha,
                                 multi.split.test=multi.split.test,
                                 se.fit=se.fit,
                                 keep.vcov=keep.vcov)
        output$score <- delong.res$score
        output$contrasts <- delong.res$contrasts
        if (keep.vcov){
            output$vcov <- delong.res$vcov
        }
        output
    }else{
        output
    }
}


AUC.survival <- function(DT,MC,se.fit,conservative,cens.model,keep.vcov=FALSE,multi.split.test,alpha,N,NT,NF,dolist,ROC,...){
    ID=model=times=risk=Cases=time=status=Controls=TPR=FPR=WTi=Wt=ipcwControls=ipcwCases=IF.AUC=lower=se=upper=AUC=NULL
    cause <- 1
    aucDT <- DT[model>0]
    ## remove null model comparisons
    dolist <- dolist[sapply(dolist,function(do){match("0",do,nomatch=0L)})==0]
    ## assign Weights before ordering
    aucDT[,ipcwControls:=1/(Wt*N)]
    aucDT[,ipcwCases:=1/(WTi*N)]
    ## order data
    data.table::setorder(aucDT,model,times,-risk)
    ## identify cases and controls
    aucDT[,Cases:=(time <= times &  status==cause)]
    aucDT[,Controls:=(time > times)]
    ## prepare Weights
    aucDT[Cases==0,ipcwCases:=0]
    aucDT[Controls==0,ipcwControls:=0]
    ## compute denominator
    aucDT[,TPR:=cumsum(ipcwCases)/sum(ipcwCases),by=list(model,times)]
    aucDT[,FPR:=(cumsum(ipcwControls))/(sum(ipcwControls)),by=list(model,times)]
    nodups <- aucDT[,c(!duplicated(risk)[-1],TRUE),by=list(model,times)]$V1
    if (ROC==TRUE) {
        output <- list(ROC=aucDT[nodups,c("model","times","risk","TPR","FPR"),with=FALSE])
    }else{
        output <- NULL
    }
    score <- aucDT[nodups,list(AUC=AireTrap(FPR,TPR)),by=list(model,times)]
    data.table::setkey(score,model,times)
    if (se.fit[[1]]==1L || multi.split.test[[1]]==TRUE){
        ## compute influence function
        ## data.table::setorder(aucDT,model,times,time,-status)
        data.table::setorder(aucDT,model,times,ID)
        aucDT[,IF.AUC:=getInfluenceCurve.AUC.survival(t=times[1],
                                                      n=N,
                                                      time=time,
                                                      risk=risk,
                                                      Cases=Cases,
                                                      Controls=Controls,
                                                      ipcwControls=ipcwControls,
                                                      ipcwCases=ipcwCases,
                                                      MC=MC), by=list(model,times)]
        se.score <- aucDT[,list(se=sd(IF.AUC)/sqrt(N)),by=list(model,times)]
        data.table::setkey(se.score,model,times)
        score <- score[se.score]
        if (se.fit==1L){
            score[,lower:=pmax(0,AUC-qnorm(1-alpha/2)*se)]
            score[,upper:=pmin(1,AUC+qnorm(1-alpha/2)*se)]
        }else{
            score[,se:=NULL]
        }
        data.table::setkey(aucDT,model,times)
        aucDT <- aucDT[score]
        if (keep.vcov){
            output <- c(output,list(vcov=getVcov(aucDT,"IF.AUC",times=TRUE)))
        }
    }
    ## add score to object
    output <- c(list(score=score),output)
    if (length(dolist)>0){
        if (se.fit[[1]]==TRUE || multi.split.test[[1]]==TRUE){
            contrasts.AUC <- aucDT[,getComparisons(data.table(x=AUC,IF=IF.AUC,model=model),
                                                   NF=NF,
                                                   N=N,
                                                   alpha=alpha,
                                                   dolist=dolist,multi.split.test=multi.split.test,se.fit=se.fit),by=list(times)]
        }else{
            contrasts.AUC <- score[,getComparisons(data.table(x=AUC,model=model),
                                                   NF=NF,
                                                   N=N,
                                                   alpha=alpha,
                                                   dolist=dolist,
                                                   multi.split.test=FALSE,
                                                   se.fit=FALSE),by=list(times)]
        }
        setnames(contrasts.AUC,"delta","delta.AUC")
        output <- c(list(score=score,contrasts=contrasts.AUC),output)
    }
    output
}

AUC.competing.risks <- function(DT,MC,se.fit,conservative,cens.model,keep.vcov=FALSE,multi.split.test,alpha,N,NT,NF,dolist,cause,states,ROC,...){
    ID=model=times=risk=Cases=time=status=event=Controls1=Controls2=TPR=FPR=WTi=Wt=ipcwControls1=ipcwControls2=ipcwCases=IF.AUC=lower=se=upper=AUC=NULL
    aucDT <- DT[model>0]
    dolist <- dolist[sapply(dolist,function(do){match("0",do,nomatch=0L)})==0]
    ## assign Weights before ordering
    aucDT[,ipcwControls1:=1/(Wt*N)]
    aucDT[,ipcwControls2:=1/(WTi*N)]
    aucDT[,ipcwCases:=1/(WTi*N)]
    aucDT[,ipcwControls2:=1/(WTi*N)]
    ## order data
    data.table::setorder(aucDT,model,times,-risk)
    ## identify cases and controls
    thecause <- match(cause,states,nomatch=0)
    if (length(thecause)==0) stop("Cannot identify cause of interest")
    aucDT[,Cases:=(time <=times &  event==thecause)]
    aucDT[,Controls1:=(time > times)]
    aucDT[,Controls2:=(time <=times &  event!=thecause & status !=0)]
    ## prepare Weights
    aucDT[Cases==0,ipcwCases:=0]
    aucDT[Controls1==0,ipcwControls1:=0]
    aucDT[Controls2==0,ipcwControls2:=0]
    ## compute denominator
    ## ROC <- aucDT[,list(TPR=c(0,cumsum(ipcwCases)),FPR=c(0,cumsum(ipcwControls1)+cumsum(ipcwControls2))),by=list(model,times)]
    aucDT[,TPR:=cumsum(ipcwCases)/sum(ipcwCases),by=list(model,times)]
    aucDT[,FPR:=(cumsum(ipcwControls1)+cumsum(ipcwControls2))/(sum(ipcwControls2)+sum(ipcwControls1)),by=list(model,times)]
    nodups <- aucDT[,c(!duplicated(risk)[-1],TRUE),by=list(model,times)]$V1
    if (ROC==TRUE) {
        output <- list(ROC=aucDT[nodups,c("model","times","risk","TPR","FPR"),with=FALSE])
    }else{
        output <- NULL
    }
    score <- aucDT[nodups,list(AUC=AireTrap(FPR,TPR)),by=list(model,times)]
    data.table::setkey(score,model,times)
    if (se.fit[[1]]==1L || multi.split.test[[1]]==TRUE){
        ## compute influence function
        ## data.table::setorder(aucDT,model,times,time,-status)
        data.table::setorder(aucDT,model,times,ID)
        aucDT[,IF.AUC:={
            if (sum(Controls2)==0){
                getInfluenceCurve.AUC.survival(t=times[1],
                                               n=N,
                                               time=time,
                                               risk=risk,
                                               Cases=Cases,
                                               Controls=Controls1,
                                               ipcwControls=ipcwControls1,
                                               ipcwCases=ipcwCases,
                                               MC=MC)
            }else{
                getInfluenceCurve.AUC.competing.risks(t=times[1],
                                                      n=N,
                                                      time=time,
                                                      risk=risk,
                                                      ipcwControls1=ipcwControls1,
                                                      ipcwControls2=ipcwControls2,
                                                      ipcwCases=ipcwCases,
                                                      Cases=Cases,
                                                      Controls1=Controls1,
                                                      Controls2=Controls2,
                                                      MC=MC)
            }
        }, by=list(model,times)]
        se.score <- aucDT[,list(se=sd(IF.AUC)/sqrt(N)),by=list(model,times)]
        data.table::setkey(se.score,model,times)
        score <- score[se.score]
        if (se.fit==1L){
            score[,lower:=pmax(0,AUC-qnorm(1-alpha/2)*se)]
            score[,upper:=pmin(1,AUC+qnorm(1-alpha/2)*se)]
        }else{
            score[,se:=NULL]
        }
        data.table::setkey(aucDT,model,times)
        aucDT <- aucDT[score]
        if (keep.vcov){
            output <- c(output,list(vcov=getVcov(aucDT,"IF.AUC",times=TRUE)))
        }
    }
    ## add score to object
    output <- c(list(score=score),output)
    if (length(dolist)>0){
        if (se.fit[[1]]==TRUE || multi.split.test[[1]]==TRUE){
            contrasts.AUC <- aucDT[,getComparisons(data.table(x=AUC,IF=IF.AUC,model=model),
                                                   NF=NF,
                                                   N=N,
                                                   alpha=alpha,
                                                   dolist=dolist,multi.split.test=multi.split.test,
                                                   se.fit=se.fit),by=list(times)]
        }else{
            contrasts.AUC <- score[,getComparisons(data.table(x=AUC,model=model),
                                                   NF=NF,
                                                   N=N,
                                                   alpha=alpha,
                                                   dolist=dolist,
                                                   multi.split.test=FALSE,
                                                   se.fit=FALSE),by=list(times)]
        }
        setnames(contrasts.AUC,"delta","delta.AUC")
        output <- c(list(score=score,contrasts=contrasts.AUC),output)
    }
    output
}

# }}}

Try the riskRegression package in your browser

Any scripts or data that you put into this service are public.

riskRegression documentation built on Jan. 13, 2021, 11:12 a.m.