Merge pull request #131 from PLN-team/fix-plnnetwork

Some Fix in plnnetwork (Issue #129 and #130)

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