Package {PLNmodels}

Title:	Poisson Lognormal Models
Version:	1.3.0
Description:	The Poisson-lognormal model and variants (Chiquet, Mariadassou and Robin, 2021 <doi:10.3389/fevo.2021.588292>) can be used for a variety of multivariate problems when count data are at play, including principal component analysis for count data, discriminant analysis, model-based clustering and network inference. Implements variational algorithms to fit such models accompanied with a set of functions for visualization and diagnostic.
URL:	https://pln-team.github.io/PLNmodels/
BugReports:	https://github.com/pln-team/PLNmodels/issues
License:	GPL (≥ 3)
Depends:	R (≥ 4.1.0)
Imports:	cli, corrplot, dplyr, ggplot2, glassoFast, grDevices, grid, gridExtra, igraph, magrittr, MASS, Matrix, methods, nloptr, parallel, purrr, R6, Rcpp, rlang, stats, tidyr, torch
Suggests:	factoextra, knitr, pheatmap, rmarkdown, spelling, quarto, testthat
LinkingTo:	nloptr, Rcpp, RcppArmadillo
VignetteBuilder:	knitr
Encoding:	UTF-8
Language:	en-US
LazyData:	true
Collate:	'PLNfit-class.R' 'PLN.R' 'PLNLDA.R' 'PLNLDAfit-S3methods.R' 'PLNLDAfit-class.R' 'PLNPCA.R' 'PLNPCAfamily-S3methods.R' 'PLNfamily-class.R' 'PLNPCAfamily-class.R' 'PLNPCAfit-S3methods.R' 'PLNPCAfit-class.R' 'PLNfamily-S3methods.R' 'PLNfit-S3methods.R' 'PLNmixture.R' 'PLNmixturefamily-S3methods.R' 'PLNmixturefamily-class.R' 'PLNmixturefit-S3methods.R' 'PLNmixturefit-class.R' 'PLNmodels-package.R' 'PLNnetwork.R' 'PLNnetworkfamily-S3methods.R' 'PLNnetworkfamily-class.R' 'PLNnetworkfit-S3methods.R' 'PLNnetworkfit-class.R' 'RcppExports.R' 'ZIPLNfit-class.R' 'ZIPLN.R' 'ZIPLNfit-S3methods.R' 'ZIPLNnetwork.R' 'barents.R' 'import_utils.R' 'microcosm.R' 'mollusk.R' 'oaks.R' 'plot_utils.R' 'scRNA.R' 'trichoptera.R' 'utils-pipe.R' 'utils-zipln.R' 'utils.R' 'zzz.R'
Config/roxygen2/version:	8.0.0
NeedsCompilation:	yes
Packaged:	2026-07-27 13:13:29 UTC; jchiquet
Author:	Julien Chiquet [aut, cre], Mahendra Mariadassou [aut], Stéphane Robin [aut], François Gindraud [aut], Julie Aubert [ctb], Bastien Batardière [ctb], Giovanni Poggiato [ctb], Cole Trapnell [ctb], Maddy Duran [ctb]
Maintainer:	Julien Chiquet <julien.chiquet@inrae.fr>
Repository:	CRAN
Date/Publication:	2026-07-27 14:10:26 UTC

PLNmodels: Poisson Lognormal Models

Description

The Poisson-lognormal model and variants (Chiquet, Mariadassou and Robin, 2021 doi:10.3389/fevo.2021.588292) can be used for a variety of multivariate problems when count data are at play, including principal component analysis for count data, discriminant analysis, model-based clustering and network inference. Implements variational algorithms to fit such models accompanied with a set of functions for visualization and diagnostic.

Author(s)

Maintainer: Julien Chiquet julien.chiquet@inrae.fr (ORCID)

Authors:

• Julien Chiquet julien.chiquet@inrae.fr (ORCID)

• Mahendra Mariadassou mahendra.mariadassou@inrae.fr (ORCID)

• Stéphane Robin stephane.robin@inrae.fr

• François Gindraud francois.gindraud@gmail.com

Other contributors:

• Julie Aubert julie.aubert@inrae.fr [contributor]

• Bastien Batardière bastien.batardiere@inrae.fr [contributor]

• Giovanni Poggiato giov.poggiato@gmail.com [contributor]

• Cole Trapnell coletrap@uw.edu [contributor]

• Maddy Duran duran@uw.edu [contributor]

See Also

Useful links:

• https://pln-team.github.io/PLNmodels/

• Report bugs at https://github.com/pln-team/PLNmodels/issues

Pipe operator

Description

See magrittr::%>% for details.

Usage

lhs %>% rhs

Akaike Information Criterion for a fitted PLN model

Description

Computes the variational AIC as loglik - nb_param (larger is better). This follows the maximization convention used throughout PLNmodels.

Usage

## S3 method for class 'PLNfit'
AIC(object, ..., k = 2)

Arguments

Value

A scalar: the variational AIC (larger is better).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- PLN(Abundance ~ 1, data = trichoptera)
AIC(model)

Akaike Information Criterion for a fitted ZIPLN model

Description

Computes the variational AIC as loglik - nb_param (larger is better). This follows the maximization convention used throughout PLNmodels.

Usage

## S3 method for class 'ZIPLNfit'
AIC(object, ..., k = 2)

Arguments

Value

A scalar: the variational AIC (larger is better).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- ZIPLN(Abundance ~ 1, data = trichoptera)
AIC(model)

Bayesian Information Criterion for a fitted PLN model

Description

Computes the variational BIC as loglik - 0.5 * log(n) * nb_param (larger is better). This follows the maximization convention used throughout PLNmodels.

Usage

## S3 method for class 'PLNfit'
BIC(object, ...)

Arguments

Value

A scalar: the variational BIC (larger is better).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- PLN(Abundance ~ 1, data = trichoptera)
BIC(model)

Bayesian Information Criterion for a fitted ZIPLN model

Description

Computes the variational BIC as loglik - 0.5 * log(n) * nb_param (larger is better). This follows the maximization convention used throughout PLNmodels.

Usage

## S3 method for class 'ZIPLNfit'
BIC(object, ...)

Arguments

Value

A scalar: the variational BIC (larger is better).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- ZIPLN(Abundance ~ 1, data = trichoptera)
BIC(model)

Integrated Classification Likelihood

Description

Generic function to compute the Integrated Classification Likelihood (ICL) of a fitted model. ICL = BIC - entropy of the variational distribution (larger is better).

ICL.PLNfit: ICL for a fitted PLNfit.

ICL.ZIPLNfit: ICL for a fitted ZIPLNfit.

Usage

ICL(object, ...)

## S3 method for class 'PLNfit'
ICL(object, ...)

## S3 method for class 'ZIPLNfit'
ICL(object, ...)

Arguments

Value

A scalar: the variational ICL (larger is better).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- PLN(Abundance ~ 1, data = trichoptera)
ICL(model)
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- ZIPLN(Abundance ~ 1, data = trichoptera)
ICL(model)

An R6 Class to virtually represent a collection of network fits

Description

The functions PLNnetwork() and ZIPLNnetwork() both produce an instance of this class, which can be thought of as a vector of PLNnetworkfits ZIPLNfit_sparses (indexed by penalty parameter)

This class comes with a set of methods mostly used to compare network fits (in terms of goodness of fit) or extract one from the family (based on penalty parameter and/or goodness of it). See the documentation for getBestModel(), getModel() and plot() for the user-facing ones.

Super class

PLNfamily -> Networkfamily

Active bindings

Methods

Public methods

• Networkfamily$new()

• Networkfamily$optimize()

• Networkfamily$coefficient_path()

• Networkfamily$getBestModel()

• Networkfamily$plot()

• Networkfamily$plot_stars()

• Networkfamily$plot_objective()

• Networkfamily$show()

• Networkfamily$clone()

Networkfamily$new()

Initialize all models in the collection

Usage

Networkfamily$new(penalties, data, control)

Arguments

Returns

Update all network fits in the family with smart starting values

Networkfamily$optimize()

Call to the C++ optimizer on all models of the collection

Usage

Networkfamily$optimize(data, config)

Arguments

Networkfamily$coefficient_path()

Extract the regularization path of a Networkfamily

Usage

Networkfamily$coefficient_path(precision = TRUE, corr = TRUE)

Arguments

Networkfamily$getBestModel()

Extract the best network in the family according to some criteria

Usage

Networkfamily$getBestModel(crit = c("BIC", "EBIC", "StARS"), stability = 0.9)

Arguments

Details

For BIC and EBIC criteria, higher is better.

Networkfamily$plot()

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (a Networkfamily)

Usage

Networkfamily$plot(
  criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
  reverse = FALSE,
  log.x = TRUE
)

Arguments

Returns

a ggplot2::ggplot graph

Networkfamily$plot_stars()

Plot stability path

Usage

Networkfamily$plot_stars(stability = 0.9, log.x = TRUE)

Arguments

Returns

a ggplot2::ggplot graph

Networkfamily$plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

Networkfamily$plot_objective()

Returns

a ggplot2::ggplot graph

Networkfamily$show()

User friendly print method

Usage

Networkfamily$show()

Networkfamily$clone()

The objects of this class are cloneable with this method.

Usage

Networkfamily$clone(deep = FALSE)

Arguments

See Also

The functions PLNnetwork(), ZIPLNnetwork() and the classes PLNnetworkfit, ZIPLNfit_sparse

Poisson lognormal model

Description

Fit the multivariate Poisson lognormal model with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets, weights).

Usage

PLN(formula, data, subset, weights, control = PLN_param())

Arguments

Value

an R6 object with class PLNfit

See Also

The class PLNfit and the configuration function PLN_param()

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1, data = trichoptera)

Poisson lognormal model towards Linear Discriminant Analysis

Description

Fit the Poisson lognormal for LDA with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets).

Usage

PLNLDA(formula, data, subset, weights, grouping, control = PLNLDA_param())

Arguments

Details

See PLNLDA_param() for a full description of the optimization parameters. Note that unlike PLN_param(), PLNLDA_param() does not expose the "fixed" covariance option or the Omega parameter, which are not meaningful in the LDA context.

Value

an R6 object with class PLNLDAfit()

See Also

The class PLNLDAfit

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ 1, grouping = Group, data = trichoptera)

Control of a PLNLDA fit

Description

Helper to define list of parameters to control the PLNLDA fit. All arguments have defaults.

Usage

PLNLDA_param(
  backend = c("builtin", "nlopt", "torch"),
  trace = 1,
  covariance = c("full", "diagonal", "spherical"),
  config_post = list(),
  config_optim = list(),
  inception = NULL
)

Arguments

Details

The list of parameters config_optim controls the optimizers. When "nlopt" is chosen the following entries are relevant

• "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)

• "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)

• "profiled" (full covariance only) if TRUE, profile both B and Omega at every nlopt evaluation instead of running an EM loop (Omega fixed for the duration of each inner nlopt solve, B profiled in closed form at every evaluation). Despite the extra O(n*p^2 + p^3) cost per evaluation, benchmarks found it consistently faster than the EM loop (and with a slightly better loglik) across a range of problem sizes. Default is TRUE; set to FALSE to recover the EM loop.

