Package 'PLNmodels'

Description

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

Usage

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

Arguments

Value

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

Examples

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

Description

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

Usage

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

Arguments

Value

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

Examples

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

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)

Description

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

Usage

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

Arguments

Value

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

Examples

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

Description

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

Usage

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

Arguments

Value

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

Examples

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

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

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)

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

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)

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

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)

Description

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

Usage

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

Arguments

Details

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

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

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

Value

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

Examples

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

## End(Not run)

Description

Fast LM-based starting point for ZIPLN: one multivariate lm.fit for the PLN component and empirical zero rates / binomial GLMs for the ZI component. Replaces the previous per-species pscl::zeroinfl loop. R is set to the empirical zero rate of each species on its zero counts, and to 0 elsewhere.

Usage

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

Arguments

Value

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

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.

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)

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)

Description

Sparse estimation of a precision matrix by the graphical Lasso, that is the minimization over positive definite matrices $\Theta$ of

$-\log\det(\Theta) + \mathrm{tr}(S\Theta) + \|\rho \circ \Theta\|_1.$

This is the solver used internally by PLNnetwork() and ZIPLNnetwork().

Usage

graphical_lasso(
  S,
  rho,
  thr = 1e-04,
  maxit = 10000L,
  w_init = NULL,
  wi_init = NULL,
  trace = FALSE,
  stall_patience = 1000L
)

Arguments

Details

The algorithm is the block coordinate descent of Friedman, Hastie and Tibshirani (2008), in the implementation of Sustik and Calderhead (2012): the code is a C++ port of the Fortran routine of the glassoFast package. It first puts the problem on a unit diagonal: for any positive diagonal matrix $D$, $\Theta$ solves the problem for $(S, \rho)$ if and only if $D^{-1}\Theta D^{-1}$ solves it for $(DSD, D\rho D)$, and $D = \mathrm{diag}(S + \rho)^{-1/2}$ is used. This change of variables is exact, and on ordinary input the result is the one of glassoFast up to the stopping rule (same support, entries within the tolerance). On a covariance matrix whose variances span orders of magnitude, as the residual covariance of PLNnetwork() does, it is what makes the descent converge, in a few sweeps where the unscaled one would cycle, fail in its inner loop, or return an indefinite precision matrix. It also departs from glassoFast on degenerate input:

• it always terminates: non-finite input, or a coordinate with $S_{ii} + \rho_{ii} \leq 0$, is rejected (the result is filled with NA and converged is FALSE), and the inner coordinate descent is bounded, where glassoFast can loop forever on a nearly collapsed covariance matrix;

• failure to converge is reported through converged rather than silently;

• it can be interrupted from R;

• when S has no off-diagonal mass, the (diagonal) solution $1 / (S_{ii} + \rho_{ii})$ is returned, where glassoFast returns $1 / \max(\rho_{ii}, \epsilon)$;

• it detects when it is cycling rather than converging, and stops.

That last point is a safeguard. Without the scaling to a unit diagonal, the sweeps could settle into a small limit cycle on an ill-conditioned S: the convergence criterion stops decreasing and oscillates just above its threshold forever, while the solution itself no longer moves. glassoFast spends its whole sweep budget on these and reports success regardless; here the cycle would be detected after stall_patience sweeps without progress, the solve stopped, and status would report "stalled". With the scaling, the covariance matrices on which this was observed (a oaks-derived one, residual covariances along PLNnetwork() paths) all converge.

The remaining statuses are "max_iter" (maxit reached while still progressing), "inner_failure" and "degenerate" (numerical trouble, the result may contain NA).

The precision matrix wi is backed out of the regression coefficients of the algorithm, so that it is the inverse of w only at the exact solution. Short of it, on an ill-conditioned S, it can come out indefinite (whatever the status), and a caller using it as a precision matrix then diverges. Its diagonal is then shifted, which keeps the estimated network, so that its smallest eigenvalue is the one of the inverse of w. The shift is reported in shift, which is 0 otherwise.

Value

a list with components

• w: the estimated covariance matrix,

• wi: the estimated precision matrix (symmetric),

• niter: the number of outer sweeps performed,

• converged: TRUE if the convergence criterion was met,

• status: how the solve ended, one of "converged", "stalled", "max_iter", "inner_failure" or "degenerate" (see Details),

• delta: the best value reached by the convergence criterion, relative to the threshold it had to cross. delta <= 1 means convergence; a stalled solve typically sits between 1 and 3, that is, just short of it,

• shift: the value added to the diagonal of wi to make it positive definite, 0 when it already was (see Details),

