R/Coxmos_sb_splsicox.R

In Coxmos: Cox MultiBlock Survival

#### ### ##
# METHODS #
#### ### ##

#' SB.sPLS-ICOX
#' @description This function performs a single-block sparse partial least squares individual Cox
#' (SB.sPLS-ICOX). The function returns a Coxmos model with the attribute model as "SB.sPLS-ICOX".
#'
#' @details
#' The `SB.sPLS-ICOX` function is designed to perform a single-block sparse partial least squares
#' individual Cox analysis. This method is particularly suited for high-dimensional datasets where
#' the number of variables (features) significantly exceeds the number of observations. The
#' "single-block" in its name indicates that while the function can handle datasets with multiple
#' blocks, it processes each block individually rather than in a multiblock manner where all blocks
#' are analyzed simultaneously.
#'
#' By analyzing one block at a time, the function ensures a focused and detailed examination of each
#' block's contribution to the survival outcome. This approach is especially beneficial when different
#' blocks represent distinct types or sources of data, allowing for a granular understanding of each
#' block's significance.
#'
#' The analysis begins by applying a penalty to select significant variables based on individual Cox
#' models. This step ensures that only the most relevant features from the current block contribute
#' to the subsequent sPLS analysis. The sPLS method then identifies latent components that capture
#' the maximum covariance between the explanatory variables (X) from the block and the response (Y),
#' which are the deviance residuals from the Cox models.
#'
#' Users have the flexibility to specify various hyperparameters, including the number of latent
#' components and the penalty for variable selection. The function also offers options for data
#' preprocessing, such as centering, scaling, and removing variables with near-zero or zero variance.
#'
#' The output provides a comprehensive overview of the analysis for the processed block, including
#' normalized data information, survival model details, and the sPLS-ICOX model. Visualization tools
#' and metrics such as AIC and BIC further aid in understanding the model's performance and
#' significance for the given block.
#'
#' In summary, the `SB.sPLS-ICOX` function offers a powerful approach for survival analysis in
#' high-dimensional settings, ensuring optimal feature selection, dimensionality reduction, and
#' predictive modeling for each individual block in the dataset.
#'
#' @param X List of numeric matrices or data.frames. Explanatory variables. Qualitative variables must be
#' transform into binary variables.
#' @param Y Numeric matrix or data.frame. Response variables. Object must have two columns named as
#' "time" and "event". For event column, accepted values are: 0/1 or FALSE/TRUE for censored and
#' event observations.
#' @param n.comp Numeric. Number of latent components to compute for the (s)PLS model (default: 10).
#' @param penalty Numeric. Penalty for variable selection for the individual cox models. Variables
#' with a lower P-Value than 1 - "penalty" in the individual cox analysis will be keep for the
#' sPLS-ICOX approach (default: 1).
#' @param x.center Logical. If x.center = TRUE, X matrix is centered to zero means (default: TRUE).
#' @param x.scale Logical. If x.scale = TRUE, X matrix is scaled to unit variances (default: FALSE).
#' @param remove_near_zero_variance Logical. If remove_near_zero_variance = TRUE, near zero variance
#' variables will be removed (default: TRUE).
#' @param remove_zero_variance Logical. If remove_zero_variance = TRUE, zero variance variables will
#' be removed (default: TRUE).
#' @param toKeep.zv Character vector. Name of variables in X to not be deleted by (near) zero variance
#' filtering (default: NULL).
#' @param remove_non_significant Logical. If remove_non_significant = TRUE, non-significant
#' variables/components in final cox model will be removed until all variables are significant by
#' forward selection (default: FALSE).
#' @param alpha Numeric. Numerical values are regarded as significant if they fall below the
#' threshold (default: 0.05).
#' @param MIN_EPV Numeric. Minimum number of Events Per Variable (EPV) you want reach for the final
#' cox model. Used to restrict the number of variables/components can be computed in final cox models.
#' If the minimum is not meet, the model cannot be computed (default: 5).
#' @param returnData Logical. Return original and normalized X and Y matrices (default: TRUE).
#' @param verbose Logical. If verbose = TRUE, extra messages could be displayed (default: FALSE).
#'
#' @return Instance of class "Coxmos" and model "sb.splsicox". The class contains the following
#' elements:
#' \code{X}: List of normalized X data information.
#' \itemize{
#'  \item \code{(data)}: normalized X matrix
#'  \item \code{(weightings)}: PLS weights
#'  \item \code{(weightings_norm)}: PLS normalize weights
#'  \item \code{(W.star)}: PLS W* vector
#'  \item \code{(scores)}: PLS scores/variates
#'  \item \code{(x.mean)}: mean values for X matrix
#'  \item \code{(x.sd)}: standard deviation for X matrix
#'  }
#' \code{Y}: List of normalized Y data information.
#' \itemize{
#'  \item \code{(deviance_residuals)}: deviance residual vector used as Y matrix in the sPLS.
#'  \item \code{(dr.mean)}: mean values for deviance residuals Y matrix
#'  \item \code{(dr.sd)}: standard deviation for deviance residuals Y matrix'
#'  \item \code{(data)}: normalized X matrix
#'  \item \code{(y.mean)}: mean values for Y matrix
#'  \item \code{(y.sd)}: standard deviation for Y matrix'
#'  }
#' \code{survival_model}: List of survival model information.
#' \itemize{
#'  \item \code{fit}: coxph object.
#'  \item \code{AIC}: AIC of cox model.
#'  \item \code{BIC}: BIC of cox model.
#'  \item \code{lp}: linear predictors for train data.
#'  \item \code{coef}: Coefficients for cox model.
#'  \item \code{YChapeau}: Y Chapeau residuals.
#'  \item \code{Yresidus}: Y residuals.
#' }
#'
#' \code{list_spls_models}: List of sPLS-ICOX models computed for each block.
#'
#' \code{n.comp}: Number of components selected.
#'
#' \code{penalty} Penalty applied.
#'
#' \code{call}: call function
#'
#' \code{X_input}: X input matrix
#'
#' \code{Y_input}: Y input matrix
#'
#' \code{nzv}: Variables removed by remove_near_zero_variance or remove_zero_variance.
#'
#' \code{nz_coeffvar}: Variables removed by coefficient variation near zero.
#'
#' \code{class}: Model class.
#'
#' \code{time}: time consumed for running the cox analysis.
#'
#' @author Pedro Salguero Garcia. Maintainer: pedsalga@upv.edu.es
#'
#' @export
#'
#' @examples
#' data("X_multiomic")
#' data("Y_multiomic")
#' X <- X_multiomic
#' X$mirna <- X$mirna[,1:50]
#' X$proteomic <- X$proteomic[,1:50]
#' Y <- Y_multiomic
#' sb.splsicox(X, Y, n.comp = 2, penalty = 0.5, x.center = TRUE, x.scale = TRUE)

