R/PLNnetwork.R

In PLNmodels: Poisson Lognormal Models

#' Poisson lognormal model towards sparse network inference
#'
#' Fit the sparse inverse covariance variant of the Poisson lognormal with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets).
#'
#' @param formula an object of class "formula": a symbolic description of the model to be fitted.
#' @param data an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which lm is called.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process.
#' @param weights an optional vector of observation weights to be used in the fitting process.
#' @param penalties an optional vector of positive real number controlling the level of sparsity of the underlying network. if NULL (the default), will be set internally. See \code{PLNnetwork_param()} for additional tuning of the penalty.
#' @param control a list-like structure for controlling the optimization, with default generated by [PLNnetwork_param()]. See the corresponding documentation for details;
#'
#' @return an R6 object with class [`PLNnetworkfamily`], which contains
#' a collection of models with class [`PLNnetworkfit`]
#'
#' @rdname PLNnetwork
#' @examples
#' data(trichoptera)
#' trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
#' fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
#' @seealso The classes [`PLNnetworkfamily`] and [`PLNnetworkfit`], and the and the configuration function [PLNnetwork_param()].
#' @importFrom stats model.frame model.matrix model.response model.offset
#' @export
PLNnetwork <- function(formula, data, subset, weights, penalties = NULL, control = PLNnetwork_param()) {

  ## Temporary test for deprecated use of list()
  if (!inherits(control, "PLNmodels_param"))
    stop("We now use the function PLNnetwork_param() to generate the list of parameters that controls the fit:
    replace 'list(my_arg = xx)' by PLN_param(my_arg = xx) and see the documentation of PLNnetwork_param().")

  ## extract the data matrices and weights
  args <- extract_model(match.call(expand.dots = FALSE), parent.frame())

  ## Instantiate the collection of models
  if (control$trace > 0) cat("\n Initialization...")
  myPLN <- PLNnetworkfamily$new(penalties, args$Y, args$X, args$O, args$w, args$formula, control)

  ## Optimization
  if (control$trace > 0) cat("\n Adjusting", length(myPLN$penalties), "PLN with sparse inverse covariance estimation\n")
  if (control$trace) cat("\tJoint optimization alternating gradient descent and graphical-lasso\n")
  myPLN$optimize(control$config_optim)

  ## Post-treatments
  if (control$trace > 0) cat("\n Post-treatments")
  config_post <- config_post_default_PLNnetwork; config_post$trace <- control$trace
  myPLN$postTreatment(config_post)

  if (control$trace > 0) cat("\n DONE!\n")
  myPLN
}

#' Control of PLNnetwork fit
#'
#' Helper to define list of parameters to control the PLN fit. All arguments have defaults.
#'
#' @param backend optimization back used, either "nlopt" or "torch". Default is "nlopt"
#' @param config_optim a list for controlling the optimizer (either "nlopt" or "torch" backend). See details
#' @param trace a integer for verbosity.
#' @param n_penalties an integer that specifies the number of values for the penalty grid when internally generated. Ignored when penalties is non `NULL`
#' @param min_ratio the penalty grid ranges from the minimal value that produces a sparse to this value multiplied by `min_ratio`. Default is 0.1.
#' @param penalize_diagonal boolean: should the diagonal terms be penalized in the graphical-Lasso? Default is \code{TRUE}
#' @param penalty_weights either a single or a list of p x p matrix of weights (default filled with 1) to adapt the amount of shrinkage to each pairs of node. Must be symmetric with positive values.
#' @param inception Set up the parameters initialization: by default, the model is initialized with a multivariate linear model applied on
#'    log-transformed data, and with the same formula as the one provided by the user. However, the user can provide a PLNfit (typically obtained from a previous fit),
#'    which sometimes speeds up the inference.
#'
#' @return list of parameters configuring the fit.
#'
#' @details The list of parameters `config_optim` controls the optimizers. When "nlopt" is chosen the following entries are relevant
#' * "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".
#' * "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000
#' * "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8
#' * "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6
#' * "ftol_out" outer solver stops when an optimization step changes the objective function by less than xtol multiply by the absolute value of the parameter. Default is 1e-6
#' * "maxit_out" outer solver stops when the number of iteration exceeds out.maxit. Default is 50
#' * "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)
#' * "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)
#' * "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)
#'
#' When "torch" backend is used, with the following entries are relevant:
#' * "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000
#' * "ftol_out" outer solver stops when an optimization step changes the objective function by less than xtol multiply by the absolute value of the parameter. Default is 1e-6
#' * "maxit_out" outer solver stops when the number of iteration exceeds out.maxit. Default is 50
#' * "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8
#' * "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6
#'
#' @export
PLNnetwork_param <- function(
    backend           = "nlopt",
    trace             = 1      ,
    n_penalties       = 30     ,
    min_ratio         = 0.1    ,
    penalize_diagonal = TRUE   ,
    penalty_weights   = NULL   ,
    config_optim  = list(),
    inception         = NULL
) {

  if (!is.null(inception)) stopifnot(isPLNfit(inception))

  ## optimization config
  backend <- match.arg(backend)
  stopifnot(backend %in% c("nlopt", "torch"))
  if (backend == "nlopt") {
    stopifnot(config_optim$algorithm %in% available_algorithms_nlopt)
    config_opt <- config_default_nlopt
  }
  if (backend == "torch") {
    stopifnot(config_optim$algorithm %in% available_algorithms_torch)
    config_opt <- config_default_torch
  }
  config_opt$trace <- trace
  config_opt$ftol_out  <- 1e-5
  config_opt$maxit_out <- 20
  config_opt[names(config_optim)] <- config_optim

  structure(list(
    backend           = backend          ,
    trace             = trace            ,
    n_penalties       = n_penalties      ,
    min_ratio         = min_ratio        ,
    penalize_diagonal = penalize_diagonal,
    penalty_weights   = penalty_weights  ,
    jackknife         = FALSE            ,
    bootstrap         = 0                ,
    variance          = TRUE             ,
    config_optim      = config_opt       ,
    inception         = inception       ), class = "PLNmodels_param")
}