mixOmics splsda - CV function for component and variable selection

CV_splsda.R

getBestFit <- function(x, y, ..., ncomp = 5, nfolds = 5, minvar = 1, maxvar = ncol(x), direction = ">", seed = 42, run = NULL){
  set.seed(seed)
  auc <- expand.grid(AUC = 0, ncomp = 1:ncomp, nvar = minvar:maxvar, stringsAsFactors = FALSE)
  folds <- sample.int(nfolds, size = nrow(x), replace = TRUE)
  for(i in minvar:maxvar){
    pred <- matrix(rep(NA,nrow(x)*ncomp), ncol = ncomp)
    for(fold in 1:nfolds){
      train <- folds != fold
      test <- folds == fold
      fit <- splsda(x[train,], y[train], keepX = rep(i,ncomp), ncomp = ncomp)
      nearzero <- fit$nzv$Position
      if(length(nearzero) != 0){
        fit <- splsda(x[train,-nearzero], y[train], keepX = rep(i,ncomp), ncomp = ncomp)
        class(fit) <- c("spls", "pls")
        pred[test,] <- predict(fit, x[test,-nearzero])$predict[,1,]
      }
      else {
        class(fit) <- c("spls", "pls")
        pred[test,] <- predict(fit, x[test,])$predict[,1,]
      }
    }
    for(comp in 1:ncomp)
      auc[auc$ncomp == comp & auc$nvar == i, "AUC"] <- roc(y, pred[,comp], direction = direction)$auc
    # message(i, " variables AUC = ", auc[i])
  }
  if(!is.null(run))
    auc <- cbind(auc,run)
  cbind(auc, nfolds = nfolds, ...)
}
library(parallel)
aucspar5comp <- mclapply(1:40, function(seed) getBestFit(x, y, ncomp = 5, minvar = 1, maxvar = 8, nfolds = 10, seed = seed, run = seed), mc.cores = 7 ) 

aucspardf5comp <- aucspar5comp %>% do.call("rbind",.)
ggplot(aucspardf5comp, aes(x = factor(nvar), y = AUC)) + geom_boxplot() + geom_jitter() + facet_grid(nfolds ~ ncomp)