sb.splsicox <- function(X, Y,
                        n.comp = 4, penalty = 0,
                        x.center = TRUE, x.scale = FALSE,
                        remove_near_zero_variance = TRUE, remove_zero_variance = TRUE, toKeep.zv = NULL,
                        remove_non_significant = FALSE, alpha = 0.05,
                        MIN_EPV = 5, returnData = TRUE, verbose = FALSE){
  # tol Numeric. Tolerance for solving: solve(t(P) %*% W) (default: 1e-15).
  tol = 1e-10

  t1 <- Sys.time()
  y.center = y.scale = FALSE
  FREQ_CUT <- 95/5

  #### Check values classes and ranges
  params_with_limits <- list("penalty" = penalty)
  check_min0_less1_variables(params_with_limits)

  params_with_limits <- list("alpha" = alpha)
  check_min0_max1_variables(params_with_limits)

  numeric_params <- list("n.comp" = n.comp,
                         "MIN_EPV" = MIN_EPV, "tol" = tol)
  check_class(numeric_params, class = "numeric")

  logical_params <- list("x.center" = unlist(x.center), "x.scale" = unlist(x.scale),
                         #"y.center" = y.center, "y.scale" = y.scale,
                         "remove_near_zero_variance" = remove_near_zero_variance, "remove_zero_variance" = remove_zero_variance,
                         "remove_non_significant" = remove_non_significant, "returnData" = returnData, "verbose" = verbose)
  check_class(logical_params, class = "logical")

  #### Check rownames
  lst_check <- checkXY.rownames.mb(X, Y, verbose = verbose)
  X <- lst_check$X
  Y <- lst_check$Y

  #### Check colnames
  X <- checkColnamesIllegalChars.mb(X)

  #### REQUIREMENTS
  checkX.colnames.mb(X)
  checkY.colnames(Y)
  lst_check <- checkXY.mb.class(X, Y, verbose = verbose)
  X <- lst_check$X
  Y <- lst_check$Y

  #### Original data
  X_original <- X
  Y_original <- Y

  time <- Y[,"time"]
  event <- Y[,"event"]

  #### ZERO VARIANCE - ALWAYS bc run/fold
  lst_dnz <- deleteZeroOrNearZeroVariance.mb(X = X,
                                             remove_near_zero_variance = remove_near_zero_variance,
                                             remove_zero_variance = remove_zero_variance,
                                             toKeep.zv = toKeep.zv,
                                             freqCut = FREQ_CUT)
  X <- lst_dnz$X
  variablesDeleted <- lst_dnz$variablesDeleted

  #### COEF VARIATION - ALWAYS bc run/fold
  lst_dnzc <- deleteNearZeroCoefficientOfVariation.mb(X = X)
  X <- lst_dnzc$X
  variablesDeleted_cvar <- lst_dnzc$variablesDeleted

  #### SCALING
  lst_scale <- XY.mb.scale(X = X, Y = Y, x.center = x.center, x.scale = x.scale, y.center = y.center, y.scale = y.scale)
  Xh <- lst_scale$Xh
  Yh <- lst_scale$Yh
  xmeans <- lst_scale$xmeans
  xsds <- lst_scale$xsds
  ymeans <- lst_scale$ymeans
  ysds <- lst_scale$ysds

  X_norm <- Xh

  #### MAX PREDICTORS
  n.comp <- check.mb.maxPredictors(X, Y, MIN_EPV, n.comp, verbose = verbose)

  # CREATE INDIVIDUAL MODELS
  lst_sb.pls <- purrr::map(names(Xh), ~splsicox(X = Xh[[.]], Y = Yh, n.comp = n.comp, penalty = penalty,
                                                x.scale = FALSE, x.center = FALSE,
                                                #y.scale = FALSE, y.center = FALSE,
                                                remove_near_zero_variance = FALSE, remove_zero_variance = FALSE,
                                                toKeep.zv = NULL, #zero_var already checked
                                                remove_non_significant = remove_non_significant, alpha = alpha,
                                                returnData = FALSE, verbose = verbose))

  ## TEST ##
  # aux <- splsicox(X = Xh$clinical, Y = Yh, n.comp = n.comp, penalty = penalty,
  #          x.scale = FALSE, x.center = FALSE,
  #          #y.scale = FALSE, y.center = FALSE,
  #          remove_near_zero_variance = FALSE, remove_zero_variance = FALSE,
  #          toKeep.zv = NULL, #zero_var already checked
  #          remove_non_significant = remove_non_significant, alpha = alpha,
  #          returnData = FALSE, verbose = verbose)

  names(lst_sb.pls) <- names(Xh)

  # CHECK - NO MODEL FOR BLOCKS
  model0comp <- NULL
  for(b in names(lst_sb.pls)){
    if(lst_sb.pls[[b]]$n.comp == 0){
      model0comp <- c(model0comp, b)
    }
  }

  lst_sb.pls[model0comp] <- NULL

  # CHECK ALL MODELS SAME COMPONENTS
  aux_ncomp <- purrr::map(lst_sb.pls, ~ifelse("n.comp" %in% names(.),.$n.comp, NA))

  # CREATE COMBINE MODEL
  data <- NULL
  cn.merge <- NULL
  for(b in names(lst_sb.pls)){

    # NA or NULL
    if(!isa(lst_sb.pls[[b]], "Coxmos")){
      next
    }

    if(!is.null(lst_sb.pls[[b]]$survival_model)){
      data <- cbind(data, lst_sb.pls[[b]]$X$scores)
      cn.merge <- c(cn.merge, paste0(colnames(lst_sb.pls[[b]]$X$scores), "_", b))
    }else{
      next
    }
  }

  colnames(data) <- cn.merge
  cox_model <- cox(X = data, Y = Yh,
                   x.center = FALSE, x.scale = FALSE,
                   #y.center = FALSE, y.scale = FALSE,
                   remove_near_zero_variance = FALSE, remove_zero_variance = FALSE,
                   remove_non_significant = remove_non_significant, FORCE = TRUE)

  # RETURN a MODEL with ALL significant Variables from complete, deleting one by one
  removed_variables <- NULL
  removed_variables_cor <- NULL
  # REMOVE NA-PVAL VARIABLES
  # p_val could be NA for some variables (if NA change to P-VAL=1)
  # DO IT ALWAYS, we do not want problems in COX models
  if(all(c("time", "event") %in% colnames(data))){
    lst_model <- removeNAorINFcoxmodel(model = cox_model$survival_model$fit, data = data, time.value = NULL, event.value = NULL)
  }else{
    lst_model <- removeNAorINFcoxmodel(model = cox_model$survival_model$fit, data = cbind(data, Yh), time.value = NULL, event.value = NULL)
  }
  cox_model$survival_model$fit <- lst_model$model
  removed_variables_cor <- c(removed_variables_cor, lst_model$removed_variables)

  # RETURN a MODEL with ALL significant Variables from complete, deleting one by one in backward method
  # already performed in cox() function
  if(remove_non_significant){
    removed_variables <- cox_model$nsv
  }else{
    removed_variables <- NULL
  }

  #### ### #
  # RETURN #
  #### ### #
  func_call <- match.call()

  if(isa(cox_model$survival_model$fit,"coxph")){
    survival_model <- getInfoCoxModel(cox_model$survival_model$fit)
  }else{
    survival_model <- NULL
  }

  if(!returnData){
    survival_model <- removeInfoSurvivalModel(survival_model)
  }else{
    survival_model <- survival_model
  }

  all_scores <- NULL
  for(b in names(lst_sb.pls)){
    aux_scores <- lst_sb.pls[[b]]$X$scores
    colnames(aux_scores) <- paste0(colnames(aux_scores), "_", b)
    all_scores <- cbind(all_scores, aux_scores)
  }

  t2 <- Sys.time()
  time <- difftime(t2,t1,units = "mins")

  # invisible(gc())
  return(sb.splsicox_class(list(X = list("data" = if(returnData) X_norm else NA,
                                         "scores" = all_scores,
                                         "x.mean" = xmeans, "x.sd" = xsds),
                                Y = list("data" = Yh,
                                         "y.mean" = ymeans, "y.sd" = ysds),
                                survival_model = survival_model,
                                list_spls_models = lst_sb.pls,
                                n.comp = aux_ncomp, #number of components used, but could be lesser than expected because not computed models
                                penalty = penalty,
                                call = if(returnData) func_call else NA,
                                X_input = if(returnData) X_original else NA,
                                Y_input = if(returnData) Y_original else NA,
                                alpha = alpha,
                                nsv = removed_variables,
                                nzv = variablesDeleted,
                                nz_coeffvar = variablesDeleted_cvar,
                                class = pkg.env$sb.splsicox,
                                time = time)))
}