• dw_trace: the per-sweep criterion when trace = TRUE.

References

J. Friedman, T. Hastie and R. Tibshirani (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3), 432–441.

M. A. Sustik and B. Calderhead (2012). GLASSOFAST: An efficient GLASSO implementation. UTCS Technical Report TR-12-29, The University of Texas at Austin.

Examples

data(trichoptera)
S <- cov(log1p(as.matrix(trichoptera$Abundance)))
fit <- graphical_lasso(S, rho = 0.1)
fit$converged
sum(fit$wi[upper.tri(fit$wi)] != 0) # number of edges

## penalty weights, leaving the diagonal unpenalized
W <- matrix(1, ncol(S), ncol(S)); diag(W) <- 0
fit <- graphical_lasso(S, rho = 0.1 * W)

Description

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

ICL.PLNfit: ICL for a fitted PLNfit.

ICL.ZIPLNfit: ICL for a fitted ZIPLNfit.

Usage

ICL(object, ...)

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

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

Arguments

Value

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

Examples

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

Description

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

Usage

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

Arguments

Value

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

Examples

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

Description

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

Usage

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

Arguments

Value

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

Examples

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

Description

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

Usage

microcosm

Format

A data frame with 6 variables:

• Abundance: A 880 samples by 259 taxa count matrix

• sample: Unique sample id

• Offset: sequencing depth

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

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

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

Source

Data from M. Mariadassou and coauthors.

References

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

See Also

prepare_data()

Examples

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

## End(Not run)

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)

Description

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

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

Super class

PLNfamily -> Networkfamily

Active bindings

Methods

Public methods

• Networkfamily$new()

• Networkfamily$optimize()

• Networkfamily$postTreatment()

• Networkfamily$coefficient_path()

• Networkfamily$getBestModel()

• Networkfamily$plot()

• Networkfamily$plot_stars()

• Networkfamily$plot_objective()

• Networkfamily$show()

• Networkfamily$clone()

Networkfamily$new()

Initialize all models in the collection

Usage

Networkfamily$new(penalties, data, control)

Arguments

Returns

Update all network fits in the family with smart starting values

Networkfamily$optimize()

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

Usage

Networkfamily$optimize(data, config)

Arguments

Networkfamily$postTreatment()

Update fields after optimization

Usage

Networkfamily$postTreatment(config_post, config_optim)

Arguments

Networkfamily$coefficient_path()

Extract the regularization path of a Networkfamily

Usage

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

Arguments

Networkfamily$getBestModel()

Extract the best network in the family according to some criteria

Usage

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

Arguments

Details

For BIC and EBIC criteria, higher is better.

Networkfamily$plot()

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

Usage

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

Arguments

Returns

a ggplot2::ggplot graph

Networkfamily$plot_stars()

Plot stability path

Usage

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

Arguments

Returns

a ggplot2::ggplot graph

Networkfamily$plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

Networkfamily$plot_objective()

Returns

a ggplot2::ggplot graph

Networkfamily$show()

User friendly print method

Usage

Networkfamily$show()

Networkfamily$clone()

The objects of this class are cloneable with this method.

Usage

Networkfamily$clone(deep = FALSE)

Arguments

See Also

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

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)

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)

Description

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

Usage

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

Arguments

Details

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• "lr" learning rate. Default is 0.1.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Value

list of parameters configuring the fit.

Description

super class for PLNPCAfamily and PLNnetworkfamily.

Public fields

Active bindings

Methods

Public methods

• PLNfamily$new()

• PLNfamily$postTreatment()

• PLNfamily$getModel()

• PLNfamily$plot()

• PLNfamily$show()

• PLNfamily$print()

• PLNfamily$clone()

PLNfamily$new()

Create a new PLNfamily object.

Usage

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

Arguments

Returns

A new PLNfamily object

PLNfamily$postTreatment()

Update fields after optimization

Usage

PLNfamily$postTreatment(config_post, config_optim)

Arguments

PLNfamily$getModel()

Extract a model from a collection of models

Usage

PLNfamily$getModel(var, index = NULL)

Arguments

Returns

A PLNfit object

PLNfamily$plot()

Lineplot of selected criteria for all models in the collection

Usage

PLNfamily$plot(criteria, reverse)

Arguments

Returns

A ggplot2::ggplot object

PLNfamily$show()

User friendly print method

Usage

PLNfamily$show()

PLNfamily$print()

User friendly print method

Usage

PLNfamily$print()

PLNfamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNfamily$clone(deep = FALSE)

Arguments

See Also

getModel()

Description

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

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

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

Active bindings

Methods

Public methods

• PLNfit$new()

• PLNfit$update()

• PLNfit$optimize()

• PLNfit$optimize_vestep()

• PLNfit$postTreatment()

• PLNfit$predict()

• PLNfit$predict_cond()

• PLNfit$show()

• PLNfit$print()

• PLNfit$clone()

PLNfit$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNfit$update()

Update a PLNfit object

Usage

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

Arguments

Returns

Update the current PLNfit object

PLNfit$optimize()

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

Usage

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

Arguments

PLNfit$optimize_vestep()

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

Usage

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

Arguments

Returns

A list with three components:

• the matrix M of variational means,

• the matrix S2 of variational variances

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

PLNfit$postTreatment()

Update R2, fisher and std_err fields after optimization

Usage

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

Arguments

Details

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

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

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

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

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

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

PLNfit$predict()

Predict position, scores or observations of new data.

Usage

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

Arguments

Details

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

Returns

A matrix with predictions scores or counts.

PLNfit$predict_cond()

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

Usage

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

Arguments

Returns

A matrix with predictions scores or counts.

PLNfit$show()

User friendly print method

Usage

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

Arguments

PLNfit$print()

User friendly print method

Usage

PLNfit$print()

PLNfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

Super class

PLNfit -> PLNfit_diagonal

Active bindings

Methods

Public methods

• PLNfit_diagonal$new()

• PLNfit_diagonal$clone()

PLNfit_diagonal$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNfit_diagonal$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_diagonal$clone(deep = FALSE)

Arguments

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)

Description

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

Super class

PLNfit -> PLNfit_fixedcov

Active bindings

Methods

Public methods

• PLNfit_fixedcov$new()

• PLNfit_fixedcov$optimize()

• PLNfit_fixedcov$clone()

PLNfit_fixedcov$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNfit_fixedcov$optimize()

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

Usage

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

Arguments

PLNfit_fixedcov$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_fixedcov$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

Super class

PLNfit -> PLNfit_genpop

Active bindings

Methods

Public methods

• PLNfit_genpop$new()

• PLNfit_genpop$clone()

PLNfit_genpop$new()

Initialize a PLNfit_genpop model

Usage

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

Arguments

PLNfit_genpop$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_genpop$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

Super class

PLNfit -> PLNfit_spherical

Active bindings

Methods

Public methods

• PLNfit_spherical$new()

• PLNfit_spherical$clone()

PLNfit_spherical$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNfit_spherical$clone()

The objects of this class are cloneable with this method.

Usage

PLNfit_spherical$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

Usage

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

Arguments

Details

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

Value

an R6 object with class PLNLDAfit()

See Also

The class PLNLDAfit

Examples

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

Description

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

Usage

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

Arguments

Details

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• "lr" learning rate. Default is 0.1.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Value

list of parameters configuring the fit.

Description

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

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

Super class

PLNfit -> PLNLDAfit

Active bindings

Methods

Public methods

• PLNLDAfit$new()

• PLNLDAfit$optimize()

• PLNLDAfit$postTreatment()

• PLNLDAfit$setVisualization()

• PLNLDAfit$plot_individual_map()

• PLNLDAfit$plot_correlation_map()

• PLNLDAfit$plot_LDA()

• PLNLDAfit$predict()

• PLNLDAfit$show()

• PLNLDAfit$clone()

PLNLDAfit$new()

Initialize a PLNLDAfit object

Usage

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

Arguments

PLNLDAfit$optimize()

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

Usage

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

Arguments

PLNLDAfit$postTreatment()

Update R2, fisher and std_err fields and visualization

Usage

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

Arguments

PLNLDAfit$setVisualization()

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

Usage

PLNLDAfit$setVisualization(scale.unit = FALSE)

Arguments

PLNLDAfit$plot_individual_map()

Plot the factorial map of the LDA

Usage

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

Arguments

Returns

a ggplot2::ggplot graphic

PLNLDAfit$plot_correlation_map()

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

Usage

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

Arguments

Returns

a ggplot2::ggplot graphic

PLNLDAfit$plot_LDA()

Plot a summary of the PLNLDAfit object

Usage

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

Arguments

Returns

a grob object

PLNLDAfit$predict()

Predict group of new samples

Usage

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

Arguments

PLNLDAfit$show()

User friendly print method

Usage

PLNLDAfit$show()

PLNLDAfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit$clone(deep = FALSE)