When "torch" backend is used (only for PLN and PLNLDA for now), the following entries are relevant:

• "algorithm" the optimizer used by torch among RPROP (default), RMSPROP, ADAM and ADAGRAD

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10 000

• "numepoch" stop training once this number of epochs exceeds numepoch. Set to -1 to enable infinite training. Default is 1 000

• "num_batch" number of batches to use during training. Defaults to 1 (use full dataset at each epoch)

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "lr" learning rate. Default is 0.1.

• "momentum" momentum factor. Default is 0 (no momentum). Only used in RMSPROP

• "weight_decay" Weight decay penalty. Default is 0 (no decay). Not used in RPROP

• "step_sizes" pair of minimal (default: 1e-6) and maximal (default: 50) allowed step sizes. Only used in RPROP

• "etas" pair of multiplicative increase and decrease factors. Default is (0.5, 1.2). Only used in RPROP

• "centered" if TRUE, compute the centered RMSProp where the gradient is normalized by an estimation of its variance weight_decay (L2 penalty). Default to FALSE. Only used in RMSPROP

When "builtin" backend is used, the following entries are relevant

• "maxeval" stop when the number of Newton steps in the inner loop exceeds maxeval. Default is 10000

• "ftol_in" stop the inner loop when the objective changes by less than ftol_in (relative). Default is 1e-8

• "maxit_em" stop the EM outer loop when the number of EM iterations exceeds maxit_em. Default is 50

• "ftol_em" stop the EM outer loop when the ELBO changes by less than ftol_em (relative). Default is 1e-8

The list of parameters config_post controls the post-treatment processing (for most ⁠PLN*()⁠ functions), with the following entries (defaults may vary depending on the specific function, check ⁠config_post_default_*⁠ for defaults values):

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimation should be used to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

Value

list of parameters configuring the fit.

An R6 Class to represent a PLNfit in a LDA framework

Description

The function PLNLDA() produces an instance of an object with class PLNLDAfit.

This class comes with a set of methods, some of them being useful for the user: See the documentation for the methods inherited by PLNfit(), the plot() method for LDA visualization and predict() method for prediction

Super class

PLNfit -> PLNLDAfit

Active bindings

Methods

Public methods

• PLNLDAfit$new()

• PLNLDAfit$optimize()

• PLNLDAfit$postTreatment()

• PLNLDAfit$setVisualization()

• PLNLDAfit$plot_individual_map()

• PLNLDAfit$plot_correlation_map()

• PLNLDAfit$plot_LDA()

• PLNLDAfit$predict()

• PLNLDAfit$show()

• PLNLDAfit$clone()

PLNLDAfit$new()

Initialize a PLNLDAfit object

Usage

PLNLDAfit$new(
  grouping,
  responses,
  covariates,
  offsets,
  weights,
  formula,
  control
)

Arguments

PLNLDAfit$optimize()

Compute group means and axis of the LDA (noted B in the model) in the latent space, update corresponding fields

Usage

PLNLDAfit$optimize(grouping, responses, covariates, offsets, weights, config)

Arguments

PLNLDAfit$postTreatment()

Update R2, fisher and std_err fields and visualization

Usage

PLNLDAfit$postTreatment(
  grouping,
  responses,
  covariates,
  offsets,
  weights = rep(1, nrow(responses)),
  config_post,
  config_optim
)

Arguments

PLNLDAfit$setVisualization()

Compute LDA scores in the latent space and update corresponding fields.

Usage

PLNLDAfit$setVisualization(scale.unit = FALSE)

Arguments

PLNLDAfit$plot_individual_map()

Plot the factorial map of the LDA

Usage

PLNLDAfit$plot_individual_map(
  axes = 1:min(2, self$rank),
  main = "Individual Factor Map",
  plot = TRUE
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNLDAfit$plot_correlation_map()

Plot the correlation circle of a specified axis for a PLNLDAfit object

Usage

PLNLDAfit$plot_correlation_map(
  axes = 1:min(2, self$rank),
  main = "Variable Factor Map",
  cols = "default",
  plot = TRUE
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNLDAfit$plot_LDA()

Plot a summary of the PLNLDAfit object

Usage

PLNLDAfit$plot_LDA(
  nb_axes = min(3, self$rank),
  var_cols = "default",
  plot = TRUE
)

Arguments

Returns

a grob object

PLNLDAfit$predict()

Predict group of new samples

Usage

PLNLDAfit$predict(
  newdata,
  type = c("posterior", "response", "scores"),
  scale = c("log", "prob"),
  prior = NULL,
  control = PLN_param(backend = "nlopt"),
  envir = parent.frame()
)

Arguments

PLNLDAfit$show()

User friendly print method

Usage

PLNLDAfit$show()

PLNLDAfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit$clone(deep = FALSE)

Arguments

See Also

The function PLNLDA.

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ 1, grouping = Group, data = trichoptera)
class(myPLNLDA)
print(myPLNLDA)

## End(Not run)

An R6 Class to represent a PLNfit in a LDA framework with diagonal covariance

Description

The function PLNLDA() produces an instance of an object with class PLNLDAfit.

This class comes with a set of methods, some of them being useful for the user: See the documentation for the methods inherited by PLNfit(), the plot() method for LDA visualization and predict() method for prediction

Super classes

PLNfit -> PLNLDAfit -> PLNLDAfit_diagonal

Active bindings

Methods

Public methods

• PLNLDAfit_diagonal$new()

• PLNLDAfit_diagonal$clone()

PLNLDAfit_diagonal$new()

Initialize a PLNfit model

Usage

PLNLDAfit_diagonal$new(
  grouping,
  responses,
  covariates,
  offsets,
  weights,
  formula,
  control
)

Arguments

PLNLDAfit_diagonal$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit_diagonal$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ 1, data = trichoptera, control = PLN_param(covariance = "diagonal"))
class(myPLNLDA)
print(myPLNLDA)

## End(Not run)

Poisson lognormal model towards Principal Component Analysis

Description

Fit the PCA variants of the Poisson lognormal with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets).

Usage

PLNPCA(formula, data, subset, weights, ranks = 1:5, control = PLNPCA_param())

Arguments

Value

an R6 object with class PLNPCAfamily, which contains a collection of models with class PLNPCAfit

See Also

The classes PLNPCAfamily and PLNPCAfit, and the configuration function PLNPCA_param().

Examples

## Use parallel to dispatch the computations on 2 workers
## Not run: 
options(mc.cores = 2)

## End(Not run)

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPCA <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:5)

## Not run: 
options(mc.cores = 1)

## End(Not run)

Control of PLNPCA fit

Description

Helper to define list of parameters to control the PLNPCA fit. All arguments have defaults.

Usage

PLNPCA_param(
  backend = c("nlopt", "builtin", "torch"),
  trace = 1,
  config_optim = list(),
  config_post = list(),
  inception = NULL,
  init_method = c("LM", "GLM", "EM"),
  init_niter = 5L,
  sequential = FALSE
)

Arguments

Details

The list of parameters config_optim controls the optimizers. When "nlopt" is chosen the following entries are relevant

• "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)

• "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)

• "profiled" (full covariance only) if TRUE, profile both B and Omega at every nlopt evaluation instead of running an EM loop (Omega fixed for the duration of each inner nlopt solve, B profiled in closed form at every evaluation). Despite the extra O(n*p^2 + p^3) cost per evaluation, benchmarks found it consistently faster than the EM loop (and with a slightly better loglik) across a range of problem sizes. Default is TRUE; set to FALSE to recover the EM loop.

When "torch" backend is used (only for PLN and PLNLDA for now), the following entries are relevant:

• "algorithm" the optimizer used by torch among RPROP (default), RMSPROP, ADAM and ADAGRAD

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10 000

• "numepoch" stop training once this number of epochs exceeds numepoch. Set to -1 to enable infinite training. Default is 1 000

• "num_batch" number of batches to use during training. Defaults to 1 (use full dataset at each epoch)

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "lr" learning rate. Default is 0.1.

• "momentum" momentum factor. Default is 0 (no momentum). Only used in RMSPROP

• "weight_decay" Weight decay penalty. Default is 0 (no decay). Not used in RPROP

• "step_sizes" pair of minimal (default: 1e-6) and maximal (default: 50) allowed step sizes. Only used in RPROP

• "etas" pair of multiplicative increase and decrease factors. Default is (0.5, 1.2). Only used in RPROP

• "centered" if TRUE, compute the centered RMSProp where the gradient is normalized by an estimation of its variance weight_decay (L2 penalty). Default to FALSE. Only used in RMSPROP

When "builtin" backend is used, the following entries are relevant

• "maxeval" stop when the number of Newton steps in the inner loop exceeds maxeval. Default is 10000

• "ftol_in" stop the inner loop when the objective changes by less than ftol_in (relative). Default is 1e-8

• "maxit_em" stop the EM outer loop when the number of EM iterations exceeds maxit_em. Default is 50

• "ftol_em" stop the EM outer loop when the ELBO changes by less than ftol_em (relative). Default is 1e-8

The list of parameters config_post controls the post-treatment processing (for most ⁠PLN*()⁠ functions), with the following entries (defaults may vary depending on the specific function, check ⁠config_post_default_*⁠ for defaults values):

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimation should be used to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

Value

list of parameters configuring the fit.

An R6 Class to represent a collection of PLNPCAfit

Description

The function PLNPCA() produces an instance of this class.

This class comes with a set of methods, some of them being useful for the user: See the documentation for getBestModel(), getModel() and plot().

Super class

PLNfamily -> PLNPCAfamily

Active bindings

Methods

Public methods

• PLNPCAfamily$new()

• PLNPCAfamily$optimize()

• PLNPCAfamily$getModel()

• PLNPCAfamily$getBestModel()

• PLNPCAfamily$plot()

• PLNPCAfamily$show()

• PLNPCAfamily$clone()

PLNPCAfamily$new()

Initialize all models in the collection. A single SVD of the residual matrix M - X*B is computed once and shared across all ranks. M and B come from either a user-provided PLNfit inception or a fast LM on log-transformed counts (default, controlled by init_method).

Usage

PLNPCAfamily$new(
  ranks,
  responses,
  covariates,
  offsets,
  weights,
  formula,
  control
)

Arguments

PLNPCAfamily$optimize()

Call to the C++ optimizer on all models of the collection

Usage

PLNPCAfamily$optimize(config)

Arguments

PLNPCAfamily$getModel()

Extract model from collection and add "PCA" class for compatibility with factoextra::fviz()

Usage

PLNPCAfamily$getModel(var, index = NULL)

Arguments

Returns

a PLNPCAfit object

PLNPCAfamily$getBestModel()

Extract best model in the collection

Usage

PLNPCAfamily$getBestModel(crit = c("ICL", "BIC"))

Arguments

Returns

a PLNPCAfit object

PLNPCAfamily$plot()

Lineplot of selected criteria for all models in the collection

Usage

PLNPCAfamily$plot(criteria = c("loglik", "BIC", "ICL"), reverse = FALSE)

Arguments

Returns

A ggplot2::ggplot object

PLNPCAfamily$show()

User friendly print method

Usage

PLNPCAfamily$show()

PLNPCAfamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNPCAfamily$clone(deep = FALSE)

Arguments

See Also

The function PLNPCA(), the class PLNPCAfit()

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPCAs <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:5)
class(myPCAs)

An R6 Class to represent a PLNfit in a PCA framework

Description

The function PLNPCA() produces a collection of models which are instances of object with class PLNPCAfit. This class comes with a set of methods, some of them being useful for the user: See the documentation for the methods inherited by PLNfit and the plot() methods for PCA visualization

Super class

PLNfit -> PLNPCAfit

Active bindings

Methods

Public methods

• PLNPCAfit$new()

• PLNPCAfit$warm_start_from()

• PLNPCAfit$update()

• PLNPCAfit$optimize()

• PLNPCAfit$optimize_vestep()

• PLNPCAfit$project()

• PLNPCAfit$setVisualization()

• PLNPCAfit$postTreatment()

• PLNPCAfit$plot_individual_map()

• PLNPCAfit$plot_correlation_circle()

• PLNPCAfit$plot_PCA()

• PLNPCAfit$show()

• PLNPCAfit$clone()

PLNPCAfit$new()

Initialize a PLNPCAfit object. Uses the shared SVD from control$svdM (computed once in PLNPCAfamily) to set the starting loadings C and scores M. The regression coefficients B are initialised by the parent PLNfit constructor (LM or user-provided inception).

Usage

PLNPCAfit$new(rank, responses, covariates, offsets, weights, formula, control)

Arguments

PLNPCAfit$warm_start_from()

Reinitialize parameters for sequential warm-starting from a lower-rank fit. Fitted loadings C, scores M, variances S, and regression coefficients B from prev_fit are carried over; new columns are padded using the inception SVD (C) or zeros/0.1 (M/S).

Usage

PLNPCAfit$warm_start_from(prev_fit, svdM)

Arguments

PLNPCAfit$update()

Update a PLNPCAfit object

Usage

PLNPCAfit$update(
  B = NA,
  Sigma = NA,
  Omega = NA,
  C = NA,
  M = NA,
  S2 = NA,
  Z = NA,
  A = NA,
  Ji = NA,
  R2 = NA,
  monitoring = NA
)

Arguments

Returns

Update the current PLNPCAfit object

PLNPCAfit$optimize()

Call to the C++ optimizer and update of the relevant fields

Usage

PLNPCAfit$optimize(responses, covariates, offsets, weights, config)

Arguments

PLNPCAfit$optimize_vestep()

Result of one call to the VE step of the optimization procedure: optimal variational parameters (M, S) and corresponding log likelihood values for fixed model parameters (C, B). Intended to position new data in the latent space for further use with PCA.

Usage

PLNPCAfit$optimize_vestep(
  covariates,
  offsets,
  responses,
  weights = rep(1, self$n),
  control = PLNPCA_param(backend = "nlopt")
)

Arguments

Returns

A list with three components:

• the matrix M of variational means,

• the matrix S2 of variational variances

• the vector log.lik of (variational) log-likelihood of each new observation

PLNPCAfit$project()

Project new samples into the PCA space using one VE step

Usage

PLNPCAfit$project(newdata, control = PLNPCA_param(), envir = parent.frame())

Arguments

Returns

• the named matrix of scores for the newdata, expressed in the same coordinate system as self$scores

PLNPCAfit$setVisualization()

Compute PCA scores in the latent space and update corresponding fields.

Usage

PLNPCAfit$setVisualization(scale.unit = FALSE)

Arguments

PLNPCAfit$postTreatment()

Update R2, fisher, std_err fields and set up visualization

Usage

PLNPCAfit$postTreatment(
  responses,
  covariates,
  offsets,
  weights = rep(1, nrow(responses)),
  config_post,
  config_optim,
  nullModel = NULL
)

Arguments

Details

The list of parameters config_post controls the post-treatment processing, with the following entries:

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

• trace integer for verbosity. should be > 1 to see output in post-treatments

PLNPCAfit$plot_individual_map()

Plot the factorial map of the PCA

Usage

PLNPCAfit$plot_individual_map(
  axes = 1:min(2, self$rank),
  main = "Individual Factor Map",
  plot = TRUE,
  cols = "default"
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNPCAfit$plot_correlation_circle()

Plot the correlation circle of a specified axis for a PLNLDAfit object

Usage

PLNPCAfit$plot_correlation_circle(
  axes = 1:min(2, self$rank),
  main = "Variable Factor Map",
  cols = "default",
  plot = TRUE
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNPCAfit$plot_PCA()

Plot a summary of the PLNPCAfit object

Usage

PLNPCAfit$plot_PCA(
  nb_axes = min(3, self$rank),
  ind_cols = "ind_cols",
  var_cols = "var_cols",
  plot = TRUE
)

Arguments

Returns

a grob object

PLNPCAfit$show()

User friendly print method

Usage

PLNPCAfit$show()

PLNPCAfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNPCAfit$clone(deep = FALSE)

Arguments

See Also

The function PLNPCA, the class PLNPCAfamily

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPCAs <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:5)
myPCA <- getBestModel(myPCAs)
class(myPCA)
print(myPCA)

Control of a PLN fit

Description

Helper to define list of parameters to control the PLN fit. All arguments have defaults.

Usage

PLN_param(
  backend = c("nlopt", "builtin", "torch"),
  trace = 1,
  covariance = c("full", "diagonal", "spherical", "fixed", "genpop"),
  Omega = NULL,
  C = NULL,
  config_post = list(),
  config_optim = list(),
  inception = NULL,
  init_method = c("LM", "GLM")
)

Arguments

Details

The list of parameters config_optim controls the optimizers. When "nlopt" is chosen the following entries are relevant

• "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)

• "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)

• "profiled" (full covariance only) if TRUE, profile both B and Omega at every nlopt evaluation instead of running an EM loop (Omega fixed for the duration of each inner nlopt solve, B profiled in closed form at every evaluation). Despite the extra O(n*p^2 + p^3) cost per evaluation, benchmarks found it consistently faster than the EM loop (and with a slightly better loglik) across a range of problem sizes. Default is TRUE; set to FALSE to recover the EM loop.

When "torch" backend is used (only for PLN and PLNLDA for now), the following entries are relevant:

• "algorithm" the optimizer used by torch among RPROP (default), RMSPROP, ADAM and ADAGRAD

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10 000

• "numepoch" stop training once this number of epochs exceeds numepoch. Set to -1 to enable infinite training. Default is 1 000

• "num_batch" number of batches to use during training. Defaults to 1 (use full dataset at each epoch)

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "lr" learning rate. Default is 0.1.

• "momentum" momentum factor. Default is 0 (no momentum). Only used in RMSPROP

• "weight_decay" Weight decay penalty. Default is 0 (no decay). Not used in RPROP

• "step_sizes" pair of minimal (default: 1e-6) and maximal (default: 50) allowed step sizes. Only used in RPROP

• "etas" pair of multiplicative increase and decrease factors. Default is (0.5, 1.2). Only used in RPROP

• "centered" if TRUE, compute the centered RMSProp where the gradient is normalized by an estimation of its variance weight_decay (L2 penalty). Default to FALSE. Only used in RMSPROP

When "builtin" backend is used, the following entries are relevant

• "maxeval" stop when the number of Newton steps in the inner loop exceeds maxeval. Default is 10000

• "ftol_in" stop the inner loop when the objective changes by less than ftol_in (relative). Default is 1e-8

• "maxit_em" stop the EM outer loop when the number of EM iterations exceeds maxit_em. Default is 50

• "ftol_em" stop the EM outer loop when the ELBO changes by less than ftol_em (relative). Default is 1e-8

The list of parameters config_post controls the post-treatment processing (for most ⁠PLN*()⁠ functions), with the following entries (defaults may vary depending on the specific function, check ⁠config_post_default_*⁠ for defaults values):

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimation should be used to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

Value

list of parameters configuring the fit.

An R6 Class to represent a collection of PLNfit

Description

super class for PLNPCAfamily and PLNnetworkfamily.

Public fields

Active bindings

Methods

Public methods

• PLNfamily$new()

• PLNfamily$postTreatment()

• PLNfamily$getModel()

• PLNfamily$plot()

• PLNfamily$show()

• PLNfamily$print()

• PLNfamily$clone()

PLNfamily$new()

Create a new PLNfamily object.

Usage

PLNfamily$new(responses, covariates, offsets, weights, control)

Arguments

Returns

A new PLNfamily object

PLNfamily$postTreatment()

Update fields after optimization

Usage

PLNfamily$postTreatment(config_post, config_optim)

Arguments

PLNfamily$getModel()

Extract a model from a collection of models

Usage

PLNfamily$getModel(var, index = NULL)

Arguments

Returns

A PLNfit object

PLNfamily$plot()

Lineplot of selected criteria for all models in the collection

Usage

PLNfamily$plot(criteria, reverse)

Arguments

Returns

A ggplot2::ggplot object

PLNfamily$show()

User friendly print method

Usage

PLNfamily$show()

PLNfamily$print()

User friendly print method

Usage

PLNfamily$print()

PLNfamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNfamily$clone(deep = FALSE)

Arguments

See Also

getModel()

An R6 Class to represent a PLNfit in a standard, general framework

Description

The function PLN() fit a model which is an instance of a object with class PLNfit. Objects produced by the functions PLNnetwork(), PLNPCA(), PLNmixture() and PLNLDA() also enjoy the methods of PLNfit() by inheritance.

This class comes with a set of R6 methods, some of them being useful for the user and exported as S3 methods. See the documentation for coef(), sigma(), predict(), vcov() and standard_error().

Fields are accessed via active binding and cannot be changed by the user.

Active bindings

Methods

Public methods

• PLNfit$new()

• PLNfit$update()

• PLNfit$optimize()

• PLNfit$optimize_vestep()

• PLNfit$postTreatment()

• PLNfit$predict()

• PLNfit$predict_cond()

• PLNfit$show()

• PLNfit$print()

• PLNfit$clone()

PLNfit$new()

Initialize a PLNfit model

Usage

PLNfit$new(responses, covariates, offsets, weights, formula, control)

Arguments

PLNfit$update()

Update a PLNfit object

Usage

PLNfit$update(
  B = NA,
  Sigma = NA,
  Omega = NA,
  M = NA,
  S2 = NA,
  Ji = NA,
  R2 = NA,
  Z = NA,
  A = NA,
  monitoring = NA
)

Arguments

Returns

Update the current PLNfit object

PLNfit$optimize()

Call to the NLopt or TORCH optimizer and update of the relevant fields

Usage

PLNfit$optimize(responses, covariates, offsets, weights, config)

Arguments

PLNfit$optimize_vestep()

Result of one call to the VE step of the optimization procedure: optimal variational parameters (M, S2) and corresponding log likelihood values for fixed model parameters (Sigma, B). Intended to position new data in the latent space.

Usage

PLNfit$optimize_vestep(
  covariates,
  offsets,
  responses,
  weights,
  B = self$model_par$B,
  Omega = self$model_par$Omega,
  control = PLN_param(backend = "nlopt")
)

Arguments

Returns

A list with three components:

• the matrix M of variational means,

• the matrix S2 of variational variances

• the vector log.lik of (variational) log-likelihood of each new observation

PLNfit$postTreatment()

Update R2, fisher and std_err fields after optimization

Usage

PLNfit$postTreatment(
  responses,
  covariates,
  offsets,
  weights = rep(1, nrow(responses)),
  config_post,
  config_optim,
  nullModel = NULL
)

Arguments

Details

The list of parameters config controls the post-treatment processing, with the following entries:

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimator should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• trace integer for verbosity. should be > 1 to see output in post-treatments

PLNfit$predict()

Predict position, scores or observations of new data.

Usage

PLNfit$predict(
  newdata,
  responses = NULL,
  type = c("link", "response"),
  level = 1,
  envir = parent.frame()
)

Arguments

Details

Note that level = 1 can only be used if responses are provided, as the variational parameters can't be estimated otherwise. In the absence of responses, level is ignored and the fitted values are returned

Returns

A matrix with predictions scores or counts.

PLNfit$predict_cond()

Predict position, scores or observations of new data, conditionally on the observation of a (set of) variables

Usage

PLNfit$predict_cond(
  newdata,
  cond_responses,
  type = c("link", "response"),
  var_par = FALSE,
  envir = parent.frame()
)

Arguments

Returns

A matrix with predictions scores or counts.

PLNfit$show()

User friendly print method

Usage

PLNfit$show(
  model = paste("A multivariate Poisson Lognormal fit with", self$vcov_model,
    "covariance model.\n")
)

Arguments

PLNfit$print()

User friendly print method

Usage

PLNfit$print()

PLNfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1, data = trichoptera)
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a PLNfit in a standard, general framework, with diagonal residual covariance

