#' Feature Selection Optimization for block (s)PLS method
#'
#' This function identify the number of feautures to keep per component and thus by cluster in \code{mixOmics::block.spls} by optimizing the silhouette coefficient, which assesses the quality of clustering.
#'
#' @param X list of numeric matrix (or data.frame) with features in columns and samples in rows (with samples order matching in all data sets).
#' @param Y (optional) numeric matrix (or data.frame) with features in columns and samples in rows (same rows as \code{X}).
#' @param indY integer, to supply if Y is missing, indicates the position of the matrix response in the list \code{X}.
#' @param ncomp integer, number of component to include in the model
#' @param test.list.keepX list of integers with the same size as X. Each entry corresponds to the different keepX value to test for each block of \code{X}.
#' @param test.keepY only if Y is provideid. 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::block.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]{block.spls}}, \code{\link[timeOmics]{getCluster}}, \code{\link[timeOmics]{plotLong}}
#'
#' @examples
#' demo <- suppressWarnings(get_demo_cluster())
#' X <- list(X = demo$X, Z = demo$Z)
#' Y <- demo$Y
#' test.list.keepX <- list("X" = c(5,10,15,20), "Z" = c(2,4,6,8))
#' test.keepY <- c(2:5)
#'
#' # tuning
#' tune.block.spls <- tuneCluster.block.spls(X= X, Y= Y,
#' test.list.keepX= test.list.keepX,
#' test.keepY= test.keepY,
#' mode= "canonical")
#' keepX <- tune.block.spls$choice.keepX
#' keepY <- tune.block.spls$choice.keepY
#'
#' # final model
#' block.spls.res <- mixOmics::block.spls(X= X, Y= Y, keepX = keepX,
#' keepY = keepY, ncomp = 2, mode = "canonical")
#' # get clusters and plot longitudinal profile by cluster
#' block.spls.cluster <- getCluster(block.spls.res)
#' @import mixOmics
#' @importFrom dplyr left_join
#' @importFrom dplyr mutate
#' @importFrom dplyr filter
#' @export
tuneCluster.block.spls <- function(X, Y = NULL, indY = NULL, ncomp = 2,
test.list.keepX = NULL, test.keepY = NULL, ...){
#-- checking input parameters ---------------------------------------------#
#--------------------------------------------------------------------------#
#-- X
X <- validate_list_matrix_X(X)
#-- Y
if(!is.null(Y)){
Y <- validate_matrix_Y(Y)
} else {
indY <- validate_indY(indY, X)
}
#-- ncomp
ncomp <- validate_ncomp(ncomp = ncomp, X =X)
#-- keepX
test.list.keepX <- validate_test_list_keepX(test.keepX = test.list.keepX, X = X, ncomp = ncomp)
#-- keepY
test.keepY <- validate_test_keepY(test.keepY = test.keepY, Y = Y)
list.keepX.keepY <- test.list.keepX
if(!is.null(Y)) {
list.keepX.keepY$Y <- test.keepY
}
list.keepX.keepY <- expand.grid(list.keepX.keepY, stringsAsFactors = FALSE,
KEEP.OUT.ATTRS = FALSE)
#-- launch tuning --------------------------------------------------------#
#--------------------------------------------------------------------------#
#-- 0. set output object
result <- matrix(ncol = 3 + length(X),nrow = nrow(list.keepX.keepY)*ncomp) %>%
as.data.frame() %>%
purrr::set_names(c("comp", names(X), "pos", "neg"))
if(!is.null(Y)){
result <- result %>%
mutate("Y" = NA)
}
result.index <- 1
#--1. compute dissimilarity matrix for silhouette coef. (once and for all)
all_data <- X
all_data_name <- as.character(unlist(imap(X, ~rep(.y,ncol(.x)))))
if(is.null(all_data$Y) & !is.null(Y)){
all_data$Y <- Y
all_data_name <- c(all_data_name, rep("Y", ncol(Y)))
}
all_data <- do.call("cbind", all_data)
dmatrix <- dmatrix.spearman.dissimilarity(all_data)
cluster <- as.data.frame(list("feature" = rownames(dmatrix),
"block" = all_data_name))
#--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 <- lapply(list.keepX.keepY, function(x) rep(min(x), ncomp))
for(block in names(list.keepX.keepY)){
kX[[block]][comp] <- list.keepX.keepY[index.list.kX.kY, block]
result[result.index, block] = list.keepX.keepY[index.list.kX.kY, block]
}
if(!is.null(Y)){
kY <- kX[["Y"]]
kX[["Y"]] <- NULL
#--3. run spls
block.spls.res <- mixOmics::block.spls(X = X, Y = Y, ncomp = ncomp,
keepX = kX, keepY = kY, ...)
} else {
#--3. run spls
block.spls.res <- mixOmics::block.spls(X = X, ncomp = ncomp,
keepX = kX, indY = indY, ...)
}
#--4. extract clusters
tmp.cluster <- getCluster(block.spls.res)
tmp.cluster <- suppressWarnings(dplyr::left_join(cluster, tmp.cluster[c(1,4)],
by = c("feature"="molecule"))) %>%
mutate(cluster = as.numeric(as.character(cluster))) %>%
mutate(cluster = ifelse(is.na(cluster), 0, cluster))
#--5. compute silhouette
sil <- silhouette(dmatrix, tmp.cluster$cluster)
#--6. store
result[result.index, "comp"] <- comp
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.index <- result.index + 1
}
}
result <- list("silhouette" = result)
result[["ncomp"]] <- ncomp
result[["test.keepX"]] <- test.list.keepX
result[["test.keepY"]] <- test.keepY
result[["block"]] <- names(list.keepX.keepY)
class(result) <- "block.spls.tune.silhouette"
#-- 7. choice.keepX / choice.keepY
result[["slopes"]] <- tune.silhouette.get_slopes(result)
tmp <- tune.silhouette.get_choice_keepX(result) %>%
do.call(what = "rbind")
if(!is.null(Y)){ # choice keepY
result[["choice.keepY"]] <- tmp$Y
tmp <- dplyr::select(tmp, -Y)
}
result[["choice.keepX"]] <- as.list(tmp)
return(result)
}
#' @importFrom purrr imap set_names map_dfr
#' @importFrom dplyr filter mutate group_by summarise left_join n
#' @importFrom stringr str_split
tune.silhouette.get_slopes <- function(object){
stopifnot(class(object) %in% c("block.spls.tune.silhouette", "spls.tune.silhouette", "spca.tune.silhouette"))
# tune.silhouette is a data.frame (comp, X, Y, bock ..., pos, neg)
# tune.silhouette <- tune.block.spls$silhouette
block <- object$block
ncomp <- object$ncomp
if(is(object, "block.spls.tune.silhouette")){
coord <- object$test.keepX
if(!is.null(object[["test.keepY"]])){
coord[["Y"]] <- object$test.keepY
coord <- lapply(coord, sort)
}
}else if(is(object, "spls.tune.silhouette")){
coord <- list(X= sort(object$test.keepX),
Y= sort(object$test.keepY))
}else if(is(object, "spca.tune.silhouette")){
coord <- list(X= sort(object$test.keepX))
}
# get all points
all.points <- unique(object$silhouette[block])
# define neighbours
neighbourhood <- map_dfr(1:nrow(all.points), ~tune.silhouette.get_neighbours(coord = coord, all.points[.x,,drop=FALSE]))
# extract pos and neg and set it as named list for better performance
split_by_comp <- split(object$silhouette, f= as.factor(object$silhouette$comp))
names_split_by_comp <- lapply(split_by_comp,
function(x) as.vector(apply(x[block], 1, function(y) paste(y, collapse = "_"))))
POS <- purrr::imap(split_by_comp, ~ purrr::set_names(.x[["pos"]], names_split_by_comp[[.y]]))
NEG <- purrr::imap(split_by_comp, ~ purrr::set_names(.x[["neg"]], names_split_by_comp[[.y]]))
# compute slope (pos / neg by comp)
slopes <- purrr::map_dfr(as.character(1:ncomp), ~{
cbind(neighbourhood,
"origin.pos" = POS[[.x]][neighbourhood$origin],
"destination.pos" = POS[[.x]][neighbourhood$destination],
"origin.neg" = NEG[[.x]][neighbourhood$origin],
"destination.neg" = NEG[[.x]][neighbourhood$destination],
"comp" = as.numeric(.x))})
slopes <- slopes %>%
dplyr::filter(origin!=destination) %>%
dplyr::mutate("slope.pos" = tune.silhouette.get_slopes_coef(x1 = lapply(stringr::str_split(.$origin, "_"), as.numeric),
x2 = lapply(stringr::str_split(.$destination, "_"), as.numeric),
y1 = .$origin.pos,
y2 = .$destination.pos)) %>%
dplyr::mutate("slope.neg" = tune.silhouette.get_slopes_coef(x1 = lapply(stringr::str_split(.$origin, "_"), as.numeric),
x2 = lapply(stringr::str_split(.$destination, "_"), as.numeric),
y1 = .$origin.neg,
y2 = .$destination.neg)) %>%
dplyr::mutate("distance_from_origin" = tune.silhouette.distance_from_origin(x1 = lapply(stringr::str_split(.$origin, "_"), as.numeric)))
# cumpute SD by comp and direction
SD <- slopes %>%
dplyr::group_by(comp, direction) %>%
dplyr::summarise(sd.pos = sd_new(slope.pos, na.rm = TRUE),
sd.neg = sd_new(slope.neg, na.rm = TRUE),
mean.pos = mean(slope.pos, na.rm = TRUE),
mean.neg = mean(slope.neg, na.rm = TRUE), N=dplyr::n())
# add Pval for signif slopes
slopes <- slopes %>%
dplyr::left_join(SD, by = c("direction", "comp")) %>%
# pos
dplyr::mutate(Z_score.pos = ifelse(.$sd.pos == 0,0,(.$slope.pos - .$mean.pos)/.$sd.pos)) %>%
dplyr::mutate(Pval.pos = ifelse(Z_score.pos >= 0,1-pnorm(Z_score.pos), pnorm(Z_score.pos))) %>%
# neg
dplyr::mutate(Z_score.neg = ifelse(.$sd.neg == 0,0,(.$slope.neg - .$mean.neg)/.$sd.neg)) %>%
dplyr::mutate(Pval.neg = ifelse(Z_score.neg >= 0,1-pnorm(Z_score.neg), pnorm(Z_score.neg)))
return(slopes)
}
#' @importFrom purrr map2
tune.silhouette.get_neighbours <- function(coord, point){
# return forward neighbourhood of a point
# valid point in coord
stopifnot(all(names(coord) == names(point)))
neighbour.max <- point
index <- purrr::map2(coord, point, ~which(.x == .y))
# get max neighbours
for(i in names(index)){
if(!(length(coord[[i]]) == index[[i]])){
neighbour.max[[i]] <- coord[[i]][index[[i]]+1]
}
} # can be simplified
# get all close possible neighbours
neighbourhood <- as.list(rbind(as.data.frame(point), as.data.frame(neighbour.max))) %>%
lapply(unique) %>%
expand.grid(stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)
# get direction for standard deviation computation
direction <- vector(mode = "numeric", length = nrow(neighbourhood))
for(i in 1:nrow(neighbourhood)){
direction[i] <- purrr::map2(coord, neighbourhood[i,], ~which(.x == .y)) %>%
purrr::map2(index, ~{.x-.y}) %>%
paste(collapse = "_")
}
neighbourhood["direction"] <- direction
neighbourhood["origin"] <- paste(point, collapse = "_")
neighbourhood["destination"] <- apply(neighbourhood[names(coord)], 1,
function(x) paste(x, collapse = "_"))
return(neighbourhood)
}
#' @importFrom purrr map_dbl
tune.silhouette.get_slopes_coef <- function(x1,x2,y1,y2){
stopifnot(length(x1) == length(x2))
stopifnot(length(x2) == length(y1))
stopifnot(length(y1) == length(y2))
res <- purrr::map_dbl(seq_along(x1), ~{
euc.dist <- sqrt(sum((x1[[.x]] - x2[[.x]])^2));
ifelse(euc.dist == 0, 0, (y2[.x] - y1[.x])/euc.dist)})
return(res)
}
#' @importFrom purrr map_dbl
tune.silhouette.distance_from_origin <- function(x1){
euc.dist <- purrr::map_dbl(seq_along(x1), ~{ sqrt(sum((x1[[.x]])^2))})
return(euc.dist)
}
#' @importFrom dplyr select summarise left_join
#' @importFrom tidyr gather
#' @importFrom magrittr %>%
tune.silhouette.get_choice_keepX <- function(tune.block.spls){
# from slopes, keep useful columns and remove NAs.
slopes <- tune.block.spls$slopes
tmp <- slopes %>%
dplyr::select(c(tune.block.spls$block, comp, direction, Pval.pos, Pval.neg, distance_from_origin)) %>%
tidyr::gather(Pval.dir, Pval.value, -c(tune.block.spls$block, comp, direction, distance_from_origin))
# for each comp, arrange by Pvalue and distance from origin and get first result
MIN <- split(tmp, f=tmp$comp) %>%
na.omit() %>%
lapply(function(x) x %>%
dplyr::arrange(Pval.value, distance_from_origin) %>%
.[1,] %>%
dplyr::select(tune.block.spls$block))
return(MIN)
}
Datasets
Models
