Title:	Poisson Lognormal Models
Version:	1.2.2
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, future, future.apply, ggplot2, glassoFast, grDevices, grid, gridExtra, igraph, magrittr, MASS, Matrix, methods, nloptr, parallel, pscl, purrr, R6, Rcpp, rlang, scales, stats, tidyr, torch
Suggests:	factoextra, knitr, rmarkdown, spelling, testthat
LinkingTo:	nloptr, Rcpp, RcppArmadillo
VignetteBuilder:	knitr
Encoding:	UTF-8
Language:	en-US
LazyData:	true
RoxygenNote:	7.3.2
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' 'mollusk.R' 'oaks.R' 'plot_utils.R' 'scRNA.R' 'trichoptera.R' 'utils-pipe.R' 'utils-zipln.R' 'utils.R' 'zzz.R'
NeedsCompilation:	yes
Packaged:	2025-03-21 11:05:39 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:	2025-03-21 17:40:06 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:

• 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

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

PLNmodels::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()

Method 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

Method optimize()

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

Usage

Networkfamily$optimize(data, config)

Arguments

Method coefficient_path()

Extract the regularization path of a Networkfamily

Usage

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

Arguments

Method 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.

Method 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

Method plot_stars()

Plot stability path

Usage

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

Arguments

Returns

a ggplot2::ggplot graph

Method plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

Networkfamily$plot_objective()

Returns

a ggplot2::ggplot graph

Method show()

User friendly print method

Usage

Networkfamily$show()

Method 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 = PLN_param())

Arguments

Details

The parameter control is a list controlling the optimization with the following entries:

• "covariance" character setting the model for the covariance matrix. Either "full" or "spherical". Default is "full".

• "trace" integer for verbosity.

• "inception" Set up the initialization. By default, the model is initialized with a multivariate linear model applied on log-transformed data. However, the user can provide a PLNfit (typically obtained from a previous fit), which often speed up the inference.

• "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

• "ftol_abs" stop when an optimization step changes the objective function by less than ftol multiplied by the absolute value of the parameter. Default is 0

• "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

• "xtol_abs" stop when an optimization step changes every parameters by less than xtol multiplied by the absolute value of the parameter. Default is 0

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

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

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

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("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)

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

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

PLNmodels::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()

Method new()

Initialize a PLNLDAfit object

Usage

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

Arguments

Method 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

Method postTreatment()

Update R2, fisher and std_err fields and visualization

Usage

PLNLDAfit$postTreatment(
  grouping,
  responses,
  covariates,
  offsets,
  config_post,
  config_optim
)

Arguments

Method setVisualization()

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

Usage

PLNLDAfit$setVisualization(scale.unit = FALSE)

Arguments

Method 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

Method 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

Method 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

Method 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

Method show()

User friendly print method

Usage

PLNLDAfit$show()

Method 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

PLNmodels::PLNfit -> PLNmodels::PLNLDAfit -> PLNLDAfit_diagonal

Active bindings

Methods

Public methods

• PLNLDAfit_diagonal$new()

• PLNLDAfit_diagonal$clone()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method 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 future to dispatch the computations on 2 workers
## Not run: 
future::plan("multisession", workers = 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)

# Shut down parallel workers
## Not run: 
future::plan("sequential")

## 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 = "nlopt",
  trace = 1,
  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)

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

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

PLNmodels::PLNfamily -> PLNPCAfamily

Active bindings

Methods

Public methods

• PLNPCAfamily$new()

• PLNPCAfamily$optimize()

• PLNPCAfamily$getModel()

• PLNPCAfamily$getBestModel()

• PLNPCAfamily$plot()

• PLNPCAfamily$show()

• PLNPCAfamily$clone()

Method new()

Initialize all models in the collection.

Usage

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

Arguments

Method optimize()

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

Usage

PLNPCAfamily$optimize(config)

Arguments

Method 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

Method getBestModel()

Extract best model in the collection

Usage

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

Arguments

Returns

a PLNPCAfit object

Method 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

Method show()

User friendly print method

Usage

PLNPCAfamily$show()

Method 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

PLNmodels::PLNfit -> PLNPCAfit

Active bindings

Methods

Public methods

• PLNPCAfit$new()

• 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()

Method new()

Initialize a PLNPCAfit object

Usage

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

Arguments

