Continuing development of ZIPLN models

DESCRIPTION

@@ -25,7 +25,7 @@ BugReports: https://github.com/pln-team/PLNmodels/issues
 License: GPL (>= 3)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Depends: R (>= 3.4)
 LazyData: true
 biocViews:
@@ -90,6 +90,7 @@ Collate:
     'ZIPLNfit-class.R'
     'ZIPLN.R'
     'ZIPLNfit-S3methods.R'
+    'ZIPLNnetwork.R'
     'barents.R'
     'import_utils.R'
     'mollusk.R'

NAMESPACE

@@ -7,12 +7,17 @@ S3method(coef,ZIPLNfit)
 S3method(fitted,PLNfit)
 S3method(fitted,PLNmixturefit)
 S3method(fitted,ZIPLNfit)
+S3method(getBestModel,Networkfamily)
 S3method(getBestModel,PLNPCAfamily)
 S3method(getBestModel,PLNmixturefamily)
 S3method(getBestModel,PLNnetworkfamily)
+S3method(getBestModel,ZIPLNnetworkfamily)
+S3method(getModel,Networkfamily)
 S3method(getModel,PLNPCAfamily)
 S3method(getModel,PLNmixturefamily)
 S3method(getModel,PLNnetworkfamily)
+S3method(getModel,ZIPLNnetworkfamily)
+S3method(plot,Networkfamily)
 S3method(plot,PLNLDAfit)
 S3method(plot,PLNPCAfamily)
 S3method(plot,PLNPCAfit)
@@ -21,6 +26,8 @@ S3method(plot,PLNmixturefamily)
 S3method(plot,PLNmixturefit)
 S3method(plot,PLNnetworkfamily)
 S3method(plot,PLNnetworkfit)
+S3method(plot,ZIPLNfit_sparse)
+S3method(plot,ZIPLNnetworkfamily)
 S3method(predict,PLNLDAfit)
 S3method(predict,PLNfit)
 S3method(predict,PLNmixturefit)
@@ -48,6 +55,8 @@ export(PLNnetwork)
 export(PLNnetwork_param)
 export(ZIPLN)
 export(ZIPLN_param)
+export(ZIPLNnetwork)
+export(ZIPLNnetwork_param)
 export(coefficient_path)
 export(compute_PLN_starting_point)
 export(compute_offset)
@@ -77,6 +86,7 @@ importFrom(corrplot,corrplot)
 importFrom(future.apply,future_lapply)
 importFrom(future.apply,future_sapply)
 importFrom(glassoFast,glassoFast)
+importFrom(grDevices,rgb)
 importFrom(grid,nullGrob)
 importFrom(grid,textGrob)
 importFrom(gridExtra,arrangeGrob)

R/PLNnetwork.R

@@ -1,6 +1,7 @@
 #' 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).
+#' Fit the sparse inverse covariance variant of the Poisson lognormal with a variational algorithm
+#' for a collection of sparsity parameter values distributed on a log scale. 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.
@@ -12,7 +13,6 @@
 #' @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)
@@ -28,16 +28,16 @@ PLNnetwork <- function(formula, data, subset, weights, penalties = NULL, control
     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())
+  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, args$Y, args$X, args$O, args$w, args$formula, control)
+  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(control$config_optim)
+  myPLN$optimize(data_, control$config_optim)
 
   ## Post-treatments
   if (control$trace > 0) cat("\n Post-treatments")
@@ -122,5 +122,6 @@ PLNnetwork_param <- function(
     variance          = TRUE             ,
     config_post       = config_pst       ,
     config_optim      = config_opt       ,
+### TODO CHECK: Why two inceptive model (cov and not ?)
     inception         = inception       ), class = "PLNmodels_param")
 }

R/PLNnetworkfamily-S3methods.R

@@ -6,11 +6,9 @@
 
 
 ## Auxiliary functions to check the given class of an objet
-isPLNnetworkfamily <- function(Robject) {inherits(Robject, "PLNnetworkfamily")}
+isNetworkfamily <- function(Robject) {inherits(Robject, "Networkfamily")}
 
-#' Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of PLNnetwork fits (a [`PLNnetworkfamily`])
-#'
-#' @name plot.PLNnetworkfamily
+#' Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (either [`PLNnetworkfamily`] or [`ZIPLNnetworkfamily`])
 #'
 #' @inheritParams plot.PLNfamily
 #' @inherit plot.PLNfamily return details
@@ -22,7 +20,6 @@ isPLNnetworkfamily <- function(Robject) {inherits(Robject, "PLNnetworkfamily")}
 #' @param log.x logical: should the x-axis be represented in log-scale? Default is `TRUE`.
 #' @param stability scalar: the targeted level of stability in stability plot. Default is .9.
 #'
-#'
 #' @examples
 #' data(trichoptera)
 #' trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
@@ -34,14 +31,13 @@ isPLNnetworkfamily <- function(Robject) {inherits(Robject, "PLNnetworkfamily")}
 #' (with \code{type = 'stability'}) or the evolution of the criteria of the different models considered
 #' (with \code{type = 'criteria'}, the default).
 #' @export
