diff --git a/R/PLNnetworkfit-class.R b/R/PLNnetworkfit-class.R index 312c16df..7aaf7312 100644 --- a/R/PLNnetworkfit-class.R +++ b/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) diff --git a/inst/case_studies/oaks_tree.R b/inst/case_studies/oaks_tree.R index 5b7983fa..c1814b83 100644 --- a/inst/case_studies/oaks_tree.R +++ b/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) diff --git a/inst/check/objective_zipln.R b/inst/check/objective_zipln.R new file mode 100644 index 00000000..5c621afb --- /dev/null +++ b/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