merging master

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.3.1
+RoxygenNote: 7.3.2
 Depends: R (>= 3.6)
 LazyData: true
 biocViews:

NAMESPACE

@@ -103,7 +103,6 @@ importFrom(pscl,zeroinfl)
 importFrom(purrr,map)
 importFrom(purrr,map2)
 importFrom(purrr,map2_dbl)
-importFrom(purrr,map_chr)
 importFrom(purrr,map_dbl)
 importFrom(purrr,map_int)
 importFrom(purrr,reduce)

R/PLNfamily-class.R

@@ -15,7 +15,7 @@
 #' @rdname PLNfamily
 #' @include PLNfamily-class.R
 #' @importFrom R6 R6Class
-#' @importFrom purrr map reduce map_chr
+#' @importFrom purrr map reduce
 PLNfamily <-
   R6Class(
     classname = "PLNfamily",
@@ -182,7 +182,7 @@ PLNfamily <-
       #' @field convergence sends back a data frame with some convergence diagnostics associated with the optimization process (method, optimal value, etc)
       convergence = function() {
         res <- purrr::map(self$models, function(model) {
-          c(nb_param = model$nb_param, purrr::map_chr(model$optim_par, tail, 1))
+          c(nb_param = model$nb_param, purrr::map(model$optim_par, tail, 1))
         }) %>% purrr::reduce(rbind)
         data.frame(param = private$params, res, stringsAsFactors = FALSE)
       }

R/PLNfit-class.R

@@ -146,7 +146,7 @@ PLNfit <- R6Class(
 
         ## Error message if objective diverges
         if (!is.finite(loss$item())) {
-          stop(sprintf("The ELBO diverged during the optimization procedure.\nConsider using:\n* a different optimizer (current optimizer: %s)\n* a smaller learning rate (current rate: %.3f)\nwith `control = PLN_param(config_optim = list(algorithm = ..., lr = ...))`",
+          stop(sprintf("The ELBO diverged during the optimization procedure.\nConsider using:\n* a different optimizer (current optimizer: %s)\n* a smaller learning rate (current rate: %.3f)\nwith `control = PLN_param(backend = 'torch', config_optim = list(algorithm = ..., lr = ...))`",
                        config$algorithm, config$lr))
         }
 
@@ -158,7 +158,7 @@ PLNfit <- R6Class(
         ## Check for convergence
         #print (delta_f)
         if (delta_f < config$ftol_rel) status <- 3
-        if (delta_x < config$xtol_rel) status <- 4
+        #if (delta_x < config$xtol_rel) status <- 4
         if (status %in% c(3,4)) {
           objective <- objective[1:iterate + 1]
           break

R/PLNnetworkfit-class.R

@@ -60,6 +60,7 @@ PLNnetworkfit <- R6Class(
       args <- list(data   = data,
                    params = list(B = private$B, M = private$M, S = private$S),
                    config = config)
+      private$Sigma <- crossprod(private$M)/self$n + diag(colMeans(private$S**2), self$p, self$p)
       while (!cond) {
         iter <- iter + 1
         if (config$trace > 1) cat("", iter)

inst/case_studies/oaks_tree.R

@@ -96,7 +96,10 @@ factoextra::fviz_pca_biplot(
   labs(col = "distance (cm)") + scale_color_viridis_c()
 
 ## Network inference with sparce covariance estimation
-system.time(myPLNnets <- PLNnetwork(Abundance ~ 0 + tree + offset(log(Offset)), data = oaks, control = PLNnetwork_param(min_ratio = 0.05)))
+
+system.time(myZIPLNnets <- ZIPLNnetwork(Abundance ~ 0 + tree + offset(log(Offset)), zi = "single", data = oaks, control = ZIPLNnetwork_param(min_ratio = 0.1)))
+
+system.time(myPLNnets <- PLNnetwork(Abundance ~ 0 + tree + offset(log(Offset)), data = oaks, control = PLNnetwork_param(min_ratio = 0.1)))
 plot(myPLNnets)
 plot(getBestModel(myPLNnets, "EBIC"))
 # stability_selection(myPLNnets)

inst/check/objective_zipln.R

@@ -0,0 +1,13 @@
+library(PLNmodels)
+data(oaks)
+
+myPLN <- PLN(Abundance ~ 0 + tree + offset(log(Offset)), data = oaks)
+myZIPLN <- ZIPLN(Abundance ~ 0 + tree + offset(log(Offset)) | 0 + tree, data = oaks)
+
+## bricolage car quelques points explosent, sans conséquences sur la convergence finale
+obj_PLN <- -myPLN$optim_par$objective
+plot(obj_PLN, ## optimisation nlopt de la vraisemblance profilée
+     ylim = c(obj_PLN[1], # fonction objective sans les constantes
+              tail(obj_PLN,1)), log = "xy")
+
+plot(myZIPLN$optim_par$objective, type = "l", log = "xy") ## Les itérations du VEM: la vraisemblance estimée