The defaults of PLNnetwork(), ZIPLNnetwork() and ZIPLN() with a sparse covariance
have changed since version 1.3.2. Penalty paths, criteria and selected models differ.
penalty_scale = "correlation"), where it used to apply on the covariance scale. The penalties
change meaning: they are dimensionless and lie between 0 and 1, a penalty of 1 or
more giving the empty network. This concerns the grid built by default as well as
the penalties given through penalties = (or penalty = in ZIPLN_param()).latent_floor = 1e-3).penalize_diagonal = FALSE), and PLNnetwork() starts from a diagonal inception (inception_cov = "diagonal").To get the former behavior back, set the scale in the control parameters:
PLNnetwork(..., control = PLNnetwork_param(penalty_scale = "covariance"))
and, to undo the other changes as well,
PLNnetwork_param(penalty_scale = "covariance", latent_floor = NULL,
penalize_diagonal = TRUE, inception_cov = "full")
(likewise in ZIPLNnetwork_param() and ZIPLN_param(), without inception_cov, whose
default has not changed there). The fits are then close to those of 1.3.2, not
identical: the graphical Lasso and the grid of penalties have changed too (see below).
Penalties given explicitly. When penalties is given while penalty_scale is left
to its default:
ZIPLN_param(), a penalty above 1 is not converted: a warning says that the
network will be empty.Setting penalty_scale explicitly, to either value, leaves the penalties as given and
silences these messages.
Why. The entries of a precision matrix are not scale invariant, so that on the
covariance scale a species with a large latent variance has nearly free edges. In
simulations with a known network (1 960 datasets, a hundred configurations of graph,
sample size, dimension, latent variances, abundances, contamination, covariates and
zero inflation), the edges were recovered better on the correlation scale in every
configuration: the F1 score at the true network size went from 0.58 to 0.74 on
uncontaminated data (from 0.45 to 0.73 when the latent variances differ between
species), and from 0.27 to 0.78 with species absent from part of the samples, for the
same computing time. The gain carries over to model selection, less strongly (BIC:
0.47 to 0.56; StARS: 0.59 to 0.71). See inst/simus_PLNnetwork/penalty_scale/.
penalty_scale in PLNnetwork_param(), ZIPLNnetwork_param() and
ZIPLN_param(), "correlation" (default) or "covariance". The l1 penalty of the graphical Lasso bears on the entries of the
precision matrix, which are not scale invariant: a species with a large latent
variance has nearly free edges. A species that is often absent but abundant when
present, whose zeros are fitted by very negative latent means, therefore ends up
connected to most of the others (see $degenerate_species below), and this happens
well before its latent variance blows up. On the correlation scale, the penalty on
the pair (i, j) is lambda * sqrt(S_ii * S_jj), recomputed at each M step from the
residual covariance S: this is the graphical Lasso on the residual correlation
matrix, as is customary for Gaussian graphical models, and the penalties become
dimensionless, between 0 and 1.latent_floor in PLNnetwork_param(): a floor on the variational means of
the degenerate species. As soon as the latent variance of a species exceeds the
threshold of $degenerate_species (100) during the optimization, its variational
means are kept above log(floor) - O, that is exp(O + M) >= floor, in this fit and
in the following ones along the penalty path. The species are found by the
optimization itself, so that no preparation of the data is needed, and are returned
by the new field $floored_species. The floor restricts the variational family, not
the model, and has no effect on a fit without degenerate species. It stops the
latent variances from diverging, which the correlation scale alone does not prevent
when species are absent from whole groups of samples: in simulations with such
groups the variances diverged in 20 replicates out of 20 without the floor, on both
scales, and in none with it. It is on by default (1e-3); latent_floor = NULL
removes it.
A floor on every cell was tried first and dropped: on mollusk, where nothing
diverges, it constrained 39 % of the cells.ZIPLNnetwork() and to ZIPLN() with a sparse
covariance (latent_floor in ZIPLNnetwork_param() and ZIPLN_param(), 1e-3 by
default), which degenerate as PLNnetwork() does: zero inflation does not protect
from it. On oaks, the largest latent variance along the default path goes from
193 000 to 93 and the largest degree at about p edges from 111 to 21. As for
PLNnetwork(), a fit without degenerate species is left exactly as it is.oaks), and the edges still
concentrate on the degenerate species.penalty_scale = "correlation" and latent_floor = 1e-3 (see the
breaking changes above). The floor changes the fits where species degenerate, and
only there: a fit without degenerate species, as on trichoptera, is exactly the one
obtained without the floor.PLNnetwork() now also names the species
bounded by the floor, since the floor is what keeps their variance under the
threshold, and says, on the covariance scale, that their edges are artefacts.inst/simus_PLNnetwork/degenerate_species/, the account in
inst/devlog/DEVLOG_2026-09-30_10-01.md.$pen_loglik is now the criterion that the M step maximizes, loglik - n/2 * sum(abs(rho * Omega)), where rho is the penalty matrix of the graphical Lasso: the
penalty times the penalty weights and, on the correlation scale, times
sqrt(S_ii * S_jj). It was loglik - penalty * sum(abs(Omega)), which ignored the
weights, the factor n/2 and the scale, and was not the criterion being optimized.
Its values change, on both scales, in $criteria and in the plots of the criteria.?PLNnetwork_param and in the PLNnetwork vignette.$penalty_scale, $latent_floor and $floored_species,
the number of cells at the floor is in $optim_par$n_floor, and
stability_selection() refits the subsamples with the settings of the collection.graphical_lasso() now solves the problem scaled to a unit diagonal. For any
positive diagonal D, Theta solves the problem for (S, rho) if and only if
D^-1 Theta D^-1 solves it for (DSD, D rho D): the change of variables is exact,
and only the stopping rule sees it. On ordinary input the result is the one of
glassoFast up to the tolerance (same supports along penalty paths, relative
differences below 1e-4). On a covariance whose variances span orders of magnitude,
as the residual covariance of PLNnetwork() does, the unscaled descent stalled in
a limit cycle, failed in its inner loop, or returned an indefinite precision
matrix; scaled, it converges in a few sweeps. On 91 such problems met along
PLNnetwork() paths, all converged once scaled, to a lower objective; on ordinary
ones, the median number of sweeps went from 646 to 4. The detection of limit
cycles ("stalled") is kept as a safeguard.graphical_lasso() always returns a positive definite wi. It is backed out of
the regression coefficients of the algorithm, so that it is the inverse of w only
at the exact solution, and could come out indefinite short of it, whatever the
status. Its diagonal is then shifted, which keeps the network, so that its
smallest eigenvalue is the one of the inverse of w; the shift is reported in the
new shift component.PLNnetwork() no longer stops the
whole path. A non-finite precision matrix from the graphical Lasso, or a
non-finite objective, ends the optimization for that penalty on the previous
iterate, with a warning; if there is none, the fit is marked as failed (NA
criteria), and getBestModel() leaves it out of the selection, with a warning,
instead of stopping. The reason is in $optim_par$failure, and the number of
shifted precision matrices in $optim_par$glasso_indefinite (also for
ZIPLNnetwork()). ZIPLN() likewise stops on a non-finite objective, on the best
iterate so far, instead of failing with "missing value where TRUE/FALSE needed".PLNnetwork(Abundance ~ 1)),
version 1.3.2 failed with an error on 3 of them, returned NA criteria on 3 more,
and non-finite or absurd log-likelihoods (up to 1e109) on 4 others. All 30 paths now
fit, the BIC of the selected model is higher (median +7 on the 20 datasets where
1.3.2 did fit), and the whole study runs twice as fast.penalize_diagonal = FALSE in PLNnetwork_param(), ZIPLNnetwork_param() and
ZIPLN_param()). Penalizing it inflates the latent variances by the penalty
(Sigma_ii = S_ii + rho), which the VE step feeds back into the residual
covariance S: along the path, the network could then never become empty,
whatever the penalty. On oaks, the top of the path had 149 edges, for a residual
covariance 73 times the largest penalty of the grid. Set penalize_diagonal = TRUE
to get the former behavior back.crossprod(M - XB) / n, whatever its covariance model, rather than
from its fitted Sigma, and PLNnetwork() now starts from a diagonal inception
(inception_cov = "diagonal", fully converged). With an unpenalized diagonal, the
inception is then the empty network, and the top of the grid the penalty above
which it is a fixed point of the alternating optimization: the path starts from
the empty network (up to a pair on the boundary). ZIPLNnetwork() keeps a full
inception, which led to a better BIC along the path.PLNnetwork() goes from -1390 to -1300 on trichoptera, from
-5331 to -5149 on barents, from -38664 to -37566 on oaks and from -37399 to
-36630 on oaks with the tree covariate, as fast or faster. For ZIPLNnetwork(),
from -1418 to -1306, -5448 to -5206, -39786 to -37585 and -38348 to -36716, in up to
2.5 times as long on oaks, where the paths are denser.oaks without covariate, the networks along the path of
PLNnetwork() are made of little else: f_OTU_1011, present on all the trees of
one type and on no other, carries 103 of the 103 edges of the first non-empty
networks. This is not new, and does not depend on the penalization of the diagonal.
PLN(), PLNnetwork(), ZIPLN() and ZIPLNnetwork() now raise a warning (one
per call, for a whole collection) naming the species whose latent variance is above
100, which the new field $degenerate_species of a fit returns. The threshold is
set by options(PLNmodels.latent_variance_threshold = ).structural_zeros(counts, covariates) finds the species absent from all the
samples of a level of a factor covariate while present elsewhere, beyond what
chance would explain (exact hypergeometric test, Bonferroni-corrected, so that rare
species absent from a level by chance are not reported). prepare_data() runs it on
all the factor, character and logical covariates and reports the result in a
message; it does not remove anything. On oaks it finds the 8 species that live on
one or two of the three types of trees.M with B profiled (B = P_X M), but the objective was then evaluated at the B of the preceding M step, which does not match the new M: the reported objective could jump by orders of magnitude, and the optimization settle far from the optimum. B is now updated to match M after the VE step, and the M step computes Omega with the new B (the joint optimum) rather than the previous one. On data with many zeros that are not excess zeros, fits were far worse than PLN()'s, although ZIPLN nests PLN; they are now at least as good. On benign data (trichoptera, oaks) the log-likelihood is unchanged up to a few units."converged". The test is now on the absolute change, the best iterate is the one returned, and the number of increases is reported in $optim_par$objective_increases (it should be 0). The same signed test is fixed in the VE step used by predict().ZIPLNnetwork(), or ZIPLN() with a
penalty) now includes the l1 penalty of the graphical Lasso, since the M step
minimizes the penalized objective: the unpenalized one can legitimately increase.R of a structural zero is now 0 where the count is positive, where the posterior is exactly 0, instead of the zero rate of the species on every cell (#186).ZIPLN() and ZIPLNnetwork() accept a formula stored in a variable (#187), as
PLN() does: f <- Abundance ~ 1 + Wind; ZIPLN(f, data) failed with "cannot set an attribute on a 'symbol'", and predict() failed on such fits with an intercept-only formula.compute_offset() and its methods):
offset_rle() keeps only species with a positive geometric mean in its
median-ratio estimator (as DESeq2 does), so a zero-inflated species no longer
drives the per-sample offset to Inf via a counts / 0 ratio even when several zero-free species are available.offset_css() falls back to TSS only for the samples with fewer than two
positive counts (instead of silently switching every sample to rowSums), and scales the whole result by a single global median — so all offsets are on the same scale (median 1). The warning names the affected samples.compute_offset() detects empty (all-zero) samples before dispatching, warns
naming them, computes the offset on the non-empty samples and returns NA for the empty ones, instead of a 0 offset (hence log(0) = -Inf) or uninformative NA/NaN/Inf errors when several methods are called directly.config_post = list(jackknife = TRUE, bootstrap = ...)) failed with a fixed or a "genpop" covariance, the fixed
precision matrix Omega or correlation matrix C not being passed to the fits on
the resampled data. They now work, also in PLNnetwork().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).setTimeLimit()/R.utils::withTimeout(), which then raise their usual error.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.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.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".$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.1 / (S_ii + rho_ii), where glassoFast returned 1 / max(rho_ii, 1.1e-16).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").$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.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).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).
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.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.
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).{backend}_optimize_{structure} convention.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).
variance_jackknife() and variance_bootstrap() to prevent estimation recycling, results from those functions are now comparable to doing jackknife / bootstrap "by hand".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)scale argument compute_offset() to force the offsets (RLE, CSS, GMPR, Wrench) to be on the same scale as the counts, like TSS.clusters) is not of the form 1:K_maxinterface for controlling the fits now use list generated by dedicated functions
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
PLNLDA() and changing extract_model() to conform with model.frame()$VEStep() for PLN-PCA, dealing with low rank matrices$project() for PLN-PCA, used to project newdata into PCA space$latent_pos() which is equivalent to active binding $latentNEWS.md file to track changes to the package.