Method update()

Update a PLNPCAfit object

Usage

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

Arguments

Returns

Update the current PLNPCAfit object

Method optimize()

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

Usage

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

Arguments

Method 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

Method 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

Method setVisualization()

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

Usage

PLNPCAfit$setVisualization(scale.unit = FALSE)

Arguments

Method postTreatment()

Update R2, fisher, std_err fields and set up visualization

Usage

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

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

Method 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

Method 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

Method 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

Method show()

User friendly print method

Usage

PLNPCAfit$show()

Method 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", "torch"),
  trace = 1,
  covariance = c("full", "diagonal", "spherical", "fixed"),
  Omega = 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)

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

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()

Method new()

Create a new PLNfamily object.

Usage

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

Arguments

Returns

A new PLNfamily object

Method postTreatment()

Update fields after optimization

Usage

PLNfamily$postTreatment(config_post, config_optim)

Arguments

Method getModel()

Extract a model from a collection of models

Usage

PLNfamily$getModel(var, index = NULL)

Arguments

Returns

A PLNfit object

Method plot()

Lineplot of selected criteria for all models in the collection

Usage

PLNfamily$plot(criteria, reverse)

Arguments

Returns

A ggplot2::ggplot object

Method show()

User friendly print method

Usage

PLNfamily$show()

Method print()

User friendly print method

Usage

PLNfamily$print()

Method 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()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method update()

Update a PLNfit object

Usage

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

Arguments

Returns

Update the current PLNfit object

Method optimize()

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

Usage

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

Arguments