Description

The function PLNLDA() produces an instance of an object with class PLNLDAfit.

This class comes with a set of methods, some of them being useful for the user: See the documentation for the methods inherited by PLNfit(), the plot() method for LDA visualization and predict() method for prediction

Super class

PLNfit -> PLNfit_diagonal

Active bindings

Methods

Public methods

• PLNfit_diagonal$new()

• PLNfit_diagonal$clone()

PLNfit_diagonal$new()

Initialize a PLNfit model

Usage

PLNfit_diagonal$new(responses, covariates, offsets, weights, formula, control)

Arguments

PLNfit_diagonal$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_diagonal$clone(deep = FALSE)

Arguments

Super classes

PLNfit -> PLNLDAfit -> PLNLDAfit_spherical

Active bindings

Methods

Public methods

• PLNLDAfit_spherical$new()

• PLNLDAfit_spherical$clone()

PLNLDAfit_spherical$new()

Initialize a PLNfit model

Usage

PLNLDAfit_spherical$new(
  grouping,
  responses,
  covariates,
  offsets,
  weights,
  formula,
  control
)

Arguments

PLNLDAfit_spherical$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit_spherical$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1, data = trichoptera)
class(myPLN)
print(myPLN)

## End(Not run)
## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ 1, data = trichoptera, control = PLN_param(covariance = "spherical"))
class(myPLNLDA)
print(myPLNLDA)

## End(Not run)

An R6 Class to represent a PLNfit in a standard, general framework, with fixed (inverse) residual covariance

Description

An R6 Class to represent a PLNfit in a standard, general framework, with fixed (inverse) residual covariance

Super class

PLNfit -> PLNfit_fixedcov

Active bindings

Methods

Public methods

• PLNfit_fixedcov$new()

• PLNfit_fixedcov$optimize()

• PLNfit_fixedcov$clone()

PLNfit_fixedcov$new()

Initialize a PLNfit model

Usage

PLNfit_fixedcov$new(responses, covariates, offsets, weights, formula, control)

Arguments

PLNfit_fixedcov$optimize()

Call to the NLopt or TORCH optimizer and update of the relevant fields

Usage

PLNfit_fixedcov$optimize(responses, covariates, offsets, weights, config)

Arguments

PLNfit_fixedcov$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_fixedcov$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1, data = trichoptera)
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a PLNfit with a residual covariance structured by a fixed correlation matrix (e.g. a genetic relationship matrix), motivated by population genetics

Description

Sigma = sigma2 * (rho * C + (1 - rho) * I_p), where C is a fixed p x p correlation matrix supplied by the user (control$C) and (sigma2, rho) are estimated. See GeneticCovTraits in src/covariance_pln.h for the C++ side.

Super class

PLNfit -> PLNfit_genpop

Active bindings

Methods

Public methods

• PLNfit_genpop$new()

• PLNfit_genpop$optimize()

• PLNfit_genpop$clone()

PLNfit_genpop$new()

Initialize a PLNfit_genpop model

Usage

PLNfit_genpop$new(responses, covariates, offsets, weights, formula, control)

Arguments

PLNfit_genpop$optimize()

Call to the NLopt or builtin optimizer and update of the relevant fields

Usage

PLNfit_genpop$optimize(responses, covariates, offsets, weights, config)

Arguments

PLNfit_genpop$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_genpop$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
p <- ncol(trichoptera$Abundance)
C <- 0.5^abs(outer(1:p, 1:p, "-")); diag(C) <- 1
myPLN <- PLN(Abundance ~ 1, data = trichoptera, control = PLN_param(covariance = "genpop", C = C))
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a PLNfit in a standard, general framework, with spherical residual covariance

Description

An R6 Class to represent a PLNfit in a standard, general framework, with spherical residual covariance

Super class

PLNfit -> PLNfit_spherical

Active bindings

Methods

Public methods

• PLNfit_spherical$new()

• PLNfit_spherical$clone()

PLNfit_spherical$new()

Initialize a PLNfit model

Usage

PLNfit_spherical$new(responses, covariates, offsets, weights, formula, control)

Arguments

PLNfit_spherical$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_spherical$clone(deep = FALSE)

Arguments

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1, data = trichoptera)
class(myPLN)
print(myPLN)

## End(Not run)

Poisson lognormal mixture model

Description

Fit the mixture variants of the Poisson lognormal with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets).

Usage

PLNmixture(formula, data, subset, clusters = 1:5, control = PLNmixture_param())

Arguments

Value

an R6 object with class PLNmixturefamily, which contains a collection of models with class PLNmixturefit

See Also

The classes PLNmixturefamily, PLNmixturefit and PLNmixture_param()

Examples

## Use parallel to dispatch the computations on 2 workers
## Not run: 
options(mc.cores = 2)

## End(Not run)

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myMixtures <- PLNmixture(Abundance ~ 1 + offset(log(Offset)), clusters = 1:4, data = trichoptera,
                         control = PLNmixture_param(smoothing = 'none'))

## Not run: 
options(mc.cores = 1)

## End(Not run)

Control of a PLNmixture fit

Description

Helper to define list of parameters to control the PLNmixture fit. All arguments have defaults.

Usage

PLNmixture_param(
  backend = c("builtin", "nlopt", "torch"),
  trace = 1,
  covariance = "spherical",
  init_cl = "kmeans",
  smoothing = "both",
  config_optim = list(),
  config_post = list(),
  inception = NULL
)

Arguments

Details

The list of parameters config_optim controls the optimizers. When "nlopt" is chosen the following entries are relevant

• "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)

• "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)

• "profiled" (full covariance only) if TRUE, profile both B and Omega at every nlopt evaluation instead of running an EM loop (Omega fixed for the duration of each inner nlopt solve, B profiled in closed form at every evaluation). Despite the extra O(n*p^2 + p^3) cost per evaluation, benchmarks found it consistently faster than the EM loop (and with a slightly better loglik) across a range of problem sizes. Default is TRUE; set to FALSE to recover the EM loop.

When "torch" backend is used (only for PLN and PLNLDA for now), the following entries are relevant:

• "algorithm" the optimizer used by torch among RPROP (default), RMSPROP, ADAM and ADAGRAD

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10 000

• "numepoch" stop training once this number of epochs exceeds numepoch. Set to -1 to enable infinite training. Default is 1 000

• "num_batch" number of batches to use during training. Defaults to 1 (use full dataset at each epoch)

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "lr" learning rate. Default is 0.1.

• "momentum" momentum factor. Default is 0 (no momentum). Only used in RMSPROP

• "weight_decay" Weight decay penalty. Default is 0 (no decay). Not used in RPROP

• "step_sizes" pair of minimal (default: 1e-6) and maximal (default: 50) allowed step sizes. Only used in RPROP

• "etas" pair of multiplicative increase and decrease factors. Default is (0.5, 1.2). Only used in RPROP

• "centered" if TRUE, compute the centered RMSProp where the gradient is normalized by an estimation of its variance weight_decay (L2 penalty). Default to FALSE. Only used in RMSPROP

When "builtin" backend is used, the following entries are relevant

• "maxeval" stop when the number of Newton steps in the inner loop exceeds maxeval. Default is 10000

• "ftol_in" stop the inner loop when the objective changes by less than ftol_in (relative). Default is 1e-8

• "maxit_em" stop the EM outer loop when the number of EM iterations exceeds maxit_em. Default is 50

• "ftol_em" stop the EM outer loop when the ELBO changes by less than ftol_em (relative). Default is 1e-8

The list of parameters config_post controls the post-treatment processing (for most ⁠PLN*()⁠ functions), with the following entries (defaults may vary depending on the specific function, check ⁠config_post_default_*⁠ for defaults values):

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimation should be used to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

Value

list of parameters configuring the fit.

Outer-loop optimization parameters

PLNmixture_param() adds parameters controlling the EM and smoothing outer loops:

• "ftol_em" outer EM solver stops when the objective changes by less than ftol_em (relative). Default is 1e-3

• "maxit_em" outer EM solver stops when the number of iterations exceeds maxit_em. Default is 50

• "it_smooth" number of the iterations of the smoothing procedure. Default is 1.

See Also

PLN_param()

An R6 Class to represent a collection of PLNmixturefit

Description

The function PLNmixture() produces an instance of this class.

This class comes with a set of methods, some of them being useful for the user: See the documentation for getBestModel(), getModel() and plot().

Super class

PLNfamily -> PLNmixturefamily

Active bindings

Methods

Public methods

• PLNmixturefamily$new()

• PLNmixturefamily$optimize()

• PLNmixturefamily$smooth()

• PLNmixturefamily$plot()

• PLNmixturefamily$plot_objective()

• PLNmixturefamily$getBestModel()

• PLNmixturefamily$show()

• PLNmixturefamily$print()

• PLNmixturefamily$clone()

PLNmixturefamily$new()

helper function for forward smoothing: split a group

Initialize all models in the collection.

Usage

PLNmixturefamily$new(
  clusters,
  responses,
  covariates,
  offsets,
  formula,
  control
)

Arguments

PLNmixturefamily$optimize()

Call to the optimizer on all models of the collection

Usage

PLNmixturefamily$optimize(config)

Arguments

PLNmixturefamily$smooth()

function to restart clustering to avoid local minima by smoothing the loglikelihood values as a function of the number of clusters

Usage

PLNmixturefamily$smooth(control)

Arguments

PLNmixturefamily$plot()

Lineplot of selected criteria for all models in the collection

Usage

PLNmixturefamily$plot(criteria = c("loglik", "BIC", "ICL"), reverse = FALSE)