Arguments

See Also

The function PLNLDA.

Examples

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

## End(Not run)

Description

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

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

Super classes

PLNfit -> PLNLDAfit -> PLNLDAfit_diagonal

Active bindings

Methods

Public methods

• PLNLDAfit_diagonal$new()

• PLNLDAfit_diagonal$clone()

PLNLDAfit_diagonal$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNLDAfit_diagonal$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit_diagonal$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

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

Super classes

PLNfit -> PLNLDAfit -> PLNLDAfit_spherical

Active bindings

Methods

Public methods

• PLNLDAfit_spherical$new()

• PLNLDAfit_spherical$clone()

PLNLDAfit_spherical$new()

Initialize a PLNfit model

Usage

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

Arguments

PLNLDAfit_spherical$clone()

The objects of this class are cloneable with this method.

Usage

PLNLDAfit_spherical$clone(deep = FALSE)

Arguments

Examples

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

## End(Not run)

Description

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

Usage

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

Arguments

Value

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

See Also

The classes PLNmixturefamily, PLNmixturefit and PLNmixture_param()

Examples

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

## End(Not run)

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

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

## End(Not run)

Description

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

Usage

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

Arguments

Details

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• "lr" learning rate. Default is 0.1.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Value

list of parameters configuring the fit.

Outer-loop optimization parameters

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

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

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

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

See Also

PLN_param()

Description

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

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

Super class

PLNfamily -> PLNmixturefamily

Active bindings

Methods

Public methods

• PLNmixturefamily$new()

• PLNmixturefamily$optimize()

• PLNmixturefamily$smooth()

• PLNmixturefamily$plot()

• PLNmixturefamily$plot_objective()

• PLNmixturefamily$getBestModel()

• PLNmixturefamily$show()

• PLNmixturefamily$print()

• PLNmixturefamily$clone()

PLNmixturefamily$new()

Initialize all models in the collection.

Usage

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

Arguments

PLNmixturefamily$optimize()

Call to the optimizer on all models of the collection

Usage

PLNmixturefamily$optimize(config)

Arguments

PLNmixturefamily$smooth()

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

Usage

PLNmixturefamily$smooth(control)

Arguments

PLNmixturefamily$plot()

Lineplot of selected criteria for all models in the collection

Usage

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

Arguments

Returns

A ggplot2::ggplot object

PLNmixturefamily$plot_objective()

Plot objective value of the optimization problem along the penalty path

Usage

PLNmixturefamily$plot_objective()

Returns

a ggplot2::ggplot graph

PLNmixturefamily$getBestModel()

Extract best model in the collection

Usage

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

Arguments

Returns

a PLNmixturefit object

PLNmixturefamily$show()

User friendly print method

Usage

PLNmixturefamily$show()

PLNmixturefamily$print()

User friendly print method

Usage

PLNmixturefamily$print()

PLNmixturefamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNmixturefamily$clone(deep = FALSE)

Arguments

See Also

The function PLNmixture, the class PLNmixturefit

Description

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

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

Active bindings

Methods

Public methods

• PLNmixturefit$new()

• PLNmixturefit$optimize()

• PLNmixturefit$predict()

• PLNmixturefit$plot_clustering_data()

• PLNmixturefit$plot_clustering_pca()

• PLNmixturefit$postTreatment()

• PLNmixturefit$show()

• PLNmixturefit$print()

• PLNmixturefit$clone()

PLNmixturefit$new()

Optimize a the

Initialize a PLNmixturefit model

Usage

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

Arguments

PLNmixturefit$optimize()

Optimize a PLNmixturefit model

Usage

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

Arguments

PLNmixturefit$predict()

Predict group of new samples

Usage

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

Arguments

PLNmixturefit$plot_clustering_data()

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

Usage

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

Arguments

Returns

a ggplot2::ggplot graphic

PLNmixturefit$plot_clustering_pca()

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

Usage

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

Arguments

Returns

a ggplot2::ggplot graphic

PLNmixturefit$postTreatment()

Update fields after optimization

Usage

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

Arguments

PLNmixturefit$show()

User friendly print method

Usage

PLNmixturefit$show()

PLNmixturefit$print()

User friendly print method

Usage

PLNmixturefit$print()

PLNmixturefit$clone()

The objects of this class are cloneable with this method.

Usage

PLNmixturefit$clone(deep = FALSE)

Arguments

See Also

The function PLNmixture, the class PLNmixturefamily

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)

Description

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

Usage

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

Arguments

Details

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

• "lr" learning rate. Default is 0.1.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Value

