gapstat_ord: Estimate the gap statistic on an ordination result

In joey711/phyloseq: Handling and analysis of high-throughput microbiome census data

Estimate the gap statistic on an ordination result

Description

This is a wrapper for the clusGap function, expecting an ordination result as the main data argument.

Usage

gapstat_ord(ord, axes = c(1:2), type = "sites",
  FUNcluster = function(x, k) {     list(cluster = pam(x, k, cluster.only
  = TRUE)) }, K.max = 8, ...)

Arguments

ord	(Required). An ordination object. The precise class can vary. Any ordination classes supported internally by the phyloseq package should work, ultimately by passing to the scores function or its internal extensions in phyloseq.
axes	(Optional). The ordination axes that you want to include.
type	(Optional). One of "sites" (the vegan package label for samples) or "species" (the vegan package label for OTUs/taxa). Default is "sites".
FUNcluster	(Optional). This is passed to clusGap. The documentation is copied here for convenience: a function which accepts as first argument a (data) matrix like x, second argument, say (the number of desired clusters) k, where k >= 2, and returns a list with a component named (or shortened to) cluster which is a vector of length n = nrow(x) of integers in 1:k determining the clustering or grouping of the n observations. The default value is the following function, which wraps partitioning around medoids, pam: function(x, k){list(cluster = pam(x, k, cluster.only=TRUE))} Any function that has these input/output properties (performing a clustering) will suffice. The more appropriate the clustering method, the better chance your gap statistic results will be useful.
K.max	(Optional). A single positive integer value. It indicates the maximum number of clusters that will be considered. Value must be at least two. This is passed to clusGap.
...	(Optional). Additional named parameters passed on to clusGap. For example, the method argument provides for extensive options regarding the method by which the “optimal” number of clusters is computed from the gap statistics (and their standard deviations). See the clusGap documentation for more details.

Value

An object of S3 class "clusGap", basically a list with components. See the clusGap documentation for more details.

Examples

data("soilrep")
sord  = ordinate(soilrep, "PCoA", "bray")
# Evaluate axes with scree plot
plot_scree(sord)
# Gap Statistic
gs = gapstat_ord(sord, axes=1:3, verbose=FALSE)
# plot_ordination(soilrep, sord,  color="Treatment")
plot_clusgap(gs)
print(gs, method="Tibs2001SEmax")

joey711/phyloseq documentation built on Nov. 4, 2022, 1:16 a.m.