Arguments

Returns

A ggplot2::ggplot object

PLNmixturefamily$plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

PLNmixturefamily$plot_objective()

Returns

a ggplot2::ggplot graph

PLNmixturefamily$getBestModel()

Extract best model in the collection

Usage

PLNmixturefamily$getBestModel(crit = c("BIC", "ICL", "loglik"))

Arguments

Returns

a PLNmixturefit object

PLNmixturefamily$show()

User friendly print method

Usage

PLNmixturefamily$show()

PLNmixturefamily$print()

User friendly print method

Usage

PLNmixturefamily$print()

PLNmixturefamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNmixturefamily$clone(deep = FALSE)

Arguments

See Also

The function PLNmixture, the class PLNmixturefit

An R6 Class to represent a PLNfit in a mixture framework

Description

The function PLNmixture produces a collection of models which are instances of object with class PLNmixturefit. A PLNmixturefit (say, with k components) is itself a collection of k PLNfit.

This class comes with a set of methods, some of them being useful for the user: See the documentation for ...

Active bindings

Methods

Public methods

• PLNmixturefit$new()

• PLNmixturefit$optimize()

• PLNmixturefit$predict()

• PLNmixturefit$plot_clustering_data()

• PLNmixturefit$plot_clustering_pca()

• PLNmixturefit$postTreatment()

• PLNmixturefit$show()

• PLNmixturefit$print()

• PLNmixturefit$clone()

PLNmixturefit$new()

Optimize a the

Initialize a PLNmixturefit model

Usage

PLNmixturefit$new(
  responses,
  covariates,
  offsets,
  posteriorProb,
  formula,
  control
)

Arguments

PLNmixturefit$optimize()

Optimize a PLNmixturefit model

Usage

PLNmixturefit$optimize(responses, covariates, offsets, config)

Arguments

PLNmixturefit$predict()

Predict group of new samples

Usage

PLNmixturefit$predict(
  newdata,
  type = c("posterior", "response", "position"),
  prior = matrix(rep(1/self$k, self$k), nrow(newdata), self$k, byrow = TRUE),
  control = PLNmixture_param(),
  envir = parent.frame()
)

Arguments

PLNmixturefit$plot_clustering_data()

Plot the matrix of expected mean counts (without offsets, without covariate effects) reordered according the inferred clustering

Usage

PLNmixturefit$plot_clustering_data(
  main = "Expected counts reorder by clustering",
  plot = TRUE,
  log_scale = TRUE
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNmixturefit$plot_clustering_pca()

Plot the individual map of a PCA performed on the latent coordinates, where individuals are colored according to the memberships

Usage

PLNmixturefit$plot_clustering_pca(
  main = "Clustering labels in Individual Factor Map",
  plot = TRUE
)

Arguments

Returns

a ggplot2::ggplot graphic

PLNmixturefit$postTreatment()

Update fields after optimization

Usage

PLNmixturefit$postTreatment(
  responses,
  covariates,
  offsets,
  weights,
  config_post,
  config_optim,
  nullModel
)

Arguments

PLNmixturefit$show()

User friendly print method

Usage

PLNmixturefit$show()

PLNmixturefit$print()

User friendly print method

Usage

PLNmixturefit$print()

PLNmixturefit$clone()

The objects of this class are cloneable with this method.

Usage

PLNmixturefit$clone(deep = FALSE)

Arguments

See Also

The function PLNmixture, the class PLNmixturefamily

Sparse Poisson lognormal model for network inference

Description

Perform sparse inverse covariance estimation for the Zero Inflated Poisson lognormal model using a variational algorithm. Iterate over a range of logarithmically spaced sparsity parameter values. Use the (g)lm syntax to specify the model (including covariates and offsets).

Usage

PLNnetwork(
  formula,
  data,
  subset,
  weights,
  penalties = NULL,
  control = PLNnetwork_param()
)

Arguments

Value

an R6 object with class PLNnetworkfamily, which contains a collection of models with class PLNnetworkfit

See Also

The classes PLNnetworkfamily and PLNnetworkfit, and the and the configuration function PLNnetwork_param().

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)

Control of PLNnetwork fit

Description

Helper to define list of parameters to control the PLN fit. All arguments have defaults.

Usage

PLNnetwork_param(
  backend = c("builtin", "nlopt", "torch"),
  inception_cov = c("full", "spherical", "diagonal"),
  inception_backend = NULL,
  inception_niter = NULL,
  maxit_ve = NULL,
  trace = 1,
  n_penalties = 30,
  min_ratio = 0.1,
  penalize_diagonal = TRUE,
  penalty_weights = NULL,
  config_post = list(),
  config_optim = list(),
  inception = NULL
)

Arguments

Details

The list of parameters config_optim controls the optimizers. When "nlopt" is chosen the following entries are relevant

• "algorithm" the optimization method used by NLOPT among LD type, e.g. "CCSAQ", "MMA", "LBFGS". See NLOPT documentation for further details. Default is "CCSAQ".

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10000

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol_abs. Default is 0.0 (disabled)

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol_abs. Default is 0.0 (disabled)

• "maxtime" stop when the optimization time (in seconds) exceeds maxtime. Default is -1 (disabled)

• "profiled" (full covariance only) if TRUE, profile both B and Omega at every nlopt evaluation instead of running an EM loop (Omega fixed for the duration of each inner nlopt solve, B profiled in closed form at every evaluation). Despite the extra O(n*p^2 + p^3) cost per evaluation, benchmarks found it consistently faster than the EM loop (and with a slightly better loglik) across a range of problem sizes. Default is TRUE; set to FALSE to recover the EM loop.

When "torch" backend is used (only for PLN and PLNLDA for now), the following entries are relevant:

• "algorithm" the optimizer used by torch among RPROP (default), RMSPROP, ADAM and ADAGRAD

• "maxeval" stop when the number of iteration exceeds maxeval. Default is 10 000

• "numepoch" stop training once this number of epochs exceeds numepoch. Set to -1 to enable infinite training. Default is 1 000

• "num_batch" number of batches to use during training. Defaults to 1 (use full dataset at each epoch)

• "ftol_rel" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 1e-8

• "xtol_rel" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 1e-6

• "lr" learning rate. Default is 0.1.

• "momentum" momentum factor. Default is 0 (no momentum). Only used in RMSPROP

• "weight_decay" Weight decay penalty. Default is 0 (no decay). Not used in RPROP

• "step_sizes" pair of minimal (default: 1e-6) and maximal (default: 50) allowed step sizes. Only used in RPROP

• "etas" pair of multiplicative increase and decrease factors. Default is (0.5, 1.2). Only used in RPROP

• "centered" if TRUE, compute the centered RMSProp where the gradient is normalized by an estimation of its variance weight_decay (L2 penalty). Default to FALSE. Only used in RMSPROP

When "builtin" backend is used, the following entries are relevant

• "maxeval" stop when the number of Newton steps in the inner loop exceeds maxeval. Default is 10000

• "ftol_in" stop the inner loop when the objective changes by less than ftol_in (relative). Default is 1e-8

• "maxit_em" stop the EM outer loop when the number of EM iterations exceeds maxit_em. Default is 50

• "ftol_em" stop the EM outer loop when the ELBO changes by less than ftol_em (relative). Default is 1e-8

The list of parameters config_post controls the post-treatment processing (for most ⁠PLN*()⁠ functions), with the following entries (defaults may vary depending on the specific function, check ⁠config_post_default_*⁠ for defaults values):

• jackknife boolean indicating whether jackknife should be performed to evaluate bias and variance of the model parameters. Default is FALSE.

• bootstrap integer indicating the number of bootstrap resamples generated to evaluate the variance of the model parameters. Default is 0 (inactivated).

• variational_var boolean indicating whether variational Fisher information matrix should be computed to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• sandwich_var boolean indicating whether sandwich estimation should be used to estimate the variance of the model parameters (highly underestimated). Default is FALSE.

• rsquared boolean indicating whether approximation of R2 based on deviance should be computed. Default is TRUE

Value

list of parameters configuring the fit.

Outer-loop optimization parameters

PLNnetwork_param() adds two parameters controlling the alternating GLASSO/VEM loop:

• "ftol_em" outer alternating solver stops when the objective changes by less than ftol_em (relative). Default is 1e-5

• "maxit_em" outer alternating solver stops when the number of iterations exceeds maxit_em. Default is 20

See Also

PLN_param()

An R6 Class to represent a collection of PLNnetworkfits

Description

The function PLNnetwork() produces an instance of this class.

This class comes with a set of methods mostly used to compare network fits (in terms of goodness of fit) or extract one from the family (based on penalty parameter and/or goodness of it). See the documentation for getBestModel(), getModel() and plot() for the user-facing ones.

Super classes

PLNfamily -> Networkfamily -> PLNnetworkfamily

Methods

Public methods

• PLNnetworkfamily$new()

• PLNnetworkfamily$stability_selection()

• PLNnetworkfamily$clone()

PLNnetworkfamily$new()

Initialize all models in the collection

Usage

PLNnetworkfamily$new(penalties, data, control)

Arguments

Returns

Update current PLNnetworkfit with smart starting values

PLNnetworkfamily$stability_selection()

Compute the stability path by stability selection

Usage

PLNnetworkfamily$stability_selection(
  subsamples = NULL,
  control = PLNnetwork_param()
)

Arguments

PLNnetworkfamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNnetworkfamily$clone(deep = FALSE)

Arguments

See Also

The function PLNnetwork(), the class PLNnetworkfit

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
class(fits)

An R6 Class to represent a PLNfit in a sparse inverse covariance framework

Description

The function PLNnetwork() produces a collection of models which are instances of object with class PLNnetworkfit. This class comes with a set of methods, some of them being useful for the user: See the documentation for plot() and methods inherited from PLNfit.

Super classes

PLNfit -> PLNfit_fixedcov -> PLNnetworkfit

Active bindings

Methods

Public methods

• PLNnetworkfit$new()

• PLNnetworkfit$optimize()

• PLNnetworkfit$latent_network()

• PLNnetworkfit$plot_network()

• PLNnetworkfit$show()

• PLNnetworkfit$clone()

PLNnetworkfit$new()

Initialize a PLNnetworkfit object

Usage

PLNnetworkfit$new(data, control)

Arguments

PLNnetworkfit$optimize()

Call to the C++ optimizer and update of the relevant fields

Usage

PLNnetworkfit$optimize(data, config)

Arguments

PLNnetworkfit$latent_network()

Extract interaction network in the latent space

Usage

PLNnetworkfit$latent_network(type = c("partial_cor", "support", "precision"))

Arguments

Returns

a square matrix of size PLNnetworkfit$n

PLNnetworkfit$plot_network()

plot the latent network.

Usage

PLNnetworkfit$plot_network(
  type = c("partial_cor", "support"),
  output = c("igraph", "corrplot"),
  edge.color = c("#F8766D", "#00BFC4"),
  remove.isolated = FALSE,
  node.labels = NULL,
  layout = layout_in_circle,
  plot = TRUE
)

Arguments

PLNnetworkfit$show()

User friendly print method

Usage

