NEWS.md

PLNmodels 1.3.2

Internal graphical Lasso, replacing glassoFast

• PLNnetwork() and ZIPLNnetwork() now use an internal graphical Lasso, a C++ port of the GLASSOFAST algorithm (Sustik and Calderhead, 2012) shared with the normalblockr package, instead of calling glassoFast::glassoFast(). glassoFast is no longer a dependency (it moves to Suggests, for tests only).
• The motivation is robustness: glassoFast's Fortran routine could loop forever on a nearly collapsed covariance matrix (entries of order 1e-8, as a rank-deficient residual covariance produces), in compiled code that no R-level timeout could stop. The new solver always terminates (non-finite input and zero-variance coordinates are rejected, the inner coordinate descent is bounded), reports non-convergence instead of hanging, and can be interrupted from R, by the user or by setTimeLimit()/R.utils::withTimeout(), which then raise their usual error.
• On ordinary input, results are those of glassoFast up to machine precision: identical supports along penalty paths and relative differences below 1e-15 on the precision matrix, for scalar as well as weighted penalties with an unpenalized diagonal; PLNnetwork() and ZIPLNnetwork() fits with the default "builtin" backend are unchanged (log-likelihoods within 1e-11). Fits with backend = "nlopt" can differ slightly along the path: that backend is sensitive to perturbations at the level of the last floating-point digit. Speed is the same or slightly better.
• The solver is exported as graphical_lasso(S, rho, thr, maxit, w_init, wi_init), returning w, wi, niter, converged, status and delta, with the same defaults as glassoFast and an optional warm start.
• It stops when it is cycling rather than converging. On an ill-conditioned covariance (ie rank-deficient, as PLNnetwork() produces when p ~ n) the sweeps settle into a small limit cycle: the convergence criterion stops decreasing and oscillates just above its threshold forever. glassoFast spends its whole 10000-sweep budget on these and reports success regardless. The cycle is now detected after 1000 sweeps without progress (tunable through stall_patience), and reported as status = "stalled".
• Non-convergence of the graphical Lasso is recorded in the fits' monitoring ($optim_par$glasso_nonconverged and $optim_par$glasso_stalled, counted along the alternating optimization). A warning is now raised only for the numerical failures ("degenerate", "inner_failure", "max_iter"), not for a stalled solve, which says something about the problem (too weak a penalty for a nearly rank-deficient covariance) rather than about the solution.
• Behaviour change in a degenerate case: when the covariance matrix has no off-diagonal mass, the (diagonal) precision matrix is now the correct 1 / (S_ii + rho_ii), where glassoFast returned 1 / max(rho_ii, 1.1e-16).

Model selection and display of networks

• The EBIC of PLNnetworkfit and ZIPLNfit_sparse is now the one of Foygel and Drton (2010), BIC - 2 gamma |E| log(p), which was designed for the graphical Lasso, instead of the Stirling approximation of Chen and Chen (2008)'s original criterion that was used so far. The tuning parameter is exposed as $ebic_gamma (default 0.5, 0 gives back the BIC), on a single fit or on a whole collection. EBIC values change, and so may the model selected by getBestModel("EBIC").
• The $density of a network is now |E| / (p (p - 1) / 2): it divided the edge count by p^2 rather than by the number of possible edges, and was thus understated by a factor (p - 1) / p.
• In the igraph output of plot() on a network fit, the opacity of an edge is now proportional to the strength of its partial correlation, as its width already was. Dense networks are no longer a solid blob of colour. The opacity of the weakest edge is the new edge.alpha argument (default 0.2, set it to 1 to restore uniformly opaque edges).

Bug fixes

• PLNPCA()'s variational bound carried a spurious p / 2 per observation, the entropy constant of the full-covariance models: a rank-q model's variational distribution is over the q-dimensional scores, and its q / 2 constant was already in the Kullback-Leibler term. loglik, BIC and ICL of every PLNPCAfit therefore decrease by n * p / 2. Rank selection is unchanged (the term does not depend on the rank), but a PLNPCA() fit is now comparable with a PLN() one, or with any other model. Reported by Nguyen Quang Huy (Actuarial Science Laboratory, National Economics University, Vietnam).