Method 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 (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

Method 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

Method 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.

Method 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.

Method show()

User friendly print method

Usage

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

Arguments

Method print()

User friendly print method

Usage

PLNfit$print()

Method 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

PLNmodels::PLNfit -> PLNfit_diagonal

Active bindings

Methods

Public methods

• PLNfit_diagonal$new()

• PLNfit_diagonal$clone()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_diagonal$clone(deep = FALSE)

Arguments

Super classes

PLNmodels::PLNfit -> PLNmodels::PLNLDAfit -> PLNLDAfit_spherical

Active bindings

Methods

Public methods

• PLNLDAfit_spherical$new()

• PLNLDAfit_spherical$clone()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method 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

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

Super class

PLNmodels::PLNfit -> PLNfit_fixedcov

Active bindings

Methods

Public methods

• PLNfit_fixedcov$new()

• PLNfit_fixedcov$optimize()

• PLNfit_fixedcov$clone()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method optimize()

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

Usage

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

Arguments

Method 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 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

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

Super class

PLNmodels::PLNfit -> PLNfit_spherical

Active bindings

Methods

Public methods

• PLNfit_spherical$new()

• PLNfit_spherical$clone()

Method new()

Initialize a PLNfit model

Usage

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

Arguments

Method 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 future to dispatch the computations on 2 workers
## Not run: 
future::plan("multisession", workers = 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'))

# Shut down parallel workers
## Not run: 
future::plan("sequential")

## 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 = "nlopt",
  trace = 1,
  covariance = "spherical",
  init_cl = "kmeans",
  smoothing = "both",
  config_optim = list(),
  config_post = list(),
  inception = NULL
)

Arguments

Details

See PLN_param() for a full description of the generic optimization parameters. PLNmixture_param() also has additional parameters controlling the optimization due the inner-outer loop structure of the optimizer:

• "ftol_out" outer solver stops when an optimization step changes the objective function by less than xtol 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 50

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

Value

list of parameters configuring the fit.

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

PLNmodels::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()

Method 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

Method optimize()

Call to the optimizer on all models of the collection

Usage

PLNmixturefamily$optimize(config)

Arguments

Method 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

Method 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

Method plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

PLNmixturefamily$plot_objective()

Returns

a ggplot2::ggplot graph

Method getBestModel()

Extract best model in the collection

Usage

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

Arguments

Returns

a PLNmixturefit object

Method show()

User friendly print method

Usage

PLNmixturefamily$show()

Method print()

User friendly print method

Usage

PLNmixturefamily$print()

Method 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()

Method new()

Optimize a the

Initialize a PLNmixturefit model

Usage

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

Arguments

Method optimize()

Optimize a PLNmixturefit model

Usage

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

Arguments

Method 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

Method 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

Method 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

Method postTreatment()

Update fields after optimization

Usage

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

Arguments

Method show()

User friendly print method

Usage

PLNmixturefit$show()

Method print()

User friendly print method

Usage

PLNmixturefit$print()

Method 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("nlopt", "torch"),
  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 PLN_param() for a full description of the generic optimization parameters. PLNnetwork_param() also has two additional parameters controlling the optimization due the inner-outer loop structure of the optimizer:

• "ftol_out" outer solver stops when an optimization step changes the objective function by less than ftol 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 50

Value

list of parameters configuring the fit.

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

PLNmodels::PLNfamily -> PLNmodels::Networkfamily -> PLNnetworkfamily

Methods

Public methods

• PLNnetworkfamily$new()

• PLNnetworkfamily$stability_selection()

• PLNnetworkfamily$clone()

Method new()

Initialize all models in the collection

Usage

PLNnetworkfamily$new(penalties, data, control)

Arguments

Returns

Update current PLNnetworkfit with smart starting values

Method stability_selection()

Compute the stability path by stability selection

Usage

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

Arguments

Method 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

PLNmodels::PLNfit -> PLNmodels::PLNfit_fixedcov -> PLNnetworkfit

Active bindings

Methods

Public methods

• PLNnetworkfit$new()

• PLNnetworkfit$optimize()

• PLNnetworkfit$latent_network()

• PLNnetworkfit$plot_network()

• PLNnetworkfit$show()

• PLNnetworkfit$clone()

Method new()

Initialize a PLNnetworkfit object

Usage

PLNnetworkfit$new(data, control)

Arguments

Method optimize()

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

Usage

PLNnetworkfit$optimize(data, config)

Arguments

Method 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

Method 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

Method show()

User friendly print method

Usage

PLNnetworkfit$show()

Method 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 PLN fit. All arguments have defaults.

Usage

ZIPLN_param(
  backend = c("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() and PLNnetwork_param() for a full description of the generic optimization parameters. Like PLNnetwork_param(), ZIPLN_param() has two parameters controlling the optimization due the inner-outer loop structure of the optimizer:

• "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 100 and one additional parameter controlling the form of the variational approximation of the zero inflation:

• "approx_ZI" either uses an exact or approximated conditional distribution for the zero inflation. Default is FALSE

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()

Method update()

Update a ZIPLNfit object

Usage

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

Arguments

Returns

Update the current ZIPLNfit object

Method new()

Initialize a ZIPLNfit model

Usage

ZIPLNfit$new(data, control)

Arguments

Method optimize()

Call to the Cpp optimizer and update of the relevant fields

Usage

ZIPLNfit$optimize(data, control)

Arguments

Method optimize_vestep()

Result of one call to the VE step of the optimization procedure: optimal variational parameters (M, S, 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 S of variational standard deviations

• 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

Method 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.

Method 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

Method print()

User friendly print method

Usage

ZIPLNfit$print()

Method 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

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

Super class

PLNmodels::ZIPLNfit -> ZIPLNfit_diagonal

Active bindings

Methods

Public methods

• ZIPLNfit_diagonal$new()

• ZIPLNfit_diagonal$clone()

Method new()

Initialize a ZIPLNfit_diagonal model

Usage

ZIPLNfit_diagonal$new(data, control)

Arguments

Method 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

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

Super class

PLNmodels::ZIPLNfit -> ZIPLNfit_fixed

Active bindings

Methods

Public methods

• ZIPLNfit_fixed$new()

• ZIPLNfit_fixed$clone()

Method new()

Initialize a ZIPLNfit_fixed model

Usage

ZIPLNfit_fixed$new(data, control)

Arguments

Method 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

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

Super class

PLNmodels::ZIPLNfit -> ZIPLNfit_sparse

Active bindings

Methods

Public methods

• ZIPLNfit_sparse$new()

• ZIPLNfit_sparse$latent_network()

• ZIPLNfit_sparse$plot_network()

• ZIPLNfit_sparse$clone()

Method new()

Initialize a ZIPLNfit_fixed model

Usage

ZIPLNfit_sparse$new(data, control)

Arguments

Method 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

Method 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

Method 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

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

Super class

PLNmodels::ZIPLNfit -> ZIPLNfit_spherical

Active bindings

Methods

Public methods

• ZIPLNfit_spherical$new()

• ZIPLNfit_spherical$clone()

Method new()

Initialize a ZIPLNfit_spherical model

Usage

ZIPLNfit_spherical$new(data, control)

Arguments

Method 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("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 100

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

PLNmodels::PLNfamily -> PLNmodels::Networkfamily -> ZIPLNnetworkfamily

Public fields

Methods

Public methods

• ZIPLNnetworkfamily$new()

• ZIPLNnetworkfamily$stability_selection()

• ZIPLNnetworkfamily$clone()

Method new()

Initialize all models in the collection

Usage

ZIPLNnetworkfamily$new(penalties, data, control)

Arguments

Returns

Update current PLNnetworkfit with smart starting values

Method stability_selection()

Compute the stability path by stability selection

Usage

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

Arguments

Method 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 S when fitting a PLN. Mostly intended for internal use.

Usage

compute_PLN_starting_point(Y, X, O, w, s = 0.1)

Arguments

Details

The default strategy to estimate B and M is to fit a linear model with covariates X to the response count matrix (after adding a pseudocount of 1, scaling by the offset and taking the log). The regression matrix is used to initialize B and the residuals to initialize M. S is initialized as a constant conformable matrix with value s.

Value

a named list of starting values for model parameter B and variational parameters M and S 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)

## End(Not run)

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)

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

Value

displays an individual and/or variable factor maps for the corresponding axes, and/or sends back a ggplot2::ggplot or gtable object

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLNLDA <- PLNLDA(Abundance ~ 1, grouping = Group, data = trichoptera)
## Not run: 
plot(myPLNLDA, map = "individual", nb_axes = 2)

## End(Not run)

Display the criteria associated with a collection of PLNPCA fits (a PLNPCAfamily)

Description

Display the criteria associated with a collection of PLNPCA fits (a PLNPCAfamily)

Usage

## S3 method for class 'PLNPCAfamily'
plot(x, criteria = c("loglik", "BIC", "ICL"), reverse = FALSE, ...)

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 a plot representing the evolution of the criteria of the different models considered, highlighting the best model in terms of BIC and ICL (see details).

Examples

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

## End(Not run)

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

Description

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

Usage

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

Arguments

Value

displays an individual and/or variable factor maps for the corresponding axes, and/or sends back a ggplot2::ggplot or gtable object

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)
## Not run: 
plot(myPCA, map = "individual", nb_axes=2, ind_cols = trichoptera$Group)
plot(myPCA, map = "variable", nb_axes=2)
plot(myPCA, map = "both", nb_axes=2, ind_cols = trichoptera$Group)

## End(Not run)

Display the criteria associated with a collection of PLN fits (a PLNfamily)

Description

Display the criteria associated with a collection of PLN fits (a PLNfamily)

Usage

## S3 method for class 'PLNfamily'
plot(x, criteria = c("loglik", "BIC", "ICL"), reverse = FALSE, ...)

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 a plot representing the evolution of the criteria of the different models considered, highlighting the best model in terms of BIC and ICL (see details).

See Also

plot.PLNPCAfamily() and plot.PLNnetworkfamily()

Examples

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

## End(Not run)

Display the criteria associated with a collection of PLNmixture fits (a PLNmixturefamily)

Description

Display the criteria associated with a collection of PLNmixture fits (a PLNmixturefamily)

Usage

## S3 method for class 'PLNmixturefamily'
plot(
  x,
  type = c("criteria", "diagnostic"),
  criteria = c("loglik", "BIC", "ICL"),
  reverse = FALSE,
  ...
)

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') or the evolution of the criteria of the different models considered (with type = 'criteria', the default).

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myMixtures <- PLNmixture(Abundance ~ 1 + offset(log(Offset)),
           data = trichoptera, control = PLNmixture_param(smoothing = "none"))
plot(myMixtures, reverse = TRUE)

Mixture visualization of a PLNmixturefit object

Description

Represent the result of the clustering either by coloring the individual in a two-dimension PCA factor map, or by representing the expected matrix of count reorder according to the clustering.

Usage

## S3 method for class 'PLNmixturefit'
plot(x, type = c("pca", "matrix"), main = NULL, plot = TRUE, ...)

Arguments

Value

a ggplot2::ggplot graphic

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()
## Not run: 
plot(myPLN, "pca")
plot(myPLN, "matrix")

## End(Not run)

Extract and plot the network (partial correlation, support or inverse covariance) from a PLNnetworkfit object

Description

Extract and plot the network (partial correlation, support or inverse covariance) from a PLNnetworkfit object

Usage

## S3 method for class 'PLNnetworkfit'
plot(
  x,
  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

Value

Send back an invisible object (igraph or Matrix, depending on the output chosen) and optionally displays a graph (via igraph or corrplot for large ones)

Examples

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

## End(Not run)

Extract and plot the network (partial correlation, support or inverse covariance) from a ZIPLNfit_sparse object

Description

Extract and plot the network (partial correlation, support or inverse covariance) from a ZIPLNfit_sparse object

Usage

## S3 method for class 'ZIPLNfit_sparse'
plot(
  x,
  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

Value

Send back an invisible object (igraph or Matrix, depending on the output chosen) and optionally displays a graph (via igraph or corrplot for large ones)

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
fit <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(penalty = 0.1))
## Not run: 
plot(fit)

## End(Not run)

Predict group of new samples

Description

Predict group of new samples

Usage

## S3 method for class 'PLNLDAfit'
predict(
  object,
  newdata,
  type = c("posterior", "response", "scores"),
  scale = c("log", "prob"),
  prior = NULL,
  control = PLN_param(backend = "nlopt"),
  ...
)

Arguments

Value

A matrix of posterior probabilities for each group (if type = "posterior"), a matrix of (average) scores in the latent space (if type = "scores") or a vector of predicted groups (if type = "response").

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myLDA <- PLNLDA(Abundance ~ 0 + offset(log(Offset)),
                grouping = Group,
                data = trichoptera)
## Not run: 
post_probs <- predict(myLDA, newdata = trichoptera, type = "posterior", scale = "prob")
head(round(post_probs, digits = 3))
predicted_group <- predict(myLDA, newdata = trichoptera, type = "response")
table(predicted_group, trichoptera$Group, dnn = c("predicted", "true"))

## End(Not run)

Predict counts of a new sample

Description

Predict counts of a new sample

Usage

## S3 method for class 'PLNfit'
predict(
  object,
  newdata,
  responses = NULL,
  level = 1,
  type = c("link", "response"),
  ...
)

Arguments

Value

A matrix of predicted log-counts (if type = "link") or predicted counts (if type = "response").

Prediction for a PLNmixturefit object

Description

Predict either posterior probabilities for each group or latent positions based on new samples

Usage

## S3 method for class 'PLNmixturefit'
predict(
  object,
  newdata,
  type = c("posterior", "response", "position"),
  prior = matrix(rep(1/object$k, object$k), nrow(newdata), object$k, byrow = TRUE),
  control = PLNmixture_param(),
  ...
)

Arguments

Value

A matrix of posterior probabilities for each group (if type = "posterior"), a matrix of (average) position in the latent space (if type = "position") or a vector of predicted groups (if type = "response").

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()
predict(myPLN, trichoptera, "posterior")
predict(myPLN, trichoptera, "position")
predict(myPLN, trichoptera, "response")

Predict counts of a new sample

Description

Predict counts of a new sample

Usage

## S3 method for class 'ZIPLNfit'
predict(
  object,
  newdata,
  responses = NULL,
  level = 1,
  type = c("link", "response", "deflated"),
  ...
)

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

Note also that when type = "response" corresponds to predicting values with (1 - \pi)A, where A is the average count in the PLN part of the model and \pi the probability of zero-inflation, whereas type = "deflated" corresponds to A.

Predict counts conditionally

Description

Predict counts of a new sample conditionally on a (set of) observed variables

Usage

predict_cond(
  object,
  newdata,
  cond_responses,
  type = c("link", "response"),
  var_par = FALSE
)

## S3 method for class 'PLNfit'
predict_cond(
  object,
  newdata,
  cond_responses,
  type = c("link", "response"),
  var_par = FALSE
)

Arguments

Value

A list containing:

Methods (by class)

• predict_cond(PLNfit): Predict counts of a new sample conditionally on a (set of) observed variables for a PLNfit

Examples

data(trichoptera)
trichoptera_prep <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ Temperature + Wind, trichoptera_prep)
#Condition on the set of the first two species in the dataset (Hym, Hys) at the ten first sites
Yc <- trichoptera$Abundance[1:10, c(1, 2), drop=FALSE]
newX <- cbind(1, trichoptera$Covariate[1:10, c("Temperature", "Wind")])
pred <- predict_cond(myPLN, newX, Yc, type = "response")

Prepare data for use in PLN models

Description

Prepare data in proper format for use in PLN model and its variants. The function (i) merges a count table and a covariate data frame in the most comprehensive way and (ii) computes offsets from the count table using one of several normalization schemes (TSS, CSS, RLE, GMPR, Wrench, etc). The function fails with informative messages when the heuristics used for sample matching fail.

Usage

prepare_data(
  counts,
  covariates,
  offset = "TSS",
  call = rlang::caller_env(),
  ...
)

Arguments

Value

A data.frame suited for use in PLN() and its variants with two specials components: an abundance count matrix (in component "Abundance") and an offset vector/matrix (in component "Offset", only if offset is not set to "none")

Note

User supplied offsets should be either vectors/column-matrices or have the same number of column as the original count matrix and either (i) dimension names or (ii) the same dimensions as the count matrix. Samples are trimmed in exactly the same way to remove empty samples.

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

See Also

compute_offset() for details on the different normalization schemes

Examples

data(trichoptera)
proper_data <- prepare_data(
 counts     = trichoptera$Abundance,
 covariates = trichoptera$Covariate,
 offset     = "GMPR",
 scale      = "count"
)
proper_data$Abundance
proper_data$Offset

PLN RNG

Description

Random generation for the PLN model with latent mean equal to mu, latent covariance matrix equal to Sigma and average depths (sum of counts in a sample) equal to depths

Usage

rPLN(
  n = 10,
  mu = rep(0, ncol(Sigma)),
  Sigma = diag(1, 5, 5),
  depths = rep(10000, n)
)

Arguments

Details

The default value for mu and Sigma assume equal abundances and no correlation between the different species.

Value

a n * p count matrix, with row-sums close to depths, with an attribute "offsets" corresponding to the true generated offsets (in log-scale).

Examples

## 10 samples of 5 species with equal abundances, no covariance and target depths of 10,000
rPLN()
## 2 samples of 10 highly correlated species with target depths 1,000 and 100,000
## very different abundances
mu <- rep(c(1, -1), each = 5)
Sigma <- matrix(0.8, 10, 10); diag(Sigma) <- 1
rPLN(n=2, mu = mu, Sigma = Sigma, depths = c(1e3, 1e5))

Single cell RNA-seq data

Description

A dataset containing the counts of the 500 most varying transcripts in the mixtures of 5 cell lines in human liver (obtained with standard 10x scRNAseq Chromium protocol).

Usage

scRNA

Format

A data frame named 'scRNA' with 3918 rows (the cells) and 3 variables:

counts

a 500 trancript by 3918 count matrix

cell_line

factor, the cell line of the current row (among 5)

total_counts

Total number of reads for that cell

...

Source

https://github.com/LuyiTian/sc_mixology/

Extract variance-covariance of residuals 'Sigma'

Description

Extract the variance-covariance matrix of the residuals, usually noted

\Sigma

in PLN models. This captures the correlation between the species in the latent space.

Usage

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

Arguments

Value

A semi definite positive matrix of size p, assuming there are p species in the model.

See Also

coef.PLNfit(), standard_error.PLNfit() and vcov.PLNfit() for other ways to access

\Sigma

.

Examples

data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
myPLN <- PLN(Abundance ~ 1 + offset(log(Offset)), data = trichoptera)
sigma(myPLN) ## Sigma

Extract variance-covariance of residuals 'Sigma'

Description

Extract the variance-covariance matrix of the residuals, usually noted

\Sigma

in PLN models. This captures the correlation between the species in the latent space. or PLNmixture, it is a weighted mean of the variance-covariance matrices of each component.

Usage

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

Arguments

Value

A semi definite positive matrix of size p, assuming there are p species in the model.

See Also

coef.PLNmixturefit() for other ways to access

\Sigma

.

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()
sigma(myPLN) ## Sigma

Extract variance-covariance of residuals 'Sigma'

Description

Extract the variance-covariance matrix of the residuals, usually noted \Sigma in ZIPLN models.

Usage

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

Arguments