PLNnetworkfit$show()

PLNnetworkfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNnetworkfit$clone(deep = FALSE)

Arguments

See Also

The function PLNnetwork(), the class PLNnetworkfamily

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
nets <- PLNnetwork(Abundance ~ 1, data = trichoptera)
myPLNnet <- getBestModel(nets)
class(myPLNnet)
print(myPLNnet)

## End(Not run)

Zero Inflated Poisson lognormal model

Description

Fit the multivariate Zero Inflated Poisson lognormal model with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets, subset).

Usage

ZIPLN(
  formula,
  data,
  subset,
  zi = c("single", "row", "col"),
  control = ZIPLN_param()
)

Arguments

Details

Covariates for the Zero-Inflation parameter (using a logistic regression model) can be specified in the formula RHS using the pipe (⁠~ PLN effect | ZI effect⁠) to separate covariates for the PLN part of the model from those for the Zero-Inflation part. Note that different covariates can be used for each part.

Value

an R6 object with class ZIPLNfit

See Also

The class ZIPLNfit

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
## Use different models for zero-inflation...
myZIPLN_single <- ZIPLN(Abundance ~ 1, data = trichoptera, zi = "single")
## Not run: 
myZIPLN_row    <- ZIPLN(Abundance ~ 1, data = trichoptera, zi = "row")
myZIPLN_col    <- ZIPLN(Abundance ~ 1, data = trichoptera, zi = "col")
## ...including logistic regression on covariates
myZIPLN_covar  <- ZIPLN(Abundance ~ 1 | 1 + Wind, data = trichoptera)

## End(Not run)

Control of a ZIPLN fit

Description

Helper to define list of parameters to control the ZIPLN fit. All arguments have defaults.

Usage

ZIPLN_param(
  backend = c("builtin", "nlopt"),
  trace = 1,
  covariance = c("full", "diagonal", "spherical", "fixed", "sparse"),
  Omega = NULL,
  penalty = 0,
  penalize_diagonal = TRUE,
  penalty_weights = NULL,
  config_post = list(),
  config_optim = list(),
  inception = NULL
)

Arguments

Details

See PLN_param() for a description of the generic config_optim entries (ftol_rel, xtol_rel, etc.). Like PLNnetwork_param(), ZIPLN_param() has two parameters controlling the outer EM loop:

• "ftol_out" outer solver stops when an optimization step changes the objective function by less than ftol_out multiplied by the absolute value of the parameter. Default is 1e-6

• "maxit_out" outer solver stops when the number of iteration exceeds maxit_out. Default is 200 for "builtin", 100 for "nlopt" and one additional parameter controlling the form of the variational approximation of the zero inflation:

Value

list of parameters used during the fit and post-processing steps

An R6 Class to represent a ZIPLNfit

Description

The function ZIPLN() fits a model which is an instance of an object with class ZIPLNfit.

This class comes with a set of R6 methods, some of which are useful for the end-user and exported as S3 methods. See the documentation for coef(), sigma(), predict().

Fields are accessed via active binding and cannot be changed by the user.

Details

Covariates for the Zero-Inflation parameter (using a logistic regression model) can be specified in the formula RHS using the pipe (⁠~ PLN effect | ZI effect⁠) to separate covariates for the PLN part of the model from those for the Zero-Inflation part. Note that different covariates can be used for each part.

Active bindings

Methods

Public methods

• ZIPLNfit$update()

• ZIPLNfit$new()

• ZIPLNfit$optimize()

• ZIPLNfit$optimize_vestep()

• ZIPLNfit$predict()

• ZIPLNfit$show()

• ZIPLNfit$print()

• ZIPLNfit$clone()

ZIPLNfit$update()

Update a ZIPLNfit object

Usage

ZIPLNfit$update(
  B = NA,
  B0 = NA,
  Pi = NA,
  Omega = NA,
  Sigma = NA,
  M = NA,
  S2 = NA,
  R = NA,
  Ji = NA,
  Z = NA,
  A = NA,
  monitoring = NA
)

Arguments

Returns

Update the current ZIPLNfit object

ZIPLNfit$new()

Initialize a ZIPLNfit model

Usage

ZIPLNfit$new(data, control)

Arguments

ZIPLNfit$optimize()

Call to the Cpp optimizer and update of the relevant fields

Usage

ZIPLNfit$optimize(data, control)

Arguments

ZIPLNfit$optimize_vestep()

Result of one call to the VE step of the optimization procedure: optimal variational parameters (M, S2, R) and corresponding log likelihood values for fixed model parameters (Sigma, B, B0). Intended to position new data in the latent space.

Usage

ZIPLNfit$optimize_vestep(
  data,
  B = self$model_par$B,
  B0 = self$model_par$B0,
  Omega = self$model_par$Omega,
  control = ZIPLN_param(backend = "nlopt")$config_optim
)

Arguments

Returns

A list with three components:

• the matrix M of variational means,

• the matrix S2 of variational variances

• the matrix R of variational ZI probabilities

• the vector Ji of (variational) log-likelihood of each new observation

• a list monitoring with information about convergence status

ZIPLNfit$predict()

Predict position, scores or observations of new data. See predict.ZIPLNfit() for the S3 method and additional details

Usage

ZIPLNfit$predict(
  newdata,
  responses = NULL,
  type = c("link", "response", "deflated"),
  level = 1,
  envir = parent.frame()
)

Arguments

Returns

A matrix with predictions scores or counts.

ZIPLNfit$show()

User friendly print method

Usage

ZIPLNfit$show(
  model = paste("A multivariate Zero Inflated Poisson Lognormal fit with",
    self$vcov_model, "covariance model.\n")
)

Arguments

ZIPLNfit$print()

User friendly print method

Usage

ZIPLNfit$print()

ZIPLNfit$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNfit$clone(deep = FALSE)

Arguments

Examples

## Not run: 
# See other examples in function ZIPLN
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- ZIPLN(Abundance ~ 1, data = trichoptera)
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a ZIPLNfit in a standard, general framework, with diagonal residual covariance

Description

An R6 Class to represent a ZIPLNfit in a standard, general framework, with diagonal residual covariance

Super class

ZIPLNfit -> ZIPLNfit_diagonal

Active bindings

Methods

Public methods

• ZIPLNfit_diagonal$new()

• ZIPLNfit_diagonal$clone()

ZIPLNfit_diagonal$new()

Initialize a ZIPLNfit_diagonal model

Usage

ZIPLNfit_diagonal$new(data, control)

Arguments

ZIPLNfit_diagonal$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNfit_diagonal$clone(deep = FALSE)

Arguments

Examples

## Not run: 
# See other examples in function ZIPLN
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(covariance = "diagonal"))
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a ZIPLNfit in a standard, general framework, with fixed (inverse) residual covariance

Description

An R6 Class to represent a ZIPLNfit in a standard, general framework, with fixed (inverse) residual covariance

Super class

ZIPLNfit -> ZIPLNfit_fixed

Active bindings

Methods

Public methods

• ZIPLNfit_fixed$new()

• ZIPLNfit_fixed$clone()

ZIPLNfit_fixed$new()

Initialize a ZIPLNfit_fixed model

Usage

ZIPLNfit_fixed$new(data, control)

Arguments

ZIPLNfit_fixed$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNfit_fixed$clone(deep = FALSE)

Arguments

Examples

## Not run: 
# See other examples in function ZIPLN
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- ZIPLN(Abundance ~ 1, data = trichoptera,
    control = ZIPLN_param(Omega = diag(ncol(trichoptera$Abundance))))
class(myPLN)
print(myPLN)

## End(Not run)

An R6 Class to represent a ZIPLNfit in a standard, general framework, with sparse inverse residual covariance

Description

An R6 Class to represent a ZIPLNfit in a standard, general framework, with sparse inverse residual covariance

Super class

ZIPLNfit -> ZIPLNfit_sparse

Active bindings

Methods

Public methods

• ZIPLNfit_sparse$new()

• ZIPLNfit_sparse$latent_network()

• ZIPLNfit_sparse$plot_network()

• ZIPLNfit_sparse$clone()

ZIPLNfit_sparse$new()

Initialize a ZIPLNfit_fixed model

Usage

ZIPLNfit_sparse$new(data, control)

Arguments

ZIPLNfit_sparse$latent_network()

Extract interaction network in the latent space

Usage

ZIPLNfit_sparse$latent_network(type = c("partial_cor", "support", "precision"))

Arguments

Returns

a square matrix of size ZIPLNfit_sparse$n

ZIPLNfit_sparse$plot_network()

plot the latent network.

Usage

ZIPLNfit_sparse$plot_network(
  type = c("partial_cor", "support"),
  output = c("igraph", "corrplot"),
  edge.color = c("#F8766D", "#00BFC4"),
  remove.isolated = FALSE,
  node.labels = NULL,
  layout = layout_in_circle,
  plot = TRUE
)

Arguments

ZIPLNfit_sparse$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNfit_sparse$clone(deep = FALSE)

Arguments

Examples

## Not run: 
# See other examples in function ZIPLN
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- ZIPLN(Abundance ~ 1, data = trichoptera, control=  ZIPLN_param(penalty = 1))
class(myPLN)
print(myPLN)
plot(myPLN)

## End(Not run)

An R6 Class to represent a ZIPLNfit in a standard, general framework, with spherical residual covariance

Description

An R6 Class to represent a ZIPLNfit in a standard, general framework, with spherical residual covariance

Super class

ZIPLNfit -> ZIPLNfit_spherical

Active bindings

Methods

Public methods

• ZIPLNfit_spherical$new()

• ZIPLNfit_spherical$clone()

ZIPLNfit_spherical$new()

Initialize a ZIPLNfit_spherical model

Usage

ZIPLNfit_spherical$new(data, control)

Arguments

ZIPLNfit_spherical$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNfit_spherical$clone(deep = FALSE)

Arguments

Examples

## Not run: 
# See other examples in function ZIPLN
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(covariance = "spherical"))
class(myPLN)
print(myPLN)

## End(Not run)

Zero Inflated Sparse Poisson lognormal model for network inference

Description

Perform sparse inverse covariance estimation for the Zero Inflated Poisson lognormal model using a variational algorithm. Iterate over a range of logarithmically spaced sparsity parameter values. Use the (g)lm syntax to specify the model (including covariates and offsets).

Usage

ZIPLNnetwork(
  formula,
  data,
  subset,
  weights,
  zi = c("single", "row", "col"),
  penalties = NULL,
  control = ZIPLNnetwork_param()
)

Arguments

Details

Covariates for the Zero-Inflation parameter (using a logistic regression model) can be specified in the formula RHS using the pipe (⁠~ PLN effect | ZI effect⁠) to separate covariates for the PLN part of the model from those for the Zero-Inflation part. Note that different covariates can be used for each part.

Value

an R6 object with class ZIPLNnetworkfamily

See Also

The classes ZIPLNfit and ZIPLNnetworkfamily

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myZIPLNs <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, zi = "single")

Control of ZIPLNnetwork fit

Description

Helper to define list of parameters to control the ZIPLNnetwork fit. All arguments have defaults.

Usage

