R/tuneCluster.spls.R

In abodein/timeOmics: Time-Course Multi-Omics data integration

#' Feature Selection Optimization for sPLS method 
#' 
#' This function identify the number of feautures to keep per component and thus by cluster in \code{mixOmics::spls} by optimizing the silhouette coefficient, which assesses the quality of clustering.
#' 
#' @param X numeric matrix (or data.frame) with features in columns and samples in rows
#' @param Y numeric matrix (or data.frame) with features in columns and samples in rows (same rows as \code{X})
#' @param ncomp integer, number of component to include in the model
#' @param test.keepX vector of integer containing the different value of keepX to test for block \code{X}.
#' @param test.keepY vector of integer containing the different value of keepY to test for block \code{Y}.
#' @param ... other parameters to be included in the spls model (see \code{mixOmics::spls})
#' 
#' @return 
#' \item{silhouette}{silhouette coef. computed for every combinasion of keepX/keepY}
#' \item{ncomp}{number of component included in the model}
#' \item{test.keepX}{list of tested keepX}
#' \item{test.keepY}{list of tested keepY}
#' \item{block}{names of blocks}
#' \item{slopes}{"slopes" computed from the silhouette coef. for each keepX and keepY, used to determine the best keepX and keepY}
#' \item{choice.keepX}{best \code{keepX} for each component}
#' \item{choice.keepY}{best \code{keepY} for each component}
#' 
#'
#' @details
#' For each component and for each keepX/keepY value, a spls is done from these parameters.
#' Then the clustering is performed and the silhouette coefficient is calculated for this clustering.
#'
#' We then calculate "slopes" where keepX/keepY are the coordinates and the silhouette is the intensity.
#' A z-score is assigned to each slope.
#' We then identify the most significant slope which indicates a drop in the silhouette coefficient and thus a deterioration of the clustering.
#'
#' 
#' @seealso 
#' \code{\link[mixOmics]{spls}}, \code{\link[timeOmics]{getCluster}}, \code{\link[timeOmics]{plotLong}}
#'
#' @examples
#' demo <- suppressWarnings(get_demo_cluster())
#' X <- demo$X
#' Y <- demo$Y
#' 
#' # tuning
#' tune.spls <- tuneCluster.spls(X, Y, ncomp= 2, test.keepX= c(5,10,15,20), test.keepY= c(2,4,6))
#' keepX <- tune.spls$choice.keepX
#' keepY <- tune.spls$choice.keepY
#' 
#' # final model
#' spls.res <- mixOmics::spls(X, Y, ncomp= 2, keepX= keepX, keepY= keepY)
#' 
#' # get clusters and plot longitudinal profile by cluster
#' spls.cluster <- getCluster(spls.res)
#' plotLong(spls.res)
#' 
#'
#' @export
#' @import mixOmics
#' @importFrom dplyr left_join
#' @importFrom dplyr mutate
#' @importFrom dplyr filter

tuneCluster.spls <- function(X, Y, ncomp = 2, test.keepX = rep(ncol(X), ncomp),
                             test.keepY = rep(ncol(Y), ncomp), ...){
    #-- checking input parameters ------------  ---------------------------------#
    #--------------------------------------------------------------------------#

    #-- X
    X <- validate_matrix_X(X)
    Y <- validate_matrix_Y(Y)

    #-- ncomp
    ncomp <- validate_ncomp(ncomp, list(X,Y))

    #-- test.keepX
    test.keepX <- validate_test_keepX(test.keepX = test.keepX, X = X)

    #-- test.keepY
    test.keepY <- validate_test_keepY(test.keepY = test.keepY, Y = Y)

    list.keepX.keepY <- list("keepX" = test.keepX, "keepY" = test.keepY) %>%
        expand.grid(stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)


    #-- launch tuning  --------------------------------------------------------#
    #--------------------------------------------------------------------------#

    #-- 0. set output object
    result <- as.data.frame(matrix(ncol = 5, nrow = nrow(list.keepX.keepY)*ncomp))
    colnames(result) <- c("comp", "X", "Y", "pos", "neg")
    result.index <- 0

    #--1. compute dissimilarity matrix for silhouette coef. (once and for all)
    all_data <- cbind(X, Y)
    dmatrix <- dmatrix.spearman.dissimilarity(all_data)

    cluster <- as.data.frame(list("feature" = rownames(dmatrix),
                                  "block" = c(rep("X", ncol(X)), rep("Y", ncol(Y)))))

    #--2. tuning
    for(comp in 1:ncomp){
        for(index.list.kX.kY in 1:nrow(list.keepX.keepY)){
            # foreach comp, keepX and keepY of other comp is set to minimum
            kX <- rep(min(test.keepX), ncomp)
            kY <- rep(min(test.keepY), ncomp)
            kX[comp] <- list.keepX.keepY[index.list.kX.kY,"keepX"]
            kY[comp] <- list.keepX.keepY[index.list.kX.kY,"keepY"]

            #--3. run spls
            spls.res <- mixOmics::spls(X = X, Y = Y, ncomp =  ncomp, keepX = kX, keepY = kY, ...)

            #--4. extract clusters
            tmp.cluster <- getCluster(spls.res)
            tmp.cluster <- suppressWarnings(dplyr::left_join(cluster, tmp.cluster[c(1,4)],
                                                             by = c("feature"="molecule"))) %>%
                dplyr::mutate(cluster = as.numeric(as.character(cluster))) %>%
                dplyr::mutate(cluster = ifelse(is.na(cluster), 0, cluster))

            #--5. compute silhouette
            sil <- silhouette(dmatrix, tmp.cluster$cluster)

            #--6. store
            result.index <- result.index + 1
            result[result.index, "comp"] <- comp
            result[result.index, "X"] <- kX[comp]
            result[result.index, "Y"] <- kY[comp]
            # result[result.index, "pos"] <- sil$average.cluster  %>%
            #     dplyr::filter(cluster == comp) %>% pull(silhouette.coef)
            # result[result.index, "neg"] <- sil$average.cluster  %>%
            #     dplyr::filter(cluster == -comp) %>% pull(silhouette.coef)
            pos.res <-  sil$average.cluster  %>% 
                dplyr::filter(cluster == comp) %>%
                dplyr::pull(silhouette.coef)
            result[result.index, "pos"] <- ifelse(length(pos.res) == 0, NA, pos.res)
            neg.res <-  sil$average.cluster  %>%
                dplyr::filter(cluster == -comp) %>%
                dplyr::pull(silhouette.coef)
            result[result.index, "neg"] <- ifelse(length(neg.res) == 0, NA, neg.res)
        }
    }
    result <- list("silhouette" = result)
    result[["ncomp"]] <- ncomp
    result[["test.keepX"]] <- test.keepX
    result[["test.keepY"]] <- test.keepY
    result[["block"]] <- c("X", "Y")
    class(result) <- "spls.tune.silhouette"
    
    #-- 7. choice.keepX / choice.keepY
    result[["slopes"]] <- tune.silhouette.get_slopes(result)
    tmp <- tune.silhouette.get_choice_keepX(result) # choice keepX/keepY
    result[["choice.keepX"]] <- unlist(lapply(tmp, function(x) x$X))
    result[["choice.keepY"]] <- unlist(lapply(tmp, function(x) x$Y))
    return(result)
}