-plot.PLNnetworkfamily <-
-  function(x,
+plot.Networkfamily <- function(x,
            type     = c("criteria", "stability", "diagnostic"),
            criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
            reverse = FALSE,
            log.x    = TRUE,
            stability = 0.9, ...) {
-  stopifnot(isPLNnetworkfamily(x))
+  stopifnot(isNetworkfamily(x))
   type <- match.arg(type)
   if (type == "criteria")
     p <- x$plot(criteria, reverse)
@@ -53,22 +49,45 @@ plot.PLNnetworkfamily <-
   p
 }
 
-#' @describeIn getModel Model extraction for [`PLNnetworkfamily`]
+#' @describeIn plot.Networkfamily Display various outputs associated with a collection of network fits
+#' @export
+plot.PLNnetworkfamily <- plot.Networkfamily
+
+#' @describeIn plot.Networkfamily Display various outputs associated with a collection of network fits
 #' @export
-getModel.PLNnetworkfamily <- function(Robject, var, index = NULL) {
-  stopifnot(isPLNnetworkfamily(Robject))
+plot.ZIPLNnetworkfamily <- plot.Networkfamily
+
+#' @describeIn getModel Model extraction for [`PLNnetworkfamily`] or [`ZIPLNnetworkfamily`]
+#' @export
+getModel.Networkfamily <- function(Robject, var, index = NULL) {
+  stopifnot(isNetworkfamily(Robject))
   Robject$getModel(var, index)
 }
 
-#' @describeIn getBestModel Model extraction for [`PLNnetworkfamily`]
+#' @describeIn getModel Model extraction for [`PLNnetworkfamily`]
 #' @export
-getBestModel.PLNnetworkfamily <- function(Robject, crit = c("BIC", "EBIC", "StARS"), ...) {
-  stopifnot(isPLNnetworkfamily(Robject))
+getModel.PLNnetworkfamily <- getModel.Networkfamily
+
+#' @describeIn getModel Model extraction for [`ZIPLNnetworkfamily`]
+#' @export
+getModel.ZIPLNnetworkfamily <- getModel.Networkfamily
+
+#' @describeIn getBestModel Model extraction for [`PLNnetworkfamily`] or [`ZIPLNnetworkfamily`]
+#' @export
+getBestModel.Networkfamily <- function(Robject, crit = c("BIC", "EBIC", "StARS"), ...) {
+  stopifnot(isNetworkfamily(Robject))
   stability <- list(...)[["stability"]]
   if (is.null(stability)) stability <- 0.9
   Robject$getBestModel(match.arg(crit), stability)
 }
 
+#' @describeIn getBestModel Model extraction for [`PLNnetworkfamily`]
+#' @export
+getBestModel.PLNnetworkfamily <- getBestModel.Networkfamily
+
+#' @describeIn getBestModel Model extraction for [`ZIPLNnetworkfamily`]
+#' @export
+getBestModel.ZIPLNnetworkfamily <- getBestModel.Networkfamily
 
 #' Extract the regularization path of a PLNnetwork fit
 #'
@@ -85,7 +104,7 @@ getBestModel.PLNnetworkfamily <- function(Robject, crit = c("BIC", "EBIC", "StAR
 #' head(coefficient_path(fits))
 #' @export
 coefficient_path <- function(Robject, precision = TRUE, corr = TRUE) {
-  stopifnot(isPLNnetworkfamily(Robject))
+  stopifnot(isNetworkfamily(Robject))
   Robject$coefficient_path(precision, corr)
 }
 
@@ -112,16 +131,15 @@ coefficient_path <- function(Robject, precision = TRUE, corr = TRUE) {
 #' }
 #' @export
 stability_selection <- function(Robject, subsamples = NULL, control = PLNnetwork_param(), force = FALSE) {
-  stopifnot(isPLNnetworkfamily(Robject))
+  stopifnot(isNetworkfamily(Robject))
+  if (inherits(Robject, "ZIPLNnetworkfamily")) control <- ZIPLNnetwork_param()
   if (force || anyNA(Robject$stability)) {
     Robject$stability_selection(subsamples, control)
   } else {
     message("Previous stability selection detected. Use \"force = TRUE\" to recompute it.")
   }
 }
 
-
-
 #' Extract edge selection frequency in bootstrap subsamples
 #'
 #' @description Extracts edge selection frequency in networks reconstructed from bootstrap subsamples
@@ -162,7 +180,7 @@ extract_probs <- function(Robject, penalty = NULL, index = NULL,
                           crit = c("StARS", "BIC", "EBIC"),
                           format = c("matrix", "vector"),
                           tol = 1e-5) {
-  stopifnot(isPLNnetworkfamily(Robject))
+  stopifnot(isNetworkfamily(Robject))
   ## Check if stability selection has been performed
   stab_path <- Robject$stability_path
   if (is.null(stab_path)) {