ZIPLNnetwork_param(
  backend = c("builtin", "nlopt"),
  inception_cov = c("full", "spherical", "diagonal"),
  trace = 1,
  n_penalties = 30,
  min_ratio = 0.1,
  penalize_diagonal = TRUE,
  penalty_weights = NULL,
  config_post = list(),
  config_optim = list(),
  inception = NULL
)

Arguments

Details

See PLNnetwork_param() for a full description of the optimization parameters. Note that some defaults values are different than those used in PLNnetwork_param():

• "ftol_out" (outer loop convergence tolerance the objective function) is set by default to 1e-6

• "maxit_out" (max number of iterations for the outer loop) is set by default to 50

Value

list of parameters configuring the fit.

See Also

PLNnetwork_param() and PLN_param()

An R6 Class to represent a collection of ZIPLNnetwork

Description

The function ZIPLNnetwork() produces an instance of this class.

This class comes with a set of methods, some of them being useful for the user: See the documentation for getBestModel(), getModel() and plot()

Super classes

PLNfamily -> Networkfamily -> ZIPLNnetworkfamily

Public fields

Methods

Public methods

• ZIPLNnetworkfamily$new()

• ZIPLNnetworkfamily$stability_selection()

• ZIPLNnetworkfamily$clone()

ZIPLNnetworkfamily$new()

Initialize all models in the collection

Usage

ZIPLNnetworkfamily$new(penalties, data, control)

Arguments

Returns

Update current PLNnetworkfit with smart starting values

ZIPLNnetworkfamily$stability_selection()

Compute the stability path by stability selection

Usage

ZIPLNnetworkfamily$stability_selection(
  subsamples = NULL,
  control = ZIPLNnetwork_param()
)

Arguments

ZIPLNnetworkfamily$clone()

The objects of this class are cloneable with this method.

Usage

ZIPLNnetworkfamily$clone(deep = FALSE)

Arguments

See Also

The function ZIPLNnetwork(), the class ZIPLNfit_sparse

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
class(fits)

Barents fish data set

Description

This data set gives the abundance of 30 fish species observed in 89 sites in the Barents sea. For each site, 4 additional covariates are known. Subsample of the original datasets studied by Fossheim et al, 2006.

Usage

barents

Format

A data frame with 6 variables:

• Abundance: A 30 fish species by 89 sites count matrix

• Offset: A 30 fish species by 89 samples offset matrix, measuring the sampling effort in each site

• 4 covariates for latitude, longitude, depth (in meters), temperature (in Celsius degrees).

Source

Data from M. Fossheim and coauthors.

References

Fossheim, Maria, Einar M. Nilssen, and Michaela Aschan. "Fish assemblages in the Barents Sea." Marine Biology Research 2.4 (2006). doi:10.1080/17451000600815698

Examples

data(barents)

Extracts model coefficients from objects returned by PLNLDA()

Description

The method for objects returned by PLNLDA() only returns coefficients associated to the

\Theta

part of the model (see the PLNLDA vignette for mathematical details).

Usage

## S3 method for class 'PLNLDAfit'
coef(object, ...)

Arguments

Value

Either NULL or a matrix of coefficients extracted from the PLNLDAfit model.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ Wind, grouping = Group, data = trichoptera)
coef(myPLNLDA)

Extract model coefficients

Description

Extracts model coefficients from objects returned by PLN() and its variants

Usage

## S3 method for class 'PLNfit'
coef(object, type = c("main", "covariance"), ...)

Arguments

Value

A matrix of coefficients extracted from the PLNfit model.

See Also

sigma.PLNfit(), vcov.PLNfit(), standard_error.PLNfit()

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera)
coef(myPLN) ## B
coef(myPLN, type = "covariance") ## Sigma

Extract model coefficients

Description

Extracts model coefficients from objects returned by PLN() and its variants

Usage

## S3 method for class 'PLNmixturefit'
coef(object, type = c("main", "means", "covariance", "mixture"), ...)

Arguments

Value

A matrix of coefficients extracted from the PLNfit model.

See Also

sigma.PLNmixturefit()

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLNmixture(Abundance ~ 1 + offset(log(Offset)),
           data = trichoptera, control = PLNmixture_param(smoothing = "none"))  %>% getBestModel()
coef(myPLN) ## Theta - empty here
coef(myPLN, type = "mixture") ## pi
coef(myPLN, type = "means") ## mu
coef(myPLN, type = "covariance") ## Sigma

Extract model coefficients

Description

Extracts model coefficients from objects returned by ZIPLN() and its variants

Usage

## S3 method for class 'ZIPLNfit'
coef(object, type = c("count", "zero", "precision", "covariance"), ...)

Arguments

Value

A matrix of coefficients extracted from the ZIPLNfit model.

See Also

sigma.ZIPLNfit()

Examples

data(scRNA)
# data subsample: only 100 random cell and the 50 most varying transcript
subset <- sample.int(nrow(scRNA), 100)
myPLN  <- ZIPLN(counts[, 1:50] ~ 1 + offset(log(total_counts)), subset = subset, data = scRNA)

Extract the regularization path of a PLNnetwork fit

Description

Extract the regularization path of a PLNnetwork fit

Usage

coefficient_path(Robject, precision = TRUE, corr = TRUE)

Arguments

Value

Sends back a tibble/data.frame.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
head(coefficient_path(fits))

Helper function for PLN initialization.

Description

Barebone function to compute starting points for B, M and S2 when fitting a PLN. Mostly intended for internal use.

Usage

compute_PLN_starting_point(Y, X, O, w, method = c("LM", "GLM"))

Arguments

Details

• B: estimated by weighted LM (method = "LM", default) or p independent Poisson GLMs (method = "GLM"). The GLM option gives better B estimates for factorial or unbalanced designs at the cost of p IRLS fits.

• M: initialized to log((1 + Y) / exp(O)) (M_full in the X*B + M_res parameterization).

• S: initialized element-wise to 1 / sqrt(2 + Y), the approximate VE-step optimum at Omega = I. This adapts automatically to count levels: high S for zero counts (high uncertainty), low S for large counts.

Value

a named list of starting values for model parameter B and variational parameters M and S2 used in the iterative optimization algorithm of PLN()

Examples

## Not run: 
data(barents)
Y <- barents$Abundance
X <- model.matrix(Abundance ~ Latitude + Longitude + Depth + Temperature, data = barents)
O <- log(barents$Offset)
w <- rep(1, nrow(Y))
compute_PLN_starting_point(Y, X, O, w)
compute_PLN_starting_point(Y, X, O, w, method = "GLM")

## End(Not run)

Helper function for ZIPLN initialization.

Description

Fast LM-based starting point for ZIPLN: one multivariate lm.fit for the PLN component and empirical zero rates / binomial GLMs for the ZI component. Replaces the previous per-species pscl::zeroinfl loop.

Usage

compute_ZIPLN_starting_point(Y, X, X0, O, w = NULL)

Arguments

Value

Named list: B (d × p), M (n × p), S2 (n × p), R (n × p), B0 (d0 × p)

Compute offsets from a count data using one of several normalization schemes

Description

Computes offsets from the count table using one of several normalization schemes (TSS, CSS, RLE, GMPR, Wrench, TMM, etc) described in the literature.

Usage

compute_offset(
  counts,
  offset = c("TSS", "GMPR", "RLE", "CSS", "Wrench", "TMM", "none"),
  scale = c("none", "count"),
  ...
)

Arguments

Details

RLE has additional pseudocounts and type arguments to add pseudocounts to the observed counts (defaults to 0L) and to compute offsets using only positive counts (if type == "poscounts"). This mimics the behavior of DESeq2::DESeq() when using sfType == "poscounts". CSS has an additional reference argument to choose the location function used to compute the reference quantiles (defaults to median as in the Nature publication but can be set to mean to reproduce behavior of functions cumNormStat* from metagenomeSeq). Wrench has two additional parameters: groups to specify sample groups and type to either reproduce exactly the default Wrench::wrench() behavior (type = "wrench", default) or to use simpler heuristics (type = "simple"). Note that (i) CSS normalization fails when the median absolute deviation around quantiles does not become instable for high quantiles (limited count variations both within and across samples) and/or one sample has less than two positive counts, (ii) RLE fails when there are no common species across all samples (unless type == "poscounts" has been specified) and (iii) GMPR fails if a sample does not share any species with all other samples. TMM code between two libraries is simplified and adapted from M. Robinson (edgeR:::.calcFactorTMM). The final output is however different from the one produced by edgeR:::.calcFactorTMM as they are intended to be used as such in the model (whereas they need to be multiplied by sequencing depths in edgeR)

Value

If offset = "none", NULL else a vector of length nrow(counts) with one offset per sample.

References

Chen, L., Reeve, J., Zhang, L., Huang, S., Wang, X. and Chen, J. (2018) GMPR: A robust normalization method for zero-inflated count data with application to microbiome sequencing data. PeerJ, 6, e4600 doi:10.7717/peerj.4600

Paulson, J. N., Colin Stine, O., Bravo, H. C. and Pop, M. (2013) Differential abundance analysis for microbial marker-gene surveys. Nature Methods, 10, 1200-1202 doi:10.1038/nmeth.2658

Anders, S. and Huber, W. (2010) Differential expression analysis for sequence count data. Genome Biology, 11, R106 doi:10.1186/gb-2010-11-10-r106

Kumar, M., Slud, E., Okrah, K. et al. (2018) Analysis and correction of compositional bias in sparse sequencing count data. BMC Genomics 19, 799 doi:10.1186/s12864-018-5160-5

Robinson, M.D., Oshlack, A. (2010) A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol 11, R25 doi:10.1186/gb-2010-11-3-r25

Examples

data(trichoptera)
counts <- trichoptera$Abundance
compute_offset(counts)
## Other normalization schemes
compute_offset(counts, offset = "RLE", pseudocounts = 1)
compute_offset(counts, offset = "Wrench", groups = trichoptera$Covariate$Group)
compute_offset(counts, offset = "GMPR")
compute_offset(counts, offset = "TMM")
## User supplied offsets
my_offset <- setNames(rep(1, nrow(counts)), rownames(counts))
compute_offset(counts, offset = my_offset)

Extract edge selection frequency in bootstrap subsamples

Description

Extracts edge selection frequency in networks reconstructed from bootstrap subsamples during the stars stability selection procedure, as either a matrix or a named vector. In the latter case, edge names follow igraph naming convention.

Usage

extract_probs(
  Robject,
  penalty = NULL,
  index = NULL,
  crit = c("StARS", "BIC", "EBIC"),
  format = c("matrix", "vector"),
  tol = 1e-05
)

Arguments

Value

Either a matrix or named vector of edge-wise probabilities. In the latter case, edge names follow igraph convention.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
nets <- PLNnetwork(Abundance ~ 1 + offset(log(Offset)), data = trichoptera)
## Not run: 
stability_selection(nets)
probs <- extract_probs(nets, crit = "StARS", format = "vector")
probs

## End(Not run)

## Not run: 
## Add edge attributes to graph using igraph
net_stars <- getBestModel(nets, "StARS")
g <- plot(net_stars, type = "partial_cor", plot=F)
library(igraph)
E(g)$prob <- probs[as_ids(E(g))]
g