• PLNnetwork() and ZIPLNnetwork() no longer fail outright when one model along the penalty path cannot be fitted: the collection is truncated with a warning instead, and stability_selection() treats a model with no estimated network as fully unstable. A second line of defense on top of the solver above (thanks @rfriedman22, #176).

PLNmodels 1.3.1

Bug fix

• ZIPLNnetwork() fitted its inception with the regularization path's truncated optimizer which is appropriate for models that are warm-started from one another, but the inception starts from scratch and was given the same budget. The inception is now fitted with an untruncated optimizer, as PLNnetworkfamily.

PLNmodels 1.3.0

New backends and optimizers

• New built-in Newton optimizer (backend = "builtin") for PLN, ZIPLN and PLNnetwork: envelope-theorem Newton steps with strong Wolfe line search, no dependency on NLOPT. Substantially faster and more accurate than nlopt on large datasets with full covariance.

• PLNPCA's "builtin" backend rewritten as a profiled trust-region Newton: the variational block (M, S) is profiled out with a per-observation Newton VE-step, and the loadings (B, C) are optimised with a saddle-aware trust-region Newton on the resulting objective (analytic Schur-complement Hessian-vector products, Jacobi-preconditioned Steihaug-CG) — replacing the previous spectral projected-gradient "builtin". It reliably reaches a higher variational bound than "nlopt" at comparable-to-better speed; tuning keys cg_maxit, maxit_out, ftol_out, gtol, delta0 in config_optim (see ?PLNPCA_param). "nlopt" remains the default for PLNPCA.

• Backend defaults revisited package-wide, based on extensive benchmarking: PLN and PLNPCA keep "nlopt" (PLN now consistently faster thanks to profiled = TRUE, see below); PLNnetwork and ZIPLNnetwork now default to "builtin", which finds a better optimum at a modest speed cost; ZIPLN keeps its "builtin" default. All backends remain configurable via the backend argument; see the corresponding *_param() documentation for the trade-offs. The torch backend is now clearly marked experimental everywhere.

• Quality and speed improvements: config_optim$profiled = TRUE is now the default for full-covariance nlopt fits (faster, slightly better loglik); ZIPLN's variational step now optimises (M, ψ, R) jointly via Newton instead of sequentially; PLNPCA shares a single SVD initialisation across ranks and can warm-start from a pre-fitted PLNfit for large ranks (inception/init_method, see ?PLNPCA_param); ZIPLN's starting point no longer relies on pscl::zeroinfl (now an internal LM + binomial GLM routine), which is both much faster and a better starting point — pscl is no longer a dependency.

• Fixed a critical nlopt convergence bug affecting PLN/PLNPCA: ill-conditioned covariate scaling could trigger the XTOL stopping criterion after very few iterations, well before convergence. The built-in backend was never affected; nlopt is now also fixed via better parameter scaling.

Internal refactoring (C++)

• Shared covariance abstraction: the optimization machinery for PLN's covariance structures (full, diagonal, spherical, fixed) is now expressed once via a small set of C++ traits (CovTraitsBase in covariance_pln.h) instead of being duplicated per structure. PLNPCA and ZIPLN's variational step now reuse the same machinery instead of separate hand-rolled implementations, removing a substantial amount of duplicated code and fixing minor inefficiencies along the way (e.g. ZIPLN's VE-step used to treat the precision matrix as dense even for diagonal/spherical covariance, at unnecessary O(np^2) cost).
• Consistent C++ naming: exported optimizer functions across PLN, PLNPCA and ZIPLN now follow the same {backend}_optimize_{structure} convention.

Other changes

• Uniform covariate normalization: a normalize_covariates() helper (zero mean, unit variance per column) is now applied consistently in all optimize() methods (PLN, PLNPCA, PLNnetwork, ZIPLN). This makes the nlopt XTOL criterion scale-invariant and stabilises the torch backend.

• Parallelism backend: future.apply::future_lapply is replaced by parallel::mclapply throughout (stability selection for PLNnetwork / ZIPLNnetwork). Use options(mc.cores = N) to set the number of cores.

• Bug fixes: PLNnetwork/ZIPLNnetwork's inception (warm-start) model didn't inherit ftol_em/maxit_em from the user's config_optim, silently falling back to defaults and producing a wrong penalty grid; the PLNPCA rank-model objective used A − Y where it should use A − Y ⊙ Z; various ZIPLN prediction/initialization fixes (#146, #149, #150, #152).

• microcosm data now included (#153, #154); AIC added for PLN and ZIPLN classes (#151); other fixes (#155).

• CRAN check fixes (no user-visible effect): removed an unused -fopenmp compilation flag in src/Makevars that was inadvertently turning on Armadillo's internal OpenMP parallelisation and inflating the CPU/elapsed time ratio of several examples; fixed a GCC -Wmismatched-new-delete false positive in src/packing.cpp's internal test helper by restructuring the code (no diagnostic-suppressing pragma involved).

PLNmodels 1.2.2 (2025-03-21)

• fix for #143 (remove LBFGS_NOCEDAL variant from the possible algorithms)

PLNmodels 1.2.1 (2025-03-10)

• fix NOTES in CRAN due to missing packages in \link{} (PR #142)
• Now requires R >= 4.1.0 because package code uses the pipe |> (PR #142)
• fix sandwich variance estimation (PR #140)
• fix use of native pipe to ensure compatibility with R 3.6 (merge PR #125, fix #124)

PLNmodels 1.2.0 (2024-03-05)

• new feature: ZIPLN (PLN with zero inflation) for standard PLN and PLN Network
  • ZIPLN() and ZIPLNfit-class to allow for zero-inflation in the standard PLN model (merge PR #116)
  • ZIPLNnetwork() and ZIPLNfit_sparse-class to allow for zero-inflation in the PLNnetwork model (merge PR #118)
  • Code factorization between PLNnetwork and ZIPLNnetwork (and associated classes)
• fix inconsistency between fitted and predict (merge PR #115)

PLNmodels 1.1.0 (2024-01-08)

• Update documentation of PLN*_param() functions to include torch optimization parameters
• Add (somehow) explicit error message when torch convergence fails
• Change initialization in variance_jackknife() and variance_bootstrap() to prevent estimation recycling, results from those functions are now comparable to doing jackknife / bootstrap "by hand".
• Merge PR #110 from Cole Trapnell to add:
  • bootstrap estimation of the variance of model parameter
  • improved interface for model initialization / optimisation parameters, which are now passed on to jackknife / bootstrap post-treatments
  • better support of GPU when using torch backend
• Change behavior of predict() function for PLNfit model to (i) return fitted values if newdata is missing or (ii) perform one VE step to improve fit if responses are provided (fix issue #114)

PLNmodels 1.0.4 (2023-08-24)

• changed initial value in optim for variational variance (1 -> 0.1) in VE-step of PLN and PLNPCA
• fix sign in objective of VE_step for PLN with full covariance Issue #100
• add a scale argument compute_offset() to force the offsets (RLE, CSS, GMPR, Wrench) to be on the same scale as the counts, like TSS.
• add a new "TMM" for compute_offset()
• fix nb_param for PLNLDA, which caused wrong BIC/ICL and erratic model selection
• fix minor issues #102, #103 plus some others
• fix package file documentation as suggested in r-lib/roxygen2#1491

PLNmodels 1.0.3 (2023-07-06)

• higher tolerance on a single test (among 700) that fails on the 'noLD' additional architecture on CRAN (tests without long double)

PLNmodels 1.0.2 (2023-06-21)

• changed initial value in optim for variational variance (1 -> 0.1), which caused failure in some cases
• fix bug when using inception in PLNnetwork()
• starting handling of missing data
• slightly faster (factorized) initialization for PCA

PLNmodels 1.0.1 (2023-02-12)

• fix in the use of future_lapply which used to make post-Treatments in PLNPCA last for ever with multicore in v1.0.0...
• prevent use of bootstrap/jackknife when not appropriate
• fix bug in PLNmixture() when the sequence of cluster numbers (clusters) is not of the form 1:K_max
• use bibentry to replace citEntry in CITATION

PLNmodels 1.0.0

Breaking changes

• interface for controlling the fits now use list generated by dedicated functions

  • PLN_param() for PLN
  • PLNLDA_param() for PLNLDA
  • PLNnetwork_param() for PLNnetwork
  • PLNPCA_param() for PLNPCA
  • PLNmixture_param() for PLNmixture The use of 'control = list()' is deprecated: the code stop and send an error.

• The regression coefficients are now denoted by B, not Theta, such as B = t(Theta). We keep on sending back Theta as a field of myPLN$model_par$Theta, but this will soon be deprecated

New features

• added Barents fish data set
• support for PLN when (inverse) covariance is known/fixed
• estimator of the variance of the model parameters
  • integration of sandwich estimator of the variance-covariance of Theta when Sigma is fixed
  • variational estimation of the variance-covariance based on variational approximation of the Fisher information
  • jackknife estimation of the variance of Theta and Sigma
  • bootstrap estimation of the variance of Theta and Sigma
• handle list of penalty weights in PLNnetwork
• first support for torch optimizers (for PLN and PLNLDA)

Bug fixes

• fix in objective functions of ve_step of standard PLN models
• fix in objective functions of main of standard PLN models

PLNmodels 0.11.7

• fix expression of ELBO in VEstep, related to #91
• typos and regeneration of documentation( HTML5)
• added an S3 method predict_cond to perform conditional predictions
• fix #89 bug by forcing an intercept in PLNLDA() and changing extract_model() to conform with model.frame()

PLNmodels 0.11.6

• fix wrong use of all.equal
• fix linking problem in new version of nloptr (>=2.0.0)

PLNmodels 0.11.5

• fixing #79 by using the same variational distribution to approximate the spherical case as in the fully parametrized and diagonal cases
• faster examples and build for vignettes
• additional R6 method $VEStep() for PLN-PCA, dealing with low rank matrices
• additional R6 method $project() for PLN-PCA, used to project newdata into PCA space
• use future_lapply in PLNmixture_family
• remove a NOTE due to a DESeq2 link and a failure on solaris on CRAN machines
• some bug fixes

PLNmodels 0.11.4

• use future_lapply in PLNPCA, PLNmixture and stability_selection (plan must be set by the user)
• bug fix in prediction for PLN-LDA
• bug fix in gradients of PLN-network and PLN-spherical
• suppressing method $latent_pos() which is equivalent to active binding $latent
• finalizing integration of PLNmixture (in particular faster smoothing)
• added an argument 'reverse' to the plot methods for criteria, so that users can get their "usual" BIC definition (-2 loglik)

PLNmodels 0.11.3

• support for covariates in PLNmixture (spherical, diagonal, full)
• more support for PLNmixture (S3/R6 methods, vignette)

PLNmodels 0.11.2

• Rewriting C++ by merging modern_cpp to dev, thanks to François Gindraud
• various bug fixes in offset
• less verbose about R squared when questionable
• correction in BIC/ICL for PLNPCA
• Enhanced vignettes for PLNPCA and PLNmixture

PLNmodels 0.11.1

• Add compatibility with factoextra for PLNPCA

PLNmodels 0.11.0

• Add development version of PLNmixture

PLNmodels 0.10.7

• add type = "poscounts" option to RLE normalization
• added wrench normalization to the list of available offsets
• added the oaks data set from Jakuschkin et al (2016)

PLNmodels 0.10.6

• Correction in likelihood of diagonal PLN
• amending test-pln to fulfill CRAN request (error on ATLAS variant of BLAS...)

PLNmodels 0.10.5

• Refactor code of R6 classes to benefit from Roxygen 7.0.0 R6-related new features for documentation

PLNmodels 0.10.4

• Change name of variational variance parameters to S2 (used to be S)
• use spell_check to check spelling, found many typos

PLNmodels 0.10.3

• Change in optimization for all PLN models (PLNs, PCA, LDA, networks): solving in S such that S = S² for the variational parameters, thus avoiding lower bound and constrained optimization. Slightly finer results/estimations for similar computational cost, but easier to maintain.

PLNmodels 0.10.2

• Fix bug in predict() methods when factor levels differ between train and test datasets.
• Fix bug in PLNPCAfit S3 plot() method
• Some simplification in C++ code
• correction/changes in PLN likelihoods? + added constant terms in all likelihoods of all PLN models
• VEstep now available for all model of covariance in PLN (full, diagonal, spherical)

PLNmodels 0.9.5 - minor release

• removed any use of rmarkdown::paged_table() in the vignettes
• added screenshot.force = FALSE, in knitr options in the vignettes

PLNmodels 0.9.4 - minor release

• removing dependencies to bioconductor packages, too cumbersome to maintain on CRAN

PLNmodels 0.9.3 - minor release

• correction in test to comply new class of matrix object

PLNmodels 0.9.2.9002 - development version

• added the possibility for matrix of weights for the penalty in PLNnetworks

PLNmodels 0.9.2

• various bug fixes

PLNmodels 0.9.1

• Use nloptr to prepare CRAN release

PLNmodels 0.8.2

• Enhancement in PLNLDA

PLNmodels 0.8.1

• Preparing first CRAN release

PLNmodels 0.7.0.1

• Added a NEWS.md file to track changes to the package.