#### ### ### ### ###
# CROSS-EVALUATION #
#### ### ### ### ###

#' Cross validation cv.sb.splsicox
#' @description This function performs cross-validated sparse partial least squares single block for splsicox.
#' The function returns the optimal number of components and the optimal sparsity penalty value based
#' on cross-validation. The performance could be based on multiple metrics as Area Under the Curve
#' (AUC), I. Brier Score or C-Index. Furthermore, the user could establish more than one metric
#' simultaneously.
#'
#' @details
#' The `cv.sb.splsicox` function performs cross-validation for the single-block sparse partial least
#' squares individual Cox analysis. While the function can handle datasets with multiple blocks, it
#' processes each block individually, ensuring a detailed examination of each block's contribution to
#' the survival outcome. This is distinct from multiblock methods where all blocks are analyzed
#' simultaneously.
#'
#' In the context of this function, "single-block" means that each block of data is analyzed
#' separately, one at a time. This approach is beneficial when different blocks represent distinct
#' types or sources of data, allowing for a granular understanding of each block's significance
#' without the interference of other blocks.
#'
#' The cross-validation process involves partitioning the dataset into multiple subsets (folds) and
#' then iteratively training the model on a subset of the data while validating it on the remaining
#' data. This helps in determining the optimal hyperparameters for the model, such as the number of
#' latent components and the penalty for variable selection.
#'
#' The function offers flexibility in specifying various hyperparameters and options for data
#' preprocessing. The output provides a comprehensive overview of the cross-validation results,
#' including metrics like AIC, C-Index, I. Brier Score, and AUC for each hyper-parameter combination.
#' Visualization tools are also provided to aid in understanding the model's performance across
#' different hyperparameters.
#'
#' In summary, the `cv.sb.splsicox` function offers a robust approach for determining the optimal
#' parameters for the single-block sparse partial least squares individual Cox analysis, ensuring
#' optimal feature selection, dimensionality reduction, and predictive modeling for each individual
#' block in the dataset.
#'
#' @param X List of numeric matrices or data.frames. Explanatory variables. Qualitative variables must be
#' transform into binary variables.
#' @param Y Numeric matrix or data.frame. Response variables. Object must have two columns named as
#' "time" and "event". For event column, accepted values are: 0/1 or FALSE/TRUE for censored and
#' event observations.
#' @param max.ncomp Numeric. Maximum number of PLS components to compute for the cross validation
#' (default: 8).
#' @param penalty.list Numeric vector. Penalty for variable selection for the individual cox
#' models. Variables with a lower P-Value than 1 - "penalty" in the individual cox analysis will
#' be keep for the sPLS-ICOX approach (default: seq(0.1,0.9,0.2)).
#' @param n_run Numeric. Number of runs for cross validation (default: 3).
#' @param k_folds Numeric. Number of folds for cross validation (default: 10).
#' @param x.center Logical. If x.center = TRUE, X matrix is centered to zero means (default: TRUE).
#' @param x.scale Logical. If x.scale = TRUE, X matrix is scaled to unit variances (default: FALSE).
#' @param remove_near_zero_variance Logical. If remove_near_zero_variance = TRUE, near zero variance
#' variables will be removed (default: TRUE).
#' @param remove_zero_variance Logical. If remove_zero_variance = TRUE, zero variance variables will
#' be removed (default: TRUE).
#' @param toKeep.zv Character vector. Name of variables in X to not be deleted by (near) zero variance
#' filtering (default: NULL).
#' @param remove_variance_at_fold_level Logical. If remove_variance_at_fold_level = TRUE, (near) zero
#' variance will be removed at fold level. Not recommended. (default: FALSE).
#' @param remove_non_significant_models Logical. If remove_non_significant_models = TRUE,
#' non-significant models are removed before computing the evaluation.
#' @param remove_non_significant Logical. If remove_non_significant = TRUE, non-significant
#' variables/components in final cox model will be removed until all variables are significant by
#' forward selection (default: FALSE).
#' @param alpha Numeric. Numerical values are regarded as significant if they fall below the
#' threshold (default: 0.05).
#' @param w_AIC Numeric. Weight for AIC evaluator. All weights must sum 1 (default: 0).
#' @param w_C.Index Numeric. Weight for C-Index evaluator. All weights must sum 1 (default: 0).
#' @param w_AUC Numeric. Weight for AUC evaluator. All weights must sum 1 (default: 1).
#' @param w_I.BRIER Numeric. Weight for BRIER SCORE evaluator. All weights must sum 1 (default: 0).
#' @param times Numeric vector. Time points where the AUC will be evaluated. If NULL, a maximum of
#' 'max_time_points' points will be selected equally distributed (default: NULL).
#' @param max_time_points Numeric. Maximum number of time points to use for evaluating the model
#' (default: 15).
#' @param MIN_AUC_INCREASE Numeric. Minimum improvement between different cross validation models to
#' continue evaluating higher values in the multiple tested parameters. If it is not reached for next
#' 'MIN_COMP_TO_CHECK' models and the minimum 'MIN_AUC' value is reached, the evaluation stops
#' (default: 0.01).
#' @param MIN_AUC Numeric. Minimum AUC desire to reach cross-validation models. If the minimum is
#' reached, the evaluation could stop if the improvement does not reach an AUC higher than adding
#' the 'MIN_AUC_INCREASE' value (default: 0.8).
#' @param MIN_COMP_TO_CHECK Numeric. Number of penalties/components to evaluate to check if the AUC
#' improves. If for the next 'MIN_COMP_TO_CHECK' the AUC is not better and the 'MIN_AUC' is meet,
#' the evaluation could stop (default: 3).
#' @param pred.attr Character. Way to evaluate the metric selected. Must be one of the following:
#' "mean" or "median" (default: "mean").
#' @param pred.method Character. AUC evaluation algorithm method for evaluate the model performance.
#' Must be one of the following: "risksetROC", "survivalROC", "cenROC", "nsROC", "smoothROCtime_C",
#' "smoothROCtime_I" (default: "cenROC").
#' @param fast_mode Logical. If fast_mode = TRUE, for each run, only one fold is evaluated
#' simultaneously. If fast_mode = FALSE, for each run, all linear predictors are computed for test
#' observations. Once all have their linear predictors, the evaluation is perform across all the
#' observations together (default: FALSE).
#' @param MIN_EPV Numeric. Minimum number of Events Per Variable (EPV) you want reach for the final
#' cox model. Used to restrict the number of variables/components can be computed in final cox models.
#' If the minimum is not meet, the model cannot be computed (default: 5).
#' @param return_models Logical. Return all models computed in cross validation (default: FALSE).
#' @param returnData Logical. Return original and normalized X and Y matrices (default: TRUE).
#' @param PARALLEL Logical. Run the cross validation with multicore option. As many cores as your
#' total cores - 1 will be used. It could lead to higher RAM consumption (default: FALSE).
#' @param verbose Logical. If verbose = TRUE, extra messages could be displayed (default: FALSE).
#' @param seed Number. Seed value for performing runs/folds divisions (default: 123).
#'
#' @return Instance of class "Coxmos" and model "cv.SB.sPLS-ICOX".
#' \code{best_model_info}: A data.frame with the information for the best model.
#' \code{df_results_folds}: A data.frame with fold-level information.
#' \code{df_results_runs}: A data.frame with run-level information.
#' \code{df_results_comps}: A data.frame with component-level information (for cv.coxEN, EN.alpha
#' information).
#'
#' \code{lst_models}: If return_models = TRUE, return a the list of all cross-validated models.
#' \code{pred.method}: AUC evaluation algorithm method for evaluate the model performance.
#'
#' \code{opt.comp}: Optimal component selected by the best_model.
#' \code{opt.penalty}: Optimal penalty value selected by the best_model.
#'
#' \code{plot_AIC}: AIC plot by each hyper-parameter.
#' \code{plot_C.Index}: C-Index plot by each hyper-parameter.
#' \code{plot_I.BRIER}: Integrative Brier Score plot by each hyper-parameter.
#' \code{plot_AUC}: AUC plot by each hyper-parameter.
#'
#' \code{class}: Cross-Validated model class.
#'
#' \code{lst_train_indexes}: List (of lists) of indexes for the observations used in each run/fold
#' for train the models.
#' \code{lst_test_indexes}: List (of lists) of indexes for the observations used in each run/fold
#' for test the models.
#'
#' \code{time}: time consumed for running the cross-validated function.
#'
#' @author Pedro Salguero Garcia. Maintainer: pedsalga@upv.edu.es
#'
#' @export
#'
#' @examples
#' data("X_multiomic")
#' data("Y_multiomic")
#' set.seed(123)
#' index_train <- caret::createDataPartition(Y_multiomic$event, p = .5, list = FALSE, times = 1)
#' X_train <- X_multiomic
#' X_train$mirna <- X_train$mirna[index_train,1:50]
#' X_train$proteomic <- X_train$proteomic[index_train,1:50]
#' Y_train <- Y_multiomic[index_train,]
#' cv.sb.splsicox_model <- cv.sb.splsicox(X_train, Y_train, max.ncomp = 2, penalty.list = c(0.5),
#' n_run = 1, k_folds = 2, x.center = TRUE, x.scale = TRUE)

