optim in PLNnetwork should not start by final PLN cov, but current es…

…timate (M'M + diag(S))

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