R/PLNnetwork.R

In PLNmodels: Poisson Lognormal Models

#' 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 PLN formula data subset weights
#' @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 `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`]
#'
#' @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
  data_ <- 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, data_, 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(data_, control$config_optim)

  ## Post-treatments
  if (control$trace > 0) cat("\n Post-treatments")
  myPLN$postTreatment(control$config_post, control$config_optim)

  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 backend, either `"builtin"` (Newton, default) or
#'   `"nlopt"` (CCSAQ) or `"torch"`. The default combines `"builtin"` with `maxit_ve = 1` and
#'   `inception_niter = 5` (see `maxit_ve` and `inception_backend`): this consistently finds a
#'   better ELBO than plain `"nlopt"`, at essentially the same speed. Without a good inception,
#'   `"builtin"` alone (`maxit_ve = NULL`) can converge to a poor basin on large datasets — use
#'   `"nlopt"` if you want to opt out of the whole combination.
#' @param inception_cov Covariance structure used for the inception PLN:
#'   `"full"` (default), `"diagonal"` or `"spherical"`. Non-full structures are now
#'   fully supported: when `inception_cov != "full"`, the penalty grid is built from the
#'   empirical covariance of latent residuals \eqn{M - XB} (a full-rank proxy for \eqn{\Sigma}),
#'   avoiding the broken `max_pen = 0` that previously occurred with diagonal/spherical.
#' @param inception_backend character or `NULL` (default, i.e. same as `backend`). Backend for
#'   the inception PLN only; the penalty grid models always use `backend`.
#'   Ignored when `inception` is supplied by the user.
#' @param inception_niter integer or `NULL`. Limits the inception PLN to at most
#'   this many iterations (EM iterations for `"builtin"`, function evaluations × 10 for
#'   `"nlopt"`). Default is `5L` when `backend = "builtin"` (`NULL`, i.e. full convergence,
#'   otherwise): fewer iterations keep the latent mean M from over-converging toward the
#'   unconstrained optimum, which would make it harder to warm-start the sparse penalty models.
#'   Values above ~20 typically hurt. When `inception_cov != "full"` or `inception_niter` is set,
#'   the penalty grid uses the empirical residual covariance \eqn{crossprod(M - XB) / n}
#'   for `max_pen`.
#' @param maxit_ve integer or `NULL`. Maximum number of inner VE-step iterations
#'   per outer GLASSO alternation turn. Default is `1L` when `backend = "builtin"` (`NULL`, i.e.
#'   full convergence — `maxit_em` for `"builtin"`, `maxeval` for `"nlopt"` — otherwise).
#'   `maxit_ve = 1` implements a **partial E-step** (generalized EM): one Newton step per outer
#'   turn prevents over-convergence that causes oscillations with the GLASSO M-step — see
#'   `backend` for the full default combination and its benchmark.
#' @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: all weights equal to 1) to adapt the amount of shrinkage to each pairs of node. Must be symmetric with positive values.
#' @inheritParams PLN_param trace config_optim config_post inception
#'
#' @return list of parameters configuring the fit.
#' @inherit PLN_param details
#' @section Outer-loop optimization parameters:
#' `PLNnetwork_param()` adds two parameters controlling the alternating GLASSO/VEM loop:
#' * "ftol_em" outer alternating solver stops when the objective changes by less than ftol_em (relative). Default is 1e-5
#' * "maxit_em" outer alternating solver stops when the number of iterations exceeds maxit_em. Default is 20
#'
#' @seealso [PLN_param()]
#' @export
PLNnetwork_param <- function(
    backend           = c("builtin", "nlopt", "torch"),
    inception_cov     = c("full", "spherical", "diagonal"),
    inception_backend = NULL   ,
    inception_niter   = NULL   ,
    maxit_ve          = NULL   ,
    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(isPLNfit(inception))

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

  ## optimization config
  backend <- match.arg(backend)
  inception_cov <- match.arg(inception_cov)
  ## default combination for "builtin" (partial E-step + short inception, see ?PLNnetwork_param);
  ## only fills in values the user did not set explicitly.
  if (backend == "builtin") {
    if (is.null(maxit_ve))        maxit_ve        <- 1L
    if (is.null(inception_niter)) inception_niter <- 5L
  }
  if (!is.null(maxit_ve)) config_optim$maxit_ve <- as.integer(maxit_ve)
  config_opt <- make_config_optim(backend, config_optim, trace,
                                  extra = list(ftol_em = 1e-5, maxit_em = 20))

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