R/PRIMsrc.internal.r

Defines functions onAttach zeroslope is.wholenumber is.empty cbindlist list2mat list2array lapply.mat lapply.array cv.folds updatecut endpoints prsp cv.comb.peel cv.ave.peel cv.comb.box cv.ave.box cv.null cv.pval cv.type.box cv.box cv.spca cv.ppl coef.pcqreg predict.pcqreg pcqreg loss.pcqreg cv.pcqreg cv.pcqr cv.prsp.univ cv.prsp.tune cv.prsp cv.type.presel cv.presel

#===============================================================================================================================#
# PRIMsrc
#===============================================================================================================================#

#===============================================================================================================================#
# SURVIVAL INTERNAL SUBROUTINES
# (never to be called by end-user)
#===============================================================================================================================#

#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.presel(X,
#                             y,
#                             delta,
#                             B,
#                             K,
#                             vs,
#                             vstype,
#                             vsarg,
#                             vscons,
#                             cv,
#                             onese,
#                             parallel.vs,
#                             parallel.rep,
#                             clus.vs,
#                             clus.rep,
#                             verbose,
#                             seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.presel <- function(X,
                      y,
                      delta,
                      B,
                      K,
                      vs,
                      vstype,
                      vsarg,
                      vscons,
                      cv,
                      onese,
                      parallel.vs,
                      parallel.rep,
                      clus.vs,
                      clus.rep,
                      verbose,
                      seed) {
  
  

  # Parsing and evaluating 'vsarg' parameters
  eval(parse(text = unlist(strsplit(x = vsarg, split = ","))))

  # Treatment of variable screening
  if (vs) {
    
    cat("Screening of informative covariates ... \n", sep="")
    
    if ((vscons < 1/K) || (vscons > 1)) {
      stop("\nVariable screening conservativeness must be a real number in [1/K, 1]. Exiting ... \n\n")
    }
    
    replic <- 1:B
    if (!parallel.rep) {
      if (!is.null(seed)) {
        seed <- (0:(B-1)) + seed
      }
      CV.type.presel.obj <- cv.type.presel(X=X,
                                           y=y,
                                           delta=delta,
                                           replic=replic,
                                           K=K,
                                           vsarg=vsarg,
                                           vstype=vstype,
                                           vscons=vscons,
                                           cv=cv,
                                           onese=onese,
                                           parallel.vs=parallel.vs,
                                           parallel.rep=parallel.rep,
                                           clus.vs=clus.vs,
                                           verbose=verbose,
                                           seed=seed)
    } else {
      clusterSetRNGStream(cl=clus.rep, iseed=seed)
      obj.cl <- clusterApply(cl=clus.rep, x=replic, fun=cv.type.presel,
                             X=X,
                             y=y,
                             delta=delta,
                             K=K,
                             vsarg=vsarg,
                             vstype=vstype,
                             vscons=vscons,
                             cv=cv,
                             onese=onese,
                             parallel.vs=parallel.vs,
                             parallel.rep=parallel.rep,
                             clus.vs=clus.vs,
                             verbose=verbose,
                             seed=NULL)
      CV.type.presel.obj <- list("varsel"=vector(mode="list", length=B),
                                 "varsign"=vector(mode="list", length=B),
                                 "boxstat.mean"=vector(mode="list", length=B),
                                 "success"=logical(B))
      for (b in 1:B) {
        CV.type.presel.obj$varsel[[b]] <- obj.cl[[b]]$varsel
        CV.type.presel.obj$varsign[[b]] <- obj.cl[[b]]$varsign
        CV.type.presel.obj$boxstat.mean[[b]] <- obj.cl[[b]]$boxstat.mean
        CV.type.presel.obj$success[b] <- obj.cl[[b]]$success
      }
      CV.type.presel.obj$vsarg <- obj.cl[[1]]$vsarg
    }
    # Collect the results from all replicates
    vsarg <- CV.type.presel.obj$vsarg
    varsel.list <- CV.type.presel.obj$varsel
    varsign.list <- CV.type.presel.obj$varsign
    boxstat.profile.list <- CV.type.presel.obj$boxstat.mean
    success <-  any(CV.type.presel.obj$success)
    
    # Parsing 'vsarg' to retrieve 'cvcriterion', 'M', and 'msize'
    if (vstype == "prsp") {
      cvcriterion <- vsarg$cvcriterion
    }
    msize <- vsarg$msize
    M <- length(msize)

    # Get the fitted model
    if (!success) {
      # Failed to select any covariates
      varsel <- NA
      varsign <- NA
      if (vstype == "prsp") {
        boxstat.profile <- matrix(data=NA, nrow=B, ncol=M)
        colnames(boxstat.profile) <- paste("msize=", msize, sep="")
        boxstat.mean <- rep(NA, M)
        boxstat.se <- rep(NA, M)
        names(boxstat.mean) <- paste("msize=", msize, sep="")
        names(boxstat.se) <- paste("msize=", msize, sep="")
      } else {
        boxstat.profile <- NA
        boxstat.mean <- NA
        boxstat.se <- NA
      }
      boxstat.opt <- NA
      boxstat.1se <- NA
      success <- FALSE
    } else {
      # Get the averaged CV profiles of screening criterion for each PRSP model from all replicates
      if (vstype == "prsp") {
        # Get the averaged CV profiles
        boxstat.profile <- list2mat(list=boxstat.profile.list, fill=NA, coltrunc=M)
        colnames(boxstat.profile) <- paste("msize=", msize, sep="")
        boxstat.mean <- rep(NA, M)
        boxstat.se <- rep(NA, M)
        names(boxstat.mean) <- paste("msize=", msize, sep="")
        names(boxstat.se) <- paste("msize=", msize, sep="")
        if (cvcriterion == "lhr") {
          for (m in 1:M) {
            sumlhr <- rep(x=NA, times=B)
            for (b in 1:B) {
              sumlhr[b] <- boxstat.profile.list[[b]][[m]]
            }
            boxstat.mean[m] <- mean(sumlhr, na.rm=TRUE)
            boxstat.se[m] <- sd(sumlhr, na.rm=TRUE)
          }
          if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
            boxstat.opt <- NA
            boxstat.1se <- NA
          } else {
            boxstat.opt <- which.max(boxstat.mean)
            w <- boxstat.mean >= boxstat.mean[boxstat.opt]-boxstat.se[boxstat.opt]
            if (all(is.na(w)) || is.empty(w)) {
              boxstat.1se <- NA
            } else {
              boxstat.1se <- min(which(w))
            }
          }
        } else if (cvcriterion == "lrt") {
          for (m in 1:M) {
            sumlrt <- rep(x=NA, times=B)
            for (b in 1:B) {
              sumlrt[b] <- boxstat.profile.list[[b]][[m]]
            }
            boxstat.mean[m] <- mean(sumlrt, na.rm=TRUE)
            boxstat.se[m] <- sd(sumlrt, na.rm=TRUE)
          }
          if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
            boxstat.opt <- NA
            boxstat.1se <- NA
          } else {
            boxstat.opt <- which.max(boxstat.mean)
            w <- boxstat.mean >= boxstat.mean[boxstat.opt]-boxstat.se[boxstat.opt]
            if (all(is.na(w)) || is.empty(w)) {
              boxstat.1se <- NA
            } else {
              boxstat.1se <- min(which(w))
            }
          }
        } else if (cvcriterion == "cer") {
          for (m in 1:M) {
            sumcer <- rep(x=NA, times=B)
            for (b in 1:B) {
              sumcer[b] <- boxstat.profile.list[[b]][[m]]
            }
            boxstat.mean[m] <- mean(sumcer, na.rm=TRUE)
            boxstat.se[m] <- sd(sumcer, na.rm=TRUE)
          }
          if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
            boxstat.opt <- NA
            boxstat.1se <- NA
          } else {
            boxstat.opt <- which.min(boxstat.mean)
            w <- boxstat.mean <= boxstat.mean[boxstat.opt]+boxstat.se[boxstat.opt]
            if (all(is.na(w)) || is.empty(w)) {
              boxstat.1se <- NA
            } else {
              boxstat.1se <- min(which(w))
            }
          }
        } else {
          stop("Invalid CV criterion option. Exiting ... \n\n")
        }
      } else {
        boxstat.profile <- NA
        boxstat.mean <- NA
        boxstat.se <- NA
        boxstat.opt <- NA
        boxstat.1se <- NA
      }
      # Get the selected covariates with their signs from all replicates for the selected model
      sel.l <- unlist(varsel.list)
      sign.l <- unlist(varsign.list)
      na <- (is.na(sel.l))
      nna <- length(which(na))
      sel.l <- sel.l[!na]
      sign.l <- sign.l[!na]
      if ((is.empty(sel.l)) || (is.empty(sign.l)) || (all(sign.l == 0)))  {
        varsel <- NA
        varsign <- NA
      } else {
        vartab <- table(names(sel.l), useNA="no")
        nvotes <- B - nna
        w <- names(which(vartab >= ceiling(nvotes*vscons)))
        if ((is.na(w)) || (is.empty(w))) {
          varsel <- NA
          varsign <- NA
        } else{
          varsel <- pmatch(x=w, table=colnames(X), nomatch = NA_integer_, duplicates.ok = FALSE)
          names(varsel) <- w
          signtab <- table(names(sign.l), sign.l, useNA="no")
          if (length(unique(sign.l)) == 1) {
            signfreq <- rep(x=unique(sign.l), times=nrow(signtab))
            names(signfreq) <- rownames(signtab)
            varsign <- signfreq[w]
          } else if (all(sign.l != 1)) {
            signfreq <- apply(signtab[,c("-1","0"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 0
            varsign <- signfreq[w]
          } else if (all(sign.l != -1)) {
            signfreq <- apply(signtab[,c("0","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- 0
            signfreq[signfreq==2] <- 1
            varsign <- signfreq[w]
          } else {
            signfreq <- apply(signtab[,c("-1","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 1
            varsign <- signfreq[w]
          }
          ord <- order(as.numeric(varsel))
          varsel <- varsel[ord]
          varsign <- varsign[ord]
        }
      }
      if (all(is.na(varsel))) {
        varsel <- NA
        varsign <- NA
        if (vstype == "prsp") {
          boxstat.profile <- matrix(data=NA, nrow=B, ncol=M)
          colnames(boxstat.profile) <- paste("msize=", msize, sep="")
          boxstat.mean <- rep(NA, M)
          boxstat.se <- rep(NA, M)
          names(boxstat.mean) <- paste("msize=", msize, sep="")
          names(boxstat.se) <- paste("msize=", msize, sep="")
        } else {
          boxstat.profile <- NA
          boxstat.mean <- NA
          boxstat.se <- NA
        }
        boxstat.opt <- NA
        boxstat.1se <- NA
        success <- FALSE
      } else {
        success <- TRUE
      }
    }
    
  } else {
    
    cat("No screening of covariates. \n")
    
    if (!is.null(seed)) {
      set.seed(seed)
    }
    
    n <- nrow(X)
    p <- ncol(X)
    
    # Maximum Model Size
    msize <- p
    M <- length(msize)
    
    # Returning argument `vsarg` of variable selection parameters
    vsarg <- list(NA)
    
    # Determination of directions of peeling of all covariates
    if (cv) {
      cv.fit <- glmnet::cv.glmnet(x=X,
                                  y=survival::Surv(time=y, event=delta),
                                  alpha=0,
                                  nlambda=100,
                                  nfolds=pmax(3,K),
                                  family="cox",
                                  maxit=1e+05)
      fit <- glmnet::glmnet(x=X,
                            y=survival::Surv(time=y, event=delta),
                            alpha=0,
                            family="cox",
                            maxit=1e+05)
      if (onese) {
        cv.coef <- as.numeric(coef(fit, s=cv.fit$lambda.1se))
      } else {
        cv.coef <- as.numeric(coef(fit, s=cv.fit$lambda.min))
      }
    } else {
      fit <- glmnet::glmnet(x=X,
                            y=survival::Surv(time=y, event=delta),
                            alpha=0,
                            family="cox",
                            maxit=1e+05)
      cv.coef <- as.numeric(coef(fit, s=fit$lambda[50]))
    }
    w <- which(!(is.na(cv.coef)) & (cv.coef != 0))
    if (is.empty(w)) {
      varsel <- NA
      varsign <- NA
      success <- FALSE
    } else {
      varsel <- (1:ncol(X))[w]
      varsign <- sign(cv.coef[w])
      names(varsel) <- colnames(X)[w]
      names(varsign) <- colnames(X)[w]
      success <- TRUE
    }
    boxstat.profile <- matrix(data=NA, nrow=B, ncol=M)
    boxstat.mean <- rep(NA, M)
    boxstat.se <- rep(NA, M)
    boxstat.opt <- NA
    boxstat.1se <- NA
    
  }
  
  return(list("vsarg"=vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "boxstat.profile"=boxstat.profile,
              "boxstat.mean"=boxstat.mean,
              "boxstat.se"=boxstat.se,
              "boxstat.opt"=boxstat.opt,
              "boxstat.1se"=boxstat.1se,
              "success"=success,
              "seed"=seed))
  
}
#===============================================================================================================================#





#===============================================================================================================================#
#
#===============================================================================================================================#

cv.type.presel <- function(X,
                           y,
                           delta,
                           replic,
                           K,
                           vsarg,
                           vstype,
                           vscons,
                           cv,
                           onese,
                           parallel.vs,
                           parallel.rep,
                           clus.vs,
                           verbose,
                           seed) {
  
  if (!parallel.rep) {
    
    B <- length(replic)
    success <- logical(B)
    varsel <- vector(mode="list", length=B)
    varsign <- vector(mode="list", length=B)
    boxstat.mean <- vector(mode="list", length=B)
    b <- 1
    while (b <= B) {
      cat("Replicate : ", b, "\n", sep="")
      cat("seed : ", seed[b], "\n", sep="")
      if (vstype == "prsp") {
        CVSEL <- cv.prsp(X=X,
                         y=y,
                         delta=delta,
                         K=K,
                         vsarg=vsarg,
                         vscons=vscons,
                         cv=cv,
                         onese=onese,
                         parallel=parallel.vs,
                         clus=clus.vs,
                         verbose=verbose,
                         seed=seed[b])
        boxstat.mean[[b]] <- CVSEL$boxstat.mean
      } else if (vstype == "pcqr") {
        CVSEL <- cv.pcqr(X=X,
                         y=y,
                         delta=delta,
                         K=K,
                         vsarg=vsarg,
                         vscons=vscons,
                         cv=cv,
                         onese=onese,
                         parallel=parallel.vs,
                         clus=clus.vs,
                         verbose=verbose,
                         seed=seed[b])
        boxstat.mean[[b]] <- NULL
      } else if (vstype == "ppl") {
        CVSEL <- cv.ppl(X=X,
                        y=y,
                        delta=delta,
                        K=K,
                        vsarg=vsarg,
                        vscons=vscons,
                        cv=cv,
                        onese=onese,
                        parallel=parallel.vs,
                        clus=clus.vs,
                        verbose=verbose,
                        seed=seed[b])
        boxstat.mean[[b]] <- NULL
      } else if (vstype == "spca") {
        CVSEL <- cv.spca(X=X,
                         y=y,
                         delta=delta,
                         K=K,
                         vsarg=vsarg,
                         vscons=vscons,
                         cv=cv,
                         parallel=parallel.vs,
                         clus=clus.vs,
                         verbose=verbose,
                         seed=seed[b])
        boxstat.mean[[b]] <- NULL
      } else {
        stop("Invalid variable screening type option. Exiting ... \n\n")
      }
      varsel[[b]] <- CVSEL$varsel
      varsign[[b]] <- CVSEL$varsign
      success[b] <- CVSEL$success
      b <- b + 1
    }
    
  } else {
    
    if (vstype == "prsp") {
      CVSEL <- cv.prsp(X=X,
                       y=y,
                       delta=delta,
                       K=K,
                       vsarg=vsarg,
                       vscons=vscons,
                       cv=cv,
                       onese=onese,
                       parallel=parallel.vs,
                       clus=clus.vs,
                       verbose=verbose,
                       seed=NULL)
      boxstat.mean <- CVSEL$boxstat.mean
    } else if (vstype == "pcqr") {
      CVSEL <- cv.pcqr(X=X,
                       y=y,
                       delta=delta,
                       K=K,
                       vsarg=vsarg,
                       vscons=vscons,
                       cv=cv,
                       onese=onese,
                       parallel=parallel.vs,
                       clus=clus.vs,
                       verbose=verbose,
                       seed=NULL)
      boxstat.mean <- NULL
    } else if (vstype == "ppl") {
      CVSEL <- cv.ppl(X=X,
                      y=y,
                      delta=delta,
                      K=K,
                      vsarg=vsarg,
                      vscons=vscons,
                      cv=cv,
                      onese=onese,
                      parallel=parallel.vs,
                      clus=clus.vs,
                      verbose=verbose,
                      seed=NULL)
      boxstat.mean <- NULL
    } else if (vstype == "spca") {
      CVSEL <- cv.spca(X=X,
                       y=y,
                       delta=delta,
                       K=K,
                       vsarg=vsarg,
                       vscons=vscons,
                       cv=cv,
                       parallel=parallel.vs,
                       clus=clus.vs,
                       verbose=verbose,
                       seed=NULL)
      boxstat.mean <- NULL
    } else {
      stop("Invalid variable screening type option. Exiting ... \n\n")
    }
    varsel <- CVSEL$varsel
    varsign <- CVSEL$varsign
    success <- CVSEL$success
    
  }
  
  return(list("vsarg"=CVSEL$vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "boxstat.mean"=boxstat.mean,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
# Variables screening by univariate PRSP screening (PRSP)
# CV over number of ordered statistics used to compute variable frequencies
#===============================================================================================================================#

cv.prsp <- function(X,
                    y,
                    delta,
                    K,
                    vsarg,
                    vscons,
                    cv,
                    onese,
                    parallel,
                    clus,
                    verbose,
                    seed) {
  
  n <- nrow(X)
  p <- ncol(X)
  
  # Parsing and evaluating 'vsarg' argument to evaluate all parameters
  alpha <- NULL
  beta <- NULL
  msize <- NULL
  peelcriterion <- NULL
  cvcriterion <- NULL
  eval(parse( text=unlist(strsplit(x=vsarg, split=",")) ))
  
  if (is.null(msize)) {
    cv <- TRUE
  } else {
    cv <- FALSE
  }
  
  # Maximal possible number of peeling steps
  Lmax <- ceiling(log(beta) / log(1 - alpha))
  
  # Maximum Model Size or Vector of Model Sizes
  if (!is.null(msize)) {
    if (msize < 1) {
      cat("Warning: Parameter `msize` was less than the minimal allowed value and was reset to ", 1, ".\n", sep="")
    }
    msize <- max(1,floor(msize))
    Smax <- max(1,floor(p))
    if (msize > Smax) {
      cat("Warning: Parameter `msize` was greater than the maximal allowed vale and was reset to ", Smax, ".\n", sep="")
    }
    msize <- min(Smax, msize)
  } else {
    Smax <- max(1,floor(p/5))
    msize <- unique(ceiling(seq(from=max(1,floor(Smax/100)), to=Smax, length=min(msize,floor(100*Smax/p)))))
    cat("Model sizes to be explored by cross-validation: ", msize, "\n", sep=" ")
  }
  M <- length(msize)
  
  # Creating argument `cvarg` of `cv.prsp.tune` for variable selection parameters
  cvarg <- paste("alpha=", alpha,
                 ",beta=", beta,
                 ",L=", Lmax,
                 ",peelcriterion=\"", peelcriterion, "\"",
                 ",cvcriterion=\"", cvcriterion, "\"", sep="")
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  
  varsel.list <- vector(mode="list", length=K)
  varsign.list <- vector(mode="list", length=K)
  boxstat.list <- vector(mode="list", length=K)
  k <- 1
  while (k <= K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    # Initialize training and test data
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
    }
    # Training the model for each length of ranked statistics
    if (!parallel) {
      cv.obj <- cv.prsp.tune(traindata=traindata,
                             traintime=traintime,
                             trainstatus=trainstatus,
                             testdata=testdata,
                             testtime=testtime,
                             teststatus=teststatus,
                             cvarg=cvarg,
                             msize=msize,
                             parallel=parallel,
                             verbose=verbose)
    } else {
      clusterSetRNGStream(cl=clus, iseed=seed)
      obj.cl <- clusterApply(cl=clus, x=msize, fun=cv.prsp.tune,
                             traindata=traindata,
                             traintime=traintime,
                             trainstatus=trainstatus,
                             testdata=testdata,
                             testtime=testtime,
                             teststatus=teststatus,
                             cvarg=cvarg,
                             parallel=parallel,
                             verbose=verbose)
      cv.obj <- list("varsel"=vector(mode="list", length=M),
                     "varsign"=vector(mode="list", length=M),
                     "boxstat"=vector(mode="list", length=M),
                     "success"=TRUE)
      for (m in 1:M) {
        cv.obj$varsel[[m]] <- obj.cl[[m]]$varsel
        cv.obj$varsign[[m]] <- obj.cl[[m]]$varsign
        cv.obj$boxstat[[m]] <- obj.cl[[m]]$boxstat
        cv.obj$success <- (cv.obj$success & obj.cl[[m]]$success)
      }
    }
    if (cv.obj$success) {
      # Store the test set results for each model from each fold
      boxstat.list[[k]] <- cv.obj$boxstat
      varsel.list[[k]] <- cv.obj$varsel
      varsign.list[[k]] <- cv.obj$varsign
    } else {
      boxstat.list[[k]] <- as.list(rep(NA, M))
      varsel.list[[k]] <- as.list(rep(NA, M))
      varsign.list[[k]] <- as.list(rep(NA, M))
    }
    k <- k + 1
  }
  
  # Get the fitted model
  if (all(is.na(varsel.list))) {
    # Failed to select any covariates
    varsel <- NA
    varsign <- NA
    boxstat.mean <- rep(NA, M)
    boxstat.opt <- NA
    boxstat.1se <- NA
    msize <- NA
    success <- FALSE
  } else {
    # Get the averaged CV profiles of screening criterion for each model from all folds
    boxstat.mean <- rep(NA, M)
    boxstat.se <- rep(NA, M)
    if (cvcriterion == "lhr") {
      for (m in 1:M) {
        sumlhr <- rep(x=NA, times=K)
        for (k in 1:K) {
          sumlhr[k] <- boxstat.list[[k]][[m]]
        }
        boxstat.mean[m] <- mean(sumlhr, na.rm=TRUE)
        boxstat.se[m] <- sd(sumlhr, na.rm=TRUE) / sqrt(n/K)
      }
      if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
        boxstat.opt <- NA
        boxstat.1se <- NA
      } else {
        boxstat.opt <- which.max(boxstat.mean)
        w <- boxstat.mean >= boxstat.mean[boxstat.opt]-boxstat.se[boxstat.opt]
        if (all(is.na(w)) || is.empty(w)) {
          boxstat.1se <- NA
        } else {
          boxstat.1se <- min(which(w))
        }
      }
    } else if (cvcriterion == "lrt") {
      for (m in 1:M) {
        sumlrt <- rep(x=NA, times=K)
        for (k in 1:K) {
          sumlrt[k] <- boxstat.list[[k]][[m]]
        }
        boxstat.mean[m] <- mean(sumlrt, na.rm=TRUE)
        boxstat.se[m] <- sd(sumlrt, na.rm=TRUE) / sqrt(n/K)
      }
      if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
        boxstat.opt <- NA
        boxstat.1se <- NA
      } else {
        boxstat.opt <- which.max(boxstat.mean)
        w <- boxstat.mean >= boxstat.mean[boxstat.opt]-boxstat.se[boxstat.opt]
        if (all(is.na(w)) || is.empty(w)) {
          boxstat.1se <- NA
        } else {
          boxstat.1se <- min(which(w))
        }
      }
    } else if (cvcriterion == "cer") {
      for (m in 1:M) {
        sumcer <- rep(x=NA, times=K)
        for (k in 1:K) {
          sumcer[k] <- boxstat.list[[k]][[m]]
        }
        boxstat.mean[m] <- mean(sumcer, na.rm=TRUE)
        boxstat.se[m] <- sd(sumcer, na.rm=TRUE) / sqrt(n/K)
      }
      if (all(is.na(boxstat.mean)) || is.empty(boxstat.mean)) {
        boxstat.opt <- NA
        boxstat.1se <- NA
      } else {
        boxstat.opt <- which.min(boxstat.mean)
        w <- boxstat.mean <= boxstat.mean[boxstat.opt]+boxstat.se[boxstat.opt]
        if (all(is.na(w)) || is.empty(w)) {
          boxstat.1se <- NA
        } else {
          boxstat.1se <- min(which(w))
        }
      }
    } else {
      stop("Invalid CV criterion option. Exiting ... \n\n")
    }
    # Get the selected covariates with their signs from all folds for each model
    varsel <- vector(mode="list", length=M)
    varsign <- vector(mode="list", length=M)
    for (m in 1:M) {
      sel.m <- numeric(0)
      sign.m <- numeric(0)
      for (k in 1:K) {
        sel.m <- c(sel.m, varsel.list[[k]][[m]])
        sign.m <- c(sign.m, varsign.list[[k]][[m]])
      }
      na <- (is.na(sel.m))
      nna <- sum(na)
      sel.m <- sel.m[!na]
      sign.m <- sign.m[!na]
      if ((is.empty(sel.m)) || (is.empty(sign.m)) || (all(sign.m == 0)))  {
        varsel[[m]] <- NA
        varsign[[m]] <- NA
      } else {
        vartab <- table(names(sel.m), useNA="no")
        nvotes <- K - nna
        w <- names(which(vartab >= ceiling(nvotes*vscons)))
        if ((is.na(w)) || (is.empty(w))) {
          varsel[[m]] <- NA
          varsign[[m]] <- NA
        } else{
          varsel[[m]] <- pmatch(x=w, table=colnames(X), nomatch = NA_integer_, duplicates.ok = FALSE)
          names(varsel[[m]]) <- w
          signtab <- table(names(sign.m), sign.m, useNA="no")
          if (length(unique(sign.m)) == 1) {
            signfreq <- rep(x=unique(sign.m), times=nrow(signtab))
            names(signfreq) <- rownames(signtab)
            varsign[[m]] <- signfreq[w]
          } else if (all(sign.m != 1)) {
            signfreq <- apply(signtab[,c("-1","0"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 0
            varsign[[m]] <- signfreq[w]
          } else if (all(sign.m != -1)) {
            signfreq <- apply(signtab[,c("0","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- 0
            signfreq[signfreq==2] <- 1
            varsign[[m]] <- signfreq[w]
          } else {
            signfreq <- apply(signtab[,c("-1","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 1
            varsign[[m]] <- signfreq[w]
          }
          ord <- order(as.numeric(varsel[[m]]))
          varsel[[m]] <- varsel[[m]][ord]
          varsign[[m]] <- varsign[[m]][ord]
        }
      }
    }
    if (all(is.na(varsel))) {
      varsel <- NA
      varsign <- NA
      boxstat.mean <- rep(NA, M)
      boxstat.opt <- NA
      boxstat.1se <- NA
      success <- FALSE
      msize <- NA
    } else {
      if (onese) {
        varsel <- varsel[[boxstat.1se]]
        varsign <- varsign[[boxstat.1se]]
      } else {
        varsel <- varsel[[boxstat.opt]]
        varsign <- varsign[[boxstat.opt]]
      }
      success <- TRUE
    }
  }
  
  # Returning updated argument `vsarg` of variable selection parameters
  vsarg <- list("alpha"=alpha,
                "beta"=beta,
                "peelcriterion"=peelcriterion,
                "cvcriterion"=cvcriterion,
                "msize"=msize)
  
  return(list("vsarg"=vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "boxstat.mean"=boxstat.mean,
              "boxstat.opt"=boxstat.opt,
              "boxstat.1se"=boxstat.1se,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
#
#===============================================================================================================================#

cv.prsp.tune <- function(traindata,
                         traintime,
                         trainstatus,
                         testdata,
                         testtime,
                         teststatus,
                         cvarg,
                         msize,
                         parallel,
                         verbose) {
  
  # Parsing and evaluating 'arg' parameters to evaluate 'cvcriterion'
  alpha <- NULL
  beta <- NULL
  L <- NULL
  peelcriterion <- NULL
  cvcriterion <- NULL
  eval(parse( text=unlist(strsplit(x=cvarg, split=",")) ))
  
  if (!parallel) {
    
    M <- length(msize)
    varsel <- vector(mode="list", length=M)
    varsign <- vector(mode="list", length=M)
    boxstat <- vector(mode="list", length=M)
    peelobj <- cv.prsp.univ(traindata=traindata,
                            traintime=traintime,
                            trainstatus=trainstatus,
                            cvarg=cvarg)
    m <- 1
    while (m <= M) {
      if (verbose) cat("Model Size: ", msize[m], "\n")
      if (is.empty(peelobj$varsel)) {
        success <- FALSE
        varsel <- NULL
        varsign <- NULL
        boxstat <- NULL
        m <- M + 1
      } else {
        success <- TRUE
        # Retrieving screened variables and box from the univariate screening
        varsel[[m]] <- peelobj$varsel[1:msize[m]]
        varsign[[m]] <- peelobj$varsign[1:msize[m]]
        # Extract the rule and sign as one vector
        varcut <- peelobj$varcut[1:msize[m]]
        test.cut <- t(t(testdata[,varsel[[m]], drop=FALSE]) * varsign[[m]])
        test.ind <- t(t(test.cut) >= varcut * varsign[[m]])
        # Select which test observations are `TRUE` for all screened covariates
        pred <- (rowMeans(test.ind) == 1)
        test.ind1 <- 1*pred
        if ((sum(pred, na.rm=TRUE) != length(pred[!is.na(pred)])) && (sum(pred, na.rm=TRUE) != 0)) {
          surv.formula <- (survival::Surv(testtime, teststatus) ~ 1 + test.ind1)
          if (cvcriterion == "lhr") {
            coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1)
            boxstat[[m]] <- coxobj$coef
          } else if (cvcriterion == "lrt") {
            boxstat[[m]] <- survival::survdiff(surv.formula, rho=0)$chisq
          } else if (cvcriterion == "cer") {
            coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
            predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
            boxstat[[m]] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(testtime, teststatus))['C Index']
          } else {
            stop("Invalid CV criterion option. Exiting ... \n\n")
          }
        } else {
          if (cvcriterion == "lhr") {
            boxstat[[m]] <- 0
          } else if (cvcriterion == "lrt") {
            boxstat[[m]] <- 0
          } else if (cvcriterion == "cer") {
            boxstat[[m]] <- 1
          } else {
            stop("Invalid CV criterion option. Exiting ... \n\n")
          }
        }
        m <- m + 1
      }
    }
    
  } else {
    
    peelobj <- cv.prsp.univ(traindata=traindata,
                            traintime=traintime,
                            trainstatus=trainstatus,
                            cvarg=cvarg)
    if (is.empty(peelobj$varsel)) {
      success <- FALSE
      varsel <- NULL
      varsign <- NULL
      boxstat <- NULL
    } else {
      # Retrieving screened variables and box from the univariate screening
      success <- TRUE
      varsel <- peelobj$varsel[1:msize]
      varsign <- peelobj$varsign[1:msize]
      varcut <- peelobj$varcut[1:msize]
      # Extract the rule and sign as one vector
      test.cut <- t(t(testdata[,varsel, drop=FALSE]) * varsign)
      test.ind <- t(t(test.cut) >= varcut * varsign)
      # Set as TRUE which observations are TRUE for all covariates
      pred <- (rowMeans(test.ind) == 1)
      test.ind1 <- 1*pred
      if ((sum(pred, na.rm=TRUE) != length(pred[!is.na(pred)])) && (sum(pred, na.rm=TRUE) != 0)) {
        surv.formula <- (survival::Surv(testtime, teststatus) ~ 1 + test.ind1)
        if (cvcriterion == "lhr") {
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1)
          boxstat <- coxobj$coef
        } else if (cvcriterion == "lrt") {
          boxstat <- survival::survdiff(surv.formula, rho=0)$chisq
        } else if (cvcriterion == "cer") {
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
          predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
          boxstat <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(testtime, teststatus))['C Index']
        } else {
          stop("Invalid CV criterion option. Exiting ... \n\n")
        }
      } else {
        if (cvcriterion == "lhr") {
          boxstat <- 0
        } else if (cvcriterion == "lrt") {
          boxstat <- 0
        } else if (cvcriterion == "cer") {
          boxstat <- 1
        } else {
          stop("Invalid CV criterion option. Exiting ... \n\n")
        }
      }
    }
  }
  
  return(list("varsel"=varsel,
              "varsign"=varsign,
              "boxstat"=boxstat,
              "success"=success))
}
#===============================================================================================================================#



#===============================================================================================================================#
#
#===============================================================================================================================#

cv.prsp.univ <- function(traindata,
                         traintime,
                         trainstatus,
                         cvarg) {
  
  # Parsing and evaluating 'arg' to evaluate all parameters
  alpha <- NULL
  beta <- NULL
  L <- NULL
  peelcriterion <- NULL
  cvcriterion <- NULL
  eval(parse( text=unlist(strsplit(x=cvarg, split=",")) ))
  
  # Creating argument `arg` of `cv.comb.peel` for PRSP parameters
  arg <- paste("alpha=", alpha,
               ",beta=", beta,
               ",L=", L,
               ",peelcriterion=\"", peelcriterion, "\"", sep="")
  
  # Univariate screening
  p <- ncol(traindata)
  varstat <- numeric(p)
  varsel <- 1:p
  varsign <- numeric(p)
  varcut <- numeric(p)
  names(varsel) <- colnames(traindata)
  names(varsign) <- colnames(traindata)
  names(varcut) <- colnames(traindata)
  for (j in 1:p) {
    prsp.fit <- prsp(traindata=traindata[, j, drop=FALSE],
                     traintime=traintime,
                     trainstatus=trainstatus,
                     varsign=NULL,
                     initcutpts=NULL,
                     arg=arg,
                     traingroups=NULL)
    cat("covariate: ", j, "\t (maxstep: ", prsp.fit$maxsteps, ")\n", sep="")
    varcut[j] <- prsp.fit$boxcut[1,1,drop=TRUE]
    varsign[j] <- prsp.fit$varsign
    Lm <- prsp.fit$maxsteps
    l <- 1
    boxstat <- numeric(Lm)
    while (l <= Lm) {
      # Extract the rule and sign as one vector
      boxcut <- prsp.fit$boxcut[l,1,drop=TRUE]
      train.cut <- t(t(traindata[,j,drop=FALSE]) * varsign[j])
      train.ind <- t(t(train.cut) >= boxcut * varsign[j])
      # Select which test observations are `TRUE`
      pred <- (rowMeans(train.ind) == 1)
      train.ind1 <- 1*pred
      if ((sum(pred, na.rm=TRUE) != length(pred[!is.na(pred)])) && (sum(pred, na.rm=TRUE) != 0)) {
        surv.formula <- (survival::Surv(traintime, trainstatus) ~ 1 + train.ind1)
        if (cvcriterion == "lhr") {
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1)
          boxstat[l] <- coxobj$coef
        } else if (cvcriterion == "lrt") {
          boxstat[l] <- survival::survdiff(surv.formula, rho=0)$chisq
        } else if (cvcriterion == "cer") {
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
          predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
          boxstat[l] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(traintime, trainstatus))['C Index']
        } else {
          stop("Invalid CV criterion option. Exiting ... \n\n")
        }
      } else {
        if (cvcriterion == "lhr") {
          boxstat[l] <- 0
        } else if (cvcriterion == "lrt") {
          boxstat[l] <- 0
        } else if (cvcriterion == "cer") {
          boxstat[l] <- 1
        } else {
          stop("Invalid CV criterion option. Exiting ... \n\n")
        }
      }
      l <- l + 1
    }
    if (cvcriterion == "lhr") {
      if (all(is.na(boxstat)) || is.empty(boxstat)) {
        varstat[j] <- NA
      } else {
        boxstat.opt <- which.max(boxstat)
        varstat[j] <- boxstat[boxstat.opt]
      }
    } else if (cvcriterion == "lrt") {
      if (all(is.na(boxstat)) || is.empty(boxstat)) {
        varstat[j] <- NA
      } else {
        boxstat.opt <- which.max(boxstat)
        varstat[j] <- boxstat[boxstat.opt]
      }
    } else if (cvcriterion == "cer") {
      if (all(is.na(boxstat)) || is.empty(boxstat)) {
        varstat[j] <- NA
      } else {
        boxstat.opt <- which.min(boxstat)
        varstat[j] <- boxstat[boxstat.opt]
      }
    } else {
      stop("Invalid CV criterion option. Exiting ... \n\n")
    }
  }
  
  # Order by variable statistics
  if (cvcriterion == "lhr") {
    ord <- order(varstat, decreasing=TRUE, na.last=TRUE)
  } else if (cvcriterion == "lrt") {
    ord <- order(varstat, decreasing=TRUE, na.last=TRUE)
  } else if (cvcriterion == "cer") {
    ord <- order(varstat, decreasing=FALSE, na.last=TRUE)
  } else {
    stop("Invalid CV criterion option. Exiting ... \n\n")
  }
  varsel <- varsel[ord]
  varsign <- varsign[ord]
  varcut <- varcut[ord]
  
  return(list("varsel"=varsel,
              "varsign"=varsign,
              "varcut"=varcut))
}
#===============================================================================================================================#



#===============================================================================================================================#
# Variables screening by Penalized Censored Quantile Regression model (PCQR)
# CV of variable selection using full cross-validation of both parameters alpha and lambda
#===============================================================================================================================#

cv.pcqr <- function(X,
                    y,
                    delta,
                    K,
                    vsarg,
                    vscons,
                    cv,
                    onese,
                    parallel,
                    clus,
                    verbose,
                    seed) {
  
  # Parsing and evaluating 'vsarg' argument to evaluate all parameters
  tau <- NULL
  alpha <- NULL
  nalpha <- NULL
  nlambda <- NULL
  eval(parse( text=unlist(strsplit(x=vsarg, split=",")) ))
  
  M <- 2
  msize <- numeric(M)
  
  if (is.null(alpha)) {
    enalpha <- seq(from=0, to=1, length.out=nalpha)
  } else {
    nalpha <- 1
    enalpha <- alpha
  }
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  
  varsel.list <- vector(mode="list", length=K)
  varsign.list <- vector(mode="list", length=K)
  k <- 1
  while (k <= K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    # Initialize training and test data
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
    }
    enlambda <- vector(mode="list", length=nalpha)
    cv.errmu <- vector(mode="list", length=nalpha)
    cv.errse <- vector(mode="list", length=nalpha)
    for (i in 1:nalpha) {
      cv.fit <- cv.pcqreg(X=traindata,
                          y=traintime,
                          delta=trainstatus,
                          method="cquantile",
                          gamma=IQR(traintime)/10,
                          tau=tau,
                          nfolds=pmax(3,K),
                          type.loss="deviance",
                          alpha=enalpha[i],
                          nlambda=nlambda,
                          lambda.min=ifelse(nrow(X)>ncol(X), 0.001, 0.05),
                          preprocess="standardize",
                          screen="ASR",
                          max.iter=1e4,
                          eps=1e-7,
                          dfmax=ncol(X)+1,
                          penalty.factor=rep(x=1, times=ncol(X)),
                          parallel=parallel,
                          verbose=verbose,
                          seed=seed)
      cv.errmu[[i]] <- cv.fit$cve
      cv.errse[[i]] <- cv.fit$cvse
      enlambda[[i]] <- cv.fit$lambda
    }
    cv.errmu <- list2mat(list=cv.errmu, coltrunc="max", fill=NA)
    cv.errse <- list2mat(list=cv.errse, coltrunc="max", fill=NA)
    cv.errmu.index <- as.matrix(which(x=(cv.errmu == as.numeric(cv.errmu)[which.min(cv.errmu)]), arr.ind=TRUE, useNames=FALSE))
    w <- which.min(cv.errmu.index[, 1, drop=TRUE])
    index.min <- as.numeric(cv.errmu.index[w,,drop=TRUE])
    ww <- as.matrix(which(x=cv.errmu < cv.errmu[index.min[1],index.min[2]] + cv.errse[index.min[1],index.min[2]], arr.ind=TRUE, useNames=TRUE))
    index.1se <- c(min(ww[,1]), min(ww[ww[,1]==min(ww[,1]),2]))
    if (is.empty(index.min) || is.empty(index.1se)) {
      varsel.list[[k]] <- list(NA, NA)
      varsign.list[[k]] <- list(NA, NA)
    } else {
      enalpha.min <- enalpha[index.min[1]]
      enalpha.1se <- enalpha[index.1se[1]]
      enlambda.min <- enlambda[[index.min[1]]][index.min[2]]
      enlambda.1se <- enlambda[[index.1se[1]]][index.1se[2]]
      fit <- pcqreg(X=traindata,
                    y=traintime,
                    delta=trainstatus,
                    method="cquantile",
                    gamma=IQR(traintime)/10,
                    tau=tau,
                    alpha=ifelse(test=onese, yes=enalpha.1se, no=enalpha.min),
                    nlambda=nlambda,
                    lambda.min=ifelse(nrow(X)>ncol(X), 0.001, 0.05),
                    lambda=c(enlambda.1se, enlambda.min),
                    preprocess="standardize",
                    screen="ASR",
                    max.iter=1e4,
                    eps=1e-7,
                    dfmax=ncol(X)+1,
                    penalty.factor=rep(x=1, times=ncol(X)),
                    verbose=verbose)
      w <- apply(X=fit$beta, MARGIN=2, FUN=function(x) {sum(!(is.na(x)) & (x != 0))})
      if (all(w == 0)) {
        varsel.list[[k]] <- list(NA, NA)
        varsign.list[[k]] <- list(NA, NA)
      } else {
        if (length(fit$lambda) <= 2) {
          cv.coef <- -coef.pcqreg(fit, lambda=c(enlambda.1se, enlambda.min), exact=TRUE)[-1,]
          cv.coef <- list(cv.coef[,1], cv.coef[,2])
          varsel <- lapply(cv.coef, FUN=function(x) {which(!(is.na(x)) & (x != 0))})
          if (all(is.na(varsel))) {
            varsel.list[[k]] <- list(NA, NA)
            varsign.list[[k]] <- list(NA, NA)
          } else {
            varsel.list[[k]] <- varsel
            varsign.list[[k]] <- lapply(cv.coef, function(x) sign(x))
          }
        } else {
          cv.coef <- -coef.pcqreg(fit, lambda=c(enlambda.1se, enlambda.min), exact=FALSE)[-1,]
          cv.coef <- list(cv.coef[,1], cv.coef[,2])
          varsel <- lapply(cv.coef, FUN=function(x) {which(!(is.na(x)) & (x != 0))})
          if (all(is.na(varsel))) {
            varsel.list[[k]] <- list(NA, NA)
            varsign.list[[k]] <- list(NA, NA)
          } else {
            varsel.list[[k]] <- varsel
            varsign.list[[k]] <- lapply(cv.coef, function(x) sign(x))
          }
        }
      }
    }
    k <- k + 1
  }
  
  # Get the fitted model
  if (all(is.na(varsel.list))) {
    # Failed to select any covariates
    varsel <- NA
    varsign <- NA
    enalpha.min <- NA
    enalpha.1se <- NA
    enlambda.min <- NA
    enlambda.1se <- NA
    msize <- rep(NA, M)
    success <- FALSE
  } else {
    # Get the selected covariates with their signs from all folds for each optimal lambda value
    varsel <- vector(mode="list", length=M)
    varsign <- vector(mode="list", length=M)
    for (l in 1:M) {
      sel.l <- numeric(0)
      sign.l <- numeric(0)
      for (k in 1:K) {
        sel.l <- c(sel.l, varsel.list[[k]][[l]])
        sign.l <- c(sign.l, varsign.list[[k]][[l]])
      }
      na <- (is.na(sel.l))
      nna <- length(which(na))
      sel.l <- sel.l[!na]
      sign.l <- sign.l[!na]
      if ((is.empty(sel.l)) || (is.empty(sign.l)) || (all(sign.l == 0)))  {
        varsel[[l]] <- NA
        varsign[[l]] <- NA
      } else {
        vartab <- table(names(sel.l), useNA="no")
        nvotes <- K - nna
        w <- names(which(vartab >= ceiling(nvotes*vscons)))
        if ((is.na(w)) || (is.empty(w))) {
          varsel[[l]] <- NA
          varsign[[l]] <- NA
        } else {
          varsel[[l]] <- pmatch(x=w, table=colnames(X), nomatch = NA_integer_, duplicates.ok = FALSE)
          names(varsel[[l]]) <- w
          signtab <- table(names(sign.l), sign.l, useNA="no")
          if (length(unique(sign.l)) == 1) {
            signfreq <- rep(x=unique(sign.l), times=nrow(signtab))
            names(signfreq) <- rownames(signtab)
            varsign[[l]] <- signfreq[w]
          } else if (all(sign.l != 1)) {
            signfreq <- apply(signtab[,c("-1","0"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 0
            varsign[[l]] <- signfreq[w]
          } else if (all(sign.l != -1)) {
            signfreq <- apply(signtab[,c("0","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- 0
            signfreq[signfreq==2] <- 1
            varsign[[l]] <- signfreq[w]
          } else {
            signfreq <- apply(signtab[,c("-1","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 1
            varsign[[l]] <- signfreq[w]
          }
          ord <- order(as.numeric(varsel[[l]]))
          varsel[[l]] <- varsel[[l]][ord]
          varsign[[l]] <- varsign[[l]][ord]
        }
      }
    }
    if (all(is.na(varsel))) {
      varsel <- NA
      varsign <- NA
      enalpha.min <- NA
      enalpha.1se <- NA
      enlambda.min <- NA
      enlambda.1se <- NA
      msize <- rep(NA, M)
      success <- FALSE
    } else {
      if (onese) {
        varsel <- varsel[[1]]
        varsign <- varsign[[1]]
        msize[1] <- length(varsel)
      } else {
        varsel <- varsel[[2]]
        varsign <- varsign[[2]]
        msize[2] <- length(varsel)
      }
      success <- TRUE
    }
  }
  
  # Returning updated argument `vsarg` of variable selection parameters
  vsarg <- list("tau"=tau,
                "alpha"=alpha,
                "nalpha"=nalpha,
                "nlambda"=nlambda,
                "msize"=msize)
  
  return(list("vsarg"=vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "enalpha.min"=enalpha.min,
              "enalpha.1se"=enalpha.1se,
              "enlambda.min"=enlambda.min,
              "enlambda.1se"=enlambda.1se,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#



#===============================================================================================================================#
#
#===============================================================================================================================#

cv.pcqreg <- function(X,
                      y,
                      delta,
                      method,
                      gamma,
                      tau,
                      nfolds,
                      type.loss,
                      alpha,
                      nlambda,
                      lambda.min,
                      preprocess,
                      screen,
                      max.iter,
                      eps,
                      dfmax,
                      penalty.factor,
                      parallel,
                      verbose,
                      seed) {
  
  cvf <- function(i, XX, y, delta, folds, cv.args, loss.args) {
    
    cv.args$X <- XX[folds$perm[(folds$which != i)], , drop=FALSE]
    cv.args$y <- y[folds$perm[(folds$which != i)]]
    cv.args$delta <- delta[folds$perm[(folds$which != i)]]
    X2 <- XX[folds$perm[(folds$which == i)], , drop=FALSE]
    y2 <- y[folds$perm[(folds$which == i)]]
    fit.i <- do.call("pcqreg", cv.args) # ensures the cross validation uses the same gamma for huber loss
    yhat <- matrix(data=predict.pcqreg(object=fit.i, X=X2, type="response", exact=FALSE),
                   nrow=length(y2))
    
    return(list(pe=loss.pcqreg(y=y2, yhat=yhat, args=loss.args),
                nl=length(fit.i$lambda)))
  }
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  n <- length(y)
  
  fit <- pcqreg(X=X,
                y=y,
                delta=delta,
                method=method,
                gamma=gamma,
                tau=tau,
                alpha=alpha,
                nlambda=nlambda,
                lambda.min=lambda.min,
                preprocess=preprocess,
                screen=screen,
                max.iter=max.iter,
                eps=eps,
                dfmax=dfmax,
                penalty.factor=penalty.factor,
                verbose=verbose)
  
  cv.args <- list("method"=method,
                  "gamma"=gamma,
                  "tau"=tau,
                  "alpha"=alpha,
                  "lambda"=fit$lambda,
                  "nlambda"=nlambda,
                  "lambda.min"=lambda.min,
                  "preprocess"=preprocess,
                  "screen"=screen,
                  "max.iter"=max.iter,
                  "eps"=eps,
                  "dfmax"=dfmax,
                  "penalty.factor"=penalty.factor,
                  "verbose"=verbose)
  
  loss.args <- list("method"=method,
                    "gamma"=gamma,
                    "tau"=tau,
                    "type.loss"=type.loss)
  
  if (is.null(delta)) {
    folds <- cv.folds(y=y, K=nfolds, seed=seed)
  } else {
    folds <- cv.folds(y=delta, K=nfolds, seed=seed)
  }
  
  E <- matrix(NA, nrow=n, ncol=length(cv.args$lambda))
  if (parallel) {
    cluster <- makeCluster(detectCores())
    clusterSetRNGStream(cl=cluster, iseed=seed)
    fold.results <- parLapply(cl=cluster,
                              X=1:nfolds,
                              fun=cvf,
                              XX=X,
                              y=y,
                              delta=delta,
                              folds=folds,
                              cv.args=cv.args,
                              loss.args=loss.args)
    stopCluster(cluster)
    for (i in 1:nfolds) {
      fit.i <- fold.results[[i]]
      E[folds$perm[(folds$which == i)], 1:fit.i$nl] <- fit.i$pe
    }
  } else {
    if (!is.null(seed)) {
      set.seed(seed)
    }
    for (i in 1:nfolds) {
      fit.i <- cvf(i=i,
                   XX=X,
                   y=y,
                   delta=delta,
                   folds=folds,
                   cv.args=cv.args,
                   loss.args=loss.args)
      E[folds$perm[(folds$which == i)], 1:fit.i$nl] <- fit.i$pe
    }
  }
  
  ## Eliminate saturated lambda values
  ind <- which(apply(is.finite(E), 2, all))
  E <- E[,ind]
  lambda <- cv.args$lambda[ind]
  
  ## Results
  cve <- apply(E, 2, mean)
  cvse <- apply(E, 2, sd) / sqrt(n/nfolds)
  index.min <- which.min(cve)
  
  # adjust the selection using 1-SD method
  index.1se <- min(which(cve < cve[index.min]+cvse[index.min]))
  
  # Output
  return(list(cve=cve,
              cvse=cvse,
              type.loss=type.loss,
              lambda=lambda,
              fit=fit,
              lambda.1se=lambda[index.1se],
              lambda.min=lambda[index.min]))
}
#===============================================================================================================================#




#===============================================================================================================================#
#
#===============================================================================================================================#

loss.pcqreg <- function(y,
                        yhat,
                        args) {
  
  hloss <- function(r, gamma) {
    rr <- abs(r)
    ifelse(rr <= gamma, rr^2/(2*gamma), rr-gamma/2)
  }
  
  qloss <- function(r, tau) {
    ifelse(r <= 0, (tau-1)*r, tau*r)
  }
  
  r <- y - yhat
  type.loss <- args$type.loss
  if (type.loss == "deviance") {
    method <- args$method
    if (method == "huber") {
      gamma <- args$gamma
      val <- hloss(r, gamma)
    } else if ((method == "quantile") || (method == "cquantile")) {
      tau <- args$tau
      val <- qloss(r, tau)
    } else {
      val <- r^2
    }
  } else if (type.loss == "mse") {
    val <- r^2
  } else {
    val <- abs(r)
  }
  return(val)
}
#===============================================================================================================================#




#===============================================================================================================================#
#
#===============================================================================================================================#

pcqreg <- function (X,
                    y,
                    delta,
                    method,
                    gamma,
                    tau,
                    alpha,
                    nlambda,
                    lambda.min,
                    lambda,
                    preprocess,
                    screen,
                    max.iter,
                    eps,
                    dfmax,
                    penalty.factor,
                    verbose) {
  
  # Error checking
  if (missing(lambda) && nlambda < 2)
    stop("nlambda should be at least 2. Exiting ... \n\n")
  if (alpha < 0 || alpha > 1)
    stop("alpha should be between 0 and 1. Exiting ... \n\n")
  if (method == "huber" && !missing(gamma) && gamma <= 0)
    stop("gamma should be positive for Huber loss. Exiting ... \n\n")
  if ((method == "quantile" || method == "cquantile")  && (tau < 0 || tau > 1))
    stop("tau should be between 0 and 1 for quantile loss. Exiting ... \n\n")
  if (length(penalty.factor) != ncol(X))
    stop("the length of penalty.factor should equal the number of columns of X. Exiting ... \n\n")
  
  # Flag for user-supplied lambda
  if (missing(lambda)) {
    lambda <- double(nlambda)
    user <- 0
  } else {
    nlambda <- length(lambda)
    user <- 1
  }
  
  # Saving function call for output
  call <- match.call()
  
  # Include a column for intercept
  n <- nrow(X)
  XX <- cbind(rep(1,n), X)
  penalty.factor <- c(0, penalty.factor) # no penalty for intercept term
  p <- ncol(XX)
  
  shift <- 0
  if (method == "huber") {
    shift <- if(gamma > sd(y))
      mean(y)
    else
      median(y)
  } else if (method == "ls") {
    shift <- mean(y)
  } else if (method == "quantile") {
    shift <- quantile(y, tau)
  } else if (method == "cquantile") {
    obj <- tryCatch({coef(object=quantreg::crq(formula=survival::Surv(time=y, event=delta, type="right") ~ 1, taus=tau, method="Portnoy"),taus=tau)},
                    error=function(){NULL})
    if ((is.na(obj)) || (is.empty(obj))) {
      shift <- quantile(y, tau)
    } else {
      shift <- coef(object=quantreg::crq(formula=survival::Surv(time=y, event=delta, type="right") ~ 1, taus=tau, method="Portnoy"),taus=tau)
    }
  }
  yy <- y - shift
  
  # Flags for preprocessing and screening
  ppflag = switch(preprocess, standardize = 1L, rescale = 2L)
  scrflag = switch(screen, ASR = 1L, SR = 2L, none = 0L)
  
  # Fitting
  if (alpha > 0) {
    if (method == "huber") {
      fit <- .C("C_huber", double(p*nlambda), integer(nlambda), as.double(lambda), integer(1), integer(1), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(gamma), as.double(alpha), as.double(eps), as.double(lambda.min), as.integer(nlambda),
                as.integer(n), as.integer(p), as.integer(ppflag), as.integer(scrflag), 1L, as.integer(dfmax), as.integer(max.iter), as.integer(user), as.integer(verbose))
    } else if ((method == "quantile") || (method == "cquantile")) {
      fit <- .C("C_quantile", double(p*nlambda), integer(nlambda), as.double(lambda), integer(1), integer(1), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(tau), as.double(alpha), as.double(eps), as.double(lambda.min), as.integer(nlambda),
                as.integer(n), as.integer(p), as.integer(ppflag), as.integer(scrflag), 1L, as.integer(dfmax), as.integer(max.iter), as.integer(user), as.integer(verbose))
    } else {
      fit <- .C("C_squared", double(p*nlambda), integer(nlambda), as.double(lambda), integer(1), integer(1), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(alpha), as.double(eps), as.double(lambda.min), as.integer(nlambda),
                as.integer(n), as.integer(p), as.integer(ppflag), as.integer(scrflag), 1L, as.integer(dfmax), as.integer(max.iter), as.integer(user), as.integer(verbose))
    }
    beta <- matrix(fit[[1]],nrow = p)
    iter <- fit[[2]]
    lambda <- fit[[3]]
    saturated <- fit[[4]]
    nv <- fit[[5]]
    # Eliminate saturated lambda values
    ind <- !is.na(iter)
    beta <- beta[, ind]
    iter <- iter[ind]
    lambda <- lambda[ind]
  } else {
    if (method == "huber") {
      fit <- .C("C_huber_l2", double(p*nlambda), integer(nlambda), as.double(lambda), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(gamma), as.double(eps), as.double(lambda.min), as.integer(nlambda),
                as.integer(n), as.integer(p), as.integer(ppflag), 1L, as.integer(max.iter), as.integer(user), as.integer(verbose))
    } else if ((method == "quantile") || (method == "cquantile")) {
      fit <- .C("C_quantile_l2", double(p*nlambda), integer(nlambda), as.double(lambda), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(tau), as.double(eps), as.double(lambda.min), as.integer(nlambda),
                as.integer(n), as.integer(p), as.integer(ppflag), 1L, as.integer(max.iter), as.integer(user), as.integer(verbose))
    } else {
      fit <- .C("C_squared_l2", double(p*nlambda), integer(nlambda), as.double(lambda), as.double(XX), as.double(yy),
                as.double(penalty.factor), as.double(eps), as.double(lambda.min), as.integer(nlambda), as.integer(n),
                as.integer(p), as.integer(ppflag), 1L, as.integer(max.iter), as.integer(user), as.integer(verbose))
    }
    beta <- matrix(fit[[1]],nrow = p)
    iter <- fit[[2]]
    lambda <- fit[[3]]
    saturated <- 0
    nv <- 0
  }
  
  # Intercept
  beta[1,] <- beta[1,] + shift
  
  # Names
  vnames <- colnames(X)
  if (is.null(vnames)) vnames=paste0("V", seq(p-1))
  vnames <- c("(Intercept)", vnames)
  dimnames(beta) <- list(vnames, paste0("L", 1:length(lambda)))
  
  # Output
  return(list(call=call,
              beta=beta,
              iter=iter,
              saturated=saturated,
              lambda=lambda,
              alpha=alpha,
              gamma=if (method == "huber") gamma else NULL,
              tau=if ((method == "quantile") || (method == "cquantile")) tau else NULL,
              penalty.factor=penalty.factor[-1],
              method=method,
              nv=nv))
}
#===============================================================================================================================#




#===============================================================================================================================#
#
#===============================================================================================================================#

predict.pcqreg <- function(object,
                           X,
                           lambda,
                           type=c("response","coefficients","nvars"),
                           exact) {
  
  type <- match.arg(type)
  
  if (missing(X) && type == "response")
    stop("Need to supply 'X'")
  
  beta <- coef.pcqreg(object=object, lambda=lambda, exact=exact)
  if (type == "coefficients")
    return(beta)
  
  if (is.matrix(beta)) {
    b0 <- beta[1,]
    b <- beta[-1,]
  } else {
    b0 <- beta[1]
    b <- beta[-1]
  }
  if (type == "nvars") {
    if (is.matrix(b))
      return(apply(b!=0, 2, sum))
    else
      return(sum(b!=0))
  }
  if (type == "response")
    return(sweep(X %*% b, 2, b0, "+"))
  
  return(NULL)
}
#===============================================================================================================================#




#===============================================================================================================================#
#
#===============================================================================================================================#

coef.pcqreg <- function(object,
                        lambda,
                        exact) {
  
  if (missing(lambda)) {
    beta <- object$beta
  } else if (exact) {
    # augment the lambda sequence with the new values, and refit the model
    ls <- object$lambda
    ind <- match(lambda, ls, 0)
    if (any(ind == 0)) {
      ls <- unique(rev(sort(c(lambda,ls))))
      object <- update(object, lambda=ls)
      ind <- match(lambda, ls)
    }
    beta <- object$beta[, ind]
  } else {
    # use linear interpolation to estimate coefficients for supplied lambda
    ls <- object$lambda
    lambda[lambda>max(ls)] <- max(ls)
    lambda[lambda<min(ls)] <- min(ls)
    ind <- approx(x=ls, y=seq(ls), xout=lambda)$y
    left <- floor(ind)
    right <- ceiling(ind)
    weight <- ind %% 1
    beta <- (1-weight)*object$beta[,left] + weight*object$beta[,right]
    if (length(lambda) > 1) colnames(beta) <- round(lambda, 4)
  }
  return(beta)
}
#===============================================================================================================================#




#===============================================================================================================================#
# Variables screening by Penalized Partial-Likelihood (PPL)
# CV of variable selection using full cross-validation of both parameters alpha and lambda
#===============================================================================================================================#

cv.ppl <- function(X,
                   y,
                   delta,
                   K,
                   vsarg,
                   vscons,
                   cv,
                   onese,
                   parallel,
                   clus,
                   verbose,
                   seed) {
  
  # Parsing and evaluating 'vsarg' argument to evaluate all parameters
  alpha <- NULL
  nalpha <- NULL
  nlambda <- NULL
  eval(parse( text=unlist(strsplit(x=vsarg, split=",")) ))
  
  M <- 2
  msize <- numeric(M)
  
  if (is.null(alpha)) {
    enalpha <- seq(from=0, to=1, length.out=nalpha)
  } else {
    nalpha <- 1
    enalpha <- alpha
  }
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  
  varsel.list <- vector(mode="list", length=K)
  varsign.list <- vector(mode="list", length=K)
  k <- 1
  while (k <= K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    # Initialize training and test data
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
    }
    enlambda <- vector(mode="list", length=nalpha)
    cv.errmu <- vector(mode="list", length=nalpha)
    cv.errse <- vector(mode="list", length=nalpha)
    for (i in 1:nalpha) {
      cv.fit <- glmnet::cv.glmnet(x=traindata,
                                  y=survival::Surv(time=traintime, event=trainstatus),
                                  alpha=enalpha[i],
                                  nlambda=nlambda,
                                  nfolds=pmax(3,K),
                                  family="cox",
                                  parallel=parallel,
                                  maxit=1e5)
      cv.errmu[[i]] <- cv.fit$cvm
      cv.errse[[i]] <- cv.fit$cvsd
      enlambda[[i]] <- cv.fit$lambda
    }
    cv.errmu <- list2mat(list=cv.errmu, coltrunc="max", fill=NA)
    cv.errse <- list2mat(list=cv.errse, coltrunc="max", fill=NA)
    cv.errmu.index <- as.matrix(which(x=(cv.errmu == as.numeric(cv.errmu)[which.min(cv.errmu)]), arr.ind=TRUE, useNames=FALSE))
    w <- which.min(cv.errmu.index[, 1, drop=TRUE])
    index.min <- as.numeric(cv.errmu.index[w,,drop=TRUE])
    ww <- as.matrix(which(x=cv.errmu < cv.errmu[index.min[1],index.min[2]] + cv.errse[index.min[1],index.min[2]], arr.ind=TRUE, useNames=TRUE))
    index.1se <- c(min(ww[,1]), min(ww[ww[,1]==min(ww[,1]),2]))
    if (is.empty(index.min) || is.empty(index.1se)) {
      varsel.list[[k]] <- list(NA, NA)
      varsign.list[[k]] <- list(NA, NA)
    } else {
      enalpha.min <- enalpha[index.min[1]]
      enalpha.1se <- enalpha[index.1se[1]]
      enlambda.min <- enlambda[[index.min[1]]][index.min[2]]
      enlambda.1se <- enlambda[[index.1se[1]]][index.1se[2]]
      fit <- glmnet::glmnet(x=traindata,
                            y=survival::Surv(time=traintime, event=trainstatus),
                            alpha=ifelse(test=onese, yes=enalpha.1se, no=enalpha.min),
                            family="cox",
                            maxit=1e5)
      w <- apply(X=fit$beta, MARGIN=2, FUN=function(x) {sum(!(is.na(x)) & (x != 0))})
      if (all(w == 0)) {
        varsel.list[[k]] <- list(NA, NA)
        varsign.list[[k]] <- list(NA, NA)
      } else {
        cv.coef <- coef(object=fit, s=c(enlambda.1se, enlambda.min))
        cv.coef <- list(cv.coef[,1], cv.coef[,2])
        varsel <- lapply(cv.coef, FUN=function(x) {which(!(is.na(x)) & (x != 0))})
        if (all(is.na(varsel))) {
          varsel.list[[k]] <- list(NA, NA)
          varsign.list[[k]] <- list(NA, NA)
        } else {
          varsel.list[[k]] <- varsel
          varsign.list[[k]] <- lapply(cv.coef, function(x) sign(x))
        }
      }
    }
    k <- k + 1
  }
  
  # Get the fitted model
  if (all(is.na(varsel.list))) {
    # Failed to select any covariates
    varsel <- NA
    varsign <- NA
    enalpha.min <- NA
    enalpha.1se <- NA
    enlambda.min <- NA
    enlambda.1se <- NA
    msize <- rep(NA, M)
    success <- FALSE
  } else {
    # Get the selected covariates with their signs from all folds for each optimal lambda value
    varsel <- vector(mode="list", length=M)
    varsign <- vector(mode="list", length=M)
    for (l in 1:M) {
      sel.l <- numeric(0)
      sign.l <- numeric(0)
      for (k in 1:K) {
        sel.l <- c(sel.l, varsel.list[[k]][[l]])
        sign.l <- c(sign.l, varsign.list[[k]][[l]])
      }
      na <- (is.na(sel.l))
      nna <- length(which(na))
      sel.l <- sel.l[!na]
      sign.l <- sign.l[!na]
      if ((is.empty(sel.l)) || (is.empty(sign.l)) || (all(sign.l == 0)))  {
        varsel[[l]] <- NA
        varsign[[l]] <- NA
      } else {
        vartab <- table(names(sel.l), useNA="no")
        nvotes <- K - nna
        w <- names(which(vartab >= ceiling(nvotes*vscons)))
        if ((is.na(w)) || (is.empty(w))) {
          varsel[[l]] <- NA
          varsign[[l]] <- NA
        } else {
          varsel[[l]] <- pmatch(x=w, table=colnames(X), nomatch = NA_integer_, duplicates.ok = FALSE)
          names(varsel[[l]]) <- w
          signtab <- table(names(sign.l), sign.l, useNA="no")
          if (length(unique(sign.l)) == 1) {
            signfreq <- rep(x=unique(sign.l), times=nrow(signtab))
            names(signfreq) <- rownames(signtab)
            varsign[[l]] <- signfreq[w]
          } else if (all(sign.l != 1)) {
            signfreq <- apply(signtab[,c("-1","0"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 0
            varsign[[l]] <- signfreq[w]
          } else if (all(sign.l != -1)) {
            signfreq <- apply(signtab[,c("0","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- 0
            signfreq[signfreq==2] <- 1
            varsign[[l]] <- signfreq[w]
          } else {
            signfreq <- apply(signtab[,c("-1","1"),drop=FALSE], 1, which.max)
            signfreq[signfreq==1] <- -1
            signfreq[signfreq==2] <- 1
            varsign[[l]] <- signfreq[w]
          }
          ord <- order(as.numeric(varsel[[l]]))
          varsel[[l]] <- varsel[[l]][ord]
          varsign[[l]] <- varsign[[l]][ord]
        }
      }
    }
    if (all(is.na(varsel))) {
      varsel <- NA
      varsign <- NA
      enalpha.min <- NA
      enalpha.1se <- NA
      enlambda.min <- NA
      enlambda.1se <- NA
      msize <- rep(NA, M)
      success <- FALSE
    } else {
      if (onese) {
        varsel <- varsel[[1]]
        varsign <- varsign[[1]]
        msize[1] <- length(varsel)
      } else {
        varsel <- varsel[[2]]
        varsign <- varsign[[2]]
        msize[2] <- length(varsel)
      }
      success <- TRUE
    }
  }
  
  # Returning updated argument `vsarg` of variable selection parameters
  vsarg <- list("alpha"=alpha,
                "nalpha"=nalpha,
                "nlambda"=nlambda,
                "msize"=msize)
  
  return(list("vsarg"=vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "enalpha.min"=enalpha.min,
              "enalpha.1se"=enalpha.1se,
              "enlambda.min"=enlambda.min,
              "enlambda.1se"=enlambda.1se,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
# Variables screening by SPCA (SPCA)
# CV over number of SPCs and ordered CPH coefficients
#===============================================================================================================================#

cv.spca <- function(X,
                    y,
                    delta,
                    K,
                    vsarg,
                    vscons,
                    cv,
                    parallel,
                    clus,
                    verbose,
                    seed) {
  
  # Parsing and evaluating 'vsarg' argument to evaluate all parameters
  n.thres <- NULL
  n.pcs <- NULL
  n.var <- NULL
  eval(parse( text=unlist(strsplit(x=vsarg, split=",")) ))
  
  M <- 1
  msize <- numeric(M)
  
  if (ncol(X) <= n.var) {
    n.var <- pmin(ncol(X), n.var) - 1
    cat("Warning: Parameter `n.var` was greater than the dimensionality and was reset to ", ncol(X) - 1, ".\n\n", sep="")
  }
  if (ncol(X) <= n.pcs) {
    n.pcs <- pmin(ncol(X), n.pcs) - 1
    cat("Warning: Parameter `n.pcs` was greater than the dimensionality and was reset to ", ncol(X) - 1, ".\n\n", sep="")
  }
  
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  
  varsel.list <- vector(mode="list", length=K)
  varsign.list <- vector(mode="list", length=K)
  k <- 1
  while (k <= K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    # Initialize training and test data
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
    }
    
    # Structure the data
    pca.train.data <- list(x=t(traindata), y=traintime, censoring.status=trainstatus)
    pca.test.data  <- list(x=t(testdata),  y=testtime,  censoring.status=teststatus)
    
    # Compute trained Wald scores for each variable
    ws.train <- superpc::superpc.train(data=pca.train.data, type="survival", s0.perc=NULL)
    lower <- quantile(x=abs(ws.train$feature.scores), probs = 0)
    upper <- quantile(x=abs(ws.train$feature.scores), probs = 1 - (n.var/ncol(traindata)))
    var.thres <- seq(from=lower, to=upper, length=n.thres)
    
    if (cv) {
      # Finding the optimal trained PCR model by internal CV of the Likelihood Ratio Statistic (LRS)
      # Return LRS from full CV w.r.t. threshold values and #PCs
      superpccv <- superpc::superpc.cv(fit=ws.train,
                                       data=pca.train.data,
                                       n.threshold=n.thres,
                                       n.fold=K,
                                       n.components=n.pcs,
                                       min.features=n.var,
                                       max.features=ncol(traindata))
      lrs <- superpccv$scor
      max.mat <- which(lrs == max(lrs, na.rm=TRUE), arr.ind=TRUE)
      max.mat <- max.mat[1,,drop=FALSE]
      max.npcs <- max.mat[1]
      max.thr <- max.mat[2]
    } else {
      max.npcs <- n.pcs[length(n.pcs)]
      max.thr <- n.thres[length(n.thres)]
    }
    
    # Trained model used to predict the survival times on the test set
    superpcfit <- superpc::superpc.predict(object=ws.train,
                                           data=pca.train.data,
                                           newdata=pca.test.data,
                                           prediction.type="continuous",
                                           threshold=var.thres[max.thr],
                                           n.components=max.npcs)
    
    # List of selected variables (exceeding threshold) used in each fold
    varsel <- which(superpcfit$which.features)
    if (is.empty(varsel)) {
      varsel.list[[k]] <- NA
      varsign.list[[k]] <- NA
    } else {
      varsel.list[[k]] <- varsel
      # Feature selection for supervised principal components
      # Forms reduced model to approximate the supervised principal component predictor.
      ft <- superpc::superpc.predict.red(fit=ws.train,
                                         data=pca.train.data,
                                         data.test=pca.test.data,
                                         threshold=var.thres[max.thr],
                                         n.components=max.npcs,
                                         n.shrinkage=n.thres,
                                         prediction.type="continuous")
      # Importance scores of selected features for PC#1
      scores <- ft$import[varsel,1]
      # Cox scores for selected features in order of decreasing absolute value of importance score for PC#1
      lt <- superpc::superpc.listfeatures(data=pca.train.data,
                                          train.obj=ws.train,
                                          fit.red=ft,
                                          component.number=1)
      ordscor <- order(abs(scores), decreasing = TRUE)
      rawscor <- as.numeric(lt[,"Raw-score"])
      names(rawscor) <- names(scores[ordscor])
      varsign.list[[k]] <- sign(rawscor)
      # Determining the separating threshold of median survival time and formation of the two prognostic groups
      pred.pcr <- superpcfit$v.pred.1df
      pred.ind <- (pred.pcr <= median(pred.pcr))
    }
    k <- k + 1
  }
  
  # Get the fitted model
  if (all(is.na(varsel.list))) {
    # Failed to select any covariates
    varsel <- NA
    varsign <- NA
    msize <- NA
    success <- FALSE
  } else {
    # Get the selected covariates with their signs from all folds
    sel.l <- unlist(varsel.list)
    sign.l <- unlist(varsign.list)
    na <- (is.na(sel.l))
    nna <- length(which(na))
    sel.l <- sel.l[!na]
    sign.l <- sign.l[!na]
    if ((is.empty(sel.l)) || (is.empty(sign.l)) || (all(sign.l == 0)))  {
      varsel <- NA
      varsign <- NA
    } else {
      vartab <- table(names(sel.l), useNA="no")
      nvotes <- K - nna
      w <- names(which(vartab >= ceiling(nvotes*vscons)))
      if ((is.na(w)) || (is.empty(w))) {
        varsel <- NA
        varsign <- NA
      } else{
        varsel <- pmatch(x=w, table=colnames(X), nomatch = NA_integer_, duplicates.ok = FALSE)
        names(varsel) <- w
        signtab <- table(names(sign.l), sign.l, useNA="no")
        if (length(unique(sign.l)) == 1) {
          signfreq <- rep(x=unique(sign.l), times=nrow(signtab))
          names(signfreq) <- rownames(signtab)
          varsign <- signfreq[w]
        } else if (all(sign.l != 1)) {
          signfreq <- apply(signtab[,c("-1","0"),drop=FALSE], 1, which.max)
          signfreq[signfreq==1] <- -1
          signfreq[signfreq==2] <- 0
          varsign <- signfreq[w]
        } else if (all(sign.l != -1)) {
          signfreq <- apply(signtab[,c("0","1"),drop=FALSE], 1, which.max)
          signfreq[signfreq==1] <- 0
          signfreq[signfreq==2] <- 1
          varsign <- signfreq[w]
        } else {
          signfreq <- apply(signtab[,c("-1","1"),drop=FALSE], 1, which.max)
          signfreq[signfreq==1] <- -1
          signfreq[signfreq==2] <- 1
          varsign <- signfreq[w]
        }
        ord <- order(as.numeric(varsel))
        varsel <- varsel[ord]
        varsign <- varsign[ord]
      }
    }
    if (all(is.na(varsel))) {
      varsel <- NA
      varsign <- NA
      msize <- NA
      success <- FALSE
    } else {
      msize <- length(varsel)
      success <- TRUE
    }
  }
  
  # Returning updated argument `vsarg` of variable selection parameters
  vsarg <- list("n.thres"=n.thres,
                "n.pcs"=n.pcs,
                "n.var"=n.var,
                "msize"=msize)
  
  return(list("vsarg"=vsarg,
              "varsel"=varsel,
              "varsign"=varsign,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.box(X,
#                          y,
#                          delta,
#                          B,
#                          K,
#                          cv,
#                          cvtype,
#                          cvarg,
#                          varsign,
#                          initcutpts,
#                          groups,
#                          decimals,
#                          probval,
#                          timeval,
#                          parallel.rep,
#                          clus.rep,
#                          verbose,
#                          seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.box <- function(X,
                   y,
                   delta,
                   B,
                   K,
                   cv,
                   cvtype,
                   cvarg,
                   varsign,
                   initcutpts,
                   groups,
                   decimals,
                   probval,
                   timeval,
                   parallel.rep,
                   clus.rep,
                   verbose,
                   seed) {
  
  seed <- seed[1]
  replic <- 1:B
  if (!parallel.rep) {
    if (!is.null(seed)) {
      seed <- (0:(B-1)) + seed
    }
    CV.type.box.obj <- cv.type.box(X=X,
                                   y=y,
                                   delta=delta,
                                   replic=replic,
                                   K=K,
                                   cvarg=cvarg,
                                   cv=cv,
                                   cvtype=cvtype,
                                   varsign=varsign,
                                   initcutpts=initcutpts,
                                   groups=groups,
                                   decimals=decimals,
                                   probval=probval,
                                   timeval=timeval,
                                   parallel.rep=parallel.rep,
                                   verbose=verbose,
                                   seed=seed)
  } else {
    clusterSetRNGStream(cl=clus.rep, iseed=seed)
    obj.cl <- clusterApply(cl=clus.rep, x=replic, fun=cv.type.box,
                           X=X,
                           y=y,
                           delta=delta,
                           K=K,
                           cvarg=cvarg,
                           cv=cv,
                           cvtype=cvtype,
                           varsign=varsign,
                           initcutpts=initcutpts,
                           groups=groups,
                           decimals=decimals,
                           probval=probval,
                           timeval=timeval,
                           parallel.rep=parallel.rep,
                           verbose=verbose,
                           seed=NULL)
    CV.type.box.obj <- list("cv.maxsteps"=numeric(0),
                            "cv.trace"=vector(mode="list", length=B),
                            "cv.boxind"=vector(mode="list", length=B),
                            "cv.boxcut"=vector(mode="list", length=B),
                            "cv.support"=vector(mode="list", length=B),
                            "cv.size"=vector(mode="list", length=B),
                            "cv.lhr"=vector(mode="list", length=B),
                            "cv.lrt"=vector(mode="list", length=B),
                            "cv.cer"=vector(mode="list", length=B),
                            "cv.grp.lhr"=vector(mode="list", length=B),
                            "cv.grp.lrt"=vector(mode="list", length=B),
                            "cv.grp.cer"=vector(mode="list", length=B),
                            "cv.time.bar"=vector(mode="list", length=B),
                            "cv.prob.bar"=vector(mode="list", length=B),
                            "cv.max.time.bar"=vector(mode="list", length=B),
                            "cv.min.prob.bar"=vector(mode="list", length=B),
                            "success"=logical(B))
    # Collect the peeling statistics for each step from all the replicates
    for (b in 1:B) {
      CV.type.box.obj$cv.maxsteps <- c(CV.type.box.obj$cv.maxsteps, obj.cl[[b]]$cv.maxsteps)
      CV.type.box.obj$cv.trace[[b]] <- obj.cl[[b]]$cv.trace
      CV.type.box.obj$cv.boxind[[b]] <- obj.cl[[b]]$cv.boxind
      CV.type.box.obj$cv.boxcut[[b]] <- obj.cl[[b]]$cv.boxcut
      CV.type.box.obj$cv.support[[b]] <- obj.cl[[b]]$cv.support
      CV.type.box.obj$cv.size[[b]] <- obj.cl[[b]]$cv.size
      CV.type.box.obj$cv.lhr[[b]] <- obj.cl[[b]]$cv.lhr
      CV.type.box.obj$cv.lrt[[b]] <- obj.cl[[b]]$cv.lrt
      CV.type.box.obj$cv.cer[[b]] <- obj.cl[[b]]$cv.cer
      CV.type.box.obj$cv.grp.lhr[[b]] <- obj.cl[[b]]$cv.grp.lhr
      CV.type.box.obj$cv.grp.lrt[[b]] <- obj.cl[[b]]$cv.grp.lrt
      CV.type.box.obj$cv.grp.cer[[b]] <- obj.cl[[b]]$cv.grp.cer
      CV.type.box.obj$cv.time.bar[[b]] <- obj.cl[[b]]$cv.time.bar
      CV.type.box.obj$cv.prob.bar[[b]] <- obj.cl[[b]]$cv.prob.bar
      CV.type.box.obj$cv.max.time.bar[[b]] <- obj.cl[[b]]$cv.max.time.bar
      CV.type.box.obj$cv.min.prob.bar[[b]] <- obj.cl[[b]]$cv.min.prob.bar
      CV.type.box.obj$success[b] <- obj.cl[[b]]$success
    }
  }
  
  return(list("cv.maxsteps"=CV.type.box.obj$cv.maxsteps,
              "cv.trace"=CV.type.box.obj$cv.trace,
              "cv.boxind"=CV.type.box.obj$cv.boxind,
              "cv.boxcut"=CV.type.box.obj$cv.boxcut,
              "cv.support"=CV.type.box.obj$cv.support,
              "cv.size"=CV.type.box.obj$cv.size,
              "cv.lhr"=CV.type.box.obj$cv.lhr,
              "cv.lrt"=CV.type.box.obj$cv.lrt,
              "cv.cer"=CV.type.box.obj$cv.cer,
              "cv.grp.lhr"=CV.type.box.obj$cv.grp.lhr,
              "cv.grp.lrt"=CV.type.box.obj$cv.grp.lrt,
              "cv.grp.cer"=CV.type.box.obj$cv.grp.cer,
              "cv.time.bar"=CV.type.box.obj$cv.time.bar,
              "cv.prob.bar"=CV.type.box.obj$cv.prob.bar,
              "cv.max.time.bar"=CV.type.box.obj$cv.max.time.bar,
              "cv.min.prob.bar"=CV.type.box.obj$cv.min.prob.bar,
              "success"=all(CV.type.box.obj$success),
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.type.box(X,
#                               y,
#                               delta,
#                               replic,
#                               K,
#                               cv,
#                               cvtype,
#                               cvarg,
#                               varsign,
#                               initcutpts,
#                               decimals,
#                               groups,
#                               probval,
#                               timeval,
#                               parallel.rep,
#                               verbose,
#                               seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.type.box <- function(X,
                        y,
                        delta,
                        replic,
                        K,
                        cv,
                        cvtype,
                        cvarg,
                        varsign,
                        initcutpts,
                        groups,
                        decimals,
                        probval,
                        timeval,
                        parallel.rep,
                        verbose,
                        seed) {
  
  if (!parallel.rep) {
    
    B <- length(replic)
    success <- logical(B)
    CV.maxsteps <- numeric(B)
    CV.boxind <- vector(mode="list", length=B)
    CV.boxcut <- vector(mode="list", length=B)
    CV.support <- vector(mode="list", length=B)
    CV.size<- vector(mode="list", length=B)
    CV.lhr <- vector(mode="list", length=B)
    CV.lrt <- vector(mode="list", length=B)
    CV.cer <- vector(mode="list", length=B)
    CV.grp.lhr <- vector(mode="list", length=B)
    CV.grp.lrt <- vector(mode="list", length=B)
    CV.grp.cer <- vector(mode="list", length=B)
    CV.trace <- vector(mode="list", length=B)
    CV.time.bar <- vector(mode="list", length=B)
    CV.prob.bar <- vector(mode="list", length=B)
    CV.max.time.bar <- vector(mode="list", length=B)
    CV.min.prob.bar <- vector(mode="list", length=B)
    
    b <- 1
    while (b <= B) {
      if (verbose) {
        cat("Replicate : ", b, "\n", sep="")
        cat("seed : ", seed[b], "\n", sep="")
      }
      if (cvtype == "averaged") {
        CVBOX <- cv.ave.box(X=X,
                            y=y,
                            delta=delta,
                            K=K,
                            cv=cv,
                            cvarg=cvarg,
                            groups=groups,
                            decimals=decimals,
                            probval=probval,
                            timeval=timeval,
                            varsign=varsign,
                            initcutpts=initcutpts,
                            verbose=verbose,
                            seed=seed[b])
      } else if (cvtype == "combined") {
        CVBOX <- cv.comb.box(X=X,
                             y=y,
                             delta=delta,
                             K=K,
                             cv=cv,
                             cvarg=cvarg,
                             groups=groups,
                             decimals=decimals,
                             probval=probval,
                             timeval=timeval,
                             varsign=varsign,
                             initcutpts=initcutpts,
                             verbose=verbose,
                             seed=seed[b])
      } else {
        stop("Invalid CV type option. Exiting ... \n\n")
      }
      success[b] <- CVBOX$success
      CV.maxsteps[b] <- CVBOX$cvfit$cv.maxsteps
      CV.trace[[b]] <- CVBOX$cvfit$cv.trace
      CV.boxind[[b]] <- CVBOX$cvfit$cv.boxind
      CV.boxcut[[b]] <- CVBOX$cvfit$cv.boxcut
      CV.support[[b]] <- CVBOX$cvfit$cv.stats$cv.support
      CV.size[[b]] <- CVBOX$cvfit$cv.stats$cv.size
      CV.lhr[[b]] <- CVBOX$cvfit$cv.stats$cv.lhr
      CV.lrt[[b]] <- CVBOX$cvfit$cv.stats$cv.lrt
      CV.cer[[b]] <- CVBOX$cvfit$cv.stats$cv.cer
      CV.grp.lhr[[b]] <- CVBOX$cvfit$cv.stats$cv.grp.lhr
      CV.grp.lrt[[b]] <- CVBOX$cvfit$cv.stats$cv.grp.lrt
      CV.grp.cer[[b]] <- CVBOX$cvfit$cv.stats$cv.grp.cer
      CV.time.bar[[b]] <- CVBOX$cvfit$cv.stats$cv.time.bar
      CV.prob.bar[[b]] <- CVBOX$cvfit$cv.stats$cv.prob.bar
      CV.max.time.bar[[b]] <- CVBOX$cvfit$cv.stats$cv.max.time.bar
      CV.min.prob.bar[[b]] <- CVBOX$cvfit$cv.stats$cv.min.prob.bar
      b <- b + 1
    }
    
  } else {
    
    if (cvtype == "averaged") {
      CVBOX <- cv.ave.box(X=X,
                          y=y,
                          delta=delta,
                          K=K,
                          cv=cv,
                          cvarg=cvarg,
                          groups=groups,
                          decimals=decimals,
                          probval=probval,
                          timeval=timeval,
                          varsign=varsign,
                          initcutpts=initcutpts,
                          verbose=verbose,
                          seed=NULL)
    } else if (cvtype == "combined") {
      CVBOX <- cv.comb.box(X=X,
                           y=y,
                           delta=delta,
                           K=K,
                           cv=cv,
                           cvarg=cvarg,
                           groups=groups,
                           decimals=decimals,
                           probval=probval,
                           timeval=timeval,
                           varsign=varsign,
                           initcutpts=initcutpts,
                           verbose=verbose,
                           seed=NULL)
    } else {
      stop("Invalid CV type option. Exiting ... \n\n")
    }
    success <- CVBOX$success
    CV.maxsteps <- CVBOX$cvfit$cv.maxsteps
    CV.trace <- CVBOX$cvfit$cv.trace
    CV.boxind <- CVBOX$cvfit$cv.boxind
    CV.boxcut <- CVBOX$cvfit$cv.boxcut
    CV.support <- CVBOX$cvfit$cv.stats$cv.support
    CV.size <- CVBOX$cvfit$cv.stats$cv.size
    CV.lhr <- CVBOX$cvfit$cv.stats$cv.lhr
    CV.lrt <- CVBOX$cvfit$cv.stats$cv.lrt
    CV.cer <- CVBOX$cvfit$cv.stats$cv.cer
    CV.grp.lhr <- CVBOX$cvfit$cv.stats$cv.grp.lhr
    CV.grp.lrt <- CVBOX$cvfit$cv.stats$cv.grp.lrt
    CV.grp.cer <- CVBOX$cvfit$cv.stats$cv.grp.cer
    CV.time.bar <- CVBOX$cvfit$cv.stats$cv.time.bar
    CV.prob.bar <- CVBOX$cvfit$cv.stats$cv.prob.bar
    CV.max.time.bar <- CVBOX$cvfit$cv.stats$cv.max.time.bar
    CV.min.prob.bar <- CVBOX$cvfit$cv.stats$cv.min.prob.bar
    
  }
  
  return(list("cv.maxsteps"=CV.maxsteps,
              "cv.trace"=CV.trace,
              "cv.boxind"=CV.boxind,
              "cv.boxcut"=CV.boxcut,
              "cv.support"=CV.support,
              "cv.size"=CV.size,
              "cv.lhr"=CV.lhr,
              "cv.lrt"=CV.lrt,
              "cv.cer"=CV.cer,
              "cv.grp.lhr"=CV.grp.lhr,
              "cv.grp.lrt"=CV.grp.lrt,
              "cv.grp.cer"=CV.grp.cer,
              "cv.time.bar"=CV.time.bar,
              "cv.prob.bar"=CV.prob.bar,
              "cv.max.time.bar"=CV.max.time.bar,
              "cv.min.prob.bar"=CV.min.prob.bar,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                    cv.pval (X,
#                             y,
#                             delta,
#                             K,
#                             A,
#                             pv,
#                             cv,
#                             cvtype,
#                             cvarg,
#                             groups,
#                             decimals,
#                             varsign,
#                             initcutpts,
#                             obs.chisq,
#                             parallel.pv,
#                             clus.pv,
#                             verbose,
#                             seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.pval <- function(X,
                    y,
                    delta,
                    K,
                    A,
                    pv,
                    cv,
                    cvtype,
                    cvarg,
                    groups,
                    decimals,
                    varsign,
                    initcutpts,
                    obs.chisq,
                    parallel.pv,
                    clus.pv,
                    verbose,
                    seed) {
  
  if (pv) {
    
    cat("Computation of p-values ... \n")
    
    # Parsing and evaluating 'cvarg' parameters to evaluate 'L' and 'peelcriterion'
    alpha <- NULL
    beta <- NULL
    L <- NULL
    peelcriterion <- NULL
    cvcriterion <- NULL
    eval(parse( text=unlist(strsplit(x=cvarg, split=",")) ))
    
    seed <- seed[1]
    
    if (cv) {
      # Computation of log-rank permutation p-values
      # using the permutation distribution for the null distribution
      
      nperm <- 1:A
      if (!parallel.pv) {
        if (!is.null(seed)) {
          seed <- (0:(A-1)) + seed
        }
        null.chisq <- cv.null(X=X,
                              y=y,
                              delta=delta,
                              nperm=nperm,
                              K=K,
                              cv=cv,
                              cvtype=cvtype,
                              cvarg=cvarg,
                              groups=groups,
                              peelcriterion=peelcriterion,
                              decimals=decimals,
                              varsign=varsign,
                              initcutpts=initcutpts,
                              parallel.pv=parallel.pv,
                              verbose=verbose,
                              seed=seed)
        null.chisq <- t(list2mat(list=null.chisq, coltrunc=L, fill=NA))
      } else {
        clusterSetRNGStream(cl=clus.pv, iseed=seed)
        null.cl <- clusterApply(cl=clus.pv, x=nperm, fun=cv.null,
                                X=X,
                                y=y,
                                delta=delta,
                                K=K,
                                cv=cv,
                                cvtype=cvtype,
                                cvarg=cvarg,
                                groups=groups,
                                peelcriterion=peelcriterion,
                                decimals=decimals,
                                varsign=varsign,
                                initcutpts=initcutpts,
                                parallel.pv=parallel.pv,
                                verbose=verbose,
                                seed=NULL)
        null.chisq <- t(list2mat(list=null.cl, coltrunc=L, fill=NA))
      }
      pval <- numeric(L)
      for (l in 1:L) {
        pval[l] <- mean((null.chisq[l,] >= obs.chisq[l]), na.rm=TRUE)
      }
      pval <- round(pval, digits=floor(log(base=10, A)))
      names(pval) <- paste("step", 0:(L-1), sep="")
      
      return(list("pval"=pval, "seed"=seed))
      
    } else {
      # Computation of log-rank Mantel-Haenszel p-values
      # using the Chi-Squared distribution (with 1 df) for the null distribution
      
      pval <- numeric(L)
      for (l in 1:(L)) {
        pval[l] <- 1 - pchisq(q=obs.chisq[l], df=1)
      }
      pval <- round(pval, digits=decimals)
      names(pval) <- paste("step", 0:(L-1), sep="")
      
      return(list("pval"=pval, "seed"=seed))
      
    }
    
  } else {
    
    cat("No computation of p-values. \n")
    
    return(NULL)
    
  }
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.null (X,
#                            y,
#                            delta,
#                            nperm,
#                            K,
#                            cv,
#                            cvtype,
#                            cvarg,
#                            groups,
#                            peelcriterion,
#                            decimals,
#                            varsign,
#                            initcutpts,
#                            parallel.pv,
#                            verbose,
#                            seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.null <- function(X,
                    y,
                    delta,
                    nperm,
                    K,
                    cv,
                    cvtype,
                    cvarg,
                    groups,
                    peelcriterion,
                    decimals,
                    varsign,
                    initcutpts,
                    parallel.pv,
                    verbose,
                    seed) {
  
  if (!parallel.pv) {
    
    A <- length(nperm)
    n <- nrow(X)
    null.chisq <- vector(mode="list", length=A)
    
    a <- 1
    while (a <= A) {
      if (verbose) {
        cat("Permutation : ", a, "\n")
        cat("seed : ", seed[a], "\n", sep="")
      }
      perm.ind <- sample(x = 1:n, size = n, replace = FALSE, prob = NULL)
      perm.y <- y[perm.ind]
      perm.delta <- delta[perm.ind]
      if (cvtype == "averaged") {
        obj <- tryCatch({cv.ave.box(X=X,
                                    y=perm.y,
                                    delta=perm.delta,
                                    K=K,
                                    cv=cv,
                                    cvarg=cvarg,
                                    varsign=varsign,
                                    initcutpts=initcutpts,
                                    groups=groups,
                                    decimals=decimals,
                                    probval=NULL,
                                    timeval=NULL,
                                    verbose=verbose,
                                    seed=seed[a])},
                        error=function(){NULL})
        if (is.list(obj)) {
          if (peelcriterion != "grp") {
            null.chisq[[a]] <- obj$cvfit$cv.stats$cv.lrt
          } else if (peelcriterion == "grp") {
            null.chisq[[a]] <- obj$cvfit$cv.stats$cv.grp.lrt
          } else {
            stop("Invalid peeling criterion. Exiting ...\n\n")
          }
          a <- a + 1
        } else {
          if (verbose) cat("Permutation sample dropped... \n")
        }
      } else if (cvtype == "combined") {
        obj <- tryCatch({cv.comb.box(X=X,
                                     y=perm.y,
                                     delta=perm.delta,
                                     K=K,
                                     cv=cv,
                                     cvarg=cvarg,
                                     varsign=varsign,
                                     initcutpts=initcutpts,
                                     groups=groups,
                                     decimals=decimals,
                                     probval=NULL,
                                     timeval=NULL,
                                     verbose=verbose,
                                     seed=seed[a])},
                        error=function(){NULL})
        if (is.list(obj)) {
          if (peelcriterion != "grp") {
            null.chisq[[a]] <- obj$cvfit$cv.stats$cv.lrt
          } else if (peelcriterion == "grp") {
            null.chisq[[a]] <- obj$cvfit$cv.stats$cv.grp.lrt
          } else {
            stop("Invalid peeling criterion. Exiting ...\n\n")
          }
          a <- a + 1
        } else {
          if (verbose) cat("Permutation sample dropped... \n")
        }
      } else {
        stop("Invalid CV type option. Exiting ... \n\n")
      }
    }
    
  } else {
    
    n <- nrow(X)
    perm.ind <- sample(x = 1:n, size = n, replace = FALSE, prob = NULL)
    perm.y <- y[perm.ind]
    perm.delta <- delta[perm.ind]
    if (cvtype == "averaged") {
      obj <- tryCatch({cv.ave.box(X=X,
                                  y=perm.y,
                                  delta=perm.delta,
                                  K=K,
                                  cv=cv,
                                  cvarg=cvarg,
                                  decimals=decimals,
                                  varsign=varsign,
                                  initcutpts=initcutpts,
                                  groups=groups,
                                  probval=NULL,
                                  timeval=NULL,
                                  verbose=verbose,
                                  seed=NULL)},
                      error=function(w){NULL})
      if (is.list(obj)) {
        if (peelcriterion != "grp") {
          null.chisq <- obj$cvfit$cv.stats$cv.lrt
        } else if (peelcriterion == "grp") {
          null.chisq <- obj$cvfit$cv.stats$cv.grp.lrt
        } else {
          stop("Invalid peeling criterion. Exiting ...\n\n")
        }
      } else {
        null.chisq <- NA
        if (verbose) cat("Permutation sample dropped... \n")
      }
    } else if (cvtype == "combined") {
      obj <- tryCatch({cv.comb.box(X=X,
                                   y=perm.y,
                                   delta=perm.delta,
                                   K=K,
                                   cv=cv,
                                   cvarg=cvarg,
                                   decimals=decimals,
                                   varsign=varsign,
                                   initcutpts=initcutpts,
                                   groups=groups,
                                   probval=NULL,
                                   timeval=NULL,
                                   verbose=verbose,
                                   seed=NULL)},
                      error=function(w){NULL})
      if (is.list(obj)) {
        if (peelcriterion != "grp") {
          null.chisq <- obj$cvfit$cv.stats$cv.lrt
        } else if (peelcriterion == "grp") {
          null.chisq <- obj$cvfit$cv.stats$cv.grp.lrt
        } else {
          stop("Invalid peeling criterion. Exiting ...\n\n")
        }
      } else {
        null.chisq <- NA
        if (verbose) cat("Permutation sample dropped... \n")
      }
    } else {
      stop("Invalid CV type option. Exiting ... \n\n")
    }
  }
  
  return(null.chisq)
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.ave.box (X,
#                               y,
#                               delta,
#                               K,
#                               cv,
#                               cvarg,
#                               varsign,
#                               initcutpts,
#                               groups,
#                               decimals,
#                               probval,
#                               timeval,
#                               verbose,
#                               seed)
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.ave.box <- function(X,
                       y,
                       delta,
                       K,
                       cv,
                       cvarg,
                       varsign,
                       initcutpts,
                       groups,
                       decimals,
                       probval,
                       timeval,
                       verbose,
                       seed) {
  
  n <- nrow(X)
  p <- ncol(X)
  
  # Parsing and evaluating 'cvarg' parameters to evaluate 'beta'
  alpha <- NULL
  beta <- NULL
  L <- NULL
  peelcriterion <- NULL
  cvcriterion <- NULL
  eval(parse( text=unlist(strsplit(x=cvarg, split=",")) ))
  
  # Creating argument `arg` of `cv.ave.peel` for PRSP parameters
  arg <- paste("alpha=", alpha,
               ",beta=", beta,
               ",L=", L,
               ",peelcriterion=\"", peelcriterion, "\"", sep="")
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  
  boxstat.list <- vector(mode="list", length=K)
  boxcut.list <- vector(mode="list", length=K)
  trace.list <- vector(mode="list", length=K)
  maxsteps <- numeric(K)
  
  for (k in 1:K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    # Initialize training and test data
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
      traingroups <- testgroups <- groups[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      traingroups <- groups[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
      testgroups <- groups[folds$perm[(folds$which == k)]]
    }
    peelobj <- cv.ave.peel(traindata=traindata, trainstatus=trainstatus, traintime=traintime,
                           testdata=testdata, teststatus=teststatus, testtime=testtime,
                           probval=probval, timeval=timeval,
                           varsign=varsign, initcutpts=initcutpts, arg=arg, 
                           traingroups=traingroups, testgroups=testgroups)
    # Store the test set results from each fold
    maxsteps[k] <- peelobj$maxsteps
    boxstat.list[[k]] <- peelobj$boxstat
    boxcut.list[[k]] <- peelobj$boxcut
    trace.list[[k]] <- peelobj$trace
  }
  
  # Cross-validated maximum peeling length from all folds
  CV.Lm <- min(maxsteps)
  
  # Truncate the cross-validated quantities from all folds to the same cross-validated length
  for (k in 1:K) {
    boxcut.list[[k]] <- boxcut.list[[k]][1:CV.Lm,,drop=FALSE]
  }
  
  # Compute the averaged box statistics for each step from all the folds
  # Each entry or row signifies a step
  CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p, dimnames=list(paste("step", 0:(CV.Lm-1), sep=""), colnames(X)))
  for (l in 1:CV.Lm) {
    summincut <- matrix(NA, K, p)
    for (k in 1:K) {
      summincut[k,] <- boxcut.list[[k]][l,]
    }
    CV.boxcut[l, ] <- colMeans(summincut, na.rm=TRUE)
  }
  rownames(CV.boxcut) <- paste("step", 0:(CV.Lm-1), sep="")
  colnames(CV.boxcut) <- colnames(X)
  
  # Get the box membership indicator vector of all observations for each step from all the folds
  # Based on the average box over the folds
  CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
  for (l in 1:CV.Lm) {
    boxcut <- CV.boxcut[l, ] * varsign
    X.cut <- t(t(X) * varsign)
    X.ind <- t(t(X.cut) >= boxcut)
    CV.boxind[l,] <- (rowMeans(X.ind) == 1)  # Set as TRUE which observations are inside the box boudaries for all axes directions
  }
  rownames(CV.boxind) <- paste("step", 0:(CV.Lm-1), sep="")
  colnames(CV.boxind) <- rownames(X)
  
  # Get the adjusted cross-validated maximum peeling length, thresholded by minimal box support
  CV.Lm <- max(which(apply(CV.boxind, 1, function(x) {length(which(x))/n >= max(1/n, beta)})))
  
  # Get the adjusted cross-validated quantities from all folds to the same cross-validated length
  for (k in 1:K) {
    boxstat.list[[k]] <- boxstat.list[[k]][1:CV.Lm]
  }
  
  # Get the adjusted test box defnition and membership indicator vector of all observations for each step from all the folds
  CV.boxind <- CV.boxind[1:CV.Lm,,drop=FALSE]
  CV.boxcut <- CV.boxcut[1:CV.Lm,,drop=FALSE]
  
  # Get the box sample size and support based on `CV.boxind`
  CV.size <- apply(CV.boxind, 1, sum)
  names(CV.size) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.support <- CV.size/n
  names(CV.support) <- paste("step", 0:(CV.Lm-1), sep="")
  
  # Get the variable traces
  # Variable traces are first stacked and truncated in a matrix where folds are by rows and steps by columns
  CV.trace <- lapply.mat(X=trace.list,
                         coltrunc=CV.Lm,
                         fill=NA,
                         MARGIN=2,
                         FUN=function(x) {
                           vote <- table(x, useNA="no")
                           w <- as.numeric(names(which.max(vote)))
                           if(is.empty(w)) {
                             return(NA)
                           } else {
                             return(w)
                           }
                         })
  names(CV.trace) <- paste("step", 0:(CV.Lm-1), sep="")
  
  # Box peeling rules for each step
  CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p, dimnames=list(paste("step", 0:(CV.Lm-1), sep=""), colnames(X))))
  for (j in 1:p) {
    if (varsign[j] > 0) {
      ss <- ">="
    } else {
      ss <- "<="
    }
    CV.rules[, j] <- paste(colnames(X)[j], ss, format(x=CV.boxcut[, j], digits=decimals, nsmall=decimals), sep="")
  }
  
  # Compute the averaged box statistics for each step from all the folds
  # Each entry or row signifies a step
  CV.lhr <- rep(x=NA, times=CV.Lm)
  names(CV.lhr) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.lrt <- rep(x=NA, times=CV.Lm)
  names(CV.lrt) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.cer <- rep(x=NA, times=CV.Lm)
  names(CV.cer) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.lhr <- rep(x=NA, times=CV.Lm)
  names(CV.grp.lhr) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.lrt <- rep(x=NA, times=CV.Lm)
  names(CV.grp.lrt) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.cer <- rep(x=NA, times=CV.Lm)
  names(CV.grp.cer) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.time.bar <- rep(x=NA, times=CV.Lm)
  names(CV.time.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.prob.bar <- rep(x=NA, times=CV.Lm)
  names(CV.prob.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.max.time.bar <- rep(x=NA, times=CV.Lm)
  names(CV.max.time.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
  names(CV.min.prob.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  for (l in 1:CV.Lm) {
    sumlhr <- rep(x=NA, times=K)
    sumlrt <- rep(x=NA, times=K)
    sumcer <- rep(x=NA, times=K)
    sumgrplrt <- rep(x=NA, times=K)
    sumgrplrt <- rep(x=NA, times=K)
    sumgrpcer <- rep(x=NA, times=K)
    sumtime <- rep(x=NA, times=K)
    sumprob <- rep(x=NA, times=K)
    summaxtime <- rep(x=NA, times=K)
    summinprob <- rep(x=NA, times=K)
    for (k in 1:K) {
      outbounds <- boxstat.list[[k]][[l]]
      if (!is.null(outbounds)) {
        sumlhr[k] <- boxstat.list[[k]][[l]][[1]]
        sumlrt[k] <- boxstat.list[[k]][[l]][[2]]
        sumcer[k] <- boxstat.list[[k]][[l]][[3]]
        sumgrplhr[k] <- boxstat.list[[k]][[l]][[4]]
        sumgrplrt[k] <- boxstat.list[[k]][[l]][[5]]
        sumgrpcer[k] <- boxstat.list[[k]][[l]][[6]]
        sumtime[k] <- boxstat.list[[k]][[l]][[7]]
        sumprob[k] <- boxstat.list[[k]][[l]][[8]]
        summaxtime[k] <- boxstat.list[[k]][[l]][[9]]
        summinprob[k] <- boxstat.list[[k]][[l]][[10]]
      }
    }
    CV.lhr[l] <- mean(sumlhr, na.rm=TRUE)
    CV.lrt[l] <- mean(sumlrt, na.rm=TRUE)
    CV.cer[l] <- mean(sumcer, na.rm=TRUE)
    CV.grp.lhr[l] <- mean(sumgrplhr, na.rm=TRUE)
    CV.grp.lrt[l] <- mean(sumgrplrt, na.rm=TRUE)
    CV.grp.cer[l] <- mean(sumgrpcer, na.rm=TRUE)
    CV.time.bar[l] <- mean(sumtime, na.rm=TRUE)
    CV.prob.bar[l] <- mean(sumprob, na.rm=TRUE)
    CV.max.time.bar[l] <- mean(summaxtime, na.rm=TRUE)
    CV.min.prob.bar[l] <- mean(summinprob, na.rm=TRUE)
  }
  CV.maxsteps <- CV.Lm

  # Formating the results depending on successful cross-validated PRSP algorithm
  if (cvcriterion == "lhr") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.lhr)) || all(is.nan(CV.lhr))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by maximization of the LHR (between inbox and outbox test samples)
        success <- TRUE
      }
    } else if (peelcriterion == "grp") {
      if (all(is.na(CV.lhr)) || all(is.nan(CV.lhr)) || all(is.na(CV.grp.lhr)) || all(is.nan(CV.grp.lhr))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else {
      # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
      success <- TRUE
    }
  } else if (cvcriterion == "lrt") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.lrt)) || all(is.nan(CV.lrt))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by maximization of the LRT (between inbox and outbox test samples)
        success <- TRUE
      }
    } else if (peelcriterion == "grp") {
      if (all(is.na(CV.lrt)) || all(is.nan(CV.lrt)) || all(is.na(CV.grp.lrt)) || all(is.nan(CV.grp.lrt))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else {
      # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
      success <- TRUE
    }
  } else if (cvcriterion == "cer") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.cer)) || all(is.nan(CV.cer))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else if (peelcriterion == "grp") {
      if (all(is.na(CV.cer)) || all(is.nan(CV.cer)) || all(is.na(CV.grp.cer)) || all(is.nan(CV.grp.cer))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else {
      stop("Invalid peeling criterion. Exiting ...\n\n")
    }
  } else {
    stop("Invalid CV criterion option. Exiting ... \n\n")
  }

  # Box statistics for each step
  CV.stats <-  data.frame("cv.support"=CV.support,
                          "cv.size"=CV.size,
                          "cv.lhr"=CV.lhr,
                          "cv.lrt"=CV.lrt,
                          "cv.cer"=CV.cer,
                          "cv.grp.lhr"=CV.grp.lhr,
                          "cv.grp.lrt"=CV.grp.lrt,
                          "cv.grp.cer"=CV.grp.cer,
                          "cv.time.bar"=CV.time.bar,
                          "cv.prob.bar"=CV.prob.bar,
                          "cv.max.time.bar"=CV.max.time.bar,
                          "cv.min.prob.bar"=CV.min.prob.bar)
  rownames(CV.stats) <- paste("step", 0:(CV.Lm-1), sep="")
  
  # Create the return object 'CV.fit'
  CV.fit <- list("cv.maxsteps"=CV.maxsteps,
                 "cv.boxcut"=CV.boxcut,
                 "cv.boxind"=CV.boxind,
                 "cv.trace"=CV.trace,
                 "cv.rules"=CV.rules,
                 "cv.stats"=CV.stats)
  
  return(list("X"=X,
              "y"=y,
              "delta"=delta,
              "cvfit"=CV.fit,
              "success"=success,
              "seed"=seed))
}
#===============================================================================================================================#




#===============================================================================================================================#
#===============#
# Usage         :
#===============#
#                   cv.comb.box (X,
#                                y,
#                                delta,
#                                K,
#                                cv,
#                                cvarg,
#                                varsign,
#                                initcutpts,
#                                groups,
#                                probval,
#                                timeval,
#                                decimals,
#                                verbose,
#                                seed)
#
#
#===============#
# Description   :
#===============#
#
#===============#
# Arguments     :
#===============#
#
#===============#
# Values        :
#===============#
#
#===============================================================================================================================#

cv.comb.box <- function(X,
                        y,
                        delta,
                        K,
                        cv,
                        cvarg,
                        varsign,
                        initcutpts,
                        groups,
                        decimals,
                        probval,
                        timeval,
                        verbose,
                        seed) {
  n <- nrow(X)
  p <- ncol(X)
  
  # Parsing and evaluating 'cvarg' parameters to evaluate 'beta' and 'peelcriterion'
  alpha <- NULL
  beta <- NULL
  L <- NULL
  peelcriterion <- NULL
  cvcriterion <- NULL
  eval(parse( text=unlist(strsplit(x=cvarg, split=",")) ))
  
  # Creating argument `arg` of `cv.comb.peel` for PRSP parameters
  arg <- paste("alpha=", alpha,
               ",beta=", beta,
               ",L=", L,
               ",peelcriterion=\"", peelcriterion, "\"", sep="")
  
  folds <- cv.folds(y=delta, K=K, seed=seed)
  ord <- folds$foldkey
  
  times.list <- vector(mode="list", length=K)
  status.list <- vector(mode="list", length=K)
  boxind.list <- vector(mode="list", length=K)
  boxcut.list <- vector(mode="list", length=K)
  trace.list <- vector(mode="list", length=K)
  maxsteps <- numeric(K)
  
  for (k in 1:K) {
    if ((verbose) && (cv)) cat("Fold : ", k, "\n", sep="")
    if (K == 1) {
      traindata <- testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      traintime <- testtime <- y[folds$perm[(folds$which == k)]]
      trainstatus <- teststatus <- delta[folds$perm[(folds$which == k)]]
      traingroups <- testgroups <- groups[folds$perm[(folds$which == k)]]
    } else {
      traindata <- X[folds$perm[(folds$which != k)], , drop=FALSE]
      traintime <- y[folds$perm[(folds$which != k)]]
      trainstatus <- delta[folds$perm[(folds$which != k)]]
      traingroups <- groups[folds$perm[(folds$which != k)]]
      testdata <- X[folds$perm[(folds$which == k)], , drop=FALSE]
      testtime <- y[folds$perm[(folds$which == k)]]
      teststatus <- delta[folds$perm[(folds$which == k)]]
      testgroups <- groups[folds$perm[(folds$which == k)]]
    }
    # Store the test set data from each fold (Note: the order of observations of times and status from each fold is kept in the list)
    times.list[[k]] <- testtime
    status.list[[k]] <- teststatus
    peelobj <- cv.comb.peel(traindata=traindata, trainstatus=trainstatus, traintime=traintime,
                            testdata=testdata, teststatus=teststatus, testtime=testtime,
                            varsign=varsign, initcutpts=initcutpts, 
                            arg=arg, traingroups=traingroups, testgroups=testgroups)
    # Store the test set results from each fold
    maxsteps[k] <- peelobj$maxsteps
    boxind.list[[k]] <- peelobj$boxind
    boxcut.list[[k]] <- peelobj$boxcut
    trace.list[[k]] <- peelobj$trace
  }
  
  # Cross-validated maximum peeling length from all folds
  CV.Lm <- min(maxsteps)
  
  # Get the test box membership indicator vector of all observations for each step from all the folds
  # Based on the combined membership indicator vectors over the folds
  # Re-ordered by initial order of observations
  CV.boxind <- cbindlist(boxind.list, trunc=CV.Lm)[,ord,drop=FALSE]
  rownames(CV.boxind) <- paste("step", 0:(CV.Lm-1), sep="")
  colnames(CV.boxind) <- rownames(X)
  
  # Get the adjusted cross-validated maximum peeling length, thresholded by minimal box support
  CV.Lm <- max(which(apply(CV.boxind, 1, function(x) {length(which(x))/n >= max(1/n, beta)})))
  
  # Get the adjusted test box membership indicator vector of all observations for each step from all the folds
  CV.boxind <- CV.boxind[1:CV.Lm,,drop=FALSE]
  
  # Get the box sample size and support based on `CV.boxind`
  CV.size <- apply(CV.boxind, 1, sum)
  names(CV.size) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.support <- CV.size/n
  names(CV.support) <- paste("step", 0:(CV.Lm-1), sep="")
  
  # Concatenates the observations of test times and status from all folds
  # Re-ordered by initial order of observations
  CV.times <- unlist(times.list)[ord]
  CV.status <- unlist(status.list)[ord]
  
  # Get the combined boxcut (truncated to the same cross-validated length) for each step from all the folds
  # using the circumscribing box to the conmbined test set in-box samples over all the folds
  CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p, dimnames=list(paste("step", 0:(CV.Lm-1), sep=""), colnames(X)))
  tmparray <- list2array(list=boxcut.list, rowtrunc=CV.Lm)
  for (l in 1:CV.Lm) {
    for (j in 1:p) {
      if (varsign[j] > 0) {
        CV.boxcut[l,j] <- min(X[CV.boxind[l,],j])
      } else {
        CV.boxcut[l,j] <- max(X[CV.boxind[l,],j])
      }
    }
  }
  
  # Get the variable traces
  # Variable traces are first stacked and truncated in a matrix where folds are by rows and steps by columns
  CV.trace <- lapply.mat(X=trace.list,
                         coltrunc=CV.Lm,
                         fill=NA,
                         MARGIN=2,
                         FUN=function(x) {
                           vote <- table(x, useNA="no")
                           w <- as.numeric(names(which.max(vote)))
                           if(is.empty(w)) {
                             return(NA)
                           } else {
                             return(w)
                           }
                         })
  names(CV.trace) <- paste("step", 0:(CV.Lm-1), sep="")
  
  # Box peeling rules for each step
  CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p, dimnames=list(paste("step", 0:(CV.Lm-1), sep=""), colnames(X))))
  for (j in 1:p) {
    if (varsign[j] > 0) {
      ss <- ">="
    } else {
      ss <- "<="
    }
    CV.rules[, j] <- paste(colnames(X)[j], ss, format(x=CV.boxcut[, j], digits=decimals, nsmall=decimals), sep="")
  }
  
  # Compute the combined test box statistics from all folds for all steps, each entry or row signifies a step
  CV.lhr <- rep(x=NA, times=CV.Lm)
  names(CV.lhr) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.lrt <- rep(x=NA, times=CV.Lm)
  names(CV.lrt) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.cer <- rep(x=NA, times=CV.Lm)
  names(CV.cer) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.lhr <- rep(x=NA, times=CV.Lm)
  names(CV.grp.lhr) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.lrt <- rep(x=NA, times=CV.Lm)
  names(CV.grp.lrt) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.grp.cer <- rep(x=NA, times=CV.Lm)
  names(CV.grp.cer) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.time.bar <- rep(x=NA, times=CV.Lm)
  names(CV.time.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.prob.bar <- rep(x=NA, times=CV.Lm)
  names(CV.prob.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.max.time.bar <- rep(x=NA, times=CV.Lm)
  names(CV.max.time.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
  names(CV.min.prob.bar) <- paste("step", 0:(CV.Lm-1), sep="")
  timemat <- matrix(NA, nrow=CV.Lm, ncol=n)
  probmat <- matrix(NA, nrow=CV.Lm, ncol=n)
  ind.rem <- numeric(0)
  grpind <- (groups == levels(groups)[1])
  grpind1 <- 1*grpind
  for (l in 1:CV.Lm) {
    boxind <- CV.boxind[l,]
    boxind1 <- 1*boxind
    wb <- which(boxind1 == 1)
    if (peelcriterion != "grp") {
      if (l == 1) {
        surv.fit <- survival::survfit(survival::Surv(CV.times[boxind], CV.status[boxind]) ~ 1, na.action=na.exclude)
        timemat[l, (1:length(surv.fit$time))] <- surv.fit$time
        probmat[l, (1:length(surv.fit$surv))] <- surv.fit$surv
        CV.lhr[l] <- 0
        CV.lrt[l] <- 0
        CV.cer[l] <- 1
        CV.grp.lhr[l] <- NA
        CV.grp.lrt[l] <- NA
        CV.grp.cer[l] <- NA
      } else if ((sum(boxind, na.rm=TRUE) != length(boxind[!is.na(boxind)])) && (sum(boxind, na.rm=TRUE) != 0)) {
        surv.fit <- survival::survfit(survival::Surv(CV.times[boxind], CV.status[boxind]) ~ 1, na.action=na.exclude)
        timemat[l, (1:length(surv.fit$time))] <- surv.fit$time
        probmat[l, (1:length(surv.fit$surv))] <- surv.fit$surv
        surv.formula <- (survival::Surv(CV.times, CV.status) ~ 1 + boxind1)
        coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
        CV.lhr[l] <- coxobj$coef
        CV.lrt[l] <- survival::survdiff(surv.formula, na.action=na.exclude, rho=0)$chisq
        predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
        CV.cer[l] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(CV.times, CV.status))['C Index']
        CV.grp.lhr[l] <- NA
        CV.grp.lrt[l] <- NA
        CV.grp.cer[l] <- NA
      } else {
        timemat[l, ] <- NA
        probmat[l, ] <- NA
        CV.lhr[l] <- NA
        CV.lrt[l] <- NA
        CV.cer[l] <- NA
        CV.grp.lhr[l] <- NA
        CV.grp.lrt[l] <- NA
        CV.grp.cer[l] <- NA
        ind.rem <- c(ind.rem, l)
      }
    } else if (peelcriterion == "grp") {
      if (l == 1) {
        surv.fit <- survival::survfit(survival::Surv(CV.times, CV.status) ~ 1 + grpind1, na.action=na.exclude)
        timemat[l, (1:length(surv.fit$time))] <- surv.fit$time
        probmat[l, (1:length(surv.fit$surv))] <- surv.fit$surv
        CV.lhr[l] <- 0
        CV.lrt[l] <- 0
        CV.cer[l] <- 1
        surv.formula <- (survival::Surv(CV.times, CV.status) ~ 1 + grpind1)
        coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
        CV.grp.lhr[l] <- coxobj$coef 
        CV.grp.lrt[l] <- survival::survdiff(surv.formula, na.action=na.exclude, rho=0)$chisq 
        predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
        CV.grp.cer[l] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(CV.times, CV.status))['C Index']
      } else  if ((sum(grpind1[wb]) != length(grpind1[wb])) && (sum(grpind1[wb]) != 0)) {
        if ((sum(boxind, na.rm=TRUE) != length(boxind[!is.na(boxind)])) && (sum(boxind, na.rm=TRUE) != 0)) {
          surv.fit <- survival::survfit(survival::Surv(CV.times[wb], CV.status[wb]) ~ 1 + grpind1[wb], na.action=na.exclude)
          timemat[l, (1:length(surv.fit$time))] <- surv.fit$time
          probmat[l, (1:length(surv.fit$surv))] <- surv.fit$surv
          surv.formula <- (survival::Surv(CV.times, CV.status) ~ 1 + boxind1)
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
          CV.lhr[l] <- coxobj$coef
          CV.lrt[l] <- survival::survdiff(surv.formula, na.action=na.exclude, rho=0)$chisq
          predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
          CV.cer[l] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(CV.times, CV.status))['C Index']
          surv.formula <- (survival::Surv(CV.times[wb], CV.status[wb]) ~ 1 + grpind1[wb])
          coxobj <- survival::coxph(surv.formula, singular.ok=TRUE, iter.max=1, na.action=na.exclude)
          CV.grp.lhr[l] <- coxobj$coef 
          CV.grp.lrt[l] <- survival::survdiff(surv.formula, na.action=na.exclude, rho=0)$chisq 
          predobj <- predict(object=coxobj, type="lp", reference="sample", na.action=na.exclude)
          CV.grp.cer[l] <- Hmisc::rcorr.cens(x=predobj, S=survival::Surv(CV.times[wb], CV.status[wb]))['C Index']
        } else {
          timemat[l, ] <- NA
          probmat[l, ] <- NA
          CV.lhr[l] <- NA
          CV.lrt[l] <- NA
          CV.cer[l] <- NA
          CV.grp.lhr[l] <- NA
          CV.grp.lrt[l] <- NA
          CV.grp.cer[l] <- NA
          ind.rem <- c(ind.rem, l)
        }
      } else {
        timemat[l, ] <- NA
        probmat[l, ] <- NA
        CV.lhr[l] <- NA
        CV.lrt[l] <- NA
        CV.cer[l] <- NA
        CV.grp.lhr[l] <- NA 
        CV.grp.lrt[l] <- NA
        CV.grp.cer[l] <- NA
        ind.rem <- c(ind.rem, l)
      }
    } else {
      stop("Invalid peeling criterion. Exiting ...\n\n")
    }
  }
  if (length(ind.rem) != CV.Lm) {
    endobj <- endpoints(ind=ind.rem, timemat=timemat, probmat=probmat, timeval=timeval, probval=probval)
    time.bar <- endobj$time.bar
    prob.bar <- endobj$prob.bar
    max.time.bar <- endobj$max.time.bar
    min.prob.bar <- endobj$min.prob.bar
    for (l in 1:CV.Lm) {
      if (!(l %in% ind.rem)) {
        CV.time.bar[l] <- time.bar[l]
        CV.prob.bar[l] <- prob.bar[l]
        CV.max.time.bar[l] <-  max.time.bar[l]
        CV.min.prob.bar[l] <- min.prob.bar[l]
      } else {
        CV.time.bar[l] <- NA
        CV.prob.bar[l] <- NA
        CV.max.time.bar[l] <- NA
        CV.min.prob.bar[l] <- NA
      }
    }
  } else {
    CV.time.bar <- rep(x=NA, times=CV.Lm)
    CV.prob.bar <- rep(x=NA, times=CV.Lm)
    CV.max.time.bar <- rep(x=NA, times=CV.Lm)
    CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
  }
  CV.maxsteps <- CV.Lm
  
  # Formating the results depending on successful cross-validated PRSP algorithm
  if (cvcriterion == "lhr") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.lhr)) || all(is.nan(CV.lhr))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by maximization of the LHR (between inbox and outbox test samples)
        success <- TRUE
      }
    } else if (peelcriterion == "grp") {
      if (all(is.na(CV.lhr)) || all(is.nan(CV.lhr)) || all(is.na(CV.grp.lhr)) || all(is.nan(CV.grp.lhr))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else {
      # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
      success <- TRUE
    }
  } else if (cvcriterion == "lrt") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.lrt)) || all(is.nan(CV.lrt))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by maximization of the LRT (between inbox and outbox test samples)
        success <- TRUE
      }
    } else if (peelcriterion == "grp") {
      if (all(is.na(CV.lrt)) || all(is.nan(CV.lrt)) || all(is.na(CV.grp.lrt)) || all(is.nan(CV.grp.lrt))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {
        # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
        success <- TRUE
      }
    } else {
      # Cross-validated optimal length from all folds by minimization of the CER (between predicted and observed inbox test samples survival times)
      success <- TRUE
    }
  } else if (cvcriterion == "cer") {
    if (peelcriterion != "grp") {
      if (all(is.na(CV.cer)) || all(is.nan(CV.cer))) {
        success <- FALSE
        CV.maxsteps <- NA
        CV.support <- rep(x=NA, times=CV.Lm)
        CV.size <- rep(x=NA, times=CV.Lm)
        CV.lhr <- rep(x=NA, times=CV.Lm)
        CV.lrt <- rep(x=NA, times=CV.Lm)
        CV.cer <- rep(x=NA, times=CV.Lm)
        CV.grp.lhr <- rep(x=NA, times=CV.Lm)
        CV.grp.lrt <- rep(x=NA, times=CV.Lm)
        CV.grp.cer <- rep(x=NA, times=CV.Lm)
        CV.time.bar <- rep(x=NA, times=CV.Lm)
        CV.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.max.time.bar <- rep(x=NA, times=CV.Lm)
        CV.min.prob.bar <- rep(x=NA, times=CV.Lm)
        CV.boxcut <- matrix(data=NA, nrow=CV.Lm, ncol=p)
        CV.boxind <- matrix(NA, nrow=CV.Lm, ncol=n)
        CV.trace <- rep(x=NA, times=CV.Lm)
        CV.rules<- as.data.frame(matrix(data=NA, nrow=CV.Lm, ncol=p))
      } else {