}}
\describe{
-\item{\code{nb_param}}{number of parameters in the current PLN model}
+\item{\code{nb_param_pln}}{number of parameters in the PLN part of the current model}
\item{\code{vcov_model}}{character: the model used for the residual covariance}
}
@@ -55,28 +55,13 @@ print(myPLN)
\subsection{Method \code{new()}}{
Initialize a \code{\link{ZIPLNfit_spherical}} model
\subsection{Usage}{
-\if{html}{\out{
}}\preformatted{ZIPLNfit_spherical$new(
- responses,
- covariates,
- offsets,
- weights,
- formula,
- control
-)}\if{html}{\out{
}}
\describe{
-\item{\code{responses}}{the matrix of responses (called Y in the model). Will usually be extracted from the corresponding field in PLNfamily-class}
-
-\item{\code{covariates}}{design matrix (called X in the model). Will usually be extracted from the corresponding field in PLNfamily-class}
-
-\item{\code{offsets}}{offset matrix (called O in the model). Will usually be extracted from the corresponding field in PLNfamily-class}
-
-\item{\code{weights}}{an optional vector of observation weights to be used in the fitting process.}
-
-\item{\code{formula}}{model formula used for fitting, extracted from the formula in the upper-level call}
+\item{\code{data}}{a named list used internally to carry the data matrices}
\item{\code{control}}{a list for controlling the optimization. See details.}
}
diff --git a/man/ZIPLNnetwork.Rd b/man/ZIPLNnetwork.Rd
new file mode 100644
index 00000000..aef11b01
--- /dev/null
+++ b/man/ZIPLNnetwork.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ZIPLNnetwork.R
+\name{ZIPLNnetwork}
+\alias{ZIPLNnetwork}
+\title{Zero Inflated Sparse Poisson lognormal model for network inference}
+\usage{
+ZIPLNnetwork(
+ formula,
+ data,
+ subset,
+ weights,
+ zi = c("single", "row", "col"),
+ penalties = NULL,
+ control = ZIPLNnetwork_param()
+)
+}
+\arguments{
+\item{formula}{an object of class "formula": a symbolic description of the model to be fitted.}
+
+\item{data}{an optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which lm is called.}
+
+\item{subset}{an optional vector specifying a subset of observations to be used in the fitting process.}
+
+\item{weights}{an optional vector of observation weights to be used in the fitting process.}
+
+\item{zi}{a character describing the model used for zero inflation, either of
+\itemize{
+\item "single" (default, one parameter shared by all counts)
+\item "col" (one parameter per variable / feature)
+\item "row" (one parameter per sample / individual).
+If covariates are specified in the formula RHS (see details) this parameter is ignored.
+}}
+
+\item{penalties}{an optional vector of positive real number controlling the level of sparsity of the underlying network. if NULL (the default), will be set internally. See \code{PLNnetwork_param()} for additional tuning of the penalty.}
+
+\item{control}{a list-like structure for controlling the optimization, with default generated by \code{\link[=ZIPLNnetwork_param]{ZIPLNnetwork_param()}}. See the associated documentation
+for details.}
+}
+\value{
+an R6 object with class \code{\link{ZIPLNnetworkfamily}}
+}
+\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).
+}
+\details{
+Covariates for the Zero-Inflation parameter (using a logistic regression model) can be specified in the formula RHS using the pipe
+(\verb{~ 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.
+}
+\examples{
+data(trichoptera)
+trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
+myZIPLNs <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, zi = "single")
+}
+\seealso{
+The classes \code{\link{ZIPLNfit}} and \code{\link{ZIPLNnetworkfamily}}
+}
diff --git a/man/ZIPLNnetwork_param.Rd b/man/ZIPLNnetwork_param.Rd
new file mode 100644
index 00000000..fe3d4daa
--- /dev/null
+++ b/man/ZIPLNnetwork_param.Rd
@@ -0,0 +1,58 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ZIPLNnetwork.R
+\name{ZIPLNnetwork_param}
+\alias{ZIPLNnetwork_param}
+\title{Control of ZIPLNnetwork fit}
+\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{
+\item{backend}{optimization back used, either "nlopt" or "torch". Default is "nlopt"}
+
+\item{inception_cov}{Covariance structure used for the inception model used to initialize the PLNfamily. Defaults to "full" and can be constrained to "diagonal" and "spherical".}
+
+\item{trace}{a integer for verbosity.}
+
+\item{n_penalties}{an integer that specifies the number of values for the penalty grid when internally generated. Ignored when penalties is non \code{NULL}}
+
+\item{min_ratio}{the penalty grid ranges from the minimal value that produces a sparse to this value multiplied by \code{min_ratio}. Default is 0.1.}
+
+\item{penalize_diagonal}{boolean: should the diagonal terms be penalized in the graphical-Lasso? Default is \code{TRUE}}
+
+\item{penalty_weights}{either a single or a list of p x p matrix of weights (default: all weights equal to 1) to adapt the amount of shrinkage to each pairs of node. Must be symmetric with positive values.}
+
+\item{config_post}{a list for controlling the post-treatment (optional bootstrap, jackknife, R2, etc).}
+
+\item{config_optim}{a list for controlling the optimizer (either "nlopt" or "torch" backend). See details}
+
+\item{inception}{Set up the parameters initialization: by default, the model is initialized with a multivariate linear model applied on
+log-transformed data, and with the same formula as the one provided by the user. However, the user can provide a PLNfit (typically obtained from a previous fit),
+which sometimes speeds up the inference.}
+}
+\value{
+list of parameters configuring the fit.
+}
+\description{
+Helper to define list of parameters to control the ZIPLNnetwork fit. All arguments have defaults.
+}
+\details{
+See \code{\link[=PLNnetwork_param]{PLNnetwork_param()}} for a full description of the optimization parameters. Note that some defaults values are different than those used in \code{\link[=PLNnetwork_param]{PLNnetwork_param()}}:
+\itemize{
+\item "ftol_out" (outer loop convergence tolerance the objective function) is set by default to 1e-6
+\item "maxit_out" (max number of iterations for the outer loop) is set by default to 100
+}
+}
+\seealso{
+\code{\link[=PLNnetwork_param]{PLNnetwork_param()}} and \code{\link[=PLN_param]{PLN_param()}}
+}
diff --git a/man/ZIPLNnetworkfamily.Rd b/man/ZIPLNnetworkfamily.Rd
new file mode 100644
index 00000000..0bf90623
--- /dev/null
+++ b/man/ZIPLNnetworkfamily.Rd
@@ -0,0 +1,119 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PLNnetworkfamily-class.R
+\name{ZIPLNnetworkfamily}
+\alias{ZIPLNnetworkfamily}
+\title{An R6 Class to represent a collection of ZIPLNnetwork}
+\description{
+The function \code{\link[=ZIPLNnetwork]{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 \code{\link[=getBestModel]{getBestModel()}},
+\code{\link[=getModel]{getModel()}} and \link[=plot.ZIPLNnetworkfamily]{plot()}
+}
+\examples{
+data(trichoptera)
+trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
+fits <- PLNnetwork(Abundance ~ 1, data = trichoptera)
+class(fits)
+}
+\seealso{
+The function \code{\link[=ZIPLNnetwork]{ZIPLNnetwork()}}, the class \code{\link{ZIPLNfit_sparse}}
+}
+\section{Super classes}{
+\code{\link[PLNmodels:PLNfamily]{PLNmodels::PLNfamily}} -> \code{\link[PLNmodels:Networkfamily]{PLNmodels::Networkfamily}} -> \code{ZIPLNnetworkfamily}
+}
+\section{Public fields}{
+\if{html}{\out{
}}
+\describe{
+\item{\code{covariates0}}{the matrix of covariates included in the ZI component}
+}
+\if{html}{\out{
}}
+}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-ZIPLNnetworkfamily-new}{\code{ZIPLNnetworkfamily$new()}}
+\item \href{#method-ZIPLNnetworkfamily-stability_selection}{\code{ZIPLNnetworkfamily$stability_selection()}}
+\item \href{#method-ZIPLNnetworkfamily-clone}{\code{ZIPLNnetworkfamily$clone()}}
+}
+}
+\if{html}{\out{
+
Inherited methods
+
+
+}}
+\if{html}{\out{
}}
+\if{html}{\out{
}}
+\if{latex}{\out{\hypertarget{method-ZIPLNnetworkfamily-new}{}}}
+\subsection{Method \code{new()}}{
+Initialize all models in the collection
+\subsection{Usage}{
+\if{html}{\out{
}}\preformatted{ZIPLNnetworkfamily$new(penalties, data, control)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{
}}
+\describe{
+\item{\code{penalties}}{a vector of positive real number controlling the level of sparsity of the underlying network.}
+
+\item{\code{data}}{a named list used internally to carry the data matrices}
+
+\item{\code{control}}{a list for controlling the optimization.}
+}
+\if{html}{\out{
}}
+}
+\subsection{Returns}{
+Update current \code{\link{PLNnetworkfit}} with smart starting values
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{
}}
+\if{latex}{\out{\hypertarget{method-ZIPLNnetworkfamily-stability_selection}{}}}
+\subsection{Method \code{stability_selection()}}{
+Compute the stability path by stability selection
+\subsection{Usage}{
+\if{html}{\out{
}}\preformatted{ZIPLNnetworkfamily$stability_selection(
+ subsamples = NULL,
+ control = ZIPLNnetwork_param()
+)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{
}}
+\describe{
+\item{\code{subsamples}}{a list of vectors describing the subsamples. The number of vectors (or list length) determines the number of subsamples used in the stability selection. Automatically set to 20 subsamples with size \code{10*sqrt(n)} if \code{n >= 144} and \code{0.8*n} otherwise following Liu et al. (2010) recommendations.}
+
+\item{\code{control}}{a list controlling the main optimization process in each call to \code{\link[=PLNnetwork]{PLNnetwork()}}. See \code{\link[=ZIPLNnetwork]{ZIPLNnetwork()}} and \code{\link[=ZIPLN_param]{ZIPLN_param()}} for details.}
+}
+\if{html}{\out{
}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{
}}
+\if{latex}{\out{\hypertarget{method-ZIPLNnetworkfamily-clone}{}}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{
}}\preformatted{ZIPLNnetworkfamily$clone(deep = FALSE)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{
}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{
}}
+}
+}
+}
diff --git a/man/coefficient_path.Rd b/man/coefficient_path.Rd
index 3a389ccb..45552555 100644
--- a/man/coefficient_path.Rd
+++ b/man/coefficient_path.Rd
@@ -7,7 +7,7 @@
coefficient_path(Robject, precision = TRUE, corr = TRUE)
}
\arguments{
-\item{Robject}{an object with class \code{\link{PLNnetworkfamily}}, i.e. an output from \code{\link[=PLNnetwork]{PLNnetwork()}}}
+\item{Robject}{an object with class \code{\link{Networkfamily}}, i.e. an output from \code{\link[=PLNnetwork]{PLNnetwork()}}}
\item{precision}{a logical, should the coefficients of the precision matrix Omega or the covariance matrix Sigma be sent back. Default is \code{TRUE}.}
diff --git a/man/getBestModel.Rd b/man/getBestModel.Rd
index aff5ac35..bce004c4 100644
--- a/man/getBestModel.Rd
+++ b/man/getBestModel.Rd
@@ -6,7 +6,9 @@
\alias{getBestModel.PLNPCAfamily}
\alias{getBestModel}
\alias{getBestModel.PLNmixturefamily}
+\alias{getBestModel.Networkfamily}
\alias{getBestModel.PLNnetworkfamily}
+\alias{getBestModel.ZIPLNnetworkfamily}
\title{Best model extraction from a collection of models}
\usage{
\method{getBestModel}{PLNPCAfamily}(Robject, crit = c("ICL", "BIC"), ...)
@@ -15,7 +17,11 @@ getBestModel(Robject, crit, ...)
\method{getBestModel}{PLNmixturefamily}(Robject, crit = c("ICL", "BIC"), ...)
+\method{getBestModel}{Networkfamily}(Robject, crit = c("BIC", "EBIC", "StARS"), ...)
+
\method{getBestModel}{PLNnetworkfamily}(Robject, crit = c("BIC", "EBIC", "StARS"), ...)
+
+\method{getBestModel}{ZIPLNnetworkfamily}(Robject, crit = c("BIC", "EBIC", "StARS"), ...)
}
\arguments{
\item{Robject}{an object with class PLNPCAfamilly ot PLNnetworkfamily}
@@ -39,8 +45,12 @@ Best model extraction from a collection of models
\item \code{getBestModel(PLNmixturefamily)}: Model extraction for \code{\link{PLNmixturefamily}}
+\item \code{getBestModel(Networkfamily)}: Model extraction for \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}}
+
\item \code{getBestModel(PLNnetworkfamily)}: Model extraction for \code{\link{PLNnetworkfamily}}
+\item \code{getBestModel(ZIPLNnetworkfamily)}: Model extraction for \code{\link{ZIPLNnetworkfamily}}
+
}}
\examples{
\dontrun{
diff --git a/man/getModel.Rd b/man/getModel.Rd
index 846dbf8b..a812d321 100644
--- a/man/getModel.Rd
+++ b/man/getModel.Rd
@@ -6,7 +6,9 @@
\alias{getModel.PLNPCAfamily}
\alias{getModel}
\alias{getModel.PLNmixturefamily}
+\alias{getModel.Networkfamily}
\alias{getModel.PLNnetworkfamily}
+\alias{getModel.ZIPLNnetworkfamily}
\title{Model extraction from a collection of models}
\usage{
\method{getModel}{PLNPCAfamily}(Robject, var, index = NULL)
@@ -15,7 +17,11 @@ getModel(Robject, var, index)
\method{getModel}{PLNmixturefamily}(Robject, var, index = NULL)
+\method{getModel}{Networkfamily}(Robject, var, index = NULL)
+
\method{getModel}{PLNnetworkfamily}(Robject, var, index = NULL)
+
+\method{getModel}{ZIPLNnetworkfamily}(Robject, var, index = NULL)
}
\arguments{
\item{Robject}{an R6 object with class \code{\link{PLNPCAfamily}} or \code{\link{PLNnetworkfamily}}}
@@ -36,8 +42,12 @@ Model extraction from a collection of models
\item \code{getModel(PLNmixturefamily)}: Model extraction for \code{\link{PLNmixturefamily}}
+\item \code{getModel(Networkfamily)}: Model extraction for \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}}
+
\item \code{getModel(PLNnetworkfamily)}: Model extraction for \code{\link{PLNnetworkfamily}}
+\item \code{getModel(ZIPLNnetworkfamily)}: Model extraction for \code{\link{ZIPLNnetworkfamily}}
+
}}
\examples{
\dontrun{
diff --git a/man/plot.PLNnetworkfamily.Rd b/man/plot.Networkfamily.Rd
similarity index 56%
rename from man/plot.PLNnetworkfamily.Rd
rename to man/plot.Networkfamily.Rd
index 64cb0f34..ab24d8b7 100644
--- a/man/plot.PLNnetworkfamily.Rd
+++ b/man/plot.Networkfamily.Rd
@@ -1,9 +1,21 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/PLNnetworkfamily-S3methods.R
-\name{plot.PLNnetworkfamily}
+\name{plot.Networkfamily}
+\alias{plot.Networkfamily}
\alias{plot.PLNnetworkfamily}
-\title{Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of PLNnetwork fits (a \code{\link{PLNnetworkfamily}})}
+\alias{plot.ZIPLNnetworkfamily}
+\title{Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (either \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}})}
\usage{
+\method{plot}{Networkfamily}(
+ x,
+ type = c("criteria", "stability", "diagnostic"),
+ criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
+ reverse = FALSE,
+ log.x = TRUE,
+ stability = 0.9,
+ ...
+)
+
\method{plot}{PLNnetworkfamily}(
x,
type = c("criteria", "stability", "diagnostic"),
@@ -13,14 +25,25 @@
stability = 0.9,
...
)
+
+\method{plot}{ZIPLNnetworkfamily}(
+ x,
+ type = c("criteria", "stability", "diagnostic"),
+ criteria = c("loglik", "pen_loglik", "BIC", "EBIC"),
+ reverse = FALSE,
+ log.x = TRUE,
+ stability = 0.9,
+ ...
+)
}
\arguments{
-\item{x}{an R6 object with class \code{\link{PLNnetworkfamily}}}
+\item{x}{an R6 object with class \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}}}
\item{type}{a character, either "criteria", "stability" or "diagnostic" for the type of plot.}
-\item{criteria}{vector of characters. The criteria to plot in c("loglik", "BIC", "ICL", "R_squared", "EBIC", "pen_loglik").
-Default is c("loglik", "pen_loglik", "BIC", "EBIC"). Only relevant when \code{type = "criteria"}.}
+\item{criteria}{Vector of criteria to plot, to be selected among "loglik" (log-likelihood),
+"BIC", "ICL", "R_squared", "EBIC" and "pen_loglik" (penalized log-likelihood).
+Default is c("loglik", "pen_loglik", "BIC", "EBIC"). Only used when \code{type = "criteria"}.}
\item{reverse}{A logical indicating whether to plot the value of the criteria in the "natural" direction
(loglik - 0.5 penalty) or in the "reverse" direction (-2 loglik + penalty). Default to FALSE, i.e use the
@@ -38,12 +61,19 @@ Produces either a diagnostic plot (with \code{type = 'diagnostic'}), a stability
(with \code{type = 'criteria'}, the default).
}
\description{
-Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of PLNnetwork fits (a \code{\link{PLNnetworkfamily}})
+Display various outputs (goodness-of-fit criteria, robustness, diagnostic) associated with a collection of network fits (either \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}})
}
\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 \code{reverse = TRUE}.
}
+\section{Functions}{
+\itemize{
+\item \code{plot(PLNnetworkfamily)}: Display various outputs associated with a collection of network fits
+
+\item \code{plot(ZIPLNnetworkfamily)}: Display various outputs associated with a collection of network fits
+
+}}
\examples{
data(trichoptera)
trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
diff --git a/man/plot.ZIPLNfit_sparse.Rd b/man/plot.ZIPLNfit_sparse.Rd
new file mode 100644
index 00000000..d358717c
--- /dev/null
+++ b/man/plot.ZIPLNfit_sparse.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ZIPLNfit-S3methods.R
+\name{plot.ZIPLNfit_sparse}
+\alias{plot.ZIPLNfit_sparse}
+\title{Extract and plot the network (partial correlation, support or inverse covariance) from a \code{\link{ZIPLNfit_sparse}} object}
+\usage{
+\method{plot}{ZIPLNfit_sparse}(
+ 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{
+\item{x}{an R6 object with class \code{\link{ZIPLNfit_sparse}}}
+
+\item{type}{character. Value of the weight of the edges in the network, either "partial_cor" (partial correlation) or "support" (binary). Default is \code{"partial_cor"}.}
+
+\item{output}{the type of output used: either 'igraph' or 'corrplot'. Default is \code{'igraph'}.}
+
+\item{edge.color}{Length 2 color vector. Color for positive/negative edges. Default is \code{c("#F8766D", "#00BFC4")}. Only relevant for igraph output.}
+
+\item{remove.isolated}{if \code{TRUE}, isolated node are remove before plotting. Only relevant for igraph output.}
+
+\item{node.labels}{vector of character. The labels of the nodes. The default will use the column names ot the response matrix.}
+
+\item{layout}{an optional igraph layout. Only relevant for igraph output.}
+
+\item{plot}{logical. Should the final network be displayed or only sent back to the user. Default is \code{TRUE}.}
+
+\item{...}{Not used (S3 compatibility).}
+}
+\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)
+}
+\description{
+Extract and plot the network (partial correlation, support or inverse covariance) from a \code{\link{ZIPLNfit_sparse}} object
+}
+\examples{
+data(trichoptera)
+trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
+fit <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(penalty = 0.1))
+\dontrun{
+plot(fit)
+}
+}
diff --git a/man/stability_selection.Rd b/man/stability_selection.Rd
index 2fb5384a..ae378f80 100644
--- a/man/stability_selection.Rd
+++ b/man/stability_selection.Rd
@@ -12,16 +12,16 @@ stability_selection(
)
}
\arguments{
-\item{Robject}{an object with class \code{\link{PLNnetworkfamily}}, i.e. an output from \code{\link[=PLNnetwork]{PLNnetwork()}}}
+\item{Robject}{an object with class \code{\link{PLNnetworkfamily}} or \code{\link{ZIPLNnetworkfamily}}, i.e. an output from \code{\link[=PLNnetwork]{PLNnetwork()}} or \code{\link[=ZIPLNnetwork]{ZIPLNnetwork()}}}
\item{subsamples}{a list of vectors describing the subsamples. The number of vectors (or list length) determines th number of subsamples used in the stability selection. Automatically set to 20 subsamples with size \code{10*sqrt(n)} if \code{n >= 144} and \code{0.8*n} otherwise following Liu et al. (2010) recommendations.}
-\item{control}{a list controlling the main optimization process in each call to PLNnetwork. See \code{\link[=PLNnetwork]{PLNnetwork()}} for details.}
+\item{control}{a list controlling the main optimization process in each call to \code{\link[=PLNnetwork]{PLNnetwork()}} or \code{\link[=ZIPLNnetwork]{ZIPLNnetwork()}}. See \code{\link[=PLN_param]{PLN_param()}} or \code{\link[=ZIPLN_param]{ZIPLN_param()}} for details.}
\item{force}{force computation of the stability path, even if a previous one has been detected.}
}
\value{
-the list of subsamples. The estimated probabilities of selection of the edges are stored in the fields \code{stability_path} of the initial Robject with class \code{\link{PLNnetworkfamily}}
+the list of subsamples. The estimated probabilities of selection of the edges are stored in the fields \code{stability_path} of the initial Robject with class \code{\link{Networkfamily}}
}
\description{
This function computes the StARS stability criteria over a path of penalties. If a path has already been computed, the functions stops with a message unless \code{force = TRUE} has been specified.
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 62bb3121..01f8105d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -250,9 +250,9 @@ BEGIN_RCPP
return rcpp_result_gen;
END_RCPP
}
-// optim_zipln_R
-arma::mat optim_zipln_R(const arma::mat& Y, const arma::mat& X, const arma::mat& O, const arma::mat& M, const arma::mat& S, const arma::mat& Pi);
-RcppExport SEXP _PLNmodels_optim_zipln_R(SEXP YSEXP, SEXP XSEXP, SEXP OSEXP, SEXP MSEXP, SEXP SSEXP, SEXP PiSEXP) {
+// optim_zipln_R_var
+arma::mat optim_zipln_R_var(const arma::mat& Y, const arma::mat& X, const arma::mat& O, const arma::mat& M, const arma::mat& S, const arma::mat& Pi, const arma::mat& B);
+RcppExport SEXP _PLNmodels_optim_zipln_R_var(SEXP YSEXP, SEXP XSEXP, SEXP OSEXP, SEXP MSEXP, SEXP SSEXP, SEXP PiSEXP, SEXP BSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
@@ -262,7 +262,25 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const arma::mat& >::type M(MSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type S(SSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type Pi(PiSEXP);
- rcpp_result_gen = Rcpp::wrap(optim_zipln_R(Y, X, O, M, S, Pi));
+ Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+ rcpp_result_gen = Rcpp::wrap(optim_zipln_R_var(Y, X, O, M, S, Pi, B));
+ return rcpp_result_gen;
+END_RCPP
+}
+// optim_zipln_R_exact
+arma::mat optim_zipln_R_exact(const arma::mat& Y, const arma::mat& X, const arma::mat& O, const arma::mat& M, const arma::mat& S, const arma::mat& Pi, const arma::mat& B);
+RcppExport SEXP _PLNmodels_optim_zipln_R_exact(SEXP YSEXP, SEXP XSEXP, SEXP OSEXP, SEXP MSEXP, SEXP SSEXP, SEXP PiSEXP, SEXP BSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject rcpp_result_gen;
+ Rcpp::RNGScope rcpp_rngScope_gen;
+ Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type O(OSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type M(MSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type S(SSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type Pi(PiSEXP);
+ Rcpp::traits::input_parameter< const arma::mat& >::type B(BSEXP);
+ rcpp_result_gen = Rcpp::wrap(optim_zipln_R_exact(Y, X, O, M, S, Pi, B));
return rcpp_result_gen;
END_RCPP
}
@@ -331,7 +349,8 @@ static const R_CallMethodDef CallEntries[] = {
{"_PLNmodels_optim_zipln_Omega_diagonal", (DL_FUNC) &_PLNmodels_optim_zipln_Omega_diagonal, 4},
{"_PLNmodels_optim_zipln_B_dense", (DL_FUNC) &_PLNmodels_optim_zipln_B_dense, 2},
{"_PLNmodels_optim_zipln_zipar_covar", (DL_FUNC) &_PLNmodels_optim_zipln_zipar_covar, 4},
- {"_PLNmodels_optim_zipln_R", (DL_FUNC) &_PLNmodels_optim_zipln_R, 6},
+ {"_PLNmodels_optim_zipln_R_var", (DL_FUNC) &_PLNmodels_optim_zipln_R_var, 7},
+ {"_PLNmodels_optim_zipln_R_exact", (DL_FUNC) &_PLNmodels_optim_zipln_R_exact, 7},
{"_PLNmodels_optim_zipln_M", (DL_FUNC) &_PLNmodels_optim_zipln_M, 9},
{"_PLNmodels_optim_zipln_S", (DL_FUNC) &_PLNmodels_optim_zipln_S, 7},
{"_PLNmodels_cpp_test_packing", (DL_FUNC) &_PLNmodels_cpp_test_packing, 0},
diff --git a/src/lambertW.cpp b/src/lambertW.cpp
new file mode 100644
index 00000000..2d97dd3c
--- /dev/null
+++ b/src/lambertW.cpp
@@ -0,0 +1,138 @@
+/* lambertW.cpp
+
+Copyright (C) 2015, Avraham Adler
+All rights reserved.
+
+SPDX-License-Identifier: BSD-2-Clause
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+* Redistributions of source code must retain the above copyright notice, this
+list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright notice,
+this list of conditions and the following disclaimer in the documentation
+and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+References:
+
+Corless, R. M.; Gonnet, G. H.; Hare, D. E.; Jeffrey, D. J. & Knuth, D. E.
+ "On the Lambert W function", Advances in Computational Mathematics,
+ Springer, 1996, 5, 329-359
+
+Fritsch, F. N.; Shafer, R. E. & Crowley, W. P.
+ "Solution of the transcendental equation (we^w = x)",
+ Communications of the ACM, Association for Computing Machinery (ACM),
+ 1973, 16, 123-124
+*/
+
+#include "lambertW.h"
+
+ /* Fritsch Iteration
+ * W_{n+1} = W_n * (1 + e_n)
+ * z_n = ln(x / W_n) - W_n
+ * q_n = 2 * (1 + W_n) * (1 + W_n + 2 / 3 * z_n)
+ * e_n = z_n / (1 + W_n) * (q_n - z_n) / (q_n - 2 * z_n)
+ */
+
+double FritschIter(double x, double w){
+ int MaxEval = 5;
+ bool CONVERGED = false;
+ double k = 2.0 / 3.0;
+ int i = 0;
+ do {
+ double z = std::log(x / w) - w;
+ double w1 = w + 1.0;
+ double q = 2.0 * w1 * (w1 + k * z);
+ double qmz = q - z;
+ double e = z / w1 * qmz / (qmz - z);
+ CONVERGED = std::abs(e) <= EPS;
+ w *= (1.0 + e);
+ ++i;
+ } while (!CONVERGED && i < MaxEval);
+ return(w);
+}
+
+double lambertW0_CS(double x) {
+ if (x == R_PosInf) {
+ return(R_PosInf);
+ } else if (x < -M_1_E) {
+ return(R_NaN);
+ } else if (std::abs(x + M_1_E) <= EPS) {
+ return(-1.0);
+ } else if (std::abs(x) <= 1e-16) {
+ /* This close to 0 the W_0 branch is best estimated by its Taylor/Pade
+ expansion whose first term is the value x and remaining terms are below
+ machine double precision. See
+ https://math.stackexchange.com/questions/1700919
+ */
+ return(x);
+ } else {
+ double w;
+ if (std::abs(x) <= 6.4e-3) {
+ /* When this close to 0 the Fritsch iteration may underflow. Instead,
+ * function will use degree-6 minimax polynomial approximation of Halley
+ * iteration-based values. Should be more accurate by three orders of
+ * magnitude than Fritsch's equation (5) in this range.
+ */
+
+ // Minimax Approximation calculated using R package minimaxApprox 0.1.0
+ return((((((-1.0805085529250425e1 * x + 5.2100070265741278) * x -
+ 2.6666665063383532) * x + 1.4999999657268301) * x -
+ 1.0000000000016802) * x + 1.0000000000001752) * x +
+ 2.6020852139652106e-18);
+
+ } else if (x <= M_E) {
+ /* Use expansion in Corliss 4.22 to create (2, 2) Pade approximant.
+ * Equation with a few extra terms is:
+ * -1 + p - 1/3p^2 + 11/72p^3 - 43/540p^4 + 689453/8398080p^4 - O(p^5)
+ * This is just used to estimate a good starting point for the Fritsch
+ * iteration process itself.
+ */
+ double p = std::sqrt(2.0 * (M_E * x + 1.0));
+ double Numer = (0.2787037037037037 * p + 0.311111111111111) * p - 1.0;
+ double Denom = (0.0768518518518518 * p + 0.688888888888889) * p + 1.0;
+ w = Numer / Denom;
+ } else {
+ /* Use first five terms of Corliss et al. 4.19 */
+ w = std::log(x);
+ double L_2 = std::log(w);
+ double L_3 = L_2 / w;
+ double L_3_sq = L_3 * L_3;
+ w += -L_2 + L_3 + 0.5 * L_3_sq - L_3 / w + L_3 / (w * w) - 1.5 * L_3_sq /
+ w + L_3_sq * L_3 / 3.0;
+ }
+ return(FritschIter(x, w));
+ }
+}
+
+double lambertWm1_CS(double x){
+ if (x == 0.0) {
+ return(R_NegInf);
+ } else if (x < -M_1_E || x > 0.0) {
+ return(R_NaN);
+ } else if (std::abs(x + M_1_E) <= EPS) {
+ return(-1.0);
+ } else {
+ double w;
+ /* Use first five terms of Corliss et al. 4.19 */
+ w = std::log(-x);
+ double L_2 = std::log(-w);
+ double L_3 = L_2 / w;
+ double L_3_sq = L_3 * L_3;
+ w += -L_2 + L_3 + 0.5 * L_3_sq - L_3 / w + L_3 / (w * w) - 1.5 * L_3_sq /
+ w + L_3_sq * L_3 / 3.0;
+ return(FritschIter(x, w));
+ }
+}
+
diff --git a/src/lambertW.h b/src/lambertW.h
new file mode 100644
index 00000000..6b56a484
--- /dev/null
+++ b/src/lambertW.h
@@ -0,0 +1,15 @@
+#include
+
+#define _USE_MATH_DEFINES
+#include
+
+using namespace Rcpp;
+
+const double EPS = 2.2204460492503131e-16;
+const double M_1_E = 1.0 / M_E;
+
+double FritschIter(double x, double w) ;
+double lambertW0_CS(double x) ;
+double lambertWm1_CS(double x) ;
+
+
diff --git a/src/optim_zi-pln.cpp b/src/optim_zi-pln.cpp
index 5665ed47..559979cb 100644
--- a/src/optim_zi-pln.cpp
+++ b/src/optim_zi-pln.cpp
@@ -6,6 +6,7 @@
#include "nlopt_wrapper.h"
#include "packing.h"
#include "utils.h"
+#include "lambertW.h"
// [[Rcpp::export]]
arma::vec zipln_vloglik(
@@ -118,13 +119,14 @@ Rcpp::List optim_zipln_zipar_covar(
}
// [[Rcpp::export]]
-arma::mat optim_zipln_R(
+arma::mat optim_zipln_R_var(
const arma::mat & Y, // responses (n,p)
const arma::mat & X, // covariates (n,d)
const arma::mat & O, // offsets (n,p)
const arma::mat & M, // (n,p)
const arma::mat & S, // (n,p)
- const arma::mat & Pi // (d,p)
+ const arma::mat & Pi, // (d,p)
+ const arma::mat & B // covariates (n,d)
) {
arma::mat A = exp(O + M + 0.5 * S % S);
arma::mat R = pow(1. + exp(- (A + logit(Pi))), -1);
@@ -144,6 +146,39 @@ arma::mat optim_zipln_R(
return R;
}
+double phi (double mu, double sigma2) {
+ double W = lambertW0_CS(sigma2 * exp(mu)) ;
+ return(exp(-(pow(W, 2) + 2 * W) / (2 * sigma2)) / sqrt(1 + W)) ;
+}
+
+// [[Rcpp::export]]
+arma::mat optim_zipln_R_exact (
+ const arma::mat & Y, // covariates (n,d)
+ const arma::mat & X, // covariates (n,d)
+ const arma::mat & O, // offsets (n,p)
+ const arma::mat & M, // (n,p)
+ const arma::mat & S, // (n,p)
+ const arma::mat & Pi, // (n,p)
+ const arma::mat & B // covariates (n,d)
+) {
+
+ arma::mat XB = X * B;
+ arma::mat M_mu = M - XB;
+ arma::uword n = M.n_rows;
+ arma::uword p = M.n_cols;
+ arma::vec diag_Sigma = arma::diagvec((1./n) * (M_mu.t() * M_mu + diagmat(sum(S % S, 0)))) ;
+ arma::mat R = arma::zeros(n,p);
+ for(arma::uword i = 0; i < n; i += 1) {
+ for(arma::uword j = 0; j < p; j += 1) {
+ if(Y(i, j) < 0.5) {
+ double Phi = phi(O(i,j) + XB(i,j), diag_Sigma(j)) ;
+ R(i,j) = Pi(i,j) / (Phi * (1 - Pi(i,j)) + Pi(i,j)) ;
+ }
+ }
+ }
+ return R;
+}
+
// [[Rcpp::export]]
Rcpp::List optim_zipln_M(
const arma::mat & init_M, // (n,p)
diff --git a/tests/testthat/test-plnnetworkfamily.R b/tests/testthat/test-plnnetworkfamily.R
index fe736451..97ffd071 100644
--- a/tests/testthat/test-plnnetworkfamily.R
+++ b/tests/testthat/test-plnnetworkfamily.R
@@ -10,19 +10,23 @@ test_that("PLNnetwork: main function, fields access and methods", {
expect_equal(getBestModel(models), getBestModel(models, "BIC"))
- X <- model.matrix(Abundance ~ 1, data = trichoptera)
Y <- as.matrix(trichoptera$Abundance)
- O <- matrix(0, nrow(Y),ncol(Y))
- w <- rep(1, nrow(Y))
+ data = list(
+ Y = Y,
+ X = model.matrix(Abundance ~ 1, data = trichoptera),
+ O = matrix(0, nrow(Y),ncol(Y)),
+ w = rep(1, nrow(Y)),
+ formula = Abundance ~ 1
+ )
## extract the data matrices and weights
ctrl <- PLNmodels:::PLNnetwork_param(trace = 0)
## instantiate
- myPLN <- PLNmodels:::PLNnetworkfamily$new(NULL, Y, X, O, w, Abundance ~ 1, ctrl)
+ myPLN <- PLNmodels:::PLNnetworkfamily$new(NULL, data, ctrl)
## optimize
- myPLN$optimize(ctrl$config_optim)
+ myPLN$optimize(data, ctrl$config_optim)
## post-treatment
config_post <- PLNmodels:::config_post_default_PLNnetwork
@@ -32,7 +36,7 @@ test_that("PLNnetwork: main function, fields access and methods", {
expect_equal(myPLN$criteria$BIC, models$criteria$BIC)
## S3 methods
- expect_true(PLNmodels:::isPLNnetworkfamily(myPLN))
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
@@ -90,7 +94,7 @@ test_that("PLNnetwork: matrix of penalties work", {
myPLN <- PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W))
## S3 methods
- expect_true(PLNmodels:::isPLNnetworkfamily(myPLN))
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
@@ -141,7 +145,7 @@ test_that("PLNnetwork: list of matrices of penalties work", {
myPLN <- PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = list_W))
## S3 methods
- expect_true(PLNmodels:::isPLNnetworkfamily(myPLN))
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
diff --git a/tests/testthat/test-plnnetworkfit.R b/tests/testthat/test-plnnetworkfit.R
index 4b311fa3..a835dc96 100644
--- a/tests/testthat/test-plnnetworkfit.R
+++ b/tests/testthat/test-plnnetworkfit.R
@@ -42,7 +42,7 @@ test_that("PLNnetwork fit: check classes, getters and field access", {
expect_equal(vcov(myPLNfit, "covariance"), myPLNfit$model_par$Sigma)
expect_equal(vcov(myPLNfit, "covariance"), sigma(myPLNfit))
expect_warning(standard_error(myPLNfit))
- expect_true(igraph::is.igraph(myPLNfit$plot_network(output = "igraph", plot = FALSE)))
+ expect_true(igraph::is_igraph(myPLNfit$plot_network(output = "igraph", plot = FALSE)))
expect_true(inherits(myPLNfit$plot_network(output = "corrplot", plot = FALSE), "Matrix"))
})
diff --git a/tests/testthat/test-zipln.R b/tests/testthat/test-zipln.R
index d6361f9e..7eac24b0 100644
--- a/tests/testthat/test-zipln.R
+++ b/tests/testthat/test-zipln.R
@@ -60,6 +60,17 @@ test_that("PLN is working with unnamed data matrix", {
expect_error(ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(config_optim = list(algorithm = "nawak"))))
})
+test_that("ZIPLN is working with exact and variational inference for the conditional distribution of the ZI component", {
+
+ approx <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(config_optim = list(approx_ZI = TRUE)))
+ exact <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(config_optim = list(approx_ZI = FALSE)))
+
+ expect_equal(approx$loglik, exact$loglik, tolerance = 1e-1) ## Almost equivalent
+ expect_equal(approx$model_par$B, exact$model_par$B, tolerance = 1e-1) ## Almost equivalent
+ expect_equal(approx$model_par$Sigma, exact$model_par$Sigma, tolerance = 1e-1) ## Almost equivalent
+
+})
+
test_that("ZIPLN: Check that univariate ZIPLN models works, with matrix of numeric format", {
expect_no_error(uniZIPLN <- ZIPLN(Abundance[,1,drop=FALSE] ~ 1, data = trichoptera))
expect_no_error(uniZIPLN <- ZIPLN(Abundance[,1] ~ 1, data = trichoptera))
diff --git a/tests/testthat/test-ziplnfit.R b/tests/testthat/test-ziplnfit.R
index 9ba0ffc5..dbc5bfd5 100644
--- a/tests/testthat/test-ziplnfit.R
+++ b/tests/testthat/test-ziplnfit.R
@@ -118,3 +118,16 @@ test_that("ZIPLN fit: Check number of parameters", {
expect_equal(model$vcov_model, "fixed")
})
+
+test_that("ZIPLN fit: check sparse output and plot", {
+
+ myPLNfit <- ZIPLN(Abundance ~ 1, data = trichoptera, control = ZIPLN_param(trace = 1, penalty = 0.1))
+
+ expect_is(myPLNfit, "ZIPLNfit_sparse")
+ expect_is(myPLNfit, "ZIPLNfit")
+
+ expect_equal(myPLNfit$vcov_model, "sparse")
+ expect_true(igraph::is.igraph(myPLNfit$plot_network(output = "igraph", plot = FALSE)))
+ expect_true(inherits(myPLNfit$plot_network(output = "corrplot", plot = FALSE), "Matrix"))
+
+})
diff --git a/tests/testthat/test-ziplnnetworkfamily.R b/tests/testthat/test-ziplnnetworkfamily.R
new file mode 100644
index 00000000..b5223dd3
--- /dev/null
+++ b/tests/testthat/test-ziplnnetworkfamily.R
@@ -0,0 +1,170 @@
+context("test-plnnetworkfamily")
+
+data(trichoptera)
+## use a subset t save some time
+trichoptera <- prepare_data(trichoptera$Abundance[1:20, 1:5], trichoptera$Covariate[1:20, ])
+
+models <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera)
+
+test_that("ZIPLNnetwork: main function, fields access and methods", {
+
+ expect_equal(getBestModel(models), getBestModel(models, "BIC"))
+
+ Y <- as.matrix(trichoptera$Abundance)
+ data = list(
+ Y = Y,
+ X = model.matrix(Abundance ~ 1, data = trichoptera),
+ X0 = matrix(NA,0,0),
+ O = matrix(0, nrow(Y),ncol(Y)),
+ w = rep(1, nrow(Y)),
+ formula = Abundance ~ 1
+ )
+
+ ## extract the data matrices and weights
+ ctrl <- PLNmodels:::ZIPLNnetwork_param(trace = 0)
+ ctrl$ziparam <- "single"
+
+ ## instantiate
+ myPLN <- PLNmodels:::ZIPLNnetworkfamily$new(NULL, data, ctrl)
+
+ ## optimize
+ myPLN$optimize(data, ctrl$config_optim)
+
+ expect_equal(myPLN$criteria$BIC, models$criteria$BIC)
+
+ ## S3 methods
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
+ expect_is(plot(myPLN), "ggplot")
+ expect_is(plot(myPLN, reverse = TRUE), "ggplot")
+ expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
+ expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
+ expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
+
+ ## Field access
+ expect_true(all(myPLN$penalties > 0))
+ expect_null(myPLN$stability_path)
+ expect_true(anyNA(myPLN$stability))
+
+ ## Other R6 methods
+ expect_true(is.data.frame(myPLN$coefficient_path()))
+ subs <- replicate(2,
+ sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
+ simplify = FALSE)
+ myPLN$stability_selection(subsamples = subs)
+ expect_is(plot(myPLN, type = "stability"), "ggplot")
+ expect_true(!is.null(myPLN$stability_path))
+ expect_true(inherits(myPLN$plot(), "ggplot"))
+ expect_true(inherits(myPLN$plot_objective(), "ggplot"))
+ expect_true(inherits(myPLN$plot_stars(), "ggplot"))
+})
+
+test_that("ZIPLNnetwork computes the stability path only once.", {
+
+ ## extract_probs fails if stability selection has not been performed.
+ expect_error(extract_probs(models),
+ "Please perform stability selection using stability_selection(Robject) first", fixed = TRUE)
+ set.seed(1077)
+ subs <- replicate(2,
+ sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
+ simplify = FALSE)
+ stability_selection(models, subsamples = subs, force = TRUE)
+ models$stability_selection(subsamples = subs)
+ ## Stability_path has correct dimensions
+ p <- getModel(models, index = 1)$p
+ expect_equal(dim(models$stability_path),
+ c(length(models$penalties) * p*(p-1)/2L, 5))
+ ## try to compute it again
+ expect_message(stability_selection(models),
+ "Previous stability selection detected. Use \"force = TRUE\" to recompute it.")
+ ## extracts the inclusion frequencies
+ expect_equal(dim(extract_probs(models, index = 1, format = "matrix")),
+ c(p, p))
+ expect_length(extract_probs(models, index = 1, format = "vector"),
+ p*(p-1)/2)
+})
+
+test_that("ZIPLNnetwork: matrix of penalties work", {
+
+ p <- ncol(trichoptera$Abundance)
+ W <- diag(1, p, p)
+ W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
+ W[lower.tri(W)] <- t(W)[lower.tri(W)]
+ myPLN <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, control = ZIPLNnetwork_param(penalty_weights = W))
+
+ ## S3 methods
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
+ expect_is(plot(myPLN), "ggplot")
+ expect_is(plot(myPLN, reverse = TRUE), "ggplot")
+ expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
+ expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
+ expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
+
+ ## Field access
+ expect_true(all(myPLN$penalties > 0))
+ expect_null(myPLN$stability_path)
+ expect_true(anyNA(myPLN$stability))
+
+ ## Other R6 methods
+ expect_true(is.data.frame(myPLN$coefficient_path()))
+ subs <- replicate(2,
+ sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
+ simplify = FALSE)
+ myPLN$stability_selection(subsamples = subs, control = ZIPLNnetwork_param(penalty_weights = W))
+ expect_is(plot(myPLN, type = "stability"), "ggplot")
+ expect_true(!is.null(myPLN$stability_path))
+ expect_true(inherits(myPLN$plot(), "ggplot"))
+ expect_true(inherits(myPLN$plot_objective(), "ggplot"))
+ expect_true(inherits(myPLN$plot_stars(), "ggplot"))
+
+ ## misspecification of penalty weights should induce errors
+ ## not symmetric
+ W <- diag(1, p, p)
+ W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
+ expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
+
+ ## not square
+ W <- matrix(1, p + 1, p)
+ expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
+
+ ## nonpositive entries
+ W <- matrix(0, p, p)
+ expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
+
+})
+
+test_that("PLNnetwork: list of matrices of penalties work", {
+
+ p <- ncol(trichoptera$Abundance)
+ W <- diag(1, p, p)
+ W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
+ W[lower.tri(W)] <- t(W)[lower.tri(W)]
+ list_W <- lapply(seq(1, 1e-2, len = 30), function(rho) rho * W)
+
+ myPLN <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, control = ZIPLNnetwork_param(penalty_weights = list_W))
+
+ ## S3 methods
+ expect_true(PLNmodels:::isNetworkfamily(myPLN))
+ expect_is(plot(myPLN), "ggplot")
+ expect_is(plot(myPLN, reverse = TRUE), "ggplot")
+ expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
+ expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
+ expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
+
+ ## Field access
+ expect_true(all(myPLN$penalties > 0))
+ expect_null(myPLN$stability_path)
+ expect_true(anyNA(myPLN$stability))
+
+ ## Other R6 methods
+ expect_true(is.data.frame(myPLN$coefficient_path()))
+ subs <- replicate(2,
+ sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
+ simplify = FALSE)
+ myPLN$stability_selection(subsamples = subs, control = ZIPLNnetwork_param(penalty_weights = W))
+ expect_is(plot(myPLN, type = "stability"), "ggplot")
+ expect_true(!is.null(myPLN$stability_path))
+ expect_true(inherits(myPLN$plot(), "ggplot"))
+ expect_true(inherits(myPLN$plot_objective(), "ggplot"))
+ expect_true(inherits(myPLN$plot_stars(), "ggplot"))
+
+})