#' Get optimal number of components
#'
#' Compute the average silhouette coefficient for a given set of components on a mixOmics result.
#' Foreach given ncomp, the mixOmics method is performed with the sames arguments and the given `ncomp`.
#' Longitudinal clustering is performed and average silhouette coefficient is computed.
#'
#' @param object A mixOmics object of the class `pca`, `spca`, `mixo_pls`, `mixo_spls`, `block.pls`, `block.spls`
#'
#' @param max.ncomp integer, maximum number of component to include.
#' If no argument is given, `max.ncomp=object$ncomp`
#'
#' @param X a numeric matrix/data.frame or a list of data.frame for \code{block.pls}
#'
#' @param Y (only for \code{pls}, optional for \code{block.spls}) a numeric matrix, with the same nrow as \code{X}
#'
#' @param indY (optional and only for \code{block.pls}, if Y is not provided), an integer which indicates the position of the matrix response in the list X
#'
#' @param ... Other arguments to be passed to methods (pca, pls, block.pls)
#'
#' @return
#' \code{getNcomp} returns a list with class "ncomp.tune.silhouette" containing the following components:
#'
#' \item{ncomp}{a vector containing the tested ncomp}
#' \item{silhouette}{a vector containing the average silhouette coefficient by ncomp}
#' \item{dmatrix}{the distance matrix used to compute silhouette coefficient}
#'
#' @seealso
#' \code{\link{getCluster}}, \code{\link{silhouette}}, \code{\link[mixOmics]{pca}}, \code{\link[mixOmics]{pls}}, \code{\link[mixOmics]{block.pls}}
#'
#' @examples
#' # random input data
#' demo <- suppressWarnings(get_demo_cluster())
#'
#' # pca
#' pca.res <- mixOmics::pca(X=demo$X, ncomp = 5)
#' res.ncomp <- getNcomp(pca.res, max.ncomp = 4, X = demo$X)
#' plot(res.ncomp)
#'
#' # pls
#' pls.res <- mixOmics::pls(X=demo$X, Y=demo$Y)
#' res.ncomp <- getNcomp(pls.res, max.ncomp = 4, X = demo$X, Y=demo$Y)
#' plot(res.ncomp)
#'
#' # block.pls
#' block.pls.res <- suppressWarnings(mixOmics::block.pls(X=list(X=demo$X, Z=demo$Z), Y=demo$Y))
#' res.ncomp <- suppressWarnings(getNcomp(block.pls.res, max.ncomp = 4,
#' X=list(X=demo$X, Z=demo$Z), Y=demo$Y))
#' plot(res.ncomp)
#'
#' @export
#' @import mixOmics
getNcomp <- function(object, max.ncomp = NULL, X, Y = NULL, indY = NULL, ...){
#-- checking input parameters ---------------------------------------------#
#--------------------------------------------------------------------------#
#-- object
allowed_object = c("pca", "mixo_pls", "block.pls")
if(!any(class(object) %in% allowed_object)){
stop("invalid object, should be one of c(pca, mixo_pls, block.pls)")
}
#-- max.ncomp
if(is_almostInteger(max.ncomp)){
if (max.ncomp < 1)
stop("'max.ncomp' should be greater than 1")
if(is(object, "block.pls")){
if (max.ncomp > min(ncol(object$X[[1]]), nrow(object$X[[1]])))
stop("use smaller 'max.ncomp'")
} else {
if (max.ncomp > min(ncol(object$X), nrow(object$X)))
stop("use smaller 'max.ncomp'")
}
} else {
max.ncomp <- unique(object$ncomp)
}
#-- run #pca / pls / block.pls
ncomp.opt.res <- ncomp.silhouette(object, X, Y, max.ncomp, indY, ...)
to_return <- list()
to_return[["ncomp"]] <- c(0,1:max.ncomp)
to_return[["silhouette"]] <- c(0,ncomp.opt.res$silhouette.res)
to_return[["dmatrix"]] <- ncomp.opt.res$dmatrix
to_return[["choice.ncomp"]] <- to_return[["ncomp"]][which.max(to_return[["silhouette"]])]
class(to_return) <- "ncomp.tune.silhouette"
return(invisible(to_return))
}
ncomp.silhouette <- function(object, X = X, max.ncomp = max.ncomp, ...){
UseMethod("ncomp.silhouette")
}
#' @import mixOmics
ncomp.silhouette.pca <- function(object, X, Y, max.ncomp, indY, ...){
#-- check X
X <- validate_matrix_X(X)
#-- dmatrix
dmatrix <- dmatrix.spearman.dissimilarity(X)
silhouette.res <- vector(length = max.ncomp)
#-- iterative ncomp silhouette coef.
for(comp in 1:max.ncomp){
#-- mixo pca
mixo.res <- mixOmics::pca(X = X, ncomp = comp, ...)
#-- cluster
cluster.res <- getCluster(X = mixo.res)
# same names, same cluster
stopifnot(all(cluster.res$molecule == colnames(dmatrix)))
#-- silhouette
sil <- silhouette(dmatrix, cluster.res$cluster)
#-- store
silhouette.res[comp] <- sil$average
}
return(list(silhouette.res = silhouette.res, dmatrix = dmatrix))
}
#' @import mixOmics
ncomp.silhouette.mixo_pls <- function(object, X, Y, max.ncomp, indY, ...){
#-- check X
X <- validate_matrix_X(X)
#-- check Y
Y <- validate_matrix_Y(Y)
#-- dmatrix
dmatrix <- dmatrix.spearman.dissimilarity(cbind(X,Y))
silhouette.res <- vector(length = max.ncomp)
#-- iterative ncomp silhouette coef.
for(comp in 1:max.ncomp){
#-- mixo pls
mixo.res <- mixOmics::pls(X = X, Y=Y, ncomp = comp, ...)
#-- cluster
cluster.res <- getCluster(X = mixo.res)
# same names, same cluster
stopifnot(all(cluster.res$molecule == colnames(dmatrix)))
#-- silhouette
sil <- silhouette(dmatrix, cluster.res$cluster)
#-- store
silhouette.res[comp] <- sil$average
}
return(list(silhouette.res = silhouette.res, dmatrix = dmatrix))
}
#' @import mixOmics
ncomp.silhouette.block.pls <- function(object, X, Y, max.ncomp, indY, ...){
#-- check X
X <- validate_list_matrix_X(X)
data <- do.call("cbind", X)
#-- Y
if(!is.null(Y)){
Y <- validate_matrix_Y(Y)
dmatrix <- dmatrix.spearman.dissimilarity(cbind(data,Y))
indY <- NULL
} else {
indY <- validate_indY(indY = indY, X=X)
dmatrix <- dmatrix.spearman.dissimilarity(data)
}
silhouette.res <- vector(length = max.ncomp)
#-- iterative ncomp silhouette coef.
for(comp in 1:max.ncomp){
#-- mixo block.pls
if(is.null(indY)){
mixo.res <- mixOmics::block.pls(X = X, ncomp = comp, Y = Y, ...)
}else{
mixo.res <- mixOmics::block.pls(X = X, ncomp = comp, indY = indY, ...)
}
#-- cluster
cluster.res <- getCluster(X = mixo.res)
# same names, same cluster
stopifnot(all(cluster.res$molecule == colnames(dmatrix)))
#-- silhouette
sil <- silhouette(dmatrix, cluster.res$cluster)
#-- store
silhouette.res[comp] <- sil$average
}
return(list(silhouette.res = silhouette.res, dmatrix = dmatrix))
}
#' @export
#' @import ggplot2
plot.ncomp.tune.silhouette <- function(x, title = NULL, ...){
stopifnot(is(x, "ncomp.tune.silhouette"))
# check title
if(!is.character(title)){title = NULL}
data <- as.data.frame(list(ncomp = x$ncomp, silhouette = x$silhouette))
ggplot_df <- ggplot2::ggplot(data, aes(x=ncomp, y = silhouette)) + geom_line() + geom_point() +
geom_vline(xintercept = x$choice.ncomp, lty=2, col = "grey") +
theme_bw() +
xlab("Number of Principal Components") +
ylab("Average Silhouette Coefficient")
# add title
if(is.character(title)){
ggplot_df <- ggplot_df + ggtitle(title)
}
print(ggplot_df)
return(invisible(ggplot_df))
}
Datasets
Models