list of parameters configuring the fit.

Outer-loop optimization parameters

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

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

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

Scale of the penalty

With penalty_scale = "correlation", the penalty on the pair $(i, j)$ is $\lambda w_{ij} \sqrt{S_{ii} S_{jj}}$, recomputed at each M step from the current residual covariance $S$. Without weights, the precision matrix is $\Omega = D^{-1/2} K D^{-1/2}$, where $D$ is the diagonal of $S$ and $K$ the graphical-Lasso estimate on the correlation matrix $D^{-1/2} S D^{-1/2}$: this is the correlation-based estimator of Rothman et al. (2008), who show that it has a better rate of convergence in operator norm than the estimator on the covariance scale. The grid of penalties is built on the residual correlation of the inception, and lies between 0 and 1. Since the weights depend on $S$, the alternating optimization does not maximize a fixed penalized criterion: it looks for a fixed point.

It is the default since version 1.3.3. In simulations with a known network (a hundred configurations of graph, sample size, dimension, latent variances and abundances, with or without species absent from groups of samples), the edges were recovered better on the correlation scale in every configuration, the more so as the latent variances differed between species, and the time was the same.

To get the behavior of version 1.3.2 back, use PLNnetwork_param(penalty_scale = "covariance"); penalties are then on the covariance scale. Penalties given through penalties are otherwise taken on the correlation scale. Since values above 1 cannot be on that scale (they all give the empty network), when some are given while penalty_scale is left to its default they are taken as penalties on the covariance scale and converted, with a warning: they are divided by the ratio of the largest residual covariance of the inception to its largest residual correlation, so that the penalty giving the empty network on one scale gives it on the other. This is exact when all latent variances are equal, and only a guide otherwise: the two scales do not give the same networks. Setting penalty_scale explicitly, to either value, leaves the penalties as given.

On either scale, the latent variance of a species absent from whole groups of samples diverges along the path without latent_floor. On the covariance scale, the floor stabilizes the fit without repairing the network: the degenerate species have fewer edges, but the edges still concentrate on them.

References

Rothman, A. J., Bickel, P. J., Levina, E. and Zhu, J. (2008). Sparse permutation invariant covariance estimation. Electronic Journal of Statistics, 2, 494–515. doi:10.1214/08-EJS176

See Also

PLN_param()

An R6 Class to represent a collection of PLNnetworkfits

Description

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

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

Super classes

PLNfamily -> Networkfamily -> PLNnetworkfamily

Methods

Public methods

• PLNnetworkfamily$new()

• PLNnetworkfamily$stability_selection()

• PLNnetworkfamily$clone()

PLNnetworkfamily$new()

Initialize all models in the collection

Usage

PLNnetworkfamily$new(penalties, data, control)

Arguments

Returns

Update current PLNnetworkfit with smart starting values

PLNnetworkfamily$stability_selection()

Compute the stability path by stability selection

Usage

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

Arguments

PLNnetworkfamily$clone()

The objects of this class are cloneable with this method.

Usage

PLNnetworkfamily$clone(deep = FALSE)

Arguments

See Also

The function PLNnetwork(), the class PLNnetworkfit

Examples

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

Description

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

Super classes

PLNfit -> PLNfit_fixedcov -> PLNnetworkfit

Active bindings

Methods

Public methods

• PLNnetworkfit$new()

• PLNnetworkfit$optimize()

• PLNnetworkfit$latent_network()

• PLNnetworkfit$plot_network()

• PLNnetworkfit$show()

• PLNnetworkfit$clone()

PLNnetworkfit$new()

Initialize a PLNnetworkfit object

Usage

PLNnetworkfit$new(data, control)

Arguments

PLNnetworkfit$optimize()

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

Usage

PLNnetworkfit$optimize(data, config)

Arguments

PLNnetworkfit$latent_network()

Extract interaction network in the latent space

Usage

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

Arguments

Returns

a square matrix of size PLNnetworkfit$n

PLNnetworkfit$plot_network()

plot the latent network.

Usage

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

Arguments

PLNnetworkfit$show()

User friendly print method

Usage

PLNnetworkfit$show()

PLNnetworkfit$clone()

The objects of this class are cloneable with this method.

Usage

PLNnetworkfit$clone(deep = FALSE)

Arguments

See Also

The function PLNnetwork(), the class PLNnetworkfamily

Examples

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

## End(Not run)

Description

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

Usage

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

Arguments

Value

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

See Also

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

Examples

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

## End(Not run)

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

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

## End(Not run)