R/Diagnostics.R
In Keefe-Murphy/IMIFA: Infinite Mixtures of Infinite Factor Analysers and Related Models
Defines functions
get_IMIFA_results.IMIFA
get_IMIFA_results
Documented in
get_IMIFA_results
#' Extract results, conduct posterior inference and compute performance metrics for MCMC samples of models from the IMIFA family
#'
#' This function post-processes simulations generated by \code{\link{mcmc_IMIFA}} for any of the IMIFA family of models. This includes accounting for label switching, and accounting for rotational invariance via Procrustean methods. It can be re-ran at little computational cost in order to extract different models explored by the sampler used for \code{sims}, without having to re-run the model itself. New results objects using different numbers of clusters and different numbers of factors (if visited by the model in question), or using different model selection criteria (if necessary) can be generated with ease. Posterior predictive checking of the appropriateness of the fitted model is also facilitated.
#' @param sims An object of class \code{"IMIFA"} generated by \code{\link{mcmc_IMIFA}}.
#' @param burnin Optional additional number of iterations to discard. Defaults to 0, corresponding to no additional burnin. See \code{\link{mixfaControl}} for the default \code{burnin} settings used previously by \code{\link{mcmc_IMIFA}}.
#' @param thinning Optional interval for extra thinning to be applied. Defaults to 1, corresponding to no additional thinning. See \code{\link{mixfaControl}} for the default \code{thinning} settings used previously by \code{\link{mcmc_IMIFA}}.
#' @param G If this argument is not specified, results will be returned with the optimal number of clusters. If different numbers of clusters were explored in \code{sims} for the \code{"MFA"} or \code{"MIFA"} methods, supplying an integer value allows pulling out a specific solution with \code{G} clusters, even if the solution is sub-optimal.
#'
#' Similarly, this allows retrieval of samples corresponding to a solution, if visited, with \code{G} clusters for the \code{"OMFA"}, \code{"OMIFA"}, \code{"IMFA"} and \code{"IMIFA"} methods.
#' @param Q If this argument is not specified, results will be returned with the optimal number of factors. If different numbers of factors were explored in \code{sims} for the \code{"FA"}, \code{"MFA"}, \code{"OMFA"} or \code{"IMFA"} methods, this allows pulling out a specific solution with \code{Q} factors, even if the solution is sub-optimal.
#'
#' Similarly, this allows retrieval of samples corresponding to a solution, if visited, with \code{Q} factors for the \code{"IFA"}, \code{"MIFA"}, \code{"OMIFA"} and \code{"IMIFA"} methods if adaptation has already taken place. Can be supplied as a scalar or a vector of values for each cluster.
#'
#' If adaptation didn't take place during model-fitting, \code{adapt=TRUE} can be supplied here to determine the optimal cluster-specific numbers of non-redundant factors (see below).
#' @param criterion The criterion to use for model selection, where model selection is only required if more than one model was run under the \code{"FA"}, \code{"MFA"}, \code{"MIFA"}, \code{"OMFA"} or \code{"IMFA"} methods when \code{sims} was created via \code{\link{mcmc_IMIFA}}. Defaults to \code{bicm}, but note that these are \emph{all} calculated; this argument merely indicates which one will form the basis of the construction of the output.
#'
#' Note that the first three options here might exhibit bias in favour of zero-factor models for the finite factor \code{"FA"}, \code{"MFA"}, \code{"OMFA"} and \code{"IMFA"} methods and might exhibit bias in favour of one-cluster models for the \code{"MFA"} and \code{"MIFA"} methods. The \code{aic.mcmc} and \code{bic.mcmc} criteria will only be returned for finite factor models.
#' @param adapt A logical indicating if adaptation should be applied to the stored loadings and scores matrices to truncate the cluster-specific number(s) of non-redundant factors. This argument is only relevant if \code{adapt=FALSE} was used when initially fitting the model (otherwise, adaptation has already taken place during sampling, by default), and hence defaults to \code{FALSE} here. Relevant parameters from \code{\link{mgpControl}} (namely \code{active.crit}, \code{eps}, \code{prop}, and \code{forceQg}) can be passed via the \code{...} construct, but will default to their values under \code{\link{mgpControl}} if not specified. Note that \code{adapt=TRUE} is only invoked if the relevant parameters were stored via \code{\link{storeControl}} when running the model: loadings only for \code{active.crit="BD"}, loadings and scores for \code{active.crit="SC"} for \code{"IFA"} models.
#' @param G.meth If the object in \code{sims} arises from the \code{"OMFA"}, \code{"OMIFA"}, \code{"IMFA"} or \code{"IMIFA"} methods, this argument determines whether the optimal number of clusters is given by the mode or median of the posterior distribution of \code{G}. Defaults to \code{"mode"}. Often the mode and median will agree in any case.
#' @param Q.meth If the object in \code{sims} arises from the \code{"IFA"}, \code{"MIFA"}, \code{"OMIFA"} or \code{"IMIFA"} methods, this argument determines whether the optimal number of latent factors is given by the mode or median of the posterior distribution of \code{Q}. Defaults to \code{"mode"}. Often the mode and median will agree in any case.
#' @param conf.level The confidence level to be used throughout for credible intervals for all parameters of inferential interest, and error metrics if \code{error.metrics=TRUE}. Defaults to \code{0.95}.
#' @param error.metrics A logical activating or deactivating posterior predictive checking: i.e. controlling whether metrics quantifying a) the posterior predictive reconstruction error (PPRE) between bin counts of the data and bin counts of replicate draws from the posterior distribution & and b) the error between the empirical and estimated covariance matrices should be computed. These are computed for every \emph{valid} retained iteration (see \code{Details}). Defaults to \code{TRUE}, but can be time-consuming for models which achieve clustering. These error metrics, and the uncertainty associated with them, can be visualised via \code{\link{plot.Results_IMIFA}}. Depending on what parameters were stored when calling \code{\link{mcmc_IMIFA}}, potentially not all error metrics will be available to compute.
#'
#' The Frobenius norm is used in the computation of the PPRE, by default, but the \code{type} of \code{\link{norm}} can be changed via the \code{...} construct below. So too can the breakpoints (\code{dbreaks}) used to bin the data and the posterior predictive replicate data sets. Some caution is advised in the latter case.
#' @param vari.rot Logical indicating whether the loadings matrix/matrices template(s) should be \code{\link[stats]{varimax}} rotated first, prior to the Procrustes rotation steps. Defaults to \code{FALSE}. Not necessary at all for clustering purposes, or inference on the covariance matrix, but useful if interpretable inferences on the loadings matrix/matrices are desired. Arguments to \code{\link[stats]{varimax}} can be passed via the \code{...} construct, but note that the argument \code{normalize} here defaults to \code{FALSE}.
#' @param z.avgsim Logical (defaults to \code{FALSE}) indicating whether the clustering should also be summarised with a call to \code{\link{Zsimilarity}} by the clustering with minimum mean squared error to the similarity matrix obtained by averaging the stored adjacency matrices, in addition to the MAP estimate.
#'
#' Note that the MAP clustering is computed \emph{conditional} on the estimate of the number of clusters (whether that be the modal estimate or the estimate according to \code{criterion}) and other parameters are extracted conditional on this estimate of \code{G}: however, in contrast, the number of distinct clusters in the summarised labels obtained by specifying \code{z.avgsim=TRUE} may not necessarily coincide with the MAP estimate of \code{G}, but it may provide a useful alternative summary of the partitions explored during the chain, and the user is free to call \code{\link{get_IMIFA_results}} again with the new suggested \code{G} value.
#'
#' Please be warned that this feature requires loading the \code{\link[mcclust]{mcclust}} package. This is liable to take considerable time to compute, and may not even be possible if the number of observations &/or number of stored iterations is large and the resulting matrix isn't sufficiently sparse. When \code{z.avgsim=TRUE}, both the summarised clustering and the similarity matrix are stored: the latter can be visualised as part of a call to \code{\link{plot.Results_IMIFA}}.
#' @param zlabels For any method that performs clustering, the true labels can be supplied if they are known in order to compute clustering performance metrics. This also has the effect of ordering the MAP labels (and thus the ordering of cluster-specific parameters) to most closely correspond to the true labels if supplied.
#' @param nonempty For \code{"MFA"} and \code{"MIFA"} models ONLY: a logical indicating whether only iterations with non-empty components should be retained. Defaults to \code{TRUE}, but may lead to empty chains - conversely, \code{FALSE} may lead to empty components and related errors.
#' @param x,object,MAP,... Arguments required for the \code{print.Results_IMIFA} and \code{summary.Results_IMIFA} functions: \code{x} and \code{object} are objects of class \code{"Results_IMIFA"} resulting from a call to \code{\link{get_IMIFA_results}}. \code{MAP} is a logical which governs whether a table of the MAP classification is printed, while \code{...} gathers additional arguments to those functions.
#'
#' Users can also pass the \code{type} argument to the \code{\link{norm}} function when \code{isTRUE(error.metrics)} and the posterior predictive reconstruction error (PPRE) is calculated. By default the Frobenius norm (\code{type="F"}) is employed.
#'
#' When \code{adapt=TRUE} here, relevant arguments from \code{\link{mgpControl}} (namely \code{active.crit}, \code{eps}, \code{prop}, and \code{forceQg}) can be supplied here too, though they retain their default values from \code{\link{mgpControl}} otherwise.
#'
#' Finally, the \code{...} construct also allows arguments to \code{\link[stats]{varimax}} to be passed to \code{\link{get_IMIFA_results}} itself, when \code{isTRUE(vari.rot)}, or arguments to \code{\link[graphics]{hist}} when \code{isTRUE(error.metrics)}, in order to guide construction of the bins. Additionally, by passing the argument \code{dbreaks} via the \code{...} construct, the bins can be specified directly. However, caution is advised in doing so; in particular, the bins must be constructed on data which has been standardised in the same way as the data modelled within \code{\link{mcmc_IMIFA}}.
#'
#' @details The function also performs post-hoc corrections for label switching, as well as post-hoc Procrustes rotation of loadings matrices and scores, in order to ensure sensible posterior parameter estimates, computes error metrics, constructs credible intervals, and generally transforms the raw \code{sims} object into an object of class \code{"Results_IMIFA"} in order to prepare the results for plotting via \code{\link{plot.Results_IMIFA}}.
#'
#' For the infinite factor methods, iterations where the maximum number of factors was greater than or equal to the maximum of the estimated cluster-specific factors are retained for posterior summaries of the scores, in order to preserve the estimated dimension of the scores matrices. Similarly, these are also the \emph{valid} iterations used for the computation of the averages and credible intervals for the error metrics. For the finite factor models, \emph{all} retained iterations are used in both instances (i.e. both for the scores and the error metrics).
#'
#' In all cases, only iterations with \code{G} non-empty components are retained.
#' @return An object of class \code{"Results_IMIFA"} to be passed to \code{\link{plot.Results_IMIFA}} for visualising results. Dedicated \code{print} and \code{summary} functions also exist for objects of this class. The object, say \code{x}, is a list of lists, the most important components of which are:
#' \item{\code{Clust}}{Everything pertaining to clustering performance can be found here for all but the \code{"FA"} and \code{"IFA"} methods (or models where the estimate number of clusters is 1), in particular \code{x$Clust$MAP}, the MAP summary of the posterior clustering, the last valid sample of cluster labels \code{x$Clust$last.z}, the matrix of posterior cluster membership probabilities \code{x$Clust$post.prob}, and the posterior confusion matrix \code{x$Clust$PCM}.
#'
#' More detail is given if known \code{zlabels} are supplied: performance is always evaluated against the MAP clustering, with additional evaluation against the alternative clustering computed if \code{z.avgsim=TRUE}. Posterior summaries of the mixing proportions, and the concentration/discount parameters, if necessary, are also included here, as well as the last valid samples of each.}
#' \item{\code{Error}}{Everything pertaining the model fit assessment can be found here, incl. the distribution of the PPRE values and associated bin counts for the replicate draws, as well as average error metrics (e.g. MSE, RMSE), and credible intervals quantifying the associated uncertainty, between the empirical and estimated covariance matrix/matrices, both of which are also included.}
#' \item{\code{GQ.results}}{Everything pertaining to model choice can be found here, incl. posterior summaries for the estimated number of clusters and estimated number of factors, if applicable to the method employed. Model selection criterion values are also accessible here.}
#' \item{\code{Means}}{Posterior summaries for the means, after conditioning on \code{G}.}
#' \item{\code{Loadings}}{Posterior summaries for the factor loadings matrix/matrices, after conditioning on \code{G} and \code{Q}. Posterior mean loadings given by \code{x$Loadings$post.load} are given the \code{\link[stats]{loadings}} class for printing purposes and thus the manner in which they are displayed can be modified.
#'
#' The number of iterations retained for posterior summaries of the loadings may vary for different clusters for the infinite factor methods, corresponding to iterations where the cluster-specific number of factors was greater than or equal to the modal estimate of the cluster-specific number of factors.}
#' \item{\code{Scores}}{Posterior summaries for the latent factor scores, after conditioning on the maximum of the estimated number of cluster-specific factors. Summaries are given for the \emph{single} matrix of factor scores. See \code{\link{scores_MAP}} to decompose these summaries into sub-matrices according to the MAP partition (for models which achieve clustering).
#'
#' For the infinite factor methods, iterations where the maximum number of factors was greater than or equal to the maximum of the estimated cluster-specific factors are retained for posterior summaries of the scores, in order to preserve the estimated dimension of the scores matrices.}
#' \item{\code{Uniquenesses}}{Posterior summaries for the uniquenesses, after conditioning on \code{G}.}
#'
#' The objects \code{Means}, \code{Loadings}, \code{Scores} and \code{Uniquenesses} (if stored when calling \code{\link{mcmc_IMIFA}}!) also contain, as well as the posterior summaries, the entire chain of valid samples of each, as well as, for convenience, the last valid samples of each (after conditioning on the modal \code{G} and \code{Q} values, and accounting for label switching, and rotational invariance via Procrustes rotation).
#' @note For the \code{"IMIFA"}, \code{"IMFA"}, \code{"OMIFA"}, and \code{"OMFA"} methods, the retained mixing proportions are renormalised after conditioning on the modal \code{G}. This is especially necessary for the computation of the \code{error.metrics}, just note that the values on which posterior inference are conducted will ever so slightly differ from the actually sampled values.
#'
#' Due to the way the offline label-switching correction is performed, different runs of this function may give \emph{very slightly} different results in terms of the cluster labellings (and by extension the parameters, which are permuted in the same way), but only if the chain was run for an extremely small number of iterations, well below the number required for convergence, and samples of the cluster labels match poorly across iterations (particularly if the number of clusters suggested by those sampled labels is high).
#' @keywords IMIFA main
#' @include MainFunction.R
#' @export
#' @importFrom Rfast "colMaxs" "colTabulate" "Median" "rowMaxs" "rowOrder" "rowSort" "rowTabulate" "rowVars" "Var"
#' @importFrom mclust "classError"
#' @importFrom matrixStats "colSums2" "colMeans2" "rowAlls" "rowMeans2" "rowMedians" "rowQuantiles" "rowSums2"
#' @importFrom slam "as.simple_triplet_matrix"
#'
#' @seealso \code{\link{plot.Results_IMIFA}}, \code{\link{mcmc_IMIFA}}, \code{\link{Zsimilarity}}, \code{\link{scores_MAP}}, \code{\link{sim_IMIFA_model}}, \code{\link{Procrustes}}, \code{\link[stats]{varimax}}, \code{\link{norm}}, \code{\link{mgpControl}}, \code{\link{storeControl}}
#' @references Murphy, K., Viroli, C., and Gormley, I. C. (2020) Infinite mixtures of infinite factor analysers, \emph{Bayesian Analysis}, 15(3): 937-963. <\doi{10.1214/19-BA1179}>.
#'
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @usage
#' get_IMIFA_results(sims = NULL,
#' burnin = 0L,
#' thinning = 1L,
#' G = NULL,
#' Q = NULL,
#' criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
#' adapt = FALSE,
#' G.meth = c("mode", "median"),
#' Q.meth = c("mode", "median"),
#' conf.level = 0.95,
#' error.metrics = TRUE,
#' vari.rot = FALSE,
#' z.avgsim = FALSE,
#' zlabels = NULL,
#' nonempty = TRUE,
#' ...)
#' @examples
#' \donttest{# data(coffee)
#' # data(olive)
#'
#' # Run a MFA model on the coffee data over a range of clusters and factors.
#' # simMFAcoffee <- mcmc_IMIFA(coffee, method="MFA", range.G=2:3, range.Q=0:3, n.iters=1000)
#'
#' # Accept all defaults to extract the optimal model.
#' # resMFAcoffee <- get_IMIFA_results(simMFAcoffee)
#'
#' # Instead let's get results for a 3-cluster model, allowing Q be chosen by aic.mcmc.
#' # resMFAcoffee2 <- get_IMIFA_results(simMFAcoffee, G=3, criterion="aic.mcmc")
#'
#' # Run an IMIFA model on the olive data, accepting all defaults.
#' # simIMIFAolive <- mcmc_IMIFA(olive, method="IMIFA", n.iters=10000)
#'
#' # Extract optimum results
#' # Estimate G & Q by the median of their posterior distributions
#' # Construct 90% credible intervals and try to return the similarity matrix.
#' # resIMIFAolive <- get_IMIFA_results(simIMIFAolive, G.meth="median", Q.meth="median",
#' # conf.level=0.9, z.avgsim=TRUE)
#' # summary(resIMIFAolive)
#'
#' # Simulate new data from the above model
#' # newdata <- sim_IMIFA_model(resIMIFAolive)
#'
#' # Fit an IFA model without adaptation and examine results with and without post-hoc adaptation
#' # Use the "SC" criterion for determining the number of active factors
#' # simIFAnoadapt <- mcmc_IMIFA(olive, method="IFA", n.iters=10000, adapt=FALSE)
#' # resIFAnoadapt <- get_IMIFA_results(simIFAnoadapt)
#' # resIFApostadapt <- get_IMIFA_results(simIFAnoadapt, adapt=TRUE, active.crit="SC")
#'
#' # Compare to an IFA model with adaptive Gibbs sampling
#' # simIFAadapt <- mcmc_IMIFA(coffee, method="IFA", n.iters=10000, active.crit="BD")
#' # resIFAadapt <- get_IMIFA_results(simIFAadapt)
#' # plot(resIFAnoadapt, "GQ")
#' # plot(resIFApostadapt, "GQ")
#' # plot(resIFAadapt, "GQ")}
get_IMIFA_results <- function(sims = NULL, burnin = 0L, thinning = 1L, G = NULL, Q = NULL, criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
adapt = FALSE, G.meth = c("mode", "median"), Q.meth = c("mode", "median"), conf.level = 0.95, error.metrics = TRUE,
vari.rot = FALSE, z.avgsim = FALSE, zlabels = NULL, nonempty = TRUE, ...) {
UseMethod("get_IMIFA_results")
}
#' @method get_IMIFA_results IMIFA
#' @export
get_IMIFA_results.IMIFA <- function(sims = NULL, burnin = 0L, thinning = 1L, G = NULL, Q = NULL, criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
adapt = FALSE, G.meth = c("mode", "median"), Q.meth = c("mode", "median"), conf.level = 0.95, error.metrics = TRUE,
vari.rot = FALSE, z.avgsim = FALSE, zlabels = NULL, nonempty = TRUE, ...) {
call <- match.call()
defopt <- options()
options(warn=1)
on.exit(suppressWarnings(options(defopt)), add=TRUE)
if(missing(sims)) stop("Simulations must be supplied", call.=FALSE)
if(!inherits(sims, "IMIFA")) stop("Object of class 'IMIFA' must be supplied", call.=FALSE)
burnin <- as.integer(burnin)
thinning <- as.integer(thinning)
store <- attr(sims, "Store")
if(any(c(length(thinning),
length(burnin)) > 1)) stop("'burnin' and 'thinning' must be of length 1", call.=FALSE)
if(any(burnin < 0, burnin >= store,
thinning < 1)) stop("Invalid 'burnin' and/or 'thinning' supplied", call.=FALSE)
store <- as.integer(seq(from=burnin + 1L, to=store, by=thinning))
if(length(store) < 10) stop("Too much 'burnin' or 'thinning' applied: not enough stored samples to proceed", call.=FALSE)
n.store <- length(store)
tmp.store <- store
dat <- attr(sims, "Dataset")
data.name <- attr(sims, "Name")
equal.pro <- attr(sims, "Equal.Pi")
label.switch <- attr(sims, "Label.Switch")
method <- attr(sims, "Method")
learn.alpha <- attr(sims, "Alph.step")
learn.d <- attr(sims, "Disc.step")
inf.G <- is.element(method, c("IMIFA", "OMIFA", "IMFA", "OMFA"))
inf.Q <- is.element(method, c("IMIFA", "OMIFA", "MIFA", "IFA"))
n.fac <- attr(sims, "Factors")
n.grp <- attr(sims, "Clusters")
n.obs <- attr(sims, "Obs")
n.var <- attr(sims, "Vars")
uni <- n.var == 1
sw <- tmpsw <- attr(sims, "Switch")
uni.meth <- attr(sims, "Uni.Meth")
uni.type <- unname(uni.meth["Uni.Type"])
conf.level <- as.numeric(conf.level)
obsnames <- attr(sims, "Obsnames")
varnames <- attr(sims, "Varnames")
adapted <- attr(sims, "Adapt")
adapt <- adapt && isFALSE(adapted)
cov.emp <- attr(sims, "Cov.Emp")
cov.range <- ifelse(uni, 1L, max(cov.emp) - min(cov.emp))
cshrink <- attr(sims, "C.Shrink")
if(any(length(conf.level) != 1,
!is.numeric(conf.level),
(conf.level <= 0 ||
conf.level >= 1))) stop("'conf.level' must be a single number in the interval (0, 1)", call.=FALSE)
conf.levels <- c((1 - conf.level)/2, (1 + conf.level)/2)
choice <- length(n.grp) * length(n.fac) > 1
crit.check <- !all(is.character(criterion))
if(isTRUE(crit.check)) stop("'criterion' must be a character vector of length 1", call.=FALSE)
criterion <- match.arg(criterion)
if(all(inf.Q, is.element(criterion,
c("aic.mcmc", "bic.mcmc")))) stop(paste0(ifelse(isTRUE(choice), "Model choice is", "Though model choice isn't"), " actually required -\n'criterion' cannot be 'aic.mcmc' or 'bic.mcmc' for the ", method, " method"), call.=FALSE)
if(length(error.metrics) != 1 ||
!is.logical(error.metrics)) stop("'error.metrics' must be a single logical indicator", call.=FALSE)
if(length(vari.rot) != 1 ||
!is.logical(vari.rot)) stop("'vari.rot' must be a single logical indicator", call.=FALSE)
if(isTRUE(error.metrics) &&
anyNA(cov.emp)) { warning("Certain error metrics cannot be computed as there are missing values in the empirical covariance matrix\n", call.=FALSE, immediate.=TRUE)
cov.met <- FALSE
} else cov.met <- error.metrics
if(length(z.avgsim) != 1 ||
!is.logical(z.avgsim)) stop("'z.avgsim' must be a single logical indicator", call.=FALSE)
if((z.avgsim <- (isTRUE(z.avgsim) &&
!is.element(method, c("FA", "IFA"))))) {
if(!(has_pkg <- suppressMessages(requireNamespace("mcclust", quietly=TRUE)))) {
z.avgsim <- FALSE; warning("Forcing 'z.avgsim' to FALSE: 'mcclust' package not installed\n", call.=FALSE, immediate.=TRUE)
} else if(n.obs >= 10000) message("Consider setting 'z.avgsim' to FALSE as the number of observations is quite large\n")
}
dots <- list(...)
dots <- dots[unique(names(dots))]
if(isTRUE(adapt)) {
active.crit <- ifelse(is.null(dots$active.crit), "BD", dots$active.crit)
epsilon <- ifelse(is.null(dots$eps), 1e-01, dots$eps)
prop <- ifelse(is.null(dots$prop), ifelse(n.var <= 2, 0.5,
switch(EXPR=active.crit,
BD=0.7, SC=0.99)), dots$prop)
prop <- switch(EXPR=active.crit, BD=floor(prop * n.var)/n.var, SC=prop)
forceQg <- ifelse(is.null(dots$forceQg), FALSE, dots$forceQg)
if(method != "IFA" &&
"SC" == active.crit) stop("'active.crit'=\"SC\" is only available for the \"IFA\" method", call.=FALSE)
if(!sw["l.sw"] &&
"BD" == active.crit) stop("Loadings must be stored for adaptation to take place", call.=FALSE)
if(any(!sw[c("l.sw",
"s.sw")]) &&
"SC" == active.crit) stop("Scores and loadings must be stored for adaption using the \"SC\" critertion to take place", call.=FALSE)
}
G.T <- !missing(G)
Q.T <- !missing(Q)
G.ind <- Q.ind <- 1L
GQs <- length(sims[[G.ind]])
GQ1x <- GQs > 1 && is.element(method, c("OMFA", "IMFA"))
if(inf.G) {
G.store2 <- lapply(seq_len(GQs), function(gq) sims[[G.ind]][[gq]]$G.store[store])
G.store <- matrix(unlist(G.store2), nrow=GQs, ncol=n.store, byrow=TRUE)
if(is.element(method, c("IMFA", "IMIFA"))) {
act.store <- lapply(seq_len(GQs), function(gq) sims[[G.ind]][[gq]]$act.store[store])
}
G.mX <- !all(is.character(G.meth))
if(isTRUE(G.mX)) stop("'G.meth' must be a character vector of length 1", call.=FALSE)
G.meth <- match.arg(G.meth)
G.tab <- if(GQ1x) lapply(apply(G.store, 1L, function(x) list(table(x, dnn=NULL))), "[[", 1L) else table(G.store, dnn=NULL)
G.prob <- if(GQ1x) lapply(G.tab, prop.table) else prop.table(G.tab)
G.mode <- if(GQ1x) unlist(lapply(G.tab, function(gt) as.numeric(names(gt[gt == max(gt)])[1L]))) else as.numeric(names(G.tab[G.tab == max(G.tab)])[1L])
G.med <- if(GQ1x) ceiling(matrixStats::rowMedians(G.store,
useNames=FALSE) * 2)/2 else ceiling(Median(G.store) * 2)/2
if(!G.T) {
G <- switch(EXPR=G.meth, mode=G.mode, floor(G.med))
}
G.CI <- if(GQ1x) round(rowQuantiles(G.store, probs=conf.levels)) else round(stats::quantile(G.store, conf.levels))
}
if(G.T) {
G <- as.integer(G)
if(any(length(G) != 1,
!is.integer(G))) stop("'G' must be an integer of length 1", call.=FALSE)
if(!inf.G) {
if(!is.element(method, c("FA", "IFA"))) {
if(!is.element(G, n.grp)) stop("This 'G' value was not used during simulation", call.=FALSE)
G.ind <- which(n.grp == G)
} else if(G > 1) message(paste0("Forced G=1 for the ", method, " method\n"))
} else {
if(GQ1x) {
if(!Q.T) stop(paste0("'G' cannot be supplied without 'Q' for the ", method, " method if a range of Q values were explored"), call.=FALSE)
tmpQ <- which(n.fac == unique(Q))
} else {
tmpQ <- Q.ind
}
if(length(tmpQ > 0) && !is.element(G,
unique(G.store[tmpQ,]))) stop("This 'G' value was not visited during simulation", call.=FALSE)
}
}
G <- if(any(inf.G, all(G.T, !is.element(method, c("FA", "IFA"))))) G else 1L
if(Q.T) {
Q <- as.integer(Q)
if(!is.integer(Q)) stop("'Q' must of integer type", call.=FALSE)
if(G.T) {
if(length(Q) == 1) Q <- rep(Q, G)
if(length(Q) != G) stop(paste0("'Q' must be supplied for each cluster, as a scalar or vector of length G=", G), call.=FALSE)
} else if(length(n.grp) != 1 && all(!is.element(length(Q),
c(1, n.grp)))) stop("'Q' must be a scalar if G=1, 'G' is not suppplied, or a range of G values were explored", call.=FALSE)
if(all(is.element(method, c("FA", "MFA", "OMFA", "IMFA")))) {
if(length(unique(Q)) != 1) stop(paste0("'Q' cannot vary across clusters for the ", method, " method"), call.=FALSE)
Q <- unique(Q)
if(!is.element(Q, n.fac)) stop("This 'Q' value was not used during simulation", call.=FALSE)
Q.ind <- which(n.fac == Q)
}
if(inf.Q) {
if(any((Q != 0) + (Q *
(n.var - Q) <= 0) > 1)) stop(paste0("'Q' must be less than the number of variables ", n.var), call.=FALSE)
Qtmp <- if(inf.G) sims[[1L]][[1L]]$Q.store[seq_len(G),store, drop=FALSE] else switch(EXPR=method, MIFA=sims[[if(G.T) G == n.grp else G.ind]][[1L]]$Q.store[,store, drop=FALSE], sims[[1L]][[1L]]$Q.store[,store, drop=FALSE])
storage.mode(Qtmp) <- switch(EXPR=method, IFA="integer", "numeric")
Qtmp <- switch(EXPR=method, IFA=max(Qtmp), Rfast::rowMaxs(Qtmp, value=TRUE))
Qlen <- length(Qtmp)
if(length(Q) == 1) Q <- rep(Q, Qlen)
if(length(Q) != Qlen) stop(paste0("'Q' must be supplied for each cluster, as a scalar or vector of length G=", G), call.=FALSE)
if(any(Q * (Qtmp - Q) < 0)) stop(paste0("'Q' can't be greater than the maximum number of factors stored in ", ifelse(method == "IFA", "", "any cluster of "), deparse(match.call()$sims)), call.=FALSE)
}
}
if(inf.G) {
tmp.store <- if(GQ1x) lapply(seq_len(GQs), function(gq) store[G.store[gq,] == G[ifelse(G.T, 1L, gq)]]) else store[G.store == G]
GQ.temp1 <- list(G = G, G.Mode = G.mode, G.Median = G.med, G.CI = G.CI, G.Probs = G.prob, G.Counts = G.tab)
GQ.temp1 <- c(GQ.temp1, list(Stored.G = switch(EXPR=method, OMIFA=provideDimnames(G.store, base=list("Non-Empty", ""), unique=FALSE),
IMIFA=provideDimnames(do.call(rbind, c(G.store2, act.store)), base=list(c("Non-Empty", "Active"), ""), unique=FALSE),
OMFA=lapply(seq_len(GQs), function(g) provideDimnames(t(G.store[g,]), base=list("Non-Empty", ""), unique=FALSE)),
IMFA=lapply(seq_len(GQs), function(g) provideDimnames(rbind(G.store2[[g]], act.store[[g]]), base=list(c("Non-Empty", "Active"), ""), unique=FALSE)))))
GQ.temp1 <- c(GQ.temp1, list(G.Last = switch(EXPR=method, OMIFA=, IMIFA=GQ.temp1$Stored.G[1L,n.store],
lapply(seq_len(GQs), function(g) GQ.temp1$Stored.G[[g]][1L,n.store]))))
}
G.range <- ifelse(G.T, 1L, length(n.grp))
Q.range <- ifelse(any(Q.T, all(!is.element(method, c("OMFA", "IMFA")), inf.Q)), 1L, length(n.fac))
crit.mat <- matrix(NA, nrow=G.range, ncol=Q.range)
# Retrieve log-likelihoods and/or tune G &/or Q according to criterion
if(all(G.T, Q.T)) {
dimnames(crit.mat) <- list(paste0("G", G), if(inf.Q) "IFA" else paste0("Q", Q))
} else if(G.T) {
dimnames(crit.mat) <- list(paste0("G", G), if(inf.Q) "IFA" else paste0("Q", n.fac))
} else if(Q.T) {
dimnames(crit.mat) <- list(paste0("G", n.grp), if(inf.Q) "IFA" else paste0("Q", Q))
} else {
dimnames(crit.mat) <- list(paste0("G", n.grp), if(inf.Q) "IFA" else paste0("Q", n.fac))
}
rownames(crit.mat) <- switch(EXPR=method, IMFA=, IMIFA="IM", OMFA=, OMIFA="OM", rownames(crit.mat))
aicm <- bicm <- aicm.sd <- dic <-
aic.mcmc <- bic.mcmc <- bicm.sd <- crit.mat
log.N <- log(n.obs)
if(!inf.Q) {
Ks <- matrix(NA, nrow=G.range, ncol=Q.range)
}
for(g in seq_len(G.range)) {
gi <- ifelse(G.T, G.ind, g)
for(q in seq_len(Q.range)) {
qi <- ifelse(Q.T, Q.ind, q)
log.likes <- if(GQ1x) sims[[gi]][[qi]]$ll.store[tmp.store[[qi]]] else sims[[gi]][[qi]]$ll.store[tmp.store]
log.likes <- log.likes[stats::complete.cases(log.likes)]
S2 <- ifelse(length(log.likes) != 1, Var(log.likes), 0L)
llbar <- mean(log.likes)
d.hat <- 2 * S2
llmax2 <- 2 * llbar + d.hat
aicm[g,q] <- 4 * llbar - llmax2
bicm[g,q] <- llmax2 - d.hat * log.N
dic[g,q] <- 2 * (3 * llbar - llmax2)
ci.tmp <- S2 * (11 * S2 + 24)
aicm.sd[g,q] <- sqrt((4 * (S2 + ci.tmp))/length(log.likes))
bicm.sd[g,q] <- sqrt((2 * d.hat + (log.N - 1)^2 * ci.tmp)/length(log.likes))
if(!inf.Q) {
Ks[g,q] <- K <- switch(EXPR=method, OMFA=, IMFA=PGMM_dfree(Q=n.fac[qi], P=n.var, G=G[ifelse(G.T, 1L, qi)],
method=switch(EXPR=uni.type, unconstrained="UUU", isotropic="UUC", constrained="UCU", single="UCC"),
equal.pro=equal.pro), attr(sims[[gi]][[qi]], "K"))
llmax2k <- 2 * max(log.likes)
aic.mcmc[g,q] <- llmax2k - K * 2
bic.mcmc[g,q] <- llmax2k - K * log.N
}
}
}
crit <- get(criterion)
crit.max <- which(crit == max(crit, na.rm=TRUE), arr.ind=TRUE)
if(!inf.Q &&
nrow(crit.max) > 1) { warning(paste0("Ties for the optimal model exist according to the '", criterion, "' criterion: choosing the most parsimonious model\n"), call.=FALSE, immediate.=TRUE)
crit.max <- which(Ks == min(Ks[crit.max], na.rm=TRUE), arr.ind=TRUE)
crit.max <- crit.max[which.min(crit.max[,1L]),]
}
# Control for supplied values of G &/or Q
if(!any(Q.T, G.T)) {
G.ind <- crit.max[1L]
Q.ind <- crit.max[2L]
if(!inf.G) {
G <- n.grp[G.ind]
}
if(!inf.Q) {
Q <- n.fac[Q.ind]
}
} else if(all(G.T, !Q.T)) {
Q.ind <- which.max(crit)
if(!inf.Q) {
Q <- n.fac[Q.ind]
}
} else if(all(Q.T, !G.T)) {
G.ind <- which.max(crit)
if(!inf.G) {
G <- n.grp[G.ind]
}
}
G <- ifelse(inf.G, ifelse(G.T, G, G[Q.ind]), ifelse(length(n.grp) == 1, n.grp, G))
Gseq <- seq_len(G)
gnames <- paste0("Cluster", Gseq)
G.ind <- if(!all(inf.G, length(n.grp) > 1)) G.ind else which(n.grp == G)
GQ.temp2 <- list(AICMs = aicm, BICMs = bicm, DICs = dic)
clust.ind <- !any(is.element(method, c("FA", "IFA")),
all(is.element(method, c("MFA", "MIFA")), G == 1))
sw.mx <- (!clust.ind || sw["mu.sw"]) && sw["u.sw"]
sw.px <- !clust.ind || sw["psi.sw"]
if(!inf.Q) {
Q <- if(length(n.fac) > 1) Q else n.fac
Q.ind <- if(all(!Q.T, length(n.fac) > 1)) Q.ind else which(n.fac == Q)
Q <- stats::setNames(if(length(Q) != G) rep(Q, G) else Q, gnames)
if(all(inf.G, Q.T)) GQ.temp1$G <- rep(G, GQs)
if(GQ1x) {
GQ.temp1$G.CI <- lapply(seq_len(GQs), function(gq) GQ.temp1$G.CI[gq,])
GQ.temp1 <- lapply(GQ.temp1, "[[", Q.ind)
} else if(inf.G) {
GQ.temp1$Stored.G <- GQ.temp1$Stored.G[[1L]]
GQ.temp1$G.Last <- GQ.temp1$G.Last[[1L]]
}
GQ.temp3 <- c(GQ.temp2, list(AIC.mcmcs = aic.mcmc, BIC.mcmcs = bic.mcmc))
GQ.res <- switch(EXPR=method, OMFA=, IMFA=c(GQ.temp1, list(Q = Q), list(Criteria = GQ.temp3)), c(list(G = G, Q = Q), list(Criteria = GQ.temp3)))
if(GQ1x) {
tmp.store <- tmp.store[[Q.ind]]
}
}
if(clust.ind) {
zadj <- sims[[G.ind]][[Q.ind]]$z.store
z <- as.matrix(zadj[tmp.store,])
}
if(condition <- all(isTRUE(nonempty), is.element(method, c("MFA", "MIFA")), G > 1)) {
emptystore <- which(rowAlls(rowTabulate(sims[[G.ind]][[Q.ind]]$z.store[tmp.store,], max_number=n.grp[G.ind]) > 0, useNames=FALSE))
if(!identical(tmp.store,
emptystore)) warning(paste0("Discarding iterations with fewer than G=", G, " non-empty components\n"), call.=FALSE, immediate.=TRUE)
tmp.store <- emptystore
z <- z[tmp.store,, drop=FALSE]
}
storeG <- seq_along(tmp.store)
TN.store <- length(tmp.store)
if(TN.store < 2) stop(paste0("Not enough samples stored to proceed", ifelse(condition && missing(nonempty), ": try supplying 'nonempty'=FALSE", ifelse(choice && any(G.T, Q.T), paste0(": try supplying different G or Q values"), ""))), call.=FALSE)
if(inf.Q) {
Q.store <- sims[[G.ind]][[Q.ind]]$Q.store[,tmp.store, drop=FALSE]
Q.mX <- !all(is.character(Q.meth))
if(isTRUE(Q.mX)) stop("'Q.meth' must be a character vector of length 1", call.=FALSE)
Q.meth <- match.arg(Q.meth)
if(cshrink <- cshrink && ((method != "MIFA") || G > 1)) {
sigmas <- sims[[G.ind]][[Q.ind]]$sigma[,tmp.store, drop=FALSE]
}
}
# Manage Label Switching & retrieve cluster labels/mixing proportions
if(clust.ind) {
label.miss <- missing(zlabels)
if(!label.miss && (all(G > 1, !is.factor(zlabels), !is.numeric(zlabels)) ||
length(zlabels) != n.obs)) stop(paste0("'zlabels' must be a factor of length N=", n.obs), call.=FALSE)
if(sw["mu.sw"]) {
mus <- sims[[G.ind]][[Q.ind]]$mu[,,tmp.store, drop=FALSE]
}
if(sw["l.sw"]) {
lmats <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$load)[,,,tmp.store, drop=FALSE] else sims[[G.ind]][[Q.ind]]$load[,,,tmp.store, drop=FALSE]
}
if(sw["psi.sw"]) {
psis <- sims[[G.ind]][[Q.ind]]$psi[,,tmp.store, drop=FALSE]
}
if(sw["pi.sw"]) {
pies <- sims[[G.ind]][[Q.ind]]$pi.prop[,tmp.store, drop=FALSE]
}
if(sw["s.sw"]) {
z2 <- z
}
zadj <- zadj[store,]
if(!label.switch) {
z.temp <- tryCatch(factor(z[1L,], labels=Gseq), error=function(e) factor(z[1L,], levels=Gseq))
for(sl in storeG) {
sw.lab <- .lab_switch(z.new=z[sl,], z.old=z.temp)
z.perm <- sw.lab$z.perm
left <- as.integer(unname(z.perm))
right <- as.integer(names(z.perm))
if(!identical(left, right)) {
z[sl,] <- sw.lab$z
if(sw["mu.sw"]) {
mus[,left,sl] <- mus[,right,sl]
}
if(sw["l.sw"]) {
lmats[,,left,sl] <- lmats[,,right,sl]
}
if(sw["psi.sw"]) {
psis[,left,sl] <- psis[,right,sl]
}
if(sw["pi.sw"]) {
pies[left,sl] <- pies[right,sl]
}
if(inf.Q) {
Q.store[left,sl] <- Q.store[right,sl]
if(cshrink) {
sigmas[left,sl] <- sigmas[right,sl]
}
}
}
}
}
if(inf.Q) {
Q.store <- provideDimnames(Q.store[Gseq,, drop=FALSE], base=list(gnames, ""), unique=FALSE)
}
if(sw["mu.sw"]) mus <- tryCatch(mus[,Gseq,, drop=FALSE], error=function(e) mus)
if(sw["l.sw"]) lmats <- tryCatch(lmats[,,Gseq,, drop=FALSE], error=function(e) lmats)
if(sw["psi.sw"]) psis <- tryCatch(psis[,Gseq,, drop=FALSE], error=function(e) psis)
if(cshrink) sigmas <- tryCatch(sigmas[Gseq,, drop=FALSE], error=function(e) sigmas)
MAP <- if(G > 1) apply(z, 2L, function(x) factor(which.max(tabulate(x, nbins=G)), levels=Gseq)) else factor(rep(1L, n.obs))
post.prob <- if(G > 1) matrix(colTabulate(z, max_number=G)/TN.store, nrow=n.obs, ncol=G, byrow=TRUE) else matrix(1L, nrow=n.obs)
if(isTRUE(z.avgsim)) {
znew <- try(Zsimilarity(zs=zadj), silent=TRUE)
condit <- all(!is.element(method, c("MIFA", "MFA")), inherits(znew, "try-error"))
if(isTRUE(condit)) {
znew <- try(Zsimilarity(zs=z), silent=TRUE)
warning("Constructing the similarity matrix failed:\ntrying again using iterations corresponding to the modal number of clusters\n", call.=FALSE, immediate.=TRUE)
}
if(!inherits(znew, "try-error")) {
zadj <- znew$z.avg
zadj <- factor(zadj, labels=seq_along(unique(zadj)))
zadj <- as.integer(levels(zadj))[zadj]
zavg <- znew$z.sim
if(!label.miss) {
zadj <- .lab_switch(z.new=zadj, z.old=factor(zlabels, labels=seq_along(unique(zlabels))))$z
tab <- table(zadj, zlabels, dnn=list("Predicted", "Observed"))
tabstat <- c(.class_agreement(tab), classError(classification=MAP, class=zlabels))
if(nrow(tab) != ncol(tab)) {
tabstat <- tabstat[-seq_len(2L)]
names(tabstat)[4L] <- "error.rate"
} else {
names(tabstat)[6L] <- "error.rate"
}
if(tabstat$error.rate == 0) {
tabstat$misclassified <- NULL
}
tabstat <- c(list(confusion.matrix = tab), tabstat)
class(tabstat) <- "listof"
}
z_simavg <- list(z.avg = zadj, z.sim = zavg)
z_simavg <- c(z_simavg, if(!label.miss) list(avgsim.perf = tabstat))
attr(z_simavg, "Conditional") <- condit
class(z_simavg) <- "listof"
} else { warning("Can't compute similarity matrix or 'average' clustering: forcing 'z.avgsim' to FALSE\n", call.=FALSE, immediate.=TRUE)
z.avgsim <- FALSE
}
}
uncertain <- if(G > 1) 1 - Rfast::rowMaxs(post.prob, value=TRUE) else integer(n.obs)
sizes <- stats::setNames(tabulate(MAP, nbins=G), gnames)
if(any(sizes == 0)) warning(paste0("Empty cluster exists in MAP partition:\nexamine trace plots", ifelse(any(is.element(method, c("OMFA", "IMFA", "OMIFA", "IMIFA")), is.element(method, c("MFA", "MIFA")) && any(n.grp < G)), ", try to supply a lower G value to get_IMIFA_results(),", ""), " or re-run the model\n"), call.=FALSE, immediate.=TRUE)
if(sw["pi.sw"]) {
pi.prop <- provideDimnames(pies[Gseq,storeG, drop=FALSE], base=list(gnames, ""), unique=FALSE)
pi.prop <- if(inf.G) sweep(pi.prop, 2L, colSums2(pi.prop, useNames=FALSE), FUN="/", check.margin=FALSE) else pi.prop
var.pi <- stats::setNames(rowVars(pi.prop), gnames)
ci.pi <- rowQuantiles(pi.prop, probs=conf.levels)
ci.pi <- if(G == 1) t(ci.pi) else ci.pi
post.pi <- stats::setNames(rowMeans2(pi.prop,
useNames=FALSE), gnames)
} else if(equal.pro) {
post.pi <- stats::setNames(rep(1/G, G), gnames)
} else { message("Mixing proportions not stored: estimating posterior mean by the proportions of the MAP clustering\n")
post.pi <- stats::setNames(sizes/n.obs, gnames)
}
tab.z <- if((!label.miss && G > 1) || (clust.ind && all(error.metrics, sw["mu.sw"], sw["psi.sw"]))) rowTabulate(z, max_number=G)
if(!label.miss && G > 1) {
sw.lab <- .lab_switch(z.new=MAP, z.old=factor(zlabels, labels=seq_along(unique(zlabels))))
MAP <- factor(factor(sw.lab$z, labels=which(sizes > 0)), levels=Gseq)
l.perm <- sw.lab$z.perm
left <- as.integer(l.perm[Gseq])
z.tab <- tab.z > 0
z.tmp <- lapply(storeG, function(i) factor(z[i,], labels=left[z.tab[i,]]))
z <- do.call(rbind, lapply(z.tmp, function(x) as.integer(levels(x))[as.integer(x)]))
index <- order(left)
tab.z <- tab.z[,index, drop=FALSE]
sizes <- stats::setNames(sizes[index], gnames)
post.prob <- post.prob[,index]
if(sw["mu.sw"]) mus <- mus[,index,, drop=FALSE]
if(sw["l.sw"]) lmats <- lmats[,,index,, drop=FALSE]
if(sw["psi.sw"]) psis <- psis[,index,, drop=FALSE]
post.pi <- stats::setNames(post.pi[index], gnames)
if(sw["pi.sw"]) {
pi.prop <- provideDimnames(unname(pi.prop[index,, drop=FALSE]), base=list(gnames, ""), unique=FALSE)
var.pi <- stats::setNames(var.pi[index], gnames)
ci.pi <- provideDimnames(unname(ci.pi[index,, drop=FALSE]), base=list(gnames, colnames(ci.pi)))
}
if(inf.Q) {
Q.store <- provideDimnames(unname(Q.store[index,, drop=FALSE]), base=list(gnames, ""), unique=FALSE)
if(cshrink) sigmas <- sigmas[index,, drop=FALSE]
}
}
if(!label.miss) {
tab <- table(MAP, zlabels, dnn=list("Predicted", "Observed"))
tab.stat <- c(.class_agreement(tab), classError(classification=MAP, class=zlabels))
if(nrow(tab) != ncol(tab)) {
tab.stat <- tab.stat[-seq_len(2L)]
names(tab.stat)[4L] <- "error.rate"
} else {
names(tab.stat)[6L] <- "error.rate"
}
if(tab.stat$error.rate == 0) {
tab.stat$misclassified <- NULL
}
tab.stat <- c(list(confusion.matrix = tab), tab.stat)
class(tab.stat) <- "listof"
}
if(learn.alpha) {
alpha <- sims[[G.ind]][[Q.ind]]$alpha[store]
post.alpha <- mean(alpha)
var.alpha <- Var(alpha)
ci.alpha <- stats::quantile(alpha, conf.levels)
rate <- sims[[G.ind]][[Q.ind]]$a.rate
DP.alpha <- list(alpha = alpha, post.alpha = post.alpha, var.alpha = var.alpha, ci.alpha = ci.alpha, alpha.rate = rate, last.alpha = alpha[n.store])
DP.alpha <- c(DP.alpha, if(isTRUE(attr(sims, "TuneZeta"))) list(avg.zeta = sims[[G.ind]][[Q.ind]]$avg.zeta))
class(DP.alpha) <- "listof"
}
if(learn.d) {
discount <- drop(sims[[G.ind]][[Q.ind]]$discount[store])
post.disc <- mean(discount)
post.kappa <- sum(discount == 0)/n.store
var.disc <- Var(discount)
ci.disc <- stats::quantile(discount, conf.levels)
rate <- sims[[G.ind]][[Q.ind]]$d.rate
post.dzero <- post.disc/(1 - post.kappa)
discount <- if(sum(discount == 0)/n.store > 0.5) as.simple_triplet_matrix(discount) else discount
PY.disc <- list(discount = discount, post.disc = post.disc, post.kappa = post.kappa, var.disc = var.disc,
ci.disc = ci.disc, disc.rate = rate, last.disc = discount[n.store], post.d_nonzero = post.dzero)
class(PY.disc) <- "listof"
}
MAP <- as.integer(levels(MAP))[MAP]
uncert.obs <- which(uncertain >= 1/G)
uncertain <- if(sum(uncertain == 0)/n.obs > 0.5) as.simple_triplet_matrix(uncertain) else uncertain
attr(uncertain, "Obs") <- if(sum(uncert.obs) != 0) uncert.obs
if(!label.miss) tab.stat$uncertain <- attr(uncertain, "Obs")
cluster <- list(MAP = MAP, z = z, uncertainty = uncertain, last.z = z[TN.store,])
cluster <- c(cluster, list(post.sizes = sizes, post.ratio = sizes/n.obs, post.pi = post.pi/sum(post.pi),
post.prob = provideDimnames(post.prob, base=list("", gnames)),
PCM = provideDimnames(post_conf_mat(post.prob), base=list(gnames, gnames))),
if(sw["pi.sw"]) list(pi.prop = pi.prop, var.pi = var.pi, ci.pi = ci.pi, last.pi = pi.prop[,TN.store]),
if(!label.miss) list(perf = tab.stat),
if(learn.alpha) list(Alpha = DP.alpha),
if(learn.d) list(Discount = PY.disc),
if(z.avgsim) list(Z.avgsim = z_simavg),
if(is.element(method, c("IMFA", "IMIFA"))) list(lab.rate = sims[[G.ind]][[Q.ind]]$lab.rate))
attr(cluster, "Z.init") <- attr(sims[[G.ind]], "Z.init")
attr(cluster, "Init.Meth") <- attr(sims, "Init.Z")
attr(cluster, "Label.Sup") <- !label.miss
class(cluster) <- "listof"
z.ind <- lapply(Gseq, function(g) MAP == g)
} else {
if(sw["l.sw"]) {
lmats <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$load) else sims[[G.ind]][[Q.ind]]$load
}
z.ind <- list(seq_len(n.obs))
sizes <- n.obs
}
if(inf.Q) {
if(G1 <- G > 1) {
if(isTRUE(adapt)) {
colvec <- lapply(Gseq, function(g) vapply(seq_len(TN.store), function(i) colSums2(abs(.a_drop(lmats[,,g,i, drop=FALSE], drop=3L:4L)) < epsilon, useNames=FALSE)/n.var < prop, logical(ncol(lmats))))
Q.store[] <- t(vapply(colvec, colSums2, numeric(TN.store), useNames=FALSE))
if(isTRUE(forceQg)) {
Q.force <- t(vapply(Gseq, function(g) pmin(Q.store[g,], sizes[g] - 1), numeric(TN.store)))
Q.big <- any(Q.force < Q.store)
Q.store[] <- Q.force
colvec <- lapply(Gseq, function(g) vapply(seq_len(TN.store), function(i) { colvec[[g]][colvec[[g]][,i],i][-seq_len(Q.store[g,i])] <- FALSE; colvec[[g]][,i] }, logical(ncol(lmats))))
}
}
Q.tab <- lapply(apply(Q.store, 1L, function(x) list(table(x, dnn=NULL))), "[[", 1L)
Q.prob <- lapply(Q.tab, prop.table)
class(Q.tab) <- class(Q.prob) <- "listof"
} else {
if(isTRUE(adapt)) {
if(active.crit == "SC") {
eta.tmp <- as.array(sims[[G.ind]][[Q.ind]]$eta)
SC <- lapply(seq_len(TN.store), function(i) .SC_crit(dat, .a_drop(eta.tmp[,,i,drop=FALSE], 3L), .a_drop(lmats[,,i, drop=FALSE], 3L), prop))
colvec <- vapply(lapply(SC, "[[", "nonred"), function(x) match(seq_len(ncol(lmats)), x, nomatch=0) > 0, logical(ncol(lmats)))
Q.store[] <- Q.store[] - vapply(SC, "[[", numeric(1L), "numred")
} else {
colvec <- vapply(seq_len(TN.store), function(i) colSums2(abs(.a_drop(lmats[,,i, drop=FALSE], 3L)) < epsilon, useNames=FALSE)/n.var < prop, logical(ncol(lmats)))
Q.store[] <- colSums2(colvec, useNames=FALSE)
}
colvec <- list(colvec)
}
Q.tab <- table(Q.store, dnn=NULL)
Q.prob <- prop.table(Q.tab)
}
Q.mode <- if(G1) unlist(lapply(Q.tab, function(qt) as.numeric(names(qt[qt == max(qt)])[1L]))) else as.numeric(names(Q.tab[Q.tab == max(Q.tab)])[1L])
Q.med <- if(G1) stats::setNames(ceiling(matrixStats::rowMedians(Q.store,
useNames=FALSE) * 2)/2, gnames) else ceiling(Median(Q.store) * 2)/2
if(!Q.T) {
Q <- switch(EXPR=Q.meth, mode=Q.mode, floor(Q.med))
} else {
Q <- if(G.T) Q else stats::setNames(if(length(Q) == G) Q else rep(Q, G), gnames)
}
Q.CI <- if(G1) round(rowQuantiles(Q.store,
probs=conf.levels)) else round(stats::quantile(Q.store, conf.levels))
GQ.temp4 <- list(Q = Q, Q.Mode = Q.mode, Q.Median = Q.med,
Q.CI = Q.CI, Q.Probs = Q.prob, Q.Counts = Q.tab,
Stored.Q = if(clust.ind) Q.store else drop(Q.store),
Q.Last = Q.store[,TN.store])
GQ.res <- if(inf.G) c(GQ.temp1, GQ.temp4) else c(list(G = G), GQ.temp4)
GQ.res <- if(cshrink) c(GQ.res, list(post.sigma = stats::setNames(rowMeans2(sigmas, useNames=FALSE), gnames))) else GQ.res
GQ.res <- c(GQ.res, list(Criteria = GQ.temp2))
attr(GQ.res, "Q.big") <- ifelse(adapt && G1 && forceQg, Q.big, attr(sims[[G.ind]][[Q.ind]], "Q.big"))
}
Q0 <- Q == 0
Qmax <- max(Q)
if(all(isTRUE(choice), is.element(criterion, c("aicm", "bicm", "dic")))) {
if(all(!G.T, !is.element(method,
c("FA", "IFA")), G == 1)) warning(paste0("Chosen model has only one group:\nNote that the ", criterion, " criterion may exhibit bias toward one-group models\n"), call.=FALSE, immediate.=TRUE)
if(all(!Q.T, method == "MIFA")) {
if(any(Q0)) warning(paste0("Chosen model has ", ifelse(sum(Q0) == G, "zero factors", "a cluster with zero factors"), ":\nNote that the ", criterion, " criterion may exhibit bias toward models ", ifelse(sum(Q0) == G, "with zero factors\n", "where some clusters have zero factors\n")), call.=FALSE, immediate.=TRUE)
} else if(all(Q0)) warning(paste0("Chosen model has zero factors:\nNote that the ", criterion, " criterion may exhibit bias toward zero-factor models\n"), call.=FALSE, immediate.=TRUE)
}
leder.b <- min(n.obs - 1L, Ledermann(n.var, isotropic=is.element(uni.type, c("isotropic", "single"))))
if(any(Q > leder.b)) { warning(paste0("Estimate of Q", ifelse(G > 1, " in one or more clusters ", " "), "is greater than ", ifelse(any(unlist(Q) > n.var), paste0("the number of variables (", n.var, ")"), paste0("the suggested Ledermann upper bound (", leder.b, ")")), "\n"), call.=FALSE, immediate.=TRUE)
} else if(any(floor((n.var - 1)/2) <
Q)) warning("It is advisable that the restriction Q <= floor((P - 1)/2) be respected\n", call.=FALSE, immediate.=TRUE)
if(any(Q >= sizes) &&
G > 1 &&
!attr(sims, "ForceQg")) warning("Estimate of Q is greater than the number of observations in one or more clusters:\nConsider setting 'forceQg' to TRUE\n", call.=FALSE, immediate.=TRUE)
# Retrieve (unrotated) scores
if(no.score <- all(Q0)) {
if(sw["s.sw"]) message("Scores & loadings not stored as model has zero factors\n")
sw["s.sw"] <- FALSE
}
if(inf.Q) {
Lstore <- lapply(Gseq, function(g) storeG[Q.store[g,] >= Q[g]])
if(any(lengths(Lstore) < 2)) {
sw["s.sw"] <- tmpsw["l.sw"] <-
sw["l.sw"] <- FALSE; warning("Forcing non-storage of scores and loadings due to shortage of retained samples\n", call.=FALSE, immediate.=TRUE)
}
eta.store <- storeG[colMaxs(Q.store, value=TRUE) >= Qmax]
} else {
eta.store <- tmp.store
Lstore <- storeG
}
Qseq <- seq_len(Qmax)
store.e <- (if(inf.Q) storeG else tmp.store) %in% eta.store
e.store <- sum(store.e)
QsE <- if(inf.Q) Q.store[,store.e, drop=!clust.ind]
Q0E <- if(inf.Q) QsE == 0 else if(clust.ind) matrix(Q == 0, nrow=G, ncol=e.store) else rep(Q == 0, e.store)
Q0X <- all(Q0E)
QX0 <- any(!Q0E)
frobenius <- error.metrics && all(sw["psi.sw"], sw["mu.sw"] || !sw["u.sw"], any(Q0X, sw["l.sw"]))
if(sw["s.sw"]) {
eta <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$eta)[,,eta.store, drop=FALSE] else sims[[G.ind]][[Q.ind]]$eta[,,eta.store, drop=FALSE]
if(!sw["l.sw"]) warning("Caution advised when examining posterior factor scores: Procrustes rotation has not taken place because loadings weren't stored\n", call.=FALSE, immediate.=TRUE)
}
# Loop over g in G to extract other results
result <- list(list())
for(g in Gseq) {
Qg <- Q[g]
Q0g <- Q0[g]
Qgs <- seq_len(Qg)
sw["l.sw"] <- tmpsw["l.sw"]
if(Q0g) {
if(all(sw["l.sw"],
!no.score)) message(paste0("Loadings ", ifelse(G > 1, paste0("for cluster ", g, " not stored as it"), " not stored as model"), " has zero factors\n"))
sw["l.sw"] <- FALSE
}
# Retrieve loadings and rotate
if(sw["l.sw"]) {
tmpind <- ifelse(inf.Q, which.max(Q.store[g,] == Qg), 1L)
lmat <- if(clust.ind) .a_drop(lmats[,,g,storeG, drop=FALSE], drop=3L) else lmats[,,storeG, drop=FALSE]
l.temp <- .a_drop(lmat[,,tmpind, drop=FALSE], drop=3L)
l.temp <- if(isTRUE(vari.rot)) .vari_max(l.temp, ...)$loadings else l.temp
if(inf.Q) {
rotind <- Lstore[[g]][-1L]
rotind <- rotind[Q.store[g,rotind] > 1]
} else {
rotind <- if(Qg > 1) Lstore[-1L]
}
for(p in rotind) {
proc <- suppressWarnings(Procrustes(X=if(uni) t(lmat[,,p]) else as.matrix(lmat[,,p]), Xstar=l.temp))
lmat[,,p] <- proc$X.new
if(sw["s.sw"] &&
(p %in%
eta.store)) {
p2 <- if(inf.Q) eta.store == p else p
if(clust.ind) {
zp <- z2[p,] == g
eta[zp,,p2] <- .a_drop(eta[zp,,p2, drop=FALSE], 3L) %*% proc$R
} else {
eta[,,p2] <- .a_drop(eta[,,p2, drop=FALSE], 3L) %*% proc$R
}
}
}
if(adapt) {
lmat <- sapply(lapply(storeG, function(i) .a_drop(lmat[,colvec[[g]][,i],i, drop=FALSE], drop=3L)),
function(x) cbind(x, matrix(0L, nrow=n.var, ncol=ncol(lmat) - ncol(x))), simplify="array")
}
}
# Retrieve means, uniquenesses & empirical covariance matrix
if(clust.ind) {
if(sw["mu.sw"]) {
mu <- if(uni) t(mus[,g,storeG]) else as.matrix(mus[,g,storeG])
}
if(sw["psi.sw"]) {
psi <- if(uni) t(psis[,g,storeG]) else as.matrix(psis[,g,storeG])
}
} else {
post.mu <- as.matrix(sims[[G.ind]][[Q.ind]]$post.mu)
post.psi <- as.matrix(sims[[G.ind]][[Q.ind]]$post.psi)
if(sw["mu.sw"]) {
mu <- sims[[G.ind]][[Q.ind]]$mu[,storeG, drop=FALSE]
}
if(sw["psi.sw"]) {
psi <- sims[[G.ind]][[Q.ind]]$psi[,storeG, drop=FALSE]
}
}
if(sw["mu.sw"] &&
is.null(rownames(mu))) {
rownames(mu) <- varnames
}
if(sw["psi.sw"] &&
is.null(rownames(psi))) {
rownames(psi) <- varnames
}
if(sw["l.sw"] &&
is.null(rownames(lmat))) {
rownames(lmat) <- varnames
}
# Compute posterior means and % variation explained
if(sw["mu.sw"]) {
post.mu <- if(clust.ind) rowMeans2(mu,
useNames=FALSE) else post.mu
var.mu <- if(uni) Var(mu) else rowVars(mu)
ci.tmp <- rowQuantiles(mu, probs=conf.levels)
ci.mu <- if(uni) t(ci.tmp) else ci.tmp
}
if(sw["psi.sw"]) {
post.psi <- if(clust.ind) rowMeans2(psi,
useNames=FALSE) else post.psi
var.psi <- if(uni) Var(psi) else rowVars(psi)
ci.tmp <- rowQuantiles(psi, probs=conf.levels)
ci.psi <- if(uni) t(ci.tmp) else ci.tmp
}
if(sw["l.sw"]) {
lmat <- lmat[,Qgs,if(inf.Q) Lstore[[g]] else storeG, drop=FALSE]
post.load <- rowMeans(lmat, dims=2)
var.load <- apply(lmat, c(1L, 2L), Var)
ci.load <- apply(lmat, c(1L, 2L), stats::quantile, conf.levels)
var.exp <- rowSums2(post.load^2, useNames=FALSE)
class(post.load) <- "loadings"
} else if(sw["psi.sw"]) {
if(sizes[g] > 1) {
dat.gg <- dat[z.ind[[g]],, drop=FALSE]
cov.gg <- stats::cov(dat.gg)
}
var.exp <- ifelse(sizes[g] <= 1, 0L, pmax.int(0L, diag(cov.gg) - post.psi))
} else {
var.exp <- NULL
}
results <- list(if(sw["mu.sw"]) list(means = mu,
var.mu = var.mu,
ci.mu = ci.mu),
if(sw["l.sw"]) list(loadings = lmat,
post.load = post.load,
var.load = var.load,
ci.load = ci.load),
if(sw["psi.sw"]) list(psis = psi,
var.psi = var.psi,
ci.psi = ci.psi),
if(sw.mx) list(post.mu = post.mu),
if(sw.px) list(post.psi = post.psi),
if(any(sw["l.sw"],
sw.px)) list(var.exp = var.exp))
result[[g]] <- unlist(results, recursive=FALSE)
}
names(result) <- gnames
GQ.res$Criteria <- c(GQ.res$Criteria, list(sd.AICMs = aicm.sd, sd.BICMs = bicm.sd,
best.models = t(vapply(GQ.res$Criteria, function(x) { inds <- arrayInd(which.max(x), dim(x));
c(clusters = rownames(x)[inds[1L]], factors = colnames(x)[inds[2L]]) }, character(2L)))))
class(GQ.res$Criteria) <- "listof"
class(GQ.res) <- "listof"
attr(GQ.res, "Clusters") <- n.grp
attr(GQ.res, "Criterion") <- criterion
attr(GQ.res, "Factors") <- n.fac
attr(GQ.res, "Supplied") <- c(Q=Q.T, G=G.T)
if(sw["mu.sw"]) {
mus2 <- lapply(result, "[[", "means")
var.mu <- provideDimnames(do.call(cbind, lapply(result, "[[", "var.mu")), base=list(if(uni) "" else varnames, gnames))
ci.mu <- lapply(result, "[[", "ci.mu")
class(mus2) <- class(ci.mu) <- "listof"
means <- list(mus = mus2, var.mu = var.mu, ci.mu = ci.mu)
} else means <- NULL
if(sw["mu.sw"] || sw.mx) {
post.mu <- provideDimnames(do.call(cbind, lapply(result, "[[", "post.mu")), base=list(if(uni) "" else varnames, gnames))
means <- c(means, list(post.mu = post.mu))
if(sw["mu.sw"]) {
means <- means[c(1L, 4L, 2L, 3L)]
}
class(means) <- "listof"
}
if(sw["mu.sw"]) {
last.mu <- if(clust.ind) .a_drop(mus[,,TN.store, drop=FALSE], drop=3L) else mu[,TN.store, drop=FALSE]
colnames(last.mu) <- gnames
means <- c(means, list(last.mu = last.mu))
class(means) <- "listof"
}
if(sw["l.sw"] <- tmpsw["l.sw"] && !all(Q == 0)) {
lnames <- lapply(Q, function(q) if(q > 0) paste0("Factor", seq_len(q)))
lmats2 <- .matnames(lapply(result, "[[", "loadings"), lnames)
post.load <- .matnames(lapply(result, "[[", "post.load"), lnames)
var.load <- .matnames(lapply(result, "[[", "var.load"), lnames)
ci.load <- .matnames(lapply(result, "[[", "ci.load"), lnames, dim=3)
last.lmat <- lapply(lmats2, function(x) { if(!is.null(x)) { x <- .a_drop(x[,,dim(x)[3L], drop=FALSE], drop=3L); class(x) <- "loadings" }; x })
loads <- list(lmats = lmats2, post.load = post.load, var.load = var.load, ci.load = ci.load, last.load = last.lmat)
loads <- lapply(loads, function(x) { class(x) <- "listof"; x} )
if(cshrink) attr(loads, "Sigma") <- GQ.res$post.sigma
class(loads) <- "listof"
}
if(sw["psi.sw"]) {
psis2 <- lapply(result, "[[", "psis")
var.psi <- provideDimnames(do.call(cbind, lapply(result, "[[", "var.psi")), base=list(if(uni) "" else varnames, gnames))
ci.psi <- lapply(result, "[[", "ci.psi")
class(psis2) <- class(ci.psi) <- "listof"
uniquenesses <- list(psis = psis2, var.psi = var.psi, ci.psi = ci.psi)
} else uniquenesses <- NULL
if(sw["psi.sw"] || sw.px) {
post.psi <- provideDimnames(do.call(cbind, lapply(result, "[[", "post.psi")), base=list(if(uni) "" else varnames, gnames))
uniquenesses <- c(uniquenesses, list(post.psi = post.psi))
if(sw["psi.sw"]) {
uniquenesses <- uniquenesses[c(1L, 4L, 2L, 3L)]
}
class(uniquenesses) <- "listof"
}
if(sw["psi.sw"]) {
last.psi <- if(clust.ind) .a_drop(psis[,,TN.store, drop=FALSE], drop=3L) else psi[,TN.store, drop=FALSE]
colnames(last.psi) <- gnames
uniquenesses <- c(uniquenesses, list(last.psi = last.psi))
class(uniquenesses) <- "listof"
}
dots$type <- ifelse(is.null(dots$type), "F", dots$type[1L])
# Calculate estimated covariance matrices & compute error metrics
if(error.metrics) {
Eseq <- seq_len(e.store)
sw.pi <- (equal.pro || G == 1) || sw["pi.sw"]
Fro <- mse <- mae <- medse <-
medae <- rmse <- nrmse <- rep(NA, e.store)
if(frobenius) {
orig.data <- as.data.frame(dat)
Pseq <- seq_len(n.var)
if(is.null(dots$dbreaks)) {
xbins <- if(isFALSE(uni)) lengths(lapply(orig.data, function(x) suppressWarnings(graphics::hist(x, plot=FALSE, ...))$counts)) else length(graphics::hist(orig.data[[1L]], plot=FALSE, ...)$counts)
dbreaks <- lapply(Pseq, function(p) c(-Inf, as.numeric(sub("\\((.+),.*", "\\1", levels(cut(orig.data[[p]], xbins[p]))[-1L])), Inf))
} else {
dbreaks <- dots$dbreaks
if(!is.list(dbreaks) ||
!all(vapply(dbreaks, is.factor, logical(1L))) ||
n.var !=
length(dbreaks)) stop(paste0("'dbreaks' must be a list of length P=", n.var, " of factors as outputted by the function 'cut'"), call.=FALSE)
dbreaks <- lapply(Pseq, function(p) unique(c(-Inf, as.numeric(sub("\\((.+),.*", "\\1", levels(dbreaks[[p]])[-1L])), Inf)))
xbins <- lengths(dbreaks) - 1L
}
dcounts <- mapply(tabulate, mapply(cut, orig.data, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
nbins <- max(xbins)
dat.bins <- vapply(dcounts, .padding, integer(nbins), nbins)
datnorm <- norm(dat.bins, dots$type)
rcounts <- vector("list", e.store)
} else if(Q0X) { warning("Need the means and uniquenesses to have been stored for zero-factor models in order to compute the posterior predictive reconstruction error\n", call.=FALSE, immediate.=TRUE)
} else warning("Need the means, uniquenesses, and loadings to have been stored in order to compute the posterior predictive reconstruction error\n", call.=FALSE, immediate.=TRUE)
}
if(clust.ind) {
if(frobenius) {
if(cov.met) {
a <- Reduce("+", lapply(Gseq, function(g) post.pi[g] * (tcrossprod(post.mu[,g]) + (if(Q0[g]) 0L else tcrossprod(post.load[[g]])) + (if(uni) post.psi[,g] else diag(post.psi[,g])))))
b <- Reduce("+", lapply(Gseq, function(g) post.pi[g] * post.mu[,g]))
cov.est <- a - tcrossprod(b)
if(sw.pi) {
if(equal.pro ||
G == 1) {
pi2 <- matrix(1/G, nrow=G, ncol=e.store)
} else {
pi2 <- pi.prop[,store.e, drop=FALSE]
}
} else warning("Mixing proportions not stored: can't compute all error metrics\n", call.=FALSE, immediate.=TRUE)
}
mu2 <- mus[,,store.e, drop=FALSE]
psi2 <- psis[,,store.e, drop=FALSE]
if(QX0) lmat2 <- lmats[,,,store.e, drop=FALSE]
if(!inf.Q) QsEr <- Q
tab.z <- tab.z[store.e,, drop=FALSE]
for(r in Eseq) {
Q0Er <- Q0E[,r]
if(frobenius) {
if(inf.Q) QsEr <- QsE[,r]
rdat <- suppressWarnings(sim_IMIFA_data(N=n.obs, G=G, P=n.var, Q=QsEr, nn=tab.z[r,], mu=mu2[,,r], psi=psi2[,,r], method="marginal", forceQg=FALSE,
loadings=if(!all(Q0Er)) lapply(Gseq, function(g) matrix(lmat2[,if(inf.Q) seq_len(QsEr[g]) else Qseq,g,r], nrow=n.var, ncol=QsEr[g]))))
rbinsx <- mapply(tabulate, mapply(cut, rdat, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
rbins <- vapply(rbinsx, .padding, integer(nbins), nbins)
Frob <- norm(dat.bins - rbins, dots$type)
rFro <- norm(rbins, dots$type)
minF <- abs(datnorm - rFro)
Fro[r] <- (Frob - minF)/(datnorm + rFro - minF)
rcounts[[r]] <- rbinsx
}
if(sw.pi && cov.met) {
a <- Reduce("+", lapply(Gseq, function(g) pi2[g,r] * (tcrossprod(mu2[,g,r]) + (if(Q0Er[g]) 0L else tcrossprod(lmat2[,,g,r])) + (if(uni) psi2[,g,r] else diag(psi2[,g,r])))))
b <- Reduce("+", lapply(Gseq, function(g) pi2[g,r] * mu2[,g,r]))
Sigma <- a - tcrossprod(b)
error <- cov.emp - Sigma
sqerr <- error^2
abserr <- abs(error)
mse[r] <- mean(sqerr)
medse[r] <- Median(sqerr)
mae[r] <- mean(abserr)
medae[r] <- Median(abserr)
rmse[r] <- sqrt(mse[r])
nrmse[r] <- rmse[r]/cov.range
}
}
} else if(cov.met) {
if(!sw["mu.sw"]) { warning("Means not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0, !sw["l.sw"],
!sw["psi.sw"])) { warning("Loadings & Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0,
!sw["l.sw"])) { warning("Loadings not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else warning("Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
}
} else {
if(any(sw["l.sw"], Q0X)) {
if(cov.met) {
cov.est <- diag(post.psi[,1L, drop=!uni]) + (if(Q0X) 0L else tcrossprod(post.load[[1L]]))
}
if(error.metrics && sw["psi.sw"]) {
if(frobenius) {
mu2 <- if(!sw["u.sw"]) matrix(0L, nrow=n.var, ncol=e.store) else mu[,store.e, drop=FALSE]
}
if(QX0) {
lmat2 <- lmats[,,store.e, drop=FALSE]
}
psi2 <- psi[,store.e, drop=FALSE]
if(!inf.Q) {
QsEr <- Q
QsMs <- Qseq
}
for(r in Eseq) {
if(frobenius) {
if(inf.Q) {
QsEr <- QsE[r]
QsMs <- seq_len(QsEr)
}
rdat <- suppressWarnings(sim_IMIFA_data(N=n.obs, G=G, P=n.var, Q=QsEr, mu=mu2[,r], psi=psi2[,r], method="marginal", forceQg=FALSE,
loadings=if(!Q0E[r]) matrix(lmat2[,QsMs,r], nrow=n.var, ncol=QsEr)))
rbinsx <- mapply(tabulate, mapply(cut, rdat, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
rbins <- vapply(rbinsx, .padding, integer(nbins), nbins)
Frob <- norm(dat.bins - rbins, dots$type)
rFro <- norm(rbins, dots$type)
minF <- abs(datnorm - rFro)
Fro[r] <- (Frob - minF)/(datnorm + rFro - minF)
rcounts[[r]] <- rbinsx
}
if(cov.met) {
Sigma <- diag(psi2[,r, drop=!uni]) + (if(Q0E[r]) 0L else if(uni) crossprod(lmat2[,,r]) else tcrossprod(lmat2[,,r]))
error <- cov.emp - Sigma
sqerr <- error^2
abserr <- abs(error)
mse[r] <- mean(sqerr)
medse[r] <- Median(sqerr)
mae[r] <- mean(abserr)
medae[r] <- Median(abserr)
rmse[r] <- sqrt(mse[r])
nrmse[r] <- rmse[r]/cov.range
}
}
} else if(error.metrics) warning("Uniquenesses not stored: can't compute error metrics\n", call.=FALSE, immediate.=TRUE)
} else if(error.metrics) {
if(all(QX0, !sw["l.sw"],
!sw["psi.sw"])) { warning("Loadings & Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0,
!sw["l.sw"])) warning("Loadings not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
}
}
if(all(sw["psi.sw"] || !clust.ind, any(sw["l.sw"], Q0X))) {
var.exps <- lapply(lapply(result, "[[", "var.exp"), function(x) if(is.null(x)) NA else x)
clust.exps <- if(sum(is.na(var.exps)) == G) NULL else vapply(var.exps, sum, na.rm=TRUE, numeric(1L))/n.var
var.exps <- if(sum(is.na(var.exps)) == G) NULL else if(clust.ind) stats::setNames(rowSums2(sweep(do.call(cbind, var.exps), 2L, FUN="*", post.pi, check.margin=FALSE), na.rm=TRUE, useNames=FALSE), varnames) else stats::setNames(var.exps[[1L]], varnames)
Err <- list(Var.Exps = var.exps, Clust.Exps = clust.exps, Exp.Var = ifelse(clust.ind, sum(clust.exps * post.pi), unname(clust.exps)))
if(all(error.metrics, sw["mu.sw"] || !clust.ind)) {
if(cov.met) {
cov.est <- as.matrix(cov.est)
dimnames(cov.est) <- list(varnames, varnames)
err <- cov.emp - cov.est
sqerr <- err^2
abserr <- abs(err)
mse2 <- mean(sqerr)
rmse2 <- sqrt(mse2)
post.met <- c(MSE = mse2, MEDSE = Median(sqerr), MAE = mean(abserr), MEDAE = Median(abserr), RMSE = rmse2, NRMSE = rmse2/cov.range)
Err <- c(list(Empirical.Cov = cov.emp, Estimated.Cov = cov.est, Post = post.met), Err)
}
if(sw["psi.sw"] && any(all(cov.met, sw.pi), frobenius)) {
metrics <- rbind(MSE = mse, MEDSE = medse, MAE = mae, MEDAE = medae, RMSE = rmse, NRMSE = nrmse)
metrics <- if(frobenius) rbind(metrics, PPRE = Fro) else metrics
metrics <- metrics[!rowAlls(is.na(metrics), useNames=FALSE),, drop=FALSE]
metricCI <- rowQuantiles(metrics, probs=conf.levels)
metricCI <- if(cov.met && sw.pi) metricCI else provideDimnames(t(metricCI), base=list("PPRE", names(metricCI)))
med.met <- stats::setNames(matrixStats::rowMedians(metrics,
useNames=FALSE), rownames(metrics))
mean.met <- stats::setNames(rowMeans2(metrics, useNames=FALSE), rownames(metrics))
last.met <- c(MSE = mse[e.store], MEDSE = medse[e.store], MAE = mae[e.store], MEDAE = medae[e.store], RMSE = rmse[e.store], NRMSE = nrmse[e.store])
last.met <- if(frobenius) c(last.met, PPRE = Fro[e.store]) else last.met
last.met <- last.met[!is.na(last.met)]
Err <- if(cov.met) c(list(Median = med.met, Avg = mean.met, CIs = metricCI), Err[seq_len(4L)], if(all(cov.met, sw.pi)) c(list(Last.Cov = Sigma)), c(list(Final = last.met)), Err[5L:6L]) else c(list(Avg = mean.met, CIs = metricCI), Err[seq_len(3L)], list(Final = last.met))
if(frobenius) {
rcounts <- stats::setNames(lapply(Pseq, function(p) t(rowQuantiles(sapply(rcounts, "[[", p), probs=c(conf.levels[1L], 0.5, conf.levels[2L])))), varnames)
class(dcounts) <- "listof"
class(rcounts) <- "listof"
Err <- c(Err, list(PPRE = Fro, DatCounts = dcounts, RepCounts = rcounts, Breaks = dbreaks))
attr(Err, "Norm") <- dots$type
}
attr(Err, "ESS") <- e.store
errs <- ifelse(frobenius, ifelse(cov.met && sw.pi, "All", "PPRE"), "Covs")
} else errs <- ifelse(cov.met, "Post", "Vars")
} else errs <- "Vars"
class(Err) <- "listof"
} else errs <- "None"
error.metrics <- errs != "None"
if(sw["s.sw"]) {
if(adapt) {
eta <- if(method == "IFA") {
sapply(lapply(seq_len(e.store), function(i) .a_drop(eta[,colvec[[g]][,store.e][,i],i, drop=FALSE], drop=3L)),
function(x) cbind(x, matrix(0L, nrow=n.obs, ncol=ncol(eta) - ncol(x))), simplify="array")
} else {
sapply(seq_len(e.store), function(i) do.call(rbind, lapply(Gseq,
function(g) .a_drop(eta[z[i,] == g, c(which(colvec[[g]][,store.e][,i]), setdiff(seq_len(ncol(eta)), which(colvec[[g]][,store.e][,i]))),i, drop=FALSE], drop=3L)))[obsnames,, drop=FALSE], simplify="array")
}
}
eta <- eta[,Qseq,, drop=FALSE]
colnames(eta) <- paste0("Factor", Qseq)
scores <- list(eta = eta, post.eta = rowMeans(eta, dims=2), var.eta = apply(eta, c(1L, 2L), Var),
ci.eta = apply(eta, c(1L, 2L), stats::quantile, conf.levels), last.eta = .a_drop(eta[,,e.store, drop=FALSE], drop=3L))
attr(scores, "Eta.store") <- e.store
class(scores) <- "listof"
}
result <- c(if(exists("cluster", envir=environment())) list(Clust = cluster),
if(error.metrics) list(Error = Err), list(GQ.results = GQ.res),
if(sw["mu.sw"] || sw.mx) list(Means = means),
if(sw["l.sw"]) list(Loadings = loads),
if(sw["s.sw"]) list(Scores = scores),
if(sw["psi.sw"] || sw.px) list(Uniquenesses = uniquenesses))
attr(result, "Adapt") <- any(adapt, adapted)
attr(result, "Alph.step") <- if(is.element(method, c("IMFA", "IMIFA", "OMFA", "OMIFA"))) learn.alpha
attr(result, "Alpha") <- if(!learn.alpha) attr(sims, "Alpha")
attr(result, "C.Shrink") <- cshrink
attr(result, "Call") <- call
attr(result, "Center") <- attr(sims, "Center")
attr(result, "Choice") <- choice
attr(result, "Conf.Level") <- conf.level
attr(result, "Criterion") <- criterion
attr(result, "Disc.step") <- if(is.element(method, c("IMFA", "IMIFA"))) learn.d
attr(result, "Discount") <- if(is.element(method, c("IMFA", "IMIFA")) && !learn.d) attr(sims, "Discount")
attr(result, "Errors") <- errs
attr(result, "Equal.Pi") <- equal.pro
attr(result, "Exchange") <- attr(sims, "Exchange")
attr(result, "ForceQg") <- attr(sims, "ForceQg") || (adapt && G1 && forceQg)
attr(result, "G.init") <- if(inf.G) attr(sims, "G.init")
attr(result, "G.Mean") <- attr(sims, "G.Mean")
attr(result, "G.Scale") <- attr(sims, "G.Scale")
attr(result, "Ind.Slice") <- if(is.element(method, c("IMFA", "IMIFA"))) attr(sims, "Ind.Slice")
attr(result, "Kappa0") <- attr(sims, "Kappa0")
attr(result, "Method") <- method
attr(result, "N.Loadstore") <- if(inf.Q) lengths(Lstore) else rep(TN.store, G)
attr(result, "Name") <- data.name
attr(result, "N.Store") <- TN.store
attr(result, "Nowarn.G") <- ifelse(inf.G || length(n.grp) == 1, TRUE,
ifelse(G == min(n.grp), paste0("Best model occurs at the min of the number of components considered\n"),
ifelse(G == max(n.grp), paste0("Best model occurs at the max of the number of components considered\n"), TRUE)))
attr(result, "Nowarn.Q") <- ifelse(inf.Q || length(n.fac) == 1, TRUE,
ifelse(Q == min(n.fac), paste0("Best model occurs at the min of the number of factors considered\n"),
ifelse(Q == max(n.fac), paste0("Best model occurs at the max of the number of factors considered\n"), TRUE)))
attr(result, "Obs") <- n.obs
attr(result, "Obsnames") <- obsnames
attr(result, "Pitman") <- attr(sims, "Pitman")
attr(result, "range.G") <- attr(sims, "Clusters")
attr(result, "range.Q") <- attr(sims, "Factors")
attr(result, "Scale") <- attr(sims, "Scale")
attr(result, "Scaling") <- attr(sims, "Scaling")
attr(result, "Sd0.drop") <- attr(sims, "Sd0.drop")
attr(result, "Store") <- tmp.store
attr(result, "Switch") <- sw
attr(result, "Thresh") <- attr(sims, "Thresh")
attr(result, "Truncated") <- attr(sims, "Truncate")
attr(result, "TuneZeta") <- attr(sims, "TuneZeta")
attr(result, "Uni.Meth") <- uni.meth
attr(result, "Varnames") <- varnames
attr(result, "Vars") <- n.var
attr(result, "Z.sim") <- z.avgsim
class(result) <- "Results_IMIFA"
cat(print(result))
return(result)
}
Keefe-Murphy/IMIFA documentation built on Jan. 31, 2024, 2:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions get_IMIFA_results.IMIFA get_IMIFA_results
Documented in get_IMIFA_results
#' Extract results, conduct posterior inference and compute performance metrics for MCMC samples of models from the IMIFA family
#'
#' This function post-processes simulations generated by \code{\link{mcmc_IMIFA}} for any of the IMIFA family of models. This includes accounting for label switching, and accounting for rotational invariance via Procrustean methods. It can be re-ran at little computational cost in order to extract different models explored by the sampler used for \code{sims}, without having to re-run the model itself. New results objects using different numbers of clusters and different numbers of factors (if visited by the model in question), or using different model selection criteria (if necessary) can be generated with ease. Posterior predictive checking of the appropriateness of the fitted model is also facilitated.
#' @param sims An object of class \code{"IMIFA"} generated by \code{\link{mcmc_IMIFA}}.
#' @param burnin Optional additional number of iterations to discard. Defaults to 0, corresponding to no additional burnin. See \code{\link{mixfaControl}} for the default \code{burnin} settings used previously by \code{\link{mcmc_IMIFA}}.
#' @param thinning Optional interval for extra thinning to be applied. Defaults to 1, corresponding to no additional thinning. See \code{\link{mixfaControl}} for the default \code{thinning} settings used previously by \code{\link{mcmc_IMIFA}}.
#' @param G If this argument is not specified, results will be returned with the optimal number of clusters. If different numbers of clusters were explored in \code{sims} for the \code{"MFA"} or \code{"MIFA"} methods, supplying an integer value allows pulling out a specific solution with \code{G} clusters, even if the solution is sub-optimal.
#'
#' Similarly, this allows retrieval of samples corresponding to a solution, if visited, with \code{G} clusters for the \code{"OMFA"}, \code{"OMIFA"}, \code{"IMFA"} and \code{"IMIFA"} methods.
#' @param Q If this argument is not specified, results will be returned with the optimal number of factors. If different numbers of factors were explored in \code{sims} for the \code{"FA"}, \code{"MFA"}, \code{"OMFA"} or \code{"IMFA"} methods, this allows pulling out a specific solution with \code{Q} factors, even if the solution is sub-optimal.
#'
#' Similarly, this allows retrieval of samples corresponding to a solution, if visited, with \code{Q} factors for the \code{"IFA"}, \code{"MIFA"}, \code{"OMIFA"} and \code{"IMIFA"} methods if adaptation has already taken place. Can be supplied as a scalar or a vector of values for each cluster.
#'
#' If adaptation didn't take place during model-fitting, \code{adapt=TRUE} can be supplied here to determine the optimal cluster-specific numbers of non-redundant factors (see below).
#' @param criterion The criterion to use for model selection, where model selection is only required if more than one model was run under the \code{"FA"}, \code{"MFA"}, \code{"MIFA"}, \code{"OMFA"} or \code{"IMFA"} methods when \code{sims} was created via \code{\link{mcmc_IMIFA}}. Defaults to \code{bicm}, but note that these are \emph{all} calculated; this argument merely indicates which one will form the basis of the construction of the output.
#'
#' Note that the first three options here might exhibit bias in favour of zero-factor models for the finite factor \code{"FA"}, \code{"MFA"}, \code{"OMFA"} and \code{"IMFA"} methods and might exhibit bias in favour of one-cluster models for the \code{"MFA"} and \code{"MIFA"} methods. The \code{aic.mcmc} and \code{bic.mcmc} criteria will only be returned for finite factor models.
#' @param adapt A logical indicating if adaptation should be applied to the stored loadings and scores matrices to truncate the cluster-specific number(s) of non-redundant factors. This argument is only relevant if \code{adapt=FALSE} was used when initially fitting the model (otherwise, adaptation has already taken place during sampling, by default), and hence defaults to \code{FALSE} here. Relevant parameters from \code{\link{mgpControl}} (namely \code{active.crit}, \code{eps}, \code{prop}, and \code{forceQg}) can be passed via the \code{...} construct, but will default to their values under \code{\link{mgpControl}} if not specified. Note that \code{adapt=TRUE} is only invoked if the relevant parameters were stored via \code{\link{storeControl}} when running the model: loadings only for \code{active.crit="BD"}, loadings and scores for \code{active.crit="SC"} for \code{"IFA"} models.
#' @param G.meth If the object in \code{sims} arises from the \code{"OMFA"}, \code{"OMIFA"}, \code{"IMFA"} or \code{"IMIFA"} methods, this argument determines whether the optimal number of clusters is given by the mode or median of the posterior distribution of \code{G}. Defaults to \code{"mode"}. Often the mode and median will agree in any case.
#' @param Q.meth If the object in \code{sims} arises from the \code{"IFA"}, \code{"MIFA"}, \code{"OMIFA"} or \code{"IMIFA"} methods, this argument determines whether the optimal number of latent factors is given by the mode or median of the posterior distribution of \code{Q}. Defaults to \code{"mode"}. Often the mode and median will agree in any case.
#' @param conf.level The confidence level to be used throughout for credible intervals for all parameters of inferential interest, and error metrics if \code{error.metrics=TRUE}. Defaults to \code{0.95}.
#' @param error.metrics A logical activating or deactivating posterior predictive checking: i.e. controlling whether metrics quantifying a) the posterior predictive reconstruction error (PPRE) between bin counts of the data and bin counts of replicate draws from the posterior distribution & and b) the error between the empirical and estimated covariance matrices should be computed. These are computed for every \emph{valid} retained iteration (see \code{Details}). Defaults to \code{TRUE}, but can be time-consuming for models which achieve clustering. These error metrics, and the uncertainty associated with them, can be visualised via \code{\link{plot.Results_IMIFA}}. Depending on what parameters were stored when calling \code{\link{mcmc_IMIFA}}, potentially not all error metrics will be available to compute.
#'
#' The Frobenius norm is used in the computation of the PPRE, by default, but the \code{type} of \code{\link{norm}} can be changed via the \code{...} construct below. So too can the breakpoints (\code{dbreaks}) used to bin the data and the posterior predictive replicate data sets. Some caution is advised in the latter case.
#' @param vari.rot Logical indicating whether the loadings matrix/matrices template(s) should be \code{\link[stats]{varimax}} rotated first, prior to the Procrustes rotation steps. Defaults to \code{FALSE}. Not necessary at all for clustering purposes, or inference on the covariance matrix, but useful if interpretable inferences on the loadings matrix/matrices are desired. Arguments to \code{\link[stats]{varimax}} can be passed via the \code{...} construct, but note that the argument \code{normalize} here defaults to \code{FALSE}.
#' @param z.avgsim Logical (defaults to \code{FALSE}) indicating whether the clustering should also be summarised with a call to \code{\link{Zsimilarity}} by the clustering with minimum mean squared error to the similarity matrix obtained by averaging the stored adjacency matrices, in addition to the MAP estimate.
#'
#' Note that the MAP clustering is computed \emph{conditional} on the estimate of the number of clusters (whether that be the modal estimate or the estimate according to \code{criterion}) and other parameters are extracted conditional on this estimate of \code{G}: however, in contrast, the number of distinct clusters in the summarised labels obtained by specifying \code{z.avgsim=TRUE} may not necessarily coincide with the MAP estimate of \code{G}, but it may provide a useful alternative summary of the partitions explored during the chain, and the user is free to call \code{\link{get_IMIFA_results}} again with the new suggested \code{G} value.
#'
#' Please be warned that this feature requires loading the \code{\link[mcclust]{mcclust}} package. This is liable to take considerable time to compute, and may not even be possible if the number of observations &/or number of stored iterations is large and the resulting matrix isn't sufficiently sparse. When \code{z.avgsim=TRUE}, both the summarised clustering and the similarity matrix are stored: the latter can be visualised as part of a call to \code{\link{plot.Results_IMIFA}}.
#' @param zlabels For any method that performs clustering, the true labels can be supplied if they are known in order to compute clustering performance metrics. This also has the effect of ordering the MAP labels (and thus the ordering of cluster-specific parameters) to most closely correspond to the true labels if supplied.
#' @param nonempty For \code{"MFA"} and \code{"MIFA"} models ONLY: a logical indicating whether only iterations with non-empty components should be retained. Defaults to \code{TRUE}, but may lead to empty chains - conversely, \code{FALSE} may lead to empty components and related errors.
#' @param x,object,MAP,... Arguments required for the \code{print.Results_IMIFA} and \code{summary.Results_IMIFA} functions: \code{x} and \code{object} are objects of class \code{"Results_IMIFA"} resulting from a call to \code{\link{get_IMIFA_results}}. \code{MAP} is a logical which governs whether a table of the MAP classification is printed, while \code{...} gathers additional arguments to those functions.
#'
#' Users can also pass the \code{type} argument to the \code{\link{norm}} function when \code{isTRUE(error.metrics)} and the posterior predictive reconstruction error (PPRE) is calculated. By default the Frobenius norm (\code{type="F"}) is employed.
#'
#' When \code{adapt=TRUE} here, relevant arguments from \code{\link{mgpControl}} (namely \code{active.crit}, \code{eps}, \code{prop}, and \code{forceQg}) can be supplied here too, though they retain their default values from \code{\link{mgpControl}} otherwise.
#'
#' Finally, the \code{...} construct also allows arguments to \code{\link[stats]{varimax}} to be passed to \code{\link{get_IMIFA_results}} itself, when \code{isTRUE(vari.rot)}, or arguments to \code{\link[graphics]{hist}} when \code{isTRUE(error.metrics)}, in order to guide construction of the bins. Additionally, by passing the argument \code{dbreaks} via the \code{...} construct, the bins can be specified directly. However, caution is advised in doing so; in particular, the bins must be constructed on data which has been standardised in the same way as the data modelled within \code{\link{mcmc_IMIFA}}.
#'
#' @details The function also performs post-hoc corrections for label switching, as well as post-hoc Procrustes rotation of loadings matrices and scores, in order to ensure sensible posterior parameter estimates, computes error metrics, constructs credible intervals, and generally transforms the raw \code{sims} object into an object of class \code{"Results_IMIFA"} in order to prepare the results for plotting via \code{\link{plot.Results_IMIFA}}.
#'
#' For the infinite factor methods, iterations where the maximum number of factors was greater than or equal to the maximum of the estimated cluster-specific factors are retained for posterior summaries of the scores, in order to preserve the estimated dimension of the scores matrices. Similarly, these are also the \emph{valid} iterations used for the computation of the averages and credible intervals for the error metrics. For the finite factor models, \emph{all} retained iterations are used in both instances (i.e. both for the scores and the error metrics).
#'
#' In all cases, only iterations with \code{G} non-empty components are retained.
#' @return An object of class \code{"Results_IMIFA"} to be passed to \code{\link{plot.Results_IMIFA}} for visualising results. Dedicated \code{print} and \code{summary} functions also exist for objects of this class. The object, say \code{x}, is a list of lists, the most important components of which are:
#' \item{\code{Clust}}{Everything pertaining to clustering performance can be found here for all but the \code{"FA"} and \code{"IFA"} methods (or models where the estimate number of clusters is 1), in particular \code{x$Clust$MAP}, the MAP summary of the posterior clustering, the last valid sample of cluster labels \code{x$Clust$last.z}, the matrix of posterior cluster membership probabilities \code{x$Clust$post.prob}, and the posterior confusion matrix \code{x$Clust$PCM}.
#'
#' More detail is given if known \code{zlabels} are supplied: performance is always evaluated against the MAP clustering, with additional evaluation against the alternative clustering computed if \code{z.avgsim=TRUE}. Posterior summaries of the mixing proportions, and the concentration/discount parameters, if necessary, are also included here, as well as the last valid samples of each.}
#' \item{\code{Error}}{Everything pertaining the model fit assessment can be found here, incl. the distribution of the PPRE values and associated bin counts for the replicate draws, as well as average error metrics (e.g. MSE, RMSE), and credible intervals quantifying the associated uncertainty, between the empirical and estimated covariance matrix/matrices, both of which are also included.}
#' \item{\code{GQ.results}}{Everything pertaining to model choice can be found here, incl. posterior summaries for the estimated number of clusters and estimated number of factors, if applicable to the method employed. Model selection criterion values are also accessible here.}
#' \item{\code{Means}}{Posterior summaries for the means, after conditioning on \code{G}.}
#' \item{\code{Loadings}}{Posterior summaries for the factor loadings matrix/matrices, after conditioning on \code{G} and \code{Q}. Posterior mean loadings given by \code{x$Loadings$post.load} are given the \code{\link[stats]{loadings}} class for printing purposes and thus the manner in which they are displayed can be modified.
#'
#' The number of iterations retained for posterior summaries of the loadings may vary for different clusters for the infinite factor methods, corresponding to iterations where the cluster-specific number of factors was greater than or equal to the modal estimate of the cluster-specific number of factors.}
#' \item{\code{Scores}}{Posterior summaries for the latent factor scores, after conditioning on the maximum of the estimated number of cluster-specific factors. Summaries are given for the \emph{single} matrix of factor scores. See \code{\link{scores_MAP}} to decompose these summaries into sub-matrices according to the MAP partition (for models which achieve clustering).
#'
#' For the infinite factor methods, iterations where the maximum number of factors was greater than or equal to the maximum of the estimated cluster-specific factors are retained for posterior summaries of the scores, in order to preserve the estimated dimension of the scores matrices.}
#' \item{\code{Uniquenesses}}{Posterior summaries for the uniquenesses, after conditioning on \code{G}.}
#'
#' The objects \code{Means}, \code{Loadings}, \code{Scores} and \code{Uniquenesses} (if stored when calling \code{\link{mcmc_IMIFA}}!) also contain, as well as the posterior summaries, the entire chain of valid samples of each, as well as, for convenience, the last valid samples of each (after conditioning on the modal \code{G} and \code{Q} values, and accounting for label switching, and rotational invariance via Procrustes rotation).
#' @note For the \code{"IMIFA"}, \code{"IMFA"}, \code{"OMIFA"}, and \code{"OMFA"} methods, the retained mixing proportions are renormalised after conditioning on the modal \code{G}. This is especially necessary for the computation of the \code{error.metrics}, just note that the values on which posterior inference are conducted will ever so slightly differ from the actually sampled values.
#'
#' Due to the way the offline label-switching correction is performed, different runs of this function may give \emph{very slightly} different results in terms of the cluster labellings (and by extension the parameters, which are permuted in the same way), but only if the chain was run for an extremely small number of iterations, well below the number required for convergence, and samples of the cluster labels match poorly across iterations (particularly if the number of clusters suggested by those sampled labels is high).
#' @keywords IMIFA main
#' @include MainFunction.R
#' @export
#' @importFrom Rfast "colMaxs" "colTabulate" "Median" "rowMaxs" "rowOrder" "rowSort" "rowTabulate" "rowVars" "Var"
#' @importFrom mclust "classError"
#' @importFrom matrixStats "colSums2" "colMeans2" "rowAlls" "rowMeans2" "rowMedians" "rowQuantiles" "rowSums2"
#' @importFrom slam "as.simple_triplet_matrix"
#'
#' @seealso \code{\link{plot.Results_IMIFA}}, \code{\link{mcmc_IMIFA}}, \code{\link{Zsimilarity}}, \code{\link{scores_MAP}}, \code{\link{sim_IMIFA_model}}, \code{\link{Procrustes}}, \code{\link[stats]{varimax}}, \code{\link{norm}}, \code{\link{mgpControl}}, \code{\link{storeControl}}
#' @references Murphy, K., Viroli, C., and Gormley, I. C. (2020) Infinite mixtures of infinite factor analysers, \emph{Bayesian Analysis}, 15(3): 937-963. <\doi{10.1214/19-BA1179}>.
#'
#' @author Keefe Murphy - <\email{keefe.murphy@@mu.ie}>
#' @usage
#' get_IMIFA_results(sims = NULL,
#' burnin = 0L,
#' thinning = 1L,
#' G = NULL,
#' Q = NULL,
#' criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
#' adapt = FALSE,
#' G.meth = c("mode", "median"),
#' Q.meth = c("mode", "median"),
#' conf.level = 0.95,
#' error.metrics = TRUE,
#' vari.rot = FALSE,
#' z.avgsim = FALSE,
#' zlabels = NULL,
#' nonempty = TRUE,
#' ...)
#' @examples
#' \donttest{# data(coffee)
#' # data(olive)
#'
#' # Run a MFA model on the coffee data over a range of clusters and factors.
#' # simMFAcoffee <- mcmc_IMIFA(coffee, method="MFA", range.G=2:3, range.Q=0:3, n.iters=1000)
#'
#' # Accept all defaults to extract the optimal model.
#' # resMFAcoffee <- get_IMIFA_results(simMFAcoffee)
#'
#' # Instead let's get results for a 3-cluster model, allowing Q be chosen by aic.mcmc.
#' # resMFAcoffee2 <- get_IMIFA_results(simMFAcoffee, G=3, criterion="aic.mcmc")
#'
#' # Run an IMIFA model on the olive data, accepting all defaults.
#' # simIMIFAolive <- mcmc_IMIFA(olive, method="IMIFA", n.iters=10000)
#'
#' # Extract optimum results
#' # Estimate G & Q by the median of their posterior distributions
#' # Construct 90% credible intervals and try to return the similarity matrix.
#' # resIMIFAolive <- get_IMIFA_results(simIMIFAolive, G.meth="median", Q.meth="median",
#' # conf.level=0.9, z.avgsim=TRUE)
#' # summary(resIMIFAolive)
#'
#' # Simulate new data from the above model
#' # newdata <- sim_IMIFA_model(resIMIFAolive)
#'
#' # Fit an IFA model without adaptation and examine results with and without post-hoc adaptation
#' # Use the "SC" criterion for determining the number of active factors
#' # simIFAnoadapt <- mcmc_IMIFA(olive, method="IFA", n.iters=10000, adapt=FALSE)
#' # resIFAnoadapt <- get_IMIFA_results(simIFAnoadapt)
#' # resIFApostadapt <- get_IMIFA_results(simIFAnoadapt, adapt=TRUE, active.crit="SC")
#'
#' # Compare to an IFA model with adaptive Gibbs sampling
#' # simIFAadapt <- mcmc_IMIFA(coffee, method="IFA", n.iters=10000, active.crit="BD")
#' # resIFAadapt <- get_IMIFA_results(simIFAadapt)
#' # plot(resIFAnoadapt, "GQ")
#' # plot(resIFApostadapt, "GQ")
#' # plot(resIFAadapt, "GQ")}
get_IMIFA_results <- function(sims = NULL, burnin = 0L, thinning = 1L, G = NULL, Q = NULL, criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
adapt = FALSE, G.meth = c("mode", "median"), Q.meth = c("mode", "median"), conf.level = 0.95, error.metrics = TRUE,
vari.rot = FALSE, z.avgsim = FALSE, zlabels = NULL, nonempty = TRUE, ...) {
UseMethod("get_IMIFA_results")
}
#' @method get_IMIFA_results IMIFA
#' @export
get_IMIFA_results.IMIFA <- function(sims = NULL, burnin = 0L, thinning = 1L, G = NULL, Q = NULL, criterion = c("bicm", "aicm", "dic", "bic.mcmc", "aic.mcmc"),
adapt = FALSE, G.meth = c("mode", "median"), Q.meth = c("mode", "median"), conf.level = 0.95, error.metrics = TRUE,
vari.rot = FALSE, z.avgsim = FALSE, zlabels = NULL, nonempty = TRUE, ...) {
call <- match.call()
defopt <- options()
options(warn=1)
on.exit(suppressWarnings(options(defopt)), add=TRUE)
if(missing(sims)) stop("Simulations must be supplied", call.=FALSE)
if(!inherits(sims, "IMIFA")) stop("Object of class 'IMIFA' must be supplied", call.=FALSE)
burnin <- as.integer(burnin)
thinning <- as.integer(thinning)
store <- attr(sims, "Store")
if(any(c(length(thinning),
length(burnin)) > 1)) stop("'burnin' and 'thinning' must be of length 1", call.=FALSE)
if(any(burnin < 0, burnin >= store,
thinning < 1)) stop("Invalid 'burnin' and/or 'thinning' supplied", call.=FALSE)
store <- as.integer(seq(from=burnin + 1L, to=store, by=thinning))
if(length(store) < 10) stop("Too much 'burnin' or 'thinning' applied: not enough stored samples to proceed", call.=FALSE)
n.store <- length(store)
tmp.store <- store
dat <- attr(sims, "Dataset")
data.name <- attr(sims, "Name")
equal.pro <- attr(sims, "Equal.Pi")
label.switch <- attr(sims, "Label.Switch")
method <- attr(sims, "Method")
learn.alpha <- attr(sims, "Alph.step")
learn.d <- attr(sims, "Disc.step")
inf.G <- is.element(method, c("IMIFA", "OMIFA", "IMFA", "OMFA"))
inf.Q <- is.element(method, c("IMIFA", "OMIFA", "MIFA", "IFA"))
n.fac <- attr(sims, "Factors")
n.grp <- attr(sims, "Clusters")
n.obs <- attr(sims, "Obs")
n.var <- attr(sims, "Vars")
uni <- n.var == 1
sw <- tmpsw <- attr(sims, "Switch")
uni.meth <- attr(sims, "Uni.Meth")
uni.type <- unname(uni.meth["Uni.Type"])
conf.level <- as.numeric(conf.level)
obsnames <- attr(sims, "Obsnames")
varnames <- attr(sims, "Varnames")
adapted <- attr(sims, "Adapt")
adapt <- adapt && isFALSE(adapted)
cov.emp <- attr(sims, "Cov.Emp")
cov.range <- ifelse(uni, 1L, max(cov.emp) - min(cov.emp))
cshrink <- attr(sims, "C.Shrink")
if(any(length(conf.level) != 1,
!is.numeric(conf.level),
(conf.level <= 0 ||
conf.level >= 1))) stop("'conf.level' must be a single number in the interval (0, 1)", call.=FALSE)
conf.levels <- c((1 - conf.level)/2, (1 + conf.level)/2)
choice <- length(n.grp) * length(n.fac) > 1
crit.check <- !all(is.character(criterion))
if(isTRUE(crit.check)) stop("'criterion' must be a character vector of length 1", call.=FALSE)
criterion <- match.arg(criterion)
if(all(inf.Q, is.element(criterion,
c("aic.mcmc", "bic.mcmc")))) stop(paste0(ifelse(isTRUE(choice), "Model choice is", "Though model choice isn't"), " actually required -\n'criterion' cannot be 'aic.mcmc' or 'bic.mcmc' for the ", method, " method"), call.=FALSE)
if(length(error.metrics) != 1 ||
!is.logical(error.metrics)) stop("'error.metrics' must be a single logical indicator", call.=FALSE)
if(length(vari.rot) != 1 ||
!is.logical(vari.rot)) stop("'vari.rot' must be a single logical indicator", call.=FALSE)
if(isTRUE(error.metrics) &&
anyNA(cov.emp)) { warning("Certain error metrics cannot be computed as there are missing values in the empirical covariance matrix\n", call.=FALSE, immediate.=TRUE)
cov.met <- FALSE
} else cov.met <- error.metrics
if(length(z.avgsim) != 1 ||
!is.logical(z.avgsim)) stop("'z.avgsim' must be a single logical indicator", call.=FALSE)
if((z.avgsim <- (isTRUE(z.avgsim) &&
!is.element(method, c("FA", "IFA"))))) {
if(!(has_pkg <- suppressMessages(requireNamespace("mcclust", quietly=TRUE)))) {
z.avgsim <- FALSE; warning("Forcing 'z.avgsim' to FALSE: 'mcclust' package not installed\n", call.=FALSE, immediate.=TRUE)
} else if(n.obs >= 10000) message("Consider setting 'z.avgsim' to FALSE as the number of observations is quite large\n")
}
dots <- list(...)
dots <- dots[unique(names(dots))]
if(isTRUE(adapt)) {
active.crit <- ifelse(is.null(dots$active.crit), "BD", dots$active.crit)
epsilon <- ifelse(is.null(dots$eps), 1e-01, dots$eps)
prop <- ifelse(is.null(dots$prop), ifelse(n.var <= 2, 0.5,
switch(EXPR=active.crit,
BD=0.7, SC=0.99)), dots$prop)
prop <- switch(EXPR=active.crit, BD=floor(prop * n.var)/n.var, SC=prop)
forceQg <- ifelse(is.null(dots$forceQg), FALSE, dots$forceQg)
if(method != "IFA" &&
"SC" == active.crit) stop("'active.crit'=\"SC\" is only available for the \"IFA\" method", call.=FALSE)
if(!sw["l.sw"] &&
"BD" == active.crit) stop("Loadings must be stored for adaptation to take place", call.=FALSE)
if(any(!sw[c("l.sw",
"s.sw")]) &&
"SC" == active.crit) stop("Scores and loadings must be stored for adaption using the \"SC\" critertion to take place", call.=FALSE)
}
G.T <- !missing(G)
Q.T <- !missing(Q)
G.ind <- Q.ind <- 1L
GQs <- length(sims[[G.ind]])
GQ1x <- GQs > 1 && is.element(method, c("OMFA", "IMFA"))
if(inf.G) {
G.store2 <- lapply(seq_len(GQs), function(gq) sims[[G.ind]][[gq]]$G.store[store])
G.store <- matrix(unlist(G.store2), nrow=GQs, ncol=n.store, byrow=TRUE)
if(is.element(method, c("IMFA", "IMIFA"))) {
act.store <- lapply(seq_len(GQs), function(gq) sims[[G.ind]][[gq]]$act.store[store])
}
G.mX <- !all(is.character(G.meth))
if(isTRUE(G.mX)) stop("'G.meth' must be a character vector of length 1", call.=FALSE)
G.meth <- match.arg(G.meth)
G.tab <- if(GQ1x) lapply(apply(G.store, 1L, function(x) list(table(x, dnn=NULL))), "[[", 1L) else table(G.store, dnn=NULL)
G.prob <- if(GQ1x) lapply(G.tab, prop.table) else prop.table(G.tab)
G.mode <- if(GQ1x) unlist(lapply(G.tab, function(gt) as.numeric(names(gt[gt == max(gt)])[1L]))) else as.numeric(names(G.tab[G.tab == max(G.tab)])[1L])
G.med <- if(GQ1x) ceiling(matrixStats::rowMedians(G.store,
useNames=FALSE) * 2)/2 else ceiling(Median(G.store) * 2)/2
if(!G.T) {
G <- switch(EXPR=G.meth, mode=G.mode, floor(G.med))
}
G.CI <- if(GQ1x) round(rowQuantiles(G.store, probs=conf.levels)) else round(stats::quantile(G.store, conf.levels))
}
if(G.T) {
G <- as.integer(G)
if(any(length(G) != 1,
!is.integer(G))) stop("'G' must be an integer of length 1", call.=FALSE)
if(!inf.G) {
if(!is.element(method, c("FA", "IFA"))) {
if(!is.element(G, n.grp)) stop("This 'G' value was not used during simulation", call.=FALSE)
G.ind <- which(n.grp == G)
} else if(G > 1) message(paste0("Forced G=1 for the ", method, " method\n"))
} else {
if(GQ1x) {
if(!Q.T) stop(paste0("'G' cannot be supplied without 'Q' for the ", method, " method if a range of Q values were explored"), call.=FALSE)
tmpQ <- which(n.fac == unique(Q))
} else {
tmpQ <- Q.ind
}
if(length(tmpQ > 0) && !is.element(G,
unique(G.store[tmpQ,]))) stop("This 'G' value was not visited during simulation", call.=FALSE)
}
}
G <- if(any(inf.G, all(G.T, !is.element(method, c("FA", "IFA"))))) G else 1L
if(Q.T) {
Q <- as.integer(Q)
if(!is.integer(Q)) stop("'Q' must of integer type", call.=FALSE)
if(G.T) {
if(length(Q) == 1) Q <- rep(Q, G)
if(length(Q) != G) stop(paste0("'Q' must be supplied for each cluster, as a scalar or vector of length G=", G), call.=FALSE)
} else if(length(n.grp) != 1 && all(!is.element(length(Q),
c(1, n.grp)))) stop("'Q' must be a scalar if G=1, 'G' is not suppplied, or a range of G values were explored", call.=FALSE)
if(all(is.element(method, c("FA", "MFA", "OMFA", "IMFA")))) {
if(length(unique(Q)) != 1) stop(paste0("'Q' cannot vary across clusters for the ", method, " method"), call.=FALSE)
Q <- unique(Q)
if(!is.element(Q, n.fac)) stop("This 'Q' value was not used during simulation", call.=FALSE)
Q.ind <- which(n.fac == Q)
}
if(inf.Q) {
if(any((Q != 0) + (Q *
(n.var - Q) <= 0) > 1)) stop(paste0("'Q' must be less than the number of variables ", n.var), call.=FALSE)
Qtmp <- if(inf.G) sims[[1L]][[1L]]$Q.store[seq_len(G),store, drop=FALSE] else switch(EXPR=method, MIFA=sims[[if(G.T) G == n.grp else G.ind]][[1L]]$Q.store[,store, drop=FALSE], sims[[1L]][[1L]]$Q.store[,store, drop=FALSE])
storage.mode(Qtmp) <- switch(EXPR=method, IFA="integer", "numeric")
Qtmp <- switch(EXPR=method, IFA=max(Qtmp), Rfast::rowMaxs(Qtmp, value=TRUE))
Qlen <- length(Qtmp)
if(length(Q) == 1) Q <- rep(Q, Qlen)
if(length(Q) != Qlen) stop(paste0("'Q' must be supplied for each cluster, as a scalar or vector of length G=", G), call.=FALSE)
if(any(Q * (Qtmp - Q) < 0)) stop(paste0("'Q' can't be greater than the maximum number of factors stored in ", ifelse(method == "IFA", "", "any cluster of "), deparse(match.call()$sims)), call.=FALSE)
}
}
if(inf.G) {
tmp.store <- if(GQ1x) lapply(seq_len(GQs), function(gq) store[G.store[gq,] == G[ifelse(G.T, 1L, gq)]]) else store[G.store == G]
GQ.temp1 <- list(G = G, G.Mode = G.mode, G.Median = G.med, G.CI = G.CI, G.Probs = G.prob, G.Counts = G.tab)
GQ.temp1 <- c(GQ.temp1, list(Stored.G = switch(EXPR=method, OMIFA=provideDimnames(G.store, base=list("Non-Empty", ""), unique=FALSE),
IMIFA=provideDimnames(do.call(rbind, c(G.store2, act.store)), base=list(c("Non-Empty", "Active"), ""), unique=FALSE),
OMFA=lapply(seq_len(GQs), function(g) provideDimnames(t(G.store[g,]), base=list("Non-Empty", ""), unique=FALSE)),
IMFA=lapply(seq_len(GQs), function(g) provideDimnames(rbind(G.store2[[g]], act.store[[g]]), base=list(c("Non-Empty", "Active"), ""), unique=FALSE)))))
GQ.temp1 <- c(GQ.temp1, list(G.Last = switch(EXPR=method, OMIFA=, IMIFA=GQ.temp1$Stored.G[1L,n.store],
lapply(seq_len(GQs), function(g) GQ.temp1$Stored.G[[g]][1L,n.store]))))
}
G.range <- ifelse(G.T, 1L, length(n.grp))
Q.range <- ifelse(any(Q.T, all(!is.element(method, c("OMFA", "IMFA")), inf.Q)), 1L, length(n.fac))
crit.mat <- matrix(NA, nrow=G.range, ncol=Q.range)
# Retrieve log-likelihoods and/or tune G &/or Q according to criterion
if(all(G.T, Q.T)) {
dimnames(crit.mat) <- list(paste0("G", G), if(inf.Q) "IFA" else paste0("Q", Q))
} else if(G.T) {
dimnames(crit.mat) <- list(paste0("G", G), if(inf.Q) "IFA" else paste0("Q", n.fac))
} else if(Q.T) {
dimnames(crit.mat) <- list(paste0("G", n.grp), if(inf.Q) "IFA" else paste0("Q", Q))
} else {
dimnames(crit.mat) <- list(paste0("G", n.grp), if(inf.Q) "IFA" else paste0("Q", n.fac))
}
rownames(crit.mat) <- switch(EXPR=method, IMFA=, IMIFA="IM", OMFA=, OMIFA="OM", rownames(crit.mat))
aicm <- bicm <- aicm.sd <- dic <-
aic.mcmc <- bic.mcmc <- bicm.sd <- crit.mat
log.N <- log(n.obs)
if(!inf.Q) {
Ks <- matrix(NA, nrow=G.range, ncol=Q.range)
}
for(g in seq_len(G.range)) {
gi <- ifelse(G.T, G.ind, g)
for(q in seq_len(Q.range)) {
qi <- ifelse(Q.T, Q.ind, q)
log.likes <- if(GQ1x) sims[[gi]][[qi]]$ll.store[tmp.store[[qi]]] else sims[[gi]][[qi]]$ll.store[tmp.store]
log.likes <- log.likes[stats::complete.cases(log.likes)]
S2 <- ifelse(length(log.likes) != 1, Var(log.likes), 0L)
llbar <- mean(log.likes)
d.hat <- 2 * S2
llmax2 <- 2 * llbar + d.hat
aicm[g,q] <- 4 * llbar - llmax2
bicm[g,q] <- llmax2 - d.hat * log.N
dic[g,q] <- 2 * (3 * llbar - llmax2)
ci.tmp <- S2 * (11 * S2 + 24)
aicm.sd[g,q] <- sqrt((4 * (S2 + ci.tmp))/length(log.likes))
bicm.sd[g,q] <- sqrt((2 * d.hat + (log.N - 1)^2 * ci.tmp)/length(log.likes))
if(!inf.Q) {
Ks[g,q] <- K <- switch(EXPR=method, OMFA=, IMFA=PGMM_dfree(Q=n.fac[qi], P=n.var, G=G[ifelse(G.T, 1L, qi)],
method=switch(EXPR=uni.type, unconstrained="UUU", isotropic="UUC", constrained="UCU", single="UCC"),
equal.pro=equal.pro), attr(sims[[gi]][[qi]], "K"))
llmax2k <- 2 * max(log.likes)
aic.mcmc[g,q] <- llmax2k - K * 2
bic.mcmc[g,q] <- llmax2k - K * log.N
}
}
}
crit <- get(criterion)
crit.max <- which(crit == max(crit, na.rm=TRUE), arr.ind=TRUE)
if(!inf.Q &&
nrow(crit.max) > 1) { warning(paste0("Ties for the optimal model exist according to the '", criterion, "' criterion: choosing the most parsimonious model\n"), call.=FALSE, immediate.=TRUE)
crit.max <- which(Ks == min(Ks[crit.max], na.rm=TRUE), arr.ind=TRUE)
crit.max <- crit.max[which.min(crit.max[,1L]),]
}
# Control for supplied values of G &/or Q
if(!any(Q.T, G.T)) {
G.ind <- crit.max[1L]
Q.ind <- crit.max[2L]
if(!inf.G) {
G <- n.grp[G.ind]
}
if(!inf.Q) {
Q <- n.fac[Q.ind]
}
} else if(all(G.T, !Q.T)) {
Q.ind <- which.max(crit)
if(!inf.Q) {
Q <- n.fac[Q.ind]
}
} else if(all(Q.T, !G.T)) {
G.ind <- which.max(crit)
if(!inf.G) {
G <- n.grp[G.ind]
}
}
G <- ifelse(inf.G, ifelse(G.T, G, G[Q.ind]), ifelse(length(n.grp) == 1, n.grp, G))
Gseq <- seq_len(G)
gnames <- paste0("Cluster", Gseq)
G.ind <- if(!all(inf.G, length(n.grp) > 1)) G.ind else which(n.grp == G)
GQ.temp2 <- list(AICMs = aicm, BICMs = bicm, DICs = dic)
clust.ind <- !any(is.element(method, c("FA", "IFA")),
all(is.element(method, c("MFA", "MIFA")), G == 1))
sw.mx <- (!clust.ind || sw["mu.sw"]) && sw["u.sw"]
sw.px <- !clust.ind || sw["psi.sw"]
if(!inf.Q) {
Q <- if(length(n.fac) > 1) Q else n.fac
Q.ind <- if(all(!Q.T, length(n.fac) > 1)) Q.ind else which(n.fac == Q)
Q <- stats::setNames(if(length(Q) != G) rep(Q, G) else Q, gnames)
if(all(inf.G, Q.T)) GQ.temp1$G <- rep(G, GQs)
if(GQ1x) {
GQ.temp1$G.CI <- lapply(seq_len(GQs), function(gq) GQ.temp1$G.CI[gq,])
GQ.temp1 <- lapply(GQ.temp1, "[[", Q.ind)
} else if(inf.G) {
GQ.temp1$Stored.G <- GQ.temp1$Stored.G[[1L]]
GQ.temp1$G.Last <- GQ.temp1$G.Last[[1L]]
}
GQ.temp3 <- c(GQ.temp2, list(AIC.mcmcs = aic.mcmc, BIC.mcmcs = bic.mcmc))
GQ.res <- switch(EXPR=method, OMFA=, IMFA=c(GQ.temp1, list(Q = Q), list(Criteria = GQ.temp3)), c(list(G = G, Q = Q), list(Criteria = GQ.temp3)))
if(GQ1x) {
tmp.store <- tmp.store[[Q.ind]]
}
}
if(clust.ind) {
zadj <- sims[[G.ind]][[Q.ind]]$z.store
z <- as.matrix(zadj[tmp.store,])
}
if(condition <- all(isTRUE(nonempty), is.element(method, c("MFA", "MIFA")), G > 1)) {
emptystore <- which(rowAlls(rowTabulate(sims[[G.ind]][[Q.ind]]$z.store[tmp.store,], max_number=n.grp[G.ind]) > 0, useNames=FALSE))
if(!identical(tmp.store,
emptystore)) warning(paste0("Discarding iterations with fewer than G=", G, " non-empty components\n"), call.=FALSE, immediate.=TRUE)
tmp.store <- emptystore
z <- z[tmp.store,, drop=FALSE]
}
storeG <- seq_along(tmp.store)
TN.store <- length(tmp.store)
if(TN.store < 2) stop(paste0("Not enough samples stored to proceed", ifelse(condition && missing(nonempty), ": try supplying 'nonempty'=FALSE", ifelse(choice && any(G.T, Q.T), paste0(": try supplying different G or Q values"), ""))), call.=FALSE)
if(inf.Q) {
Q.store <- sims[[G.ind]][[Q.ind]]$Q.store[,tmp.store, drop=FALSE]
Q.mX <- !all(is.character(Q.meth))
if(isTRUE(Q.mX)) stop("'Q.meth' must be a character vector of length 1", call.=FALSE)
Q.meth <- match.arg(Q.meth)
if(cshrink <- cshrink && ((method != "MIFA") || G > 1)) {
sigmas <- sims[[G.ind]][[Q.ind]]$sigma[,tmp.store, drop=FALSE]
}
}
# Manage Label Switching & retrieve cluster labels/mixing proportions
if(clust.ind) {
label.miss <- missing(zlabels)
if(!label.miss && (all(G > 1, !is.factor(zlabels), !is.numeric(zlabels)) ||
length(zlabels) != n.obs)) stop(paste0("'zlabels' must be a factor of length N=", n.obs), call.=FALSE)
if(sw["mu.sw"]) {
mus <- sims[[G.ind]][[Q.ind]]$mu[,,tmp.store, drop=FALSE]
}
if(sw["l.sw"]) {
lmats <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$load)[,,,tmp.store, drop=FALSE] else sims[[G.ind]][[Q.ind]]$load[,,,tmp.store, drop=FALSE]
}
if(sw["psi.sw"]) {
psis <- sims[[G.ind]][[Q.ind]]$psi[,,tmp.store, drop=FALSE]
}
if(sw["pi.sw"]) {
pies <- sims[[G.ind]][[Q.ind]]$pi.prop[,tmp.store, drop=FALSE]
}
if(sw["s.sw"]) {
z2 <- z
}
zadj <- zadj[store,]
if(!label.switch) {
z.temp <- tryCatch(factor(z[1L,], labels=Gseq), error=function(e) factor(z[1L,], levels=Gseq))
for(sl in storeG) {
sw.lab <- .lab_switch(z.new=z[sl,], z.old=z.temp)
z.perm <- sw.lab$z.perm
left <- as.integer(unname(z.perm))
right <- as.integer(names(z.perm))
if(!identical(left, right)) {
z[sl,] <- sw.lab$z
if(sw["mu.sw"]) {
mus[,left,sl] <- mus[,right,sl]
}
if(sw["l.sw"]) {
lmats[,,left,sl] <- lmats[,,right,sl]
}
if(sw["psi.sw"]) {
psis[,left,sl] <- psis[,right,sl]
}
if(sw["pi.sw"]) {
pies[left,sl] <- pies[right,sl]
}
if(inf.Q) {
Q.store[left,sl] <- Q.store[right,sl]
if(cshrink) {
sigmas[left,sl] <- sigmas[right,sl]
}
}
}
}
}
if(inf.Q) {
Q.store <- provideDimnames(Q.store[Gseq,, drop=FALSE], base=list(gnames, ""), unique=FALSE)
}
if(sw["mu.sw"]) mus <- tryCatch(mus[,Gseq,, drop=FALSE], error=function(e) mus)
if(sw["l.sw"]) lmats <- tryCatch(lmats[,,Gseq,, drop=FALSE], error=function(e) lmats)
if(sw["psi.sw"]) psis <- tryCatch(psis[,Gseq,, drop=FALSE], error=function(e) psis)
if(cshrink) sigmas <- tryCatch(sigmas[Gseq,, drop=FALSE], error=function(e) sigmas)
MAP <- if(G > 1) apply(z, 2L, function(x) factor(which.max(tabulate(x, nbins=G)), levels=Gseq)) else factor(rep(1L, n.obs))
post.prob <- if(G > 1) matrix(colTabulate(z, max_number=G)/TN.store, nrow=n.obs, ncol=G, byrow=TRUE) else matrix(1L, nrow=n.obs)
if(isTRUE(z.avgsim)) {
znew <- try(Zsimilarity(zs=zadj), silent=TRUE)
condit <- all(!is.element(method, c("MIFA", "MFA")), inherits(znew, "try-error"))
if(isTRUE(condit)) {
znew <- try(Zsimilarity(zs=z), silent=TRUE)
warning("Constructing the similarity matrix failed:\ntrying again using iterations corresponding to the modal number of clusters\n", call.=FALSE, immediate.=TRUE)
}
if(!inherits(znew, "try-error")) {
zadj <- znew$z.avg
zadj <- factor(zadj, labels=seq_along(unique(zadj)))
zadj <- as.integer(levels(zadj))[zadj]
zavg <- znew$z.sim
if(!label.miss) {
zadj <- .lab_switch(z.new=zadj, z.old=factor(zlabels, labels=seq_along(unique(zlabels))))$z
tab <- table(zadj, zlabels, dnn=list("Predicted", "Observed"))
tabstat <- c(.class_agreement(tab), classError(classification=MAP, class=zlabels))
if(nrow(tab) != ncol(tab)) {
tabstat <- tabstat[-seq_len(2L)]
names(tabstat)[4L] <- "error.rate"
} else {
names(tabstat)[6L] <- "error.rate"
}
if(tabstat$error.rate == 0) {
tabstat$misclassified <- NULL
}
tabstat <- c(list(confusion.matrix = tab), tabstat)
class(tabstat) <- "listof"
}
z_simavg <- list(z.avg = zadj, z.sim = zavg)
z_simavg <- c(z_simavg, if(!label.miss) list(avgsim.perf = tabstat))
attr(z_simavg, "Conditional") <- condit
class(z_simavg) <- "listof"
} else { warning("Can't compute similarity matrix or 'average' clustering: forcing 'z.avgsim' to FALSE\n", call.=FALSE, immediate.=TRUE)
z.avgsim <- FALSE
}
}
uncertain <- if(G > 1) 1 - Rfast::rowMaxs(post.prob, value=TRUE) else integer(n.obs)
sizes <- stats::setNames(tabulate(MAP, nbins=G), gnames)
if(any(sizes == 0)) warning(paste0("Empty cluster exists in MAP partition:\nexamine trace plots", ifelse(any(is.element(method, c("OMFA", "IMFA", "OMIFA", "IMIFA")), is.element(method, c("MFA", "MIFA")) && any(n.grp < G)), ", try to supply a lower G value to get_IMIFA_results(),", ""), " or re-run the model\n"), call.=FALSE, immediate.=TRUE)
if(sw["pi.sw"]) {
pi.prop <- provideDimnames(pies[Gseq,storeG, drop=FALSE], base=list(gnames, ""), unique=FALSE)
pi.prop <- if(inf.G) sweep(pi.prop, 2L, colSums2(pi.prop, useNames=FALSE), FUN="/", check.margin=FALSE) else pi.prop
var.pi <- stats::setNames(rowVars(pi.prop), gnames)
ci.pi <- rowQuantiles(pi.prop, probs=conf.levels)
ci.pi <- if(G == 1) t(ci.pi) else ci.pi
post.pi <- stats::setNames(rowMeans2(pi.prop,
useNames=FALSE), gnames)
} else if(equal.pro) {
post.pi <- stats::setNames(rep(1/G, G), gnames)
} else { message("Mixing proportions not stored: estimating posterior mean by the proportions of the MAP clustering\n")
post.pi <- stats::setNames(sizes/n.obs, gnames)
}
tab.z <- if((!label.miss && G > 1) || (clust.ind && all(error.metrics, sw["mu.sw"], sw["psi.sw"]))) rowTabulate(z, max_number=G)
if(!label.miss && G > 1) {
sw.lab <- .lab_switch(z.new=MAP, z.old=factor(zlabels, labels=seq_along(unique(zlabels))))
MAP <- factor(factor(sw.lab$z, labels=which(sizes > 0)), levels=Gseq)
l.perm <- sw.lab$z.perm
left <- as.integer(l.perm[Gseq])
z.tab <- tab.z > 0
z.tmp <- lapply(storeG, function(i) factor(z[i,], labels=left[z.tab[i,]]))
z <- do.call(rbind, lapply(z.tmp, function(x) as.integer(levels(x))[as.integer(x)]))
index <- order(left)
tab.z <- tab.z[,index, drop=FALSE]
sizes <- stats::setNames(sizes[index], gnames)
post.prob <- post.prob[,index]
if(sw["mu.sw"]) mus <- mus[,index,, drop=FALSE]
if(sw["l.sw"]) lmats <- lmats[,,index,, drop=FALSE]
if(sw["psi.sw"]) psis <- psis[,index,, drop=FALSE]
post.pi <- stats::setNames(post.pi[index], gnames)
if(sw["pi.sw"]) {
pi.prop <- provideDimnames(unname(pi.prop[index,, drop=FALSE]), base=list(gnames, ""), unique=FALSE)
var.pi <- stats::setNames(var.pi[index], gnames)
ci.pi <- provideDimnames(unname(ci.pi[index,, drop=FALSE]), base=list(gnames, colnames(ci.pi)))
}
if(inf.Q) {
Q.store <- provideDimnames(unname(Q.store[index,, drop=FALSE]), base=list(gnames, ""), unique=FALSE)
if(cshrink) sigmas <- sigmas[index,, drop=FALSE]
}
}
if(!label.miss) {
tab <- table(MAP, zlabels, dnn=list("Predicted", "Observed"))
tab.stat <- c(.class_agreement(tab), classError(classification=MAP, class=zlabels))
if(nrow(tab) != ncol(tab)) {
tab.stat <- tab.stat[-seq_len(2L)]
names(tab.stat)[4L] <- "error.rate"
} else {
names(tab.stat)[6L] <- "error.rate"
}
if(tab.stat$error.rate == 0) {
tab.stat$misclassified <- NULL
}
tab.stat <- c(list(confusion.matrix = tab), tab.stat)
class(tab.stat) <- "listof"
}
if(learn.alpha) {
alpha <- sims[[G.ind]][[Q.ind]]$alpha[store]
post.alpha <- mean(alpha)
var.alpha <- Var(alpha)
ci.alpha <- stats::quantile(alpha, conf.levels)
rate <- sims[[G.ind]][[Q.ind]]$a.rate
DP.alpha <- list(alpha = alpha, post.alpha = post.alpha, var.alpha = var.alpha, ci.alpha = ci.alpha, alpha.rate = rate, last.alpha = alpha[n.store])
DP.alpha <- c(DP.alpha, if(isTRUE(attr(sims, "TuneZeta"))) list(avg.zeta = sims[[G.ind]][[Q.ind]]$avg.zeta))
class(DP.alpha) <- "listof"
}
if(learn.d) {
discount <- drop(sims[[G.ind]][[Q.ind]]$discount[store])
post.disc <- mean(discount)
post.kappa <- sum(discount == 0)/n.store
var.disc <- Var(discount)
ci.disc <- stats::quantile(discount, conf.levels)
rate <- sims[[G.ind]][[Q.ind]]$d.rate
post.dzero <- post.disc/(1 - post.kappa)
discount <- if(sum(discount == 0)/n.store > 0.5) as.simple_triplet_matrix(discount) else discount
PY.disc <- list(discount = discount, post.disc = post.disc, post.kappa = post.kappa, var.disc = var.disc,
ci.disc = ci.disc, disc.rate = rate, last.disc = discount[n.store], post.d_nonzero = post.dzero)
class(PY.disc) <- "listof"
}
MAP <- as.integer(levels(MAP))[MAP]
uncert.obs <- which(uncertain >= 1/G)
uncertain <- if(sum(uncertain == 0)/n.obs > 0.5) as.simple_triplet_matrix(uncertain) else uncertain
attr(uncertain, "Obs") <- if(sum(uncert.obs) != 0) uncert.obs
if(!label.miss) tab.stat$uncertain <- attr(uncertain, "Obs")
cluster <- list(MAP = MAP, z = z, uncertainty = uncertain, last.z = z[TN.store,])
cluster <- c(cluster, list(post.sizes = sizes, post.ratio = sizes/n.obs, post.pi = post.pi/sum(post.pi),
post.prob = provideDimnames(post.prob, base=list("", gnames)),
PCM = provideDimnames(post_conf_mat(post.prob), base=list(gnames, gnames))),
if(sw["pi.sw"]) list(pi.prop = pi.prop, var.pi = var.pi, ci.pi = ci.pi, last.pi = pi.prop[,TN.store]),
if(!label.miss) list(perf = tab.stat),
if(learn.alpha) list(Alpha = DP.alpha),
if(learn.d) list(Discount = PY.disc),
if(z.avgsim) list(Z.avgsim = z_simavg),
if(is.element(method, c("IMFA", "IMIFA"))) list(lab.rate = sims[[G.ind]][[Q.ind]]$lab.rate))
attr(cluster, "Z.init") <- attr(sims[[G.ind]], "Z.init")
attr(cluster, "Init.Meth") <- attr(sims, "Init.Z")
attr(cluster, "Label.Sup") <- !label.miss
class(cluster) <- "listof"
z.ind <- lapply(Gseq, function(g) MAP == g)
} else {
if(sw["l.sw"]) {
lmats <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$load) else sims[[G.ind]][[Q.ind]]$load
}
z.ind <- list(seq_len(n.obs))
sizes <- n.obs
}
if(inf.Q) {
if(G1 <- G > 1) {
if(isTRUE(adapt)) {
colvec <- lapply(Gseq, function(g) vapply(seq_len(TN.store), function(i) colSums2(abs(.a_drop(lmats[,,g,i, drop=FALSE], drop=3L:4L)) < epsilon, useNames=FALSE)/n.var < prop, logical(ncol(lmats))))
Q.store[] <- t(vapply(colvec, colSums2, numeric(TN.store), useNames=FALSE))
if(isTRUE(forceQg)) {
Q.force <- t(vapply(Gseq, function(g) pmin(Q.store[g,], sizes[g] - 1), numeric(TN.store)))
Q.big <- any(Q.force < Q.store)
Q.store[] <- Q.force
colvec <- lapply(Gseq, function(g) vapply(seq_len(TN.store), function(i) { colvec[[g]][colvec[[g]][,i],i][-seq_len(Q.store[g,i])] <- FALSE; colvec[[g]][,i] }, logical(ncol(lmats))))
}
}
Q.tab <- lapply(apply(Q.store, 1L, function(x) list(table(x, dnn=NULL))), "[[", 1L)
Q.prob <- lapply(Q.tab, prop.table)
class(Q.tab) <- class(Q.prob) <- "listof"
} else {
if(isTRUE(adapt)) {
if(active.crit == "SC") {
eta.tmp <- as.array(sims[[G.ind]][[Q.ind]]$eta)
SC <- lapply(seq_len(TN.store), function(i) .SC_crit(dat, .a_drop(eta.tmp[,,i,drop=FALSE], 3L), .a_drop(lmats[,,i, drop=FALSE], 3L), prop))
colvec <- vapply(lapply(SC, "[[", "nonred"), function(x) match(seq_len(ncol(lmats)), x, nomatch=0) > 0, logical(ncol(lmats)))
Q.store[] <- Q.store[] - vapply(SC, "[[", numeric(1L), "numred")
} else {
colvec <- vapply(seq_len(TN.store), function(i) colSums2(abs(.a_drop(lmats[,,i, drop=FALSE], 3L)) < epsilon, useNames=FALSE)/n.var < prop, logical(ncol(lmats)))
Q.store[] <- colSums2(colvec, useNames=FALSE)
}
colvec <- list(colvec)
}
Q.tab <- table(Q.store, dnn=NULL)
Q.prob <- prop.table(Q.tab)
}
Q.mode <- if(G1) unlist(lapply(Q.tab, function(qt) as.numeric(names(qt[qt == max(qt)])[1L]))) else as.numeric(names(Q.tab[Q.tab == max(Q.tab)])[1L])
Q.med <- if(G1) stats::setNames(ceiling(matrixStats::rowMedians(Q.store,
useNames=FALSE) * 2)/2, gnames) else ceiling(Median(Q.store) * 2)/2
if(!Q.T) {
Q <- switch(EXPR=Q.meth, mode=Q.mode, floor(Q.med))
} else {
Q <- if(G.T) Q else stats::setNames(if(length(Q) == G) Q else rep(Q, G), gnames)
}
Q.CI <- if(G1) round(rowQuantiles(Q.store,
probs=conf.levels)) else round(stats::quantile(Q.store, conf.levels))
GQ.temp4 <- list(Q = Q, Q.Mode = Q.mode, Q.Median = Q.med,
Q.CI = Q.CI, Q.Probs = Q.prob, Q.Counts = Q.tab,
Stored.Q = if(clust.ind) Q.store else drop(Q.store),
Q.Last = Q.store[,TN.store])
GQ.res <- if(inf.G) c(GQ.temp1, GQ.temp4) else c(list(G = G), GQ.temp4)
GQ.res <- if(cshrink) c(GQ.res, list(post.sigma = stats::setNames(rowMeans2(sigmas, useNames=FALSE), gnames))) else GQ.res
GQ.res <- c(GQ.res, list(Criteria = GQ.temp2))
attr(GQ.res, "Q.big") <- ifelse(adapt && G1 && forceQg, Q.big, attr(sims[[G.ind]][[Q.ind]], "Q.big"))
}
Q0 <- Q == 0
Qmax <- max(Q)
if(all(isTRUE(choice), is.element(criterion, c("aicm", "bicm", "dic")))) {
if(all(!G.T, !is.element(method,
c("FA", "IFA")), G == 1)) warning(paste0("Chosen model has only one group:\nNote that the ", criterion, " criterion may exhibit bias toward one-group models\n"), call.=FALSE, immediate.=TRUE)
if(all(!Q.T, method == "MIFA")) {
if(any(Q0)) warning(paste0("Chosen model has ", ifelse(sum(Q0) == G, "zero factors", "a cluster with zero factors"), ":\nNote that the ", criterion, " criterion may exhibit bias toward models ", ifelse(sum(Q0) == G, "with zero factors\n", "where some clusters have zero factors\n")), call.=FALSE, immediate.=TRUE)
} else if(all(Q0)) warning(paste0("Chosen model has zero factors:\nNote that the ", criterion, " criterion may exhibit bias toward zero-factor models\n"), call.=FALSE, immediate.=TRUE)
}
leder.b <- min(n.obs - 1L, Ledermann(n.var, isotropic=is.element(uni.type, c("isotropic", "single"))))
if(any(Q > leder.b)) { warning(paste0("Estimate of Q", ifelse(G > 1, " in one or more clusters ", " "), "is greater than ", ifelse(any(unlist(Q) > n.var), paste0("the number of variables (", n.var, ")"), paste0("the suggested Ledermann upper bound (", leder.b, ")")), "\n"), call.=FALSE, immediate.=TRUE)
} else if(any(floor((n.var - 1)/2) <
Q)) warning("It is advisable that the restriction Q <= floor((P - 1)/2) be respected\n", call.=FALSE, immediate.=TRUE)
if(any(Q >= sizes) &&
G > 1 &&
!attr(sims, "ForceQg")) warning("Estimate of Q is greater than the number of observations in one or more clusters:\nConsider setting 'forceQg' to TRUE\n", call.=FALSE, immediate.=TRUE)
# Retrieve (unrotated) scores
if(no.score <- all(Q0)) {
if(sw["s.sw"]) message("Scores & loadings not stored as model has zero factors\n")
sw["s.sw"] <- FALSE
}
if(inf.Q) {
Lstore <- lapply(Gseq, function(g) storeG[Q.store[g,] >= Q[g]])
if(any(lengths(Lstore) < 2)) {
sw["s.sw"] <- tmpsw["l.sw"] <-
sw["l.sw"] <- FALSE; warning("Forcing non-storage of scores and loadings due to shortage of retained samples\n", call.=FALSE, immediate.=TRUE)
}
eta.store <- storeG[colMaxs(Q.store, value=TRUE) >= Qmax]
} else {
eta.store <- tmp.store
Lstore <- storeG
}
Qseq <- seq_len(Qmax)
store.e <- (if(inf.Q) storeG else tmp.store) %in% eta.store
e.store <- sum(store.e)
QsE <- if(inf.Q) Q.store[,store.e, drop=!clust.ind]
Q0E <- if(inf.Q) QsE == 0 else if(clust.ind) matrix(Q == 0, nrow=G, ncol=e.store) else rep(Q == 0, e.store)
Q0X <- all(Q0E)
QX0 <- any(!Q0E)
frobenius <- error.metrics && all(sw["psi.sw"], sw["mu.sw"] || !sw["u.sw"], any(Q0X, sw["l.sw"]))
if(sw["s.sw"]) {
eta <- if(inf.Q) as.array(sims[[G.ind]][[Q.ind]]$eta)[,,eta.store, drop=FALSE] else sims[[G.ind]][[Q.ind]]$eta[,,eta.store, drop=FALSE]
if(!sw["l.sw"]) warning("Caution advised when examining posterior factor scores: Procrustes rotation has not taken place because loadings weren't stored\n", call.=FALSE, immediate.=TRUE)
}
# Loop over g in G to extract other results
result <- list(list())
for(g in Gseq) {
Qg <- Q[g]
Q0g <- Q0[g]
Qgs <- seq_len(Qg)
sw["l.sw"] <- tmpsw["l.sw"]
if(Q0g) {
if(all(sw["l.sw"],
!no.score)) message(paste0("Loadings ", ifelse(G > 1, paste0("for cluster ", g, " not stored as it"), " not stored as model"), " has zero factors\n"))
sw["l.sw"] <- FALSE
}
# Retrieve loadings and rotate
if(sw["l.sw"]) {
tmpind <- ifelse(inf.Q, which.max(Q.store[g,] == Qg), 1L)
lmat <- if(clust.ind) .a_drop(lmats[,,g,storeG, drop=FALSE], drop=3L) else lmats[,,storeG, drop=FALSE]
l.temp <- .a_drop(lmat[,,tmpind, drop=FALSE], drop=3L)
l.temp <- if(isTRUE(vari.rot)) .vari_max(l.temp, ...)$loadings else l.temp
if(inf.Q) {
rotind <- Lstore[[g]][-1L]
rotind <- rotind[Q.store[g,rotind] > 1]
} else {
rotind <- if(Qg > 1) Lstore[-1L]
}
for(p in rotind) {
proc <- suppressWarnings(Procrustes(X=if(uni) t(lmat[,,p]) else as.matrix(lmat[,,p]), Xstar=l.temp))
lmat[,,p] <- proc$X.new
if(sw["s.sw"] &&
(p %in%
eta.store)) {
p2 <- if(inf.Q) eta.store == p else p
if(clust.ind) {
zp <- z2[p,] == g
eta[zp,,p2] <- .a_drop(eta[zp,,p2, drop=FALSE], 3L) %*% proc$R
} else {
eta[,,p2] <- .a_drop(eta[,,p2, drop=FALSE], 3L) %*% proc$R
}
}
}
if(adapt) {
lmat <- sapply(lapply(storeG, function(i) .a_drop(lmat[,colvec[[g]][,i],i, drop=FALSE], drop=3L)),
function(x) cbind(x, matrix(0L, nrow=n.var, ncol=ncol(lmat) - ncol(x))), simplify="array")
}
}
# Retrieve means, uniquenesses & empirical covariance matrix
if(clust.ind) {
if(sw["mu.sw"]) {
mu <- if(uni) t(mus[,g,storeG]) else as.matrix(mus[,g,storeG])
}
if(sw["psi.sw"]) {
psi <- if(uni) t(psis[,g,storeG]) else as.matrix(psis[,g,storeG])
}
} else {
post.mu <- as.matrix(sims[[G.ind]][[Q.ind]]$post.mu)
post.psi <- as.matrix(sims[[G.ind]][[Q.ind]]$post.psi)
if(sw["mu.sw"]) {
mu <- sims[[G.ind]][[Q.ind]]$mu[,storeG, drop=FALSE]
}
if(sw["psi.sw"]) {
psi <- sims[[G.ind]][[Q.ind]]$psi[,storeG, drop=FALSE]
}
}
if(sw["mu.sw"] &&
is.null(rownames(mu))) {
rownames(mu) <- varnames
}
if(sw["psi.sw"] &&
is.null(rownames(psi))) {
rownames(psi) <- varnames
}
if(sw["l.sw"] &&
is.null(rownames(lmat))) {
rownames(lmat) <- varnames
}
# Compute posterior means and % variation explained
if(sw["mu.sw"]) {
post.mu <- if(clust.ind) rowMeans2(mu,
useNames=FALSE) else post.mu
var.mu <- if(uni) Var(mu) else rowVars(mu)
ci.tmp <- rowQuantiles(mu, probs=conf.levels)
ci.mu <- if(uni) t(ci.tmp) else ci.tmp
}
if(sw["psi.sw"]) {
post.psi <- if(clust.ind) rowMeans2(psi,
useNames=FALSE) else post.psi
var.psi <- if(uni) Var(psi) else rowVars(psi)
ci.tmp <- rowQuantiles(psi, probs=conf.levels)
ci.psi <- if(uni) t(ci.tmp) else ci.tmp
}
if(sw["l.sw"]) {
lmat <- lmat[,Qgs,if(inf.Q) Lstore[[g]] else storeG, drop=FALSE]
post.load <- rowMeans(lmat, dims=2)
var.load <- apply(lmat, c(1L, 2L), Var)
ci.load <- apply(lmat, c(1L, 2L), stats::quantile, conf.levels)
var.exp <- rowSums2(post.load^2, useNames=FALSE)
class(post.load) <- "loadings"
} else if(sw["psi.sw"]) {
if(sizes[g] > 1) {
dat.gg <- dat[z.ind[[g]],, drop=FALSE]
cov.gg <- stats::cov(dat.gg)
}
var.exp <- ifelse(sizes[g] <= 1, 0L, pmax.int(0L, diag(cov.gg) - post.psi))
} else {
var.exp <- NULL
}
results <- list(if(sw["mu.sw"]) list(means = mu,
var.mu = var.mu,
ci.mu = ci.mu),
if(sw["l.sw"]) list(loadings = lmat,
post.load = post.load,
var.load = var.load,
ci.load = ci.load),
if(sw["psi.sw"]) list(psis = psi,
var.psi = var.psi,
ci.psi = ci.psi),
if(sw.mx) list(post.mu = post.mu),
if(sw.px) list(post.psi = post.psi),
if(any(sw["l.sw"],
sw.px)) list(var.exp = var.exp))
result[[g]] <- unlist(results, recursive=FALSE)
}
names(result) <- gnames
GQ.res$Criteria <- c(GQ.res$Criteria, list(sd.AICMs = aicm.sd, sd.BICMs = bicm.sd,
best.models = t(vapply(GQ.res$Criteria, function(x) { inds <- arrayInd(which.max(x), dim(x));
c(clusters = rownames(x)[inds[1L]], factors = colnames(x)[inds[2L]]) }, character(2L)))))
class(GQ.res$Criteria) <- "listof"
class(GQ.res) <- "listof"
attr(GQ.res, "Clusters") <- n.grp
attr(GQ.res, "Criterion") <- criterion
attr(GQ.res, "Factors") <- n.fac
attr(GQ.res, "Supplied") <- c(Q=Q.T, G=G.T)
if(sw["mu.sw"]) {
mus2 <- lapply(result, "[[", "means")
var.mu <- provideDimnames(do.call(cbind, lapply(result, "[[", "var.mu")), base=list(if(uni) "" else varnames, gnames))
ci.mu <- lapply(result, "[[", "ci.mu")
class(mus2) <- class(ci.mu) <- "listof"
means <- list(mus = mus2, var.mu = var.mu, ci.mu = ci.mu)
} else means <- NULL
if(sw["mu.sw"] || sw.mx) {
post.mu <- provideDimnames(do.call(cbind, lapply(result, "[[", "post.mu")), base=list(if(uni) "" else varnames, gnames))
means <- c(means, list(post.mu = post.mu))
if(sw["mu.sw"]) {
means <- means[c(1L, 4L, 2L, 3L)]
}
class(means) <- "listof"
}
if(sw["mu.sw"]) {
last.mu <- if(clust.ind) .a_drop(mus[,,TN.store, drop=FALSE], drop=3L) else mu[,TN.store, drop=FALSE]
colnames(last.mu) <- gnames
means <- c(means, list(last.mu = last.mu))
class(means) <- "listof"
}
if(sw["l.sw"] <- tmpsw["l.sw"] && !all(Q == 0)) {
lnames <- lapply(Q, function(q) if(q > 0) paste0("Factor", seq_len(q)))
lmats2 <- .matnames(lapply(result, "[[", "loadings"), lnames)
post.load <- .matnames(lapply(result, "[[", "post.load"), lnames)
var.load <- .matnames(lapply(result, "[[", "var.load"), lnames)
ci.load <- .matnames(lapply(result, "[[", "ci.load"), lnames, dim=3)
last.lmat <- lapply(lmats2, function(x) { if(!is.null(x)) { x <- .a_drop(x[,,dim(x)[3L], drop=FALSE], drop=3L); class(x) <- "loadings" }; x })
loads <- list(lmats = lmats2, post.load = post.load, var.load = var.load, ci.load = ci.load, last.load = last.lmat)
loads <- lapply(loads, function(x) { class(x) <- "listof"; x} )
if(cshrink) attr(loads, "Sigma") <- GQ.res$post.sigma
class(loads) <- "listof"
}
if(sw["psi.sw"]) {
psis2 <- lapply(result, "[[", "psis")
var.psi <- provideDimnames(do.call(cbind, lapply(result, "[[", "var.psi")), base=list(if(uni) "" else varnames, gnames))
ci.psi <- lapply(result, "[[", "ci.psi")
class(psis2) <- class(ci.psi) <- "listof"
uniquenesses <- list(psis = psis2, var.psi = var.psi, ci.psi = ci.psi)
} else uniquenesses <- NULL
if(sw["psi.sw"] || sw.px) {
post.psi <- provideDimnames(do.call(cbind, lapply(result, "[[", "post.psi")), base=list(if(uni) "" else varnames, gnames))
uniquenesses <- c(uniquenesses, list(post.psi = post.psi))
if(sw["psi.sw"]) {
uniquenesses <- uniquenesses[c(1L, 4L, 2L, 3L)]
}
class(uniquenesses) <- "listof"
}
if(sw["psi.sw"]) {
last.psi <- if(clust.ind) .a_drop(psis[,,TN.store, drop=FALSE], drop=3L) else psi[,TN.store, drop=FALSE]
colnames(last.psi) <- gnames
uniquenesses <- c(uniquenesses, list(last.psi = last.psi))
class(uniquenesses) <- "listof"
}
dots$type <- ifelse(is.null(dots$type), "F", dots$type[1L])
# Calculate estimated covariance matrices & compute error metrics
if(error.metrics) {
Eseq <- seq_len(e.store)
sw.pi <- (equal.pro || G == 1) || sw["pi.sw"]
Fro <- mse <- mae <- medse <-
medae <- rmse <- nrmse <- rep(NA, e.store)
if(frobenius) {
orig.data <- as.data.frame(dat)
Pseq <- seq_len(n.var)
if(is.null(dots$dbreaks)) {
xbins <- if(isFALSE(uni)) lengths(lapply(orig.data, function(x) suppressWarnings(graphics::hist(x, plot=FALSE, ...))$counts)) else length(graphics::hist(orig.data[[1L]], plot=FALSE, ...)$counts)
dbreaks <- lapply(Pseq, function(p) c(-Inf, as.numeric(sub("\\((.+),.*", "\\1", levels(cut(orig.data[[p]], xbins[p]))[-1L])), Inf))
} else {
dbreaks <- dots$dbreaks
if(!is.list(dbreaks) ||
!all(vapply(dbreaks, is.factor, logical(1L))) ||
n.var !=
length(dbreaks)) stop(paste0("'dbreaks' must be a list of length P=", n.var, " of factors as outputted by the function 'cut'"), call.=FALSE)
dbreaks <- lapply(Pseq, function(p) unique(c(-Inf, as.numeric(sub("\\((.+),.*", "\\1", levels(dbreaks[[p]])[-1L])), Inf)))
xbins <- lengths(dbreaks) - 1L
}
dcounts <- mapply(tabulate, mapply(cut, orig.data, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
nbins <- max(xbins)
dat.bins <- vapply(dcounts, .padding, integer(nbins), nbins)
datnorm <- norm(dat.bins, dots$type)
rcounts <- vector("list", e.store)
} else if(Q0X) { warning("Need the means and uniquenesses to have been stored for zero-factor models in order to compute the posterior predictive reconstruction error\n", call.=FALSE, immediate.=TRUE)
} else warning("Need the means, uniquenesses, and loadings to have been stored in order to compute the posterior predictive reconstruction error\n", call.=FALSE, immediate.=TRUE)
}
if(clust.ind) {
if(frobenius) {
if(cov.met) {
a <- Reduce("+", lapply(Gseq, function(g) post.pi[g] * (tcrossprod(post.mu[,g]) + (if(Q0[g]) 0L else tcrossprod(post.load[[g]])) + (if(uni) post.psi[,g] else diag(post.psi[,g])))))
b <- Reduce("+", lapply(Gseq, function(g) post.pi[g] * post.mu[,g]))
cov.est <- a - tcrossprod(b)
if(sw.pi) {
if(equal.pro ||
G == 1) {
pi2 <- matrix(1/G, nrow=G, ncol=e.store)
} else {
pi2 <- pi.prop[,store.e, drop=FALSE]
}
} else warning("Mixing proportions not stored: can't compute all error metrics\n", call.=FALSE, immediate.=TRUE)
}
mu2 <- mus[,,store.e, drop=FALSE]
psi2 <- psis[,,store.e, drop=FALSE]
if(QX0) lmat2 <- lmats[,,,store.e, drop=FALSE]
if(!inf.Q) QsEr <- Q
tab.z <- tab.z[store.e,, drop=FALSE]
for(r in Eseq) {
Q0Er <- Q0E[,r]
if(frobenius) {
if(inf.Q) QsEr <- QsE[,r]
rdat <- suppressWarnings(sim_IMIFA_data(N=n.obs, G=G, P=n.var, Q=QsEr, nn=tab.z[r,], mu=mu2[,,r], psi=psi2[,,r], method="marginal", forceQg=FALSE,
loadings=if(!all(Q0Er)) lapply(Gseq, function(g) matrix(lmat2[,if(inf.Q) seq_len(QsEr[g]) else Qseq,g,r], nrow=n.var, ncol=QsEr[g]))))
rbinsx <- mapply(tabulate, mapply(cut, rdat, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
rbins <- vapply(rbinsx, .padding, integer(nbins), nbins)
Frob <- norm(dat.bins - rbins, dots$type)
rFro <- norm(rbins, dots$type)
minF <- abs(datnorm - rFro)
Fro[r] <- (Frob - minF)/(datnorm + rFro - minF)
rcounts[[r]] <- rbinsx
}
if(sw.pi && cov.met) {
a <- Reduce("+", lapply(Gseq, function(g) pi2[g,r] * (tcrossprod(mu2[,g,r]) + (if(Q0Er[g]) 0L else tcrossprod(lmat2[,,g,r])) + (if(uni) psi2[,g,r] else diag(psi2[,g,r])))))
b <- Reduce("+", lapply(Gseq, function(g) pi2[g,r] * mu2[,g,r]))
Sigma <- a - tcrossprod(b)
error <- cov.emp - Sigma
sqerr <- error^2
abserr <- abs(error)
mse[r] <- mean(sqerr)
medse[r] <- Median(sqerr)
mae[r] <- mean(abserr)
medae[r] <- Median(abserr)
rmse[r] <- sqrt(mse[r])
nrmse[r] <- rmse[r]/cov.range
}
}
} else if(cov.met) {
if(!sw["mu.sw"]) { warning("Means not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0, !sw["l.sw"],
!sw["psi.sw"])) { warning("Loadings & Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0,
!sw["l.sw"])) { warning("Loadings not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else warning("Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
}
} else {
if(any(sw["l.sw"], Q0X)) {
if(cov.met) {
cov.est <- diag(post.psi[,1L, drop=!uni]) + (if(Q0X) 0L else tcrossprod(post.load[[1L]]))
}
if(error.metrics && sw["psi.sw"]) {
if(frobenius) {
mu2 <- if(!sw["u.sw"]) matrix(0L, nrow=n.var, ncol=e.store) else mu[,store.e, drop=FALSE]
}
if(QX0) {
lmat2 <- lmats[,,store.e, drop=FALSE]
}
psi2 <- psi[,store.e, drop=FALSE]
if(!inf.Q) {
QsEr <- Q
QsMs <- Qseq
}
for(r in Eseq) {
if(frobenius) {
if(inf.Q) {
QsEr <- QsE[r]
QsMs <- seq_len(QsEr)
}
rdat <- suppressWarnings(sim_IMIFA_data(N=n.obs, G=G, P=n.var, Q=QsEr, mu=mu2[,r], psi=psi2[,r], method="marginal", forceQg=FALSE,
loadings=if(!Q0E[r]) matrix(lmat2[,QsMs,r], nrow=n.var, ncol=QsEr)))
rbinsx <- mapply(tabulate, mapply(cut, rdat, dbreaks, SIMPLIFY=FALSE), xbins, SIMPLIFY=FALSE)
rbins <- vapply(rbinsx, .padding, integer(nbins), nbins)
Frob <- norm(dat.bins - rbins, dots$type)
rFro <- norm(rbins, dots$type)
minF <- abs(datnorm - rFro)
Fro[r] <- (Frob - minF)/(datnorm + rFro - minF)
rcounts[[r]] <- rbinsx
}
if(cov.met) {
Sigma <- diag(psi2[,r, drop=!uni]) + (if(Q0E[r]) 0L else if(uni) crossprod(lmat2[,,r]) else tcrossprod(lmat2[,,r]))
error <- cov.emp - Sigma
sqerr <- error^2
abserr <- abs(error)
mse[r] <- mean(sqerr)
medse[r] <- Median(sqerr)
mae[r] <- mean(abserr)
medae[r] <- Median(abserr)
rmse[r] <- sqrt(mse[r])
nrmse[r] <- rmse[r]/cov.range
}
}
} else if(error.metrics) warning("Uniquenesses not stored: can't compute error metrics\n", call.=FALSE, immediate.=TRUE)
} else if(error.metrics) {
if(all(QX0, !sw["l.sw"],
!sw["psi.sw"])) { warning("Loadings & Uniquenesses not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
} else if(all(QX0,
!sw["l.sw"])) warning("Loadings not stored: can't compute error metrics or estimate posterior mean covariance matrix\n", call.=FALSE, immediate.=TRUE)
}
}
if(all(sw["psi.sw"] || !clust.ind, any(sw["l.sw"], Q0X))) {
var.exps <- lapply(lapply(result, "[[", "var.exp"), function(x) if(is.null(x)) NA else x)
clust.exps <- if(sum(is.na(var.exps)) == G) NULL else vapply(var.exps, sum, na.rm=TRUE, numeric(1L))/n.var
var.exps <- if(sum(is.na(var.exps)) == G) NULL else if(clust.ind) stats::setNames(rowSums2(sweep(do.call(cbind, var.exps), 2L, FUN="*", post.pi, check.margin=FALSE), na.rm=TRUE, useNames=FALSE), varnames) else stats::setNames(var.exps[[1L]], varnames)
Err <- list(Var.Exps = var.exps, Clust.Exps = clust.exps, Exp.Var = ifelse(clust.ind, sum(clust.exps * post.pi), unname(clust.exps)))
if(all(error.metrics, sw["mu.sw"] || !clust.ind)) {
if(cov.met) {
cov.est <- as.matrix(cov.est)
dimnames(cov.est) <- list(varnames, varnames)
err <- cov.emp - cov.est
sqerr <- err^2
abserr <- abs(err)
mse2 <- mean(sqerr)
rmse2 <- sqrt(mse2)
post.met <- c(MSE = mse2, MEDSE = Median(sqerr), MAE = mean(abserr), MEDAE = Median(abserr), RMSE = rmse2, NRMSE = rmse2/cov.range)
Err <- c(list(Empirical.Cov = cov.emp, Estimated.Cov = cov.est, Post = post.met), Err)
}
if(sw["psi.sw"] && any(all(cov.met, sw.pi), frobenius)) {
metrics <- rbind(MSE = mse, MEDSE = medse, MAE = mae, MEDAE = medae, RMSE = rmse, NRMSE = nrmse)
metrics <- if(frobenius) rbind(metrics, PPRE = Fro) else metrics
metrics <- metrics[!rowAlls(is.na(metrics), useNames=FALSE),, drop=FALSE]
metricCI <- rowQuantiles(metrics, probs=conf.levels)
metricCI <- if(cov.met && sw.pi) metricCI else provideDimnames(t(metricCI), base=list("PPRE", names(metricCI)))
med.met <- stats::setNames(matrixStats::rowMedians(metrics,
useNames=FALSE), rownames(metrics))
mean.met <- stats::setNames(rowMeans2(metrics, useNames=FALSE), rownames(metrics))
last.met <- c(MSE = mse[e.store], MEDSE = medse[e.store], MAE = mae[e.store], MEDAE = medae[e.store], RMSE = rmse[e.store], NRMSE = nrmse[e.store])
last.met <- if(frobenius) c(last.met, PPRE = Fro[e.store]) else last.met
last.met <- last.met[!is.na(last.met)]
Err <- if(cov.met) c(list(Median = med.met, Avg = mean.met, CIs = metricCI), Err[seq_len(4L)], if(all(cov.met, sw.pi)) c(list(Last.Cov = Sigma)), c(list(Final = last.met)), Err[5L:6L]) else c(list(Avg = mean.met, CIs = metricCI), Err[seq_len(3L)], list(Final = last.met))
if(frobenius) {
rcounts <- stats::setNames(lapply(Pseq, function(p) t(rowQuantiles(sapply(rcounts, "[[", p), probs=c(conf.levels[1L], 0.5, conf.levels[2L])))), varnames)
class(dcounts) <- "listof"
class(rcounts) <- "listof"
Err <- c(Err, list(PPRE = Fro, DatCounts = dcounts, RepCounts = rcounts, Breaks = dbreaks))
attr(Err, "Norm") <- dots$type
}
attr(Err, "ESS") <- e.store
errs <- ifelse(frobenius, ifelse(cov.met && sw.pi, "All", "PPRE"), "Covs")
} else errs <- ifelse(cov.met, "Post", "Vars")
} else errs <- "Vars"
class(Err) <- "listof"
} else errs <- "None"
error.metrics <- errs != "None"
if(sw["s.sw"]) {
if(adapt) {
eta <- if(method == "IFA") {
sapply(lapply(seq_len(e.store), function(i) .a_drop(eta[,colvec[[g]][,store.e][,i],i, drop=FALSE], drop=3L)),
function(x) cbind(x, matrix(0L, nrow=n.obs, ncol=ncol(eta) - ncol(x))), simplify="array")
} else {
sapply(seq_len(e.store), function(i) do.call(rbind, lapply(Gseq,
function(g) .a_drop(eta[z[i,] == g, c(which(colvec[[g]][,store.e][,i]), setdiff(seq_len(ncol(eta)), which(colvec[[g]][,store.e][,i]))),i, drop=FALSE], drop=3L)))[obsnames,, drop=FALSE], simplify="array")
}
}
eta <- eta[,Qseq,, drop=FALSE]
colnames(eta) <- paste0("Factor", Qseq)
scores <- list(eta = eta, post.eta = rowMeans(eta, dims=2), var.eta = apply(eta, c(1L, 2L), Var),
ci.eta = apply(eta, c(1L, 2L), stats::quantile, conf.levels), last.eta = .a_drop(eta[,,e.store, drop=FALSE], drop=3L))
attr(scores, "Eta.store") <- e.store
class(scores) <- "listof"
}
result <- c(if(exists("cluster", envir=environment())) list(Clust = cluster),
if(error.metrics) list(Error = Err), list(GQ.results = GQ.res),
if(sw["mu.sw"] || sw.mx) list(Means = means),
if(sw["l.sw"]) list(Loadings = loads),
if(sw["s.sw"]) list(Scores = scores),
if(sw["psi.sw"] || sw.px) list(Uniquenesses = uniquenesses))
attr(result, "Adapt") <- any(adapt, adapted)
attr(result, "Alph.step") <- if(is.element(method, c("IMFA", "IMIFA", "OMFA", "OMIFA"))) learn.alpha
attr(result, "Alpha") <- if(!learn.alpha) attr(sims, "Alpha")
attr(result, "C.Shrink") <- cshrink
attr(result, "Call") <- call
attr(result, "Center") <- attr(sims, "Center")
attr(result, "Choice") <- choice
attr(result, "Conf.Level") <- conf.level
attr(result, "Criterion") <- criterion
attr(result, "Disc.step") <- if(is.element(method, c("IMFA", "IMIFA"))) learn.d
attr(result, "Discount") <- if(is.element(method, c("IMFA", "IMIFA")) && !learn.d) attr(sims, "Discount")
attr(result, "Errors") <- errs
attr(result, "Equal.Pi") <- equal.pro
attr(result, "Exchange") <- attr(sims, "Exchange")
attr(result, "ForceQg") <- attr(sims, "ForceQg") || (adapt && G1 && forceQg)
attr(result, "G.init") <- if(inf.G) attr(sims, "G.init")
attr(result, "G.Mean") <- attr(sims, "G.Mean")
attr(result, "G.Scale") <- attr(sims, "G.Scale")
attr(result, "Ind.Slice") <- if(is.element(method, c("IMFA", "IMIFA"))) attr(sims, "Ind.Slice")
attr(result, "Kappa0") <- attr(sims, "Kappa0")
attr(result, "Method") <- method
attr(result, "N.Loadstore") <- if(inf.Q) lengths(Lstore) else rep(TN.store, G)
attr(result, "Name") <- data.name
attr(result, "N.Store") <- TN.store
attr(result, "Nowarn.G") <- ifelse(inf.G || length(n.grp) == 1, TRUE,
ifelse(G == min(n.grp), paste0("Best model occurs at the min of the number of components considered\n"),
ifelse(G == max(n.grp), paste0("Best model occurs at the max of the number of components considered\n"), TRUE)))
attr(result, "Nowarn.Q") <- ifelse(inf.Q || length(n.fac) == 1, TRUE,
ifelse(Q == min(n.fac), paste0("Best model occurs at the min of the number of factors considered\n"),
ifelse(Q == max(n.fac), paste0("Best model occurs at the max of the number of factors considered\n"), TRUE)))
attr(result, "Obs") <- n.obs
attr(result, "Obsnames") <- obsnames
attr(result, "Pitman") <- attr(sims, "Pitman")
attr(result, "range.G") <- attr(sims, "Clusters")
attr(result, "range.Q") <- attr(sims, "Factors")
attr(result, "Scale") <- attr(sims, "Scale")
attr(result, "Scaling") <- attr(sims, "Scaling")
attr(result, "Sd0.drop") <- attr(sims, "Sd0.drop")
attr(result, "Store") <- tmp.store
attr(result, "Switch") <- sw
attr(result, "Thresh") <- attr(sims, "Thresh")
attr(result, "Truncated") <- attr(sims, "Truncate")
attr(result, "TuneZeta") <- attr(sims, "TuneZeta")
attr(result, "Uni.Meth") <- uni.meth
attr(result, "Varnames") <- varnames
attr(result, "Vars") <- n.var
attr(result, "Z.sim") <- z.avgsim
class(result) <- "Results_IMIFA"
cat(print(result))
return(result)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.