cv.sb.splsicox <- function(X, Y,
                           max.ncomp = 8, penalty.list = seq(0.1,0.9,0.2),
                           n_run = 3, k_folds = 10,
                           x.center = TRUE, x.scale = FALSE,
                           remove_near_zero_variance = TRUE, remove_zero_variance = TRUE, toKeep.zv = NULL,
                           remove_variance_at_fold_level = FALSE,
                           remove_non_significant_models = FALSE, remove_non_significant = FALSE, alpha = 0.05,
                           w_AIC = 0, w_C.Index = 0, w_AUC = 1, w_I.BRIER = 0, times = NULL,
                           max_time_points = 15,
                           MIN_AUC_INCREASE = 0.01, MIN_AUC = 0.8, MIN_COMP_TO_CHECK = 3,
                           pred.attr = "mean", pred.method = "cenROC", fast_mode = FALSE,
                           MIN_EPV = 5, return_models = FALSE, returnData = FALSE,
                           PARALLEL = FALSE, verbose = FALSE, seed = 123){
  # tol Numeric. Tolerance for solving: solve(t(P) %*% W) (default: 1e-15).
  tol = 1e-10

  t1 <- Sys.time()
  y.center = y.scale = FALSE
  FREQ_CUT <- 95/5

  #### ### ###
  # WARNINGS #
  #### ### ###

  #### Check evaluator installed:
  checkLibraryEvaluator(pred.method)

  #### Check values classes and ranges
  params_with_limits <- list("penalty.list" = penalty.list)
  check_min0_less1_variables(params_with_limits)

  params_with_limits <- list("MIN_AUC_INCREASE" = MIN_AUC_INCREASE, "MIN_AUC" = MIN_AUC, "alpha" = alpha,
                 "w_AIC" = w_AIC, "w_C.Index" = w_C.Index, "w_AUC" = w_AUC, "w_I.BRIER" = w_I.BRIER)
  check_min0_max1_variables(params_with_limits)

  numeric_params <- list("max.ncomp" = max.ncomp,
                  "n_run" = n_run, "k_folds" = k_folds, "max_time_points" = max_time_points,
                  "MIN_COMP_TO_CHECK" = MIN_COMP_TO_CHECK, "MIN_EPV" = MIN_EPV, "seed" = seed, "tol" = tol)
  check_class(numeric_params, class = "numeric")

  logical_params <- list("x.center" = unlist(x.center), "x.scale" = unlist(x.scale),
                         #"y.center" = y.center, "y.scale" = y.scale,
                      "remove_near_zero_variance" = remove_near_zero_variance, "remove_zero_variance" = remove_zero_variance,
                      "remove_variance_at_fold_level" = remove_variance_at_fold_level,
                      "remove_non_significant_models" = remove_non_significant_models,
                      "remove_non_significant" = remove_non_significant,
                      "return_models" = return_models,"returnData" = returnData, "verbose" = verbose, "PARALLEL" = PARALLEL)
  check_class(logical_params, class = "logical")

  character_params <- list("pred.attr" = pred.attr, "pred.method" = pred.method)
  check_class(character_params, class = "character")

  #### FIX possible SEQ() problems
  penalty.list <- as.character(penalty.list)
  penalty.list <- as.numeric(penalty.list)

  #### Check cv-folds
  lst_checkFR <- checkFoldRuns(Y, n_run, k_folds, fast_mode)
  n_run <- lst_checkFR$n_run
  fast_mode <- lst_checkFR$fast_mode

  #### Check rownames
  lst_check <- checkXY.rownames.mb(X, Y, verbose = verbose)
  X <- lst_check$X
  Y <- lst_check$Y

  #### Illegal chars in colnames
  X <- checkColnamesIllegalChars.mb(X)

  #### REQUIREMENTS
  checkX.colnames.mb(X)
  checkY.colnames(Y)
  lst_check <- checkXY.mb.class(X, Y, verbose = verbose)
  X <- lst_check$X
  Y <- lst_check$Y

  check.cv.weights(c(w_AIC, w_C.Index, w_I.BRIER, w_AUC))
  # if(!pred.method %in% c("risksetROC", "survivalROC", "cenROC", "nsROC", "smoothROCtime_C", "smoothROCtime_I")){
  #   stop_quietly(paste0("pred.method must be one of the following: ", paste0(c("risksetROC", "survivalROC", "cenROC", "nsROC", "smoothROCtime_C", "smoothROCtime_I"), collapse = ", ")))
  # }
  if(!pred.method %in% pkg.env$AUC_evaluators){
    stop_quietly(paste0("pred.method must be one of the following: ", paste0(pkg.env$AUC_evaluators, collapse = ", ")))
  }

  #### ZERO VARIANCE - ALWAYS
  if(!remove_variance_at_fold_level & (remove_near_zero_variance | remove_zero_variance)){
    lst_dnz <- deleteZeroOrNearZeroVariance.mb(X = X,
                                               remove_near_zero_variance = remove_near_zero_variance,
                                               remove_zero_variance = remove_zero_variance,
                                               toKeep.zv = toKeep.zv,
                                               freqCut = FREQ_CUT)
    X <- lst_dnz$X
    variablesDeleted <- lst_dnz$variablesDeleted
  }else{
    variablesDeleted <- NULL
  }

  #### COEF VARIATION
  if(!remove_variance_at_fold_level & (remove_near_zero_variance | remove_zero_variance)){
    lst_dnzc <- deleteNearZeroCoefficientOfVariation.mb(X = X)
    X <- lst_dnzc$X
    variablesDeleted_cvar <- lst_dnzc$variablesDeleted
  }else{
    variablesDeleted_cvar <- NULL
  }

  #### MAX PREDICTORS
  max.ncomp <- check.mb.ncomp(X, max.ncomp)
  max.ncomp <- check.mb.maxPredictors(X, Y, MIN_EPV, max.ncomp, verbose = verbose)
  if(MIN_COMP_TO_CHECK >= max.ncomp){
    MIN_COMP_TO_CHECK = max(max.ncomp-1, 1)
  }

  #### #
  # CV #
  #### #
  # lst_data <- splitData_Iterations_Folds.mb(X, Y, n_run = n_run, k_folds = k_folds, seed = seed) #FOR TEST
  # lst_X_train <- lst_data$lst_X_train
  # lst_Y_train <- lst_data$lst_Y_train
  # lst_X_test <- lst_data$lst_X_test
  # lst_Y_test <- lst_data$lst_Y_test
  # k_folds <- lst_data$k_folds
  #
  # lst_train_indexes <- lst_data$lst_train_index
  # lst_test_indexes <- lst_data$lst_test_index

  lst_data <- splitData_Iterations_Folds_indexes(Y, n_run = n_run, k_folds = k_folds, seed = seed) #FOR TEST

  lst_train_indexes <- lst_data$lst_train_index
  lst_test_indexes <- lst_data$lst_test_index

  #### ### ### ###
  # TRAIN MODELS #
  #### ### ### ###
  #total_models <- 1 * k_folds * n_run
  total_models <- max.ncomp * k_folds * n_run * length(penalty.list)

  lst_model <- get_Coxmos_models2.0(method = pkg.env$sb.splsicox,
                                    X_train = X, Y_train = Y,
                                    lst_X_train = lst_train_indexes, lst_Y_train = lst_train_indexes,
                                    max.ncomp = max.ncomp, penalty.list = penalty.list, EN.alpha.list = NULL, max.variables = NULL, vector = NULL,
                                    n_run = n_run, k_folds = k_folds,
                                    MIN_NVAR = NULL, MAX_NVAR = NULL, MIN_AUC_INCREASE = NULL, EVAL_METHOD = NULL,
                                    n.cut_points = NULL,
                                    x.center = x.center, x.scale = x.scale,
                                    y.center = y.center, y.scale = y.scale,
                                    remove_near_zero_variance = remove_variance_at_fold_level, remove_zero_variance = FALSE, toKeep.zv = NULL,
                                    alpha = alpha, MIN_EPV = MIN_EPV,
                                    remove_non_significant = remove_non_significant, tol = tol, max.iter = NULL,
                                    returnData = returnData, total_models = total_models,
                                    PARALLEL = PARALLEL, verbose = verbose)

  # already check in Coxmos_models
  # if(all(is.na(unlist(comp_model_lst)))){
  #   message(paste0("Best model could NOT be obtained. All models computed present problems. Try to remove variance at fold level. If problem persists, try to delete manually some problematic variables."))
  #
  #   t2 <- Sys.time()
  #   time <- difftime(t2,t1,units = "mins")
  #   if(return_models){
  #     return(cv.sb.splsicox_class(list(best_model_info = NULL, df_results_folds = NULL, df_results_runs = NULL, df_results_comps = NULL, lst_models = lst_model, pred.method = pred.method, opt.comp = NULL, opt.penalty = NULL, plot_AIC = NULL, plot_C.Index = NULL, plot_I.BRIER = NULL, plot_AUC = NULL, class = pkg.env$cv.sb.splsicox, lst_train_indexes = lst_train_indexes, lst_test_indexes = lst_test_indexes, time = time)))
  #   }else{
  #     return(cv.sb.splsicox_class(list(best_model_info = NULL, df_results_folds = NULL, df_results_runs = NULL, df_results_comps = NULL, lst_models = NULL, pred.method = pred.method, opt.comp = NULL, opt.penalty = NULL, plot_AIC = NULL, plot_C.Index = NULL, plot_I.BRIER = NULL, plot_AUC = NULL, class= pkg.env$cv.sb.splsicox, lst_train_indexes = lst_train_indexes, lst_test_indexes = lst_test_indexes, time = time)))
  #   }
  # }

  #### ### ### ### ### ### #
  # BEST MODEL FOR CV DATA #
  #### ### ### ### ### ### #
  total_models <- max.ncomp * k_folds * n_run * length(penalty.list)
  df_results_evals <- get_COX_evaluation_AIC_CINDEX(comp_model_lst = lst_model$comp_model_lst, alpha = alpha,
                                                    max.ncomp = max.ncomp, penalty.list = penalty.list, n_run = n_run, k_folds = k_folds,
                                                    total_models = total_models, remove_non_significant_models = remove_non_significant_models, verbose = verbose)

  if(all(is.null(df_results_evals))){
    message(paste0("Best model could NOT be obtained. All models computed present problems."))

    t2 <- Sys.time()
    time <- difftime(t2,t1,units = "mins")
    if(return_models){
      return(cv.sb.splsicox_class(list(best_model_info = NULL, df_results_folds = NULL, df_results_runs = NULL, df_results_comps = NULL, lst_models = lst_model, pred.method = pred.method, opt.comp = NULL, opt.penalty = NULL, plot_AIC = NULL, plot_C.Index = NULL, plot_I.BRIER = NULL, plot_AUC = NULL, class = pkg.env$cv.sb.splsicox, lst_train_indexes = lst_train_indexes, lst_test_indexes = lst_test_indexes, time = time)))
    }else{
      return(cv.sb.splsicox_class(list(best_model_info = NULL, df_results_folds = NULL, df_results_runs = NULL, df_results_comps = NULL, lst_models = NULL, pred.method = pred.method, opt.comp = NULL, opt.penalty = NULL, plot_AIC = NULL, plot_C.Index = NULL, plot_I.BRIER = NULL, plot_AUC = NULL, class= pkg.env$cv.sb.splsicox, lst_train_indexes = lst_train_indexes, lst_test_indexes = lst_test_indexes, time = time)))
    }
  }

  #### ### ### ### ### ### #
  # EVALUATING BRIER SCORE #
  #### ### ### ### ### ### #
  df_results_evals_comp <- NULL
  df_results_evals_run <- NULL
  df_results_evals_fold <- NULL
  optimal_comp_index <- NULL
  optimal_comp_flag <- FALSE
  optimal_eta_index <- NULL
  optimal_eta <- NULL

  if(TRUE){ #compute always BRIER SCORE
    #calculate time vector if still NULL
    if(is.null(times)){
      times <- getTimesVector(Y, max_time_points = max_time_points)
    }

    #As we are measuring just one evaluator and one method - PARALLEL = FALSE
    lst_df <- get_COX_evaluation_BRIER_sPLS(comp_model_lst = lst_model$comp_model_lst,
                                            fast_mode = fast_mode,
                                            X_test = X, Y_test = Y,
                                            lst_X_test = lst_test_indexes, lst_Y_test = lst_test_indexes,
                                            df_results_evals = df_results_evals, times = times,
                                            pred.method = pred.method, pred.attr = pred.attr,
                                            max.ncomp = max.ncomp, penalty.list = penalty.list, n_run = n_run, k_folds = k_folds,
                                            MIN_AUC_INCREASE = MIN_AUC_INCREASE, MIN_AUC = MIN_AUC, MIN_COMP_TO_CHECK = MIN_COMP_TO_CHECK,
                                            w_I.BRIER = w_I.BRIER, method.train = pkg.env$sb.splsicox, PARALLEL = FALSE, verbose = verbose)

    df_results_evals_comp <- lst_df$df_results_evals_comp
    df_results_evals_run <- lst_df$df_results_evals_run
    df_results_evals_fold <- lst_df$df_results_evals_fold
  }

  #### ### ### ### #
  # EVALUATING AUC #
  #### ### ### ### #

  if(w_AUC!=0){
    #total_models <- ifelse(!fast_mode, n_run * max.ncomp, k_folds * n_run * max.ncomp)#inside get_COX_evaluation_AUC

    #times should be the same for all folds
    #calculate time vector if still NULL
    if(is.null(times)){
      times <- getTimesVector(Y, max_time_points = max_time_points)
    }

    lst_df <- get_COX_evaluation_AUC_sPLS(comp_model_lst = lst_model$comp_model_lst,
                                          X_test = X, Y_test = Y,
                                          lst_X_test = lst_test_indexes, lst_Y_test = lst_test_indexes,
                                          df_results_evals = df_results_evals, times = times,
                                          fast_mode = fast_mode, pred.method = pred.method, pred.attr = pred.attr,
                                          max.ncomp = max.ncomp, penalty.list = penalty.list, n_run = n_run, k_folds = k_folds,
                                          MIN_AUC_INCREASE = MIN_AUC_INCREASE, MIN_AUC = MIN_AUC, MIN_COMP_TO_CHECK = MIN_COMP_TO_CHECK,
                                          w_AUC = w_AUC, method.train = pkg.env$sb.splsicox, PARALLEL = FALSE, verbose = verbose)

    if(is.null(df_results_evals_comp)){
      df_results_evals_comp <- lst_df$df_results_evals_comp
    }else{
      df_results_evals_comp$AUC <- lst_df$df_results_evals_comp$AUC
    }

    if(is.null(df_results_evals_run)){
      df_results_evals_run <- lst_df$df_results_evals_run
    }else{
      df_results_evals_run$AUC <- lst_df$df_results_evals_run$AUC
    }

    if(is.null(df_results_evals_fold)){
      df_results_evals_fold <- lst_df$df_results_evals_fold
    }else{
      df_results_evals_fold$AUC <- lst_df$df_results_evals_fold$AUC
    }

    optimal_comp_index <- lst_df$optimal_comp_index
    optimal_comp_flag <- lst_df$optimal_comp_flag
    optimal_eta <- lst_df$optimal_eta
    optimal_eta_index <- lst_df$optimal_eta_index
  }

  #### ### ### #
  # BEST MODEL #
  #### ### ### #

  df_results_evals_comp <- cv.getScoreFromWeight(lst_cox_mean = df_results_evals_comp, w_AIC = w_AIC, w_C.Index = w_C.Index, w_I.BRIER = w_I.BRIER, w_AUC = w_AUC,
                                                 colname_AIC = "AIC", colname_c_index = "C.Index", colname_AUC = "AUC", colname_BRIER = "IBS")

  if(optimal_comp_flag){
    best_model_info <- df_results_evals_comp[df_results_evals_comp[,"n.comps"]==optimal_comp_index,, drop = FALSE][1,]
    best_model_info <- as.data.frame(best_model_info)
  }else{
    best_model_info <- df_results_evals_comp[which(df_results_evals_comp[,"score"] == max(df_results_evals_comp[,"score"], na.rm = TRUE)),, drop = FALSE][1,]
    best_model_info <- as.data.frame(best_model_info)
  }

  #### ###
  # PLOT #
  #### ###
  class = pkg.env$sb.splsicox
  lst_EVAL_PLOTS <- get_EVAL_PLOTS(fast_mode = fast_mode, best_model_info = best_model_info, w_AUC = w_AUC, w_I.BRIER = w_I.BRIER, max.ncomp = max.ncomp, penalty.list = penalty.list,
                                   df_results_evals_fold = df_results_evals_fold, df_results_evals_run = df_results_evals_run, df_results_evals_comp = df_results_evals_comp,
                                   colname_AIC = "AIC", colname_c_index = "C.Index", colname_AUC = "AUC", colname_BRIER = "IBS", x.text = "Component",
                                   class = class)

  ggp_AUC <- lst_EVAL_PLOTS$ggp_AUC
  ggp_IBS <- lst_EVAL_PLOTS$ggp_IBS
  ggp_C.Index <- lst_EVAL_PLOTS$ggp_C.Index
  ggp_AIC <- lst_EVAL_PLOTS$ggp_AIC

  df_results_evals_comp <- lst_EVAL_PLOTS$df_results_evals_comp

  #### ### #
  # RETURN #
  #### ### #
  message(paste0("Best model obtained.\n"))

  t2 <- Sys.time()
  time <- difftime(t2,t1,units = "mins")

  # invisible(gc())
  if(return_models){
    return(cv.sb.splsicox_class(list(best_model_info = best_model_info,
                                     df_results_folds = df_results_evals_fold,
                                     df_results_runs = df_results_evals_run,
                                     df_results_comps = df_results_evals_comp,
                                     lst_models = lst_model,
                                     pred.method = pred.method,
                                     opt.comp = best_model_info$n.comps,
                                     opt.penalty = best_model_info$penalty,
                                     plot_AIC = ggp_AIC,
                                     plot_C.Index = ggp_C.Index,
                                     plot_I.BRIER = ggp_IBS,
                                     plot_AUC = ggp_AUC,
                                     class = pkg.env$cv.sb.splsicox,
                                     lst_train_indexes = lst_train_indexes,
                                     lst_test_indexes = lst_test_indexes,
                                     time = time)))
  }else{
    return(cv.sb.splsicox_class(list(best_model_info = best_model_info,
                                     df_results_folds = df_results_evals_fold,
                                     df_results_runs = df_results_evals_run,
                                     df_results_comps = df_results_evals_comp,
                                     lst_models = NULL, pred.method = pred.method,
                                     opt.comp = best_model_info$n.comps,
                                     opt.penalty = best_model_info$penalty,
                                     plot_AIC = ggp_AIC,
                                     plot_C.Index = ggp_C.Index,
                                     plot_I.BRIER = ggp_IBS,
                                     plot_AUC = ggp_AUC,
                                     class = pkg.env$cv.sb.splsicox,
                                     lst_train_indexes = lst_train_indexes,
                                     lst_test_indexes = lst_test_indexes,
                                     time = time)))
  }
}

### ## ##
# CLASS #
### ## ##

sb.splsicox_class = function(pls_model, ...) {
  model = structure(pls_model, class = pkg.env$model_class,
                    model = pkg.env$sb.splsicox)
  return(model)
}

cv.sb.splsicox_class = function(pls_model, ...) {
  model = structure(pls_model, class = pkg.env$model_class,
                    model = pkg.env$cv.sb.splsicox)
  return(model)
}