R/ZIPLNnetwork.R

In PLN-team/PLNmodels: Poisson Lognormal Models

#' Zero Inflated Sparse Poisson lognormal model for network inference
#'
#' Perform sparse inverse covariance estimation for the Zero Inflated Poisson lognormal model
#' using a variational algorithm. Iterate over a range of logarithmically spaced sparsity parameter values.
#' Use the (g)lm syntax to specify the model (including covariates and offsets).
#'
#' @inheritParams PLNnetwork
#' @param control a list-like structure for controlling the optimization, with default generated by [ZIPLNnetwork_param()]. See the associated documentation
#' for details.
#' @param zi a character describing the model used for zero inflation, either of
#' - "single" (default, one parameter shared by all counts)
#' - "col" (one parameter per variable / feature)
#' - "row" (one parameter per sample / individual).
#' If covariates are specified in the formula RHS (see details) this parameter is ignored.
#'
#' @details
#' Covariates for the Zero-Inflation parameter (using a logistic regression model) can be specified in the formula RHS using the pipe
#' (`~ PLN effect | ZI effect`) to separate covariates for the PLN part of the model from those for the Zero-Inflation part.
#' Note that different covariates can be used for each part.
#'
#' @return an R6 object with class [`ZIPLNnetworkfamily`]
#'
#' @include PLNnetworkfamily-class.R
#' @examples
#' data(trichoptera)
#' trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
#' myZIPLNs <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, zi = "single")
#' @seealso The classes [`ZIPLNfit`] and [`ZIPLNnetworkfamily`]
#' @export
ZIPLNnetwork <- function(formula, data, subset, weights, zi = c("single", "row", "col"), penalties = NULL, control = ZIPLNnetwork_param()) {

  ## extract the data matrices and weights
  data_ <- extract_model_zi(match.call(expand.dots = FALSE), parent.frame())
  control$ziparam <- ifelse((data_$zicovar), "covar", match.arg(zi))

  ## initialization
  if (control$trace > 0) cat("\n Initialization...")
  myPLN <- ZIPLNnetworkfamily$new(penalties, data_, control)

  ## optimization
  if (control$trace > 0) cat("\n Adjusting",
                             length(myPLN$penalties), "ZI-PLN with sparse inverse covariance estimation and",
                             control$ziparam, "specific parameter(s) in Zero inflation component.\n")
  myPLN$optimize(data_, control$config_optim)

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

#' Control of ZIPLNnetwork fit
#'
#' Helper to define list of parameters to control the ZIPLNnetwork fit. All arguments have defaults.
#'
#' @inheritParams PLNnetwork_param
#'
#' @inherit PLN_param details return
#' @details See [PLNnetwork_param()] for a full description of the optimization parameters. Note that some defaults values are different than those used in [PLNnetwork_param()]:
#' * "ftol_out" (outer loop convergence tolerance the objective function) is set by default to 1e-6
#' * "maxit_out" (max number of iterations for the outer loop) is set by default to 100
#'
#' @seealso [PLNnetwork_param()] and [PLN_param()]
#' @export
ZIPLNnetwork_param <- function(
    backend           = c("nlopt"),
    inception_cov     = c("full", "spherical", "diagonal"),
    trace             = 1      ,
    n_penalties       = 30     ,
    min_ratio         = 0.1    ,
    penalize_diagonal = TRUE   ,
    penalty_weights   = NULL   ,
    config_post       = list(),
    config_optim      = list(),
    inception         = NULL
) {

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

  ## post-treatment config
  config_pst <- config_post_default_PLNnetwork
  config_pst[names(config_post)] <- config_post
  config_pst$trace <- trace

  ## optimization config
  stopifnot(backend %in% c("nlopt"))
  stopifnot(config_optim$algorithm %in% available_algorithms_nlopt)
  config_opt <- config_default_nlopt
  config_opt$trace <- trace
  config_opt$ftol_out  <- 1e-6
  config_opt$maxit_out <- 50
  config_opt$approx_ZI <- TRUE
  config_opt[names(config_optim)] <- config_optim
  inception_cov <- match.arg(inception_cov)

  structure(list(
    backend           = backend          ,
    trace             = trace            ,
    inception_cov     = inception_cov    ,
    n_penalties       = n_penalties      ,
    min_ratio         = min_ratio        ,
    penalize_diagonal = penalize_diagonal,
    penalty_weights   = penalty_weights  ,
    jackknife         = FALSE            ,
    bootstrap         = 0                ,
    variance          = FALSE            ,
    config_post       = config_pst       ,
    config_optim      = config_opt       ,
    inception         = inception       ), class = "ZIPLNmodels_param")
}