R/bs.R

In coglasso: Collaborative Graphical Lasso - Multi-Omics Network Reconstruction

#' Build multiple networks and select the best one from a multi-omics data set
#'
#' `bs()` wraps the two main functions of the package in a single one:
#' [coglasso()], to build multiple multi-omics networks, and [select_coglasso()]
#' to select the best one according to the chosen criterion. 
#' 
#' When using `bs()`, first, [coglasso()] estimates multiple multi-omics networks
#' with the algorithm *collaborative graphical lasso*, one for each combination
#' of input values for the hyperparameters \eqn{\lambda_w}, \eqn{\lambda_b} and 
#' \eqn{c}. Then, [select_coglasso()] selects the best combination of 
#' hyperparameters given to `coglasso()` according to the selected model 
#' selection method. The three availble options that can be set for the argument
#' `method` are  "xstars", "xestars" and "ebic". For more information on these 
#' selection methods, visit the help page of [select_coglasso()].
#'
#' @encoding UTF-8
#' @inherit coglasso
#' @inherit select_coglasso
#' @param verbose Print information regarding the network building and the 
#' network selection processes.
#' 
#' @return `bs()` returns an object of `S3` class `select_coglasso` containing 
#'   several elements. The most 
#'   important is probably `sel_adj`, the adjacency matrix of the 
#'   selected network. Some output elements depend on the chosen model selection
#'   method.\cr
#'   These elements are always returned, and they are the result of network 
#'   estimation with [coglasso()]:
#' \itemize{
#'   \item `loglik` is a numerical vector containing the \eqn{log} likelihoods of all
#'   the estimated networks.
#'   \item `density` is a numerical vector containing a measure of the density of all
#'   the estimated networks.
#'   \item `df` is an integer vector containing the degrees of freedom of all the
#'   estimated networks.
#'   \item `convergence` is a binary vector containing whether a network was
#'   successfully estimated for the given combination of hyperparameters or not.
#'   \item `path` is a list containing the adjacency matrices of all the estimated
#'   networks.
#'   \item `icov` is a list containing the inverse covariance matrices of all the
#'   estimated networks.
#'   \item `nexploded` is the number of combinations of hyperparameters for which
#'   `coglasso()` failed to converge.
#'   \item `data` is the input multi-omics data set.
#'   \item `hpars` is the ordered table of all the combinations of hyperparameters
#'   given as input to `bs()`, with \eqn{\alpha(\lambda_w+\lambda_b)}
#'   being the key to sort rows.
#'   \item `lambda_w`, `lambda_b`, and `c` are numerical vectors with, 
#'   respectively, all the \eqn{\lambda_w}, \eqn{\lambda_b}, and \eqn{c} values
#'   `bs()` used.
#'   \item `p` is the vector with the number of variables for each omic layer of the 
#'   data set.
#'   \item `D` is the number of omics layers in the data set.
#'   \item `cov` optional, returned when `cov_output` is TRUE, is a list containing
#'   the variance-covariance matrices of all the estimated networks.
#' }
#'   These elements are returned by all selection methods available:
#' \itemize{
#'   \item `sel_index_c`, `sel_index_lw` and `sel_index_lb` are the indexes of the
#'     final selected parameters \eqn{c}, \eqn{\lambda_w} and \eqn{\lambda_b}
#'     leading to the most stable sparse network.
#'   \item `sel_c`, `sel_lambda_w` and `sel_lambda_b` are the final selected
#'     parameters \eqn{c}, \eqn{\lambda_w} and \eqn{\lambda_b} leading to the most
#'     stable sparse network.
#'   \item `sel_adj` is the adjacency matrix of the final selected network.
#'   \item `sel_density` is the density of the final selected network.
#'   \item `sel_icov` is the inverse covariance matrix of the final selected network.
#'   \item `sel_cov` optional, given only when `coglasso()` was called with 
#'   `cov_output = TRUE`. It is the covariance matrix associated with the final 
#'   selected network.
#'   \item `call` is the matched call.
#'   \item `method` is the chosen model selection method.
#' }
#'   These are the additional elements returned when choosing "xestars" or "xstars":
#' \itemize{
#'   \item `merge` is the "merged" adjacency matrix, the average of all the adjacency 
#'     matrices estimated across all the different subsamples for the selected 
#'     combination of \eqn{\lambda_w}, \eqn{\lambda_b}, and \eqn{c} values in the
#'     last path explored before convergence. Each entry is a measure of how 
#'     recurrent the corresponding edge is across the subsamples.
#'   \item `variability_lw`, `variability_lb` and `variability_c` are numeric vectors
#'     of as many items as the number of \eqn{\lambda_w}, \eqn{\lambda_b}, and 
#'     \eqn{c} values explored. Each item is the variability of the network 
#'     estimated for the corresponding hyperparameter value, keeping the other two 
#'     hyperparameters fixed to their selected value.
#'   \item `sel_variability` is the variability of the final selected network.
#' }
#'   These are the additional elements returned when choosing "ebic":
#' \itemize{
#'   \item `ebic_scores` is a numerical vector containing the eBIC scores for all the
#'     hyperparameter combination.
#' }
#' @export
#'
#' @examples
#' # Suggested usage: give the input data set, set the values for `p` and the 
#' # number of hyperparameters to explore (to choose how extensively to explore 
#' # the possible hyperparameters). Then, let the default behavior do the rest:
#' 
#' sel_mo_net <- bs(multi_omics_sd_micro, p = c(4, 2), nlambda_w = 3, 
#'                  nlambda_b = 3, nc = 3, verbose = FALSE)
#' 
bs <- function(data, p = NULL, pX = lifecycle::deprecated(), 
               lambda_w = NULL, lambda_b = NULL, c = NULL, 
               nlambda_w = NULL, nlambda_b = NULL, nc = NULL, 
               lambda_w_max = NULL, lambda_b_max = NULL, c_max = NULL, 
               lambda_w_min_ratio = NULL, lambda_b_min_ratio = NULL, 
               c_min = NULL, icov_guess = NULL, cov_output = FALSE, 
               lock_lambdas = FALSE, method = "xestars", stars_thresh = 0.1, 
               stars_subsample_ratio = NULL, rep_num = 20, max_iter = 10,
               old_sampling = FALSE, ebic_gamma = 0.5, 
               verbose = TRUE){
  
  if (lifecycle::is_present(pX)) {
    lifecycle::deprecate_warn("1.1.0", "bs(pX)", "bs(p)")
    p <- pX
  }
  
  call <- match.call()
  
  cg <- coglasso(data = data, p = p, pX = pX, lambda_w = lambda_w, 
                 lambda_b = lambda_b, c = c, nlambda_w = nlambda_w, 
                 nlambda_b = nlambda_b, nc = nc, lambda_w_max = lambda_w_max, 
                 lambda_b_max = lambda_b_max, c_max = c_max, 
                 lambda_w_min_ratio = lambda_w_min_ratio, 
                 lambda_b_min_ratio = lambda_b_min_ratio, 
                 c_min = c_min, icov_guess = icov_guess, 
                 cov_output = cov_output, lock_lambdas = lock_lambdas, 
                 verbose = verbose)
  
  cg <- select_coglasso(cg, method = method, stars_thresh = stars_thresh, 
                        stars_subsample_ratio = stars_subsample_ratio, 
                        rep_num = rep_num, max_iter = max_iter, 
                        old_sampling = old_sampling, 
                        ebic_gamma = ebic_gamma, verbose = verbose)
  
  cg$call <- call
  
  cg
}