## End(Not run)

Extracts model fitted values from objects returned by PLN() and its variants

Description

Extracts model fitted values from objects returned by PLN() and its variants

Usage

## S3 method for class 'PLNfit'
fitted(object, ...)

Arguments

Value

A matrix of Fitted values extracted from the object object.

Extracts model fitted values from objects returned by PLNmixture() and its variants

Description

Extracts model fitted values from objects returned by PLNmixture() and its variants

Usage

## S3 method for class 'PLNmixturefit'
fitted(object, ...)

Arguments

Value

A matrix of Fitted values extracted from the object object.

Extracts model fitted values from objects returned by ZIPLN() and its variants

Description

Extracts model fitted values from objects returned by ZIPLN() and its variants

Usage

## S3 method for class 'ZIPLNfit'
fitted(object, ...)

Arguments

Value

A matrix of Fitted values extracted from the object object.

Best model extraction from a collection of models

Description

Best model extraction from a collection of models

Usage

## S3 method for class 'PLNPCAfamily'
getBestModel(Robject, crit = c("ICL", "BIC"), ...)

getBestModel(Robject, crit, ...)

## S3 method for class 'PLNmixturefamily'
getBestModel(Robject, crit = c("ICL", "BIC"), ...)

## S3 method for class 'Networkfamily'
getBestModel(Robject, crit = c("BIC", "EBIC", "StARS"), ...)

## S3 method for class 'PLNnetworkfamily'
getBestModel(Robject, crit = c("BIC", "EBIC", "StARS"), ...)

## S3 method for class 'ZIPLNnetworkfamily'
getBestModel(Robject, crit = c("BIC", "EBIC", "StARS"), ...)

Arguments

Value

Send back an object with class PLNPCAfit or PLNnetworkfit

Methods (by class)

• getBestModel(PLNPCAfamily): Model extraction for PLNPCAfamily

• getBestModel(PLNmixturefamily): Model extraction for PLNmixturefamily

• getBestModel(Networkfamily): Model extraction for PLNnetworkfamily or ZIPLNnetworkfamily

• getBestModel(PLNnetworkfamily): Model extraction for PLNnetworkfamily

• getBestModel(ZIPLNnetworkfamily): Model extraction for ZIPLNnetworkfamily

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPCA <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:4)
myModel <- getBestModel(myPCA)

## End(Not run)

Model extraction from a collection of models

Description

Model extraction from a collection of models

Usage

## S3 method for class 'PLNPCAfamily'
getModel(Robject, var, index = NULL)

getModel(Robject, var, index)

## S3 method for class 'PLNmixturefamily'
getModel(Robject, var, index = NULL)

## S3 method for class 'Networkfamily'
getModel(Robject, var, index = NULL)

## S3 method for class 'PLNnetworkfamily'
getModel(Robject, var, index = NULL)

## S3 method for class 'ZIPLNnetworkfamily'
getModel(Robject, var, index = NULL)

Arguments

Value

Sends back an object with class PLNPCAfit or PLNnetworkfit.

Methods (by class)

• getModel(PLNPCAfamily): Model extraction for PLNPCAfamily

• getModel(PLNmixturefamily): Model extraction for PLNmixturefamily

• getModel(Networkfamily): Model extraction for PLNnetworkfamily or ZIPLNnetworkfamily

• getModel(PLNnetworkfamily): Model extraction for PLNnetworkfamily

• getModel(ZIPLNnetworkfamily): Model extraction for ZIPLNnetworkfamily

Examples

## Not run: 
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPCA <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:5)
myModel <- getModel(myPCA, 2)

## End(Not run)

Extract log-likelihood of a fitted PLN model

Description

Returns the variational lower bound of the log-likelihood as a "logLik" object, compatible with stats::AIC() and stats::BIC().

Usage

## S3 method for class 'PLNfit'
logLik(object, ...)

Arguments

Value

An object of class "logLik". The numeric value is the variational ELBO. Attributes df and nobs hold the number of parameters and observations.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- PLN(Abundance ~ 1, data = trichoptera)
logLik(model)

Extract log-likelihood of a fitted ZIPLN model

Description

Returns the variational lower bound of the log-likelihood as a "logLik" object, compatible with stats::AIC() and stats::BIC().

Usage

## S3 method for class 'ZIPLNfit'
logLik(object, ...)

Arguments

Value

An object of class "logLik". The numeric value is the variational ELBO. Attributes df and nobs hold the number of parameters and observations.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
model <- ZIPLN(Abundance ~ 1, data = trichoptera)
logLik(model)

Cow microbiome data set

Description

This data gives the evolution of the microbiota of 45 lactating cows before and after calving in terms of counts of Amplicon Sequence Variants (ASV) in various body sites. Three body sites (vagina, mouth, nose) in addition to the milk of the 4 teats were sampled at 4 times points: 1 week before calving (except for the milk), 1 month, 3 months and 7 months after calving. We present a reduced version of the original data consisting of 880 samples with sequencing depths ranging from 1,001 to 71,881 reads per sample (see doi:10.1186/s42523-023-00252-w for details). ASVs with prevalence in the dataset lower than 5% were filtered out and samples for which the total depth (after ASV filtering were removed), resulting in a count table of n = 880 samples with p = 259 ASV and a mean proportion of zeroes of 90.3%.

Usage

microcosm

Format

A data frame with 6 variables:

• Abundance: A 880 samples by 259 taxa count matrix

• sample: Unique sample id

• Offset: sequencing depth

• site: sampling site (O: oral; N: nasal; V: vaginal; M: milk)

• time: sampling time (⁠-1W⁠: 1 week before calving; ⁠1M⁠: 1 month after calving; ⁠3M⁠: 3 months after calving; ⁠7M⁠: 7 months after calving)

• site_time: factor of possible pairs of (site, time). The combination ⁠M -1W⁠ is absent.

Source

Data from M. Mariadassou and coauthors.

References

Mariadassou, M., Nouvel, L.X., Constant, F. et al. Microbiota members from body sites of dairy cows are largely shared within individual hosts throughout lactation but sharing is limited in the herd. anim microbiome 5, 32 (2023). doi:10.1186/s42523-023-00252-w

See Also

prepare_data()

Examples

data(microcosm)
## Not run: 
my_ZIPLN <- ZIPLN(formula = Abundance ~ 0 + site + offset(log(Offset)) | 1, data = microcosm)

## End(Not run)

Mollusk data set

Description

This data set gives the abundance of 32 mollusk species in 163 samples. For each sample, 4 additional covariates are known.

Usage

mollusk

Format

A list with 2 two data frames:

Abundance

a 163 x 32 data frame of abundancies/counts (163 samples and 32 mollusk species)

Covariate

a 163 x 4 data frame of covariates:

site

a factor with 8 levels indicating the sampling site

season

a factor with 4 levels indicating the season

method

a factor with 2 levels for the method of sampling - wood or string

duration

a numeric with 3 levels for the time of exposure in week

In order to prepare the data for using formula in multivariate analysis (multiple outputs and inputs), use prepare_data(). Original data set has been extracted from ade4.

Source

Data from Richardot-Coulet, Chessel and Bournaud.

References

Richardot-Coulet, M., Chessel D. and Bournaud M. (1986) Typological value of the benthos of old beds of a large river. Methodological approach. Archiv fùr Hydrobiologie, 107, 363–383.

See Also

prepare_data()

Examples

data(mollusk)
mollusc <- prepare_data(mollusk$Abundance, mollusk$Covariate)

Oaks amplicon data set

Description

This data set gives the abundance of 114 taxa (66 bacterial OTU, 48 fungal OTUs) in 116 samples. For each sample, 11 additional covariates are known.

Usage

oaks

Format

A data frame with 13 variables:

• Abundance: A 114 taxa by 116 samples count matrix

• Offset: A 114 taxa by 116 samples offset matrix

• Sample: Unique sample id

• tree: Tree status with respect to the pathogen (susceptible, intermediate or resistant)

• branch: Unique branch id in each tree (4 branches were sampled in each tree, with 10 leaves per branch)

• leafNO: Unique leaf id in each tree (40 leaves were sampled in each tree)

• distTObase: Distance of the sampled leaf to the base of the branch

• distTOtrunk: Distance of the sampled leaf to the base of the tree trunk

• distTOground: Distance of the sampled leaf to the base of the ground

• pmInfection: Powdery mildew infection, proportion of the upper leaf area displaying mildew symptoms

• orientation: Orientation of the branch (South-West SW or North-East NE)

• readsTOTfun: Total number of ITS1 reads for that leaf

• readsTOTbac: Total number of 16S reads for that leaf

Source

Data from B. Jakuschkin and coauthors.

References

Jakuschkin, B., Fievet, V., Schwaller, L. et al. Deciphering the Pathobiome: Intra- and Interkingdom Interactions Involving the Pathogen Erysiphe alphitoides . Microb Ecol 72, 870–880 (2016). doi:10.1007/s00248-016-0777-x

See Also

prepare_data()

Examples

data(oaks)
## Not run: 
oaks_networks <- PLNnetwork(formula = Abundance ~ 1 + offset(log(Offset)), data = oaks)

## End(Not run)

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (either PLNnetworkfamily or ZIPLNnetworkfamily)

Description

Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (either PLNnetworkfamily or ZIPLNnetworkfamily)

Usage

## S3 method for class 'Networkfamily'
plot(
  x,
  type = c("criteria", "stability", "diagnostic"),
  criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
  reverse = FALSE,
  log.x = TRUE,
  stability = 0.9,
  ...
)

## S3 method for class 'PLNnetworkfamily'
plot(
  x,
  type = c("criteria", "stability", "diagnostic"),
  criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
  reverse = FALSE,
  log.x = TRUE,
  stability = 0.9,
  ...
)

## S3 method for class 'ZIPLNnetworkfamily'
plot(
  x,
  type = c("criteria", "stability", "diagnostic"),
  criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
  reverse = FALSE,
  log.x = TRUE,
  stability = 0.9,
  ...
)

Arguments

Details

The BIC and ICL criteria have the form 'loglik - 1/2 * penalty' so that they are on the same scale as the model log-likelihood. You can change this direction and use the alternate form '-2*loglik + penalty', as some authors do, by setting reverse = TRUE.

Value

Produces either a diagnostic plot (with type = 'diagnostic'), a stability plot (with type = 'stability') or the evolution of the criteria of the different models considered (with type = 'criteria', the default).

Functions

• plot(PLNnetworkfamily): Display various outputs associated with a collection of network fits

• plot(ZIPLNnetworkfamily): Display various outputs associated with a collection of network fits

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
## Not run: 
plot(fits)

## End(Not run)

LDA visualization (individual and/or variable factor map(s)) for a PLNPCAfit object

Description

LDA visualization (individual and/or variable factor map(s)) for a PLNPCAfit object

Usage

## S3 method for class 'PLNLDAfit'
plot(
  x,
  map = c("both", "individual", "variable"),
  nb_axes = min(3, x$rank),
  axes = seq.int(min(2, x$rank)),
  var_cols = "var_colors",
  plot = TRUE,
  main = NULL,
  ...
)

Arguments