#' Poisson lognormal model towards Principal Component Analysis
#'
#' Fit the PCA variants of the Poisson lognormal with a variational algorithm. Use the (g)lm syntax for model specification (covariates, offsets).
#'
#' @param formula an object of class "formula": a symbolic description of the model to be fitted.
#' @param 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.
#' @param subset an optional vector specifying a subset of observations to be used in the fitting process.
#' @param weights an optional vector of observation weights to be used in the fitting process.
#' @param ranks a vector of integer containing the successive ranks (or number of axes to be considered)
#' @param control a list-like structure for controlling the optimization, with default generated by [PLNPCA_param()]. See the associated documentation.
#' for details.
#'
#' @return an R6 object with class [`PLNPCAfamily`], which contains
#' a collection of models with class [`PLNPCAfit`]
#'
#' @rdname PLNPCA
#' @examples
#' #' ## Use future to dispatch the computations on 2 workers
#' \dontrun{
#' future::plan("multisession", workers = 2)
#' }
#'
#' data(trichoptera)
#' trichoptera <- prepare_data(trichoptera$Abundance, trichoptera$Covariate)
#' myPCA <- PLNPCA(Abundance ~ 1 + offset(log(Offset)), data = trichoptera, ranks = 1:5)
#'
#' # Shut down parallel workers
#' \dontrun{
#' future::plan("sequential")
#' }
#' @seealso The classes [`PLNPCAfamily`] and [`PLNPCAfit`], and the configuration function [PLNPCA_param()].
#' @importFrom stats model.frame model.matrix model.response model.offset
#' @export
PLNPCA <- function(formula, data, subset, weights, ranks = 1:5, control = PLNPCA_param()) {
## Temporary test for deprecated use of list()
if (!inherits(control, "PLNmodels_param"))
stop("We now use the function PLNPCA_param() to generate the list of parameters that controls the fit:
replace 'list(my_arg = xx)' by PLN_param(my_arg = xx) and see the documentation of PLNPCA_param().")
## extract the data matrices and weights
args <- extract_model(match.call(expand.dots = FALSE), parent.frame())
## Instantiate the collection of PLN models, initialized by PLN with full rank
if (control$trace > 0) cat("\n Initialization...")
myPCA <- PLNPCAfamily$new(ranks, args$Y, args$X, args$O, args$w, args$formula, control)
## Adjust the PLN models
if (control$trace > 0) cat("\n\n Adjusting", length(ranks), "PLN models for PCA analysis.\n")
myPCA$optimize(control$config_optim)
## Post-treatments: pseudo-R2, rearrange criteria and prepare PCA visualization
if (control$trace > 0) cat("\n Post-treatments")
myPCA$postTreatment(control$config_post, control$config_optim)
if (control$trace > 0) cat("\n DONE!\n")
myPCA
}
#' Control of PLNPCA fit
#'
#' Helper to define list of parameters to control the PLNPCA fit. All arguments have defaults.
#'
#' @param backend optimization back used, either "nlopt" or "torch". Default is "nlopt"
#' @param trace a integer for verbosity.
#' @param config_optim a list for controlling the optimizer (either "nlopt" or "torch" backend). See details
#' @param config_post a list for controlling the post-treatments (optional bootstrap, jackknife, R2, etc.). See details
#' @param 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.
#'
#' @return list of parameters configuring the fit.
#'
#' @inherit PLN_param details
#' @export
PLNPCA_param <- function(
backend = "nlopt",
trace = 1 ,
config_optim = list() ,
config_post = list() ,
inception = NULL # pretrained PLNfit used as initialization
) {
if (!is.null(inception)) stopifnot(isPLNfit(inception))
## post-treatment config
config_pst <- config_post_default_PLNPCA
config_pst[names(config_post)] <- config_post
config_pst$trace <- trace
## optimization config
backend <- match.arg(backend)
stopifnot(backend %in% c("nlopt", "torch"))
if (backend == "nlopt") {
stopifnot(config_optim$algorithm %in% available_algorithms_nlopt)
config_opt <- config_default_nlopt
}
if (backend == "torch") {
stopifnot(config_optim$algorithm %in% available_algorithms_torch)
config_opt <- config_default_torch
}
config_opt[names(config_optim)] <- config_optim
config_opt$trace <- trace
structure(list(
backend = backend ,
trace = trace ,
config_optim = config_opt,
config_post = config_pst,
inception = inception ), class = "PLNmodels_param")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.