plot_clusgap: Create a ggplot summary of gap statistic results

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

Description

Create a ggplot summary of gap statistic results

Usage

plot_clusgap(clusgap, title = "Gap Statistic results")

Arguments

clusgap	(Required). An object of S3 class "clusGap", basically a list with components. See the clusGap documentation for more details. In most cases this will be the output of gapstat_ord, or clusGap if you called it directly.
title	(Optional). Character string. The main title for the graphic. Default is "Gap Statistic results".

Value

A ggplot plot object. The rendered graphic should be a plot of the gap statistic score versus values for k, the number of clusters.

See Also

gapstat_ord

clusGap

ggplot

Examples

# Load and process data
data("soilrep")
soilr = rarefy_even_depth(soilrep, rngseed=888)
print(soilr)
sample_variables(soilr)
# Ordination
sord  = ordinate(soilr, "DCA")
# Gap Statistic
gs = gapstat_ord(sord, axes=1:4, verbose=FALSE)
# Evaluate results with plots, etc.
plot_scree(sord)
plot_ordination(soilr, sord,  color="Treatment")
plot_clusgap(gs)
print(gs, method="Tibs2001SEmax")
# Non-ordination example, use cluster::clusGap function directly
library("cluster")
pam1 = function(x, k){list(cluster = pam(x, k, cluster.only=TRUE))}
gs.pam.RU = clusGap(ruspini, FUN = pam1, K.max = 8, B = 60)
gs.pam.RU
plot(gs.pam.RU, main = "Gap statistic for the 'ruspini' data")
mtext("k = 4 is best .. and  k = 5  pretty close")
plot_clusgap(gs.pam.RU)

Example output

`set.seed(888)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(888); .Random.seed` for the full vector
...
5448OTUs were removed because they are no longer 
present in any sample after random subsampling

...
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 11377 taxa and 56 samples ]
sample_data() Sample Data:       [ 56 samples by 4 sample variables ]
[1] "Treatment" "warmed"    "clipped"   "Sample"   
Clustering Gap statistic ["clusGap"] from call:
clusGap(x = x[, axes], FUNcluster = FUNcluster, K.max = K.max,     verbose = FALSE)
B=100 simulated reference sets, k = 1..8; spaceH0="scaledPCA"
 --> Number of clusters (method 'Tibs2001SEmax', SE.factor=1): 2
         logW   E.logW       gap     SE.sim
[1,] 3.405192 3.658805 0.2536129 0.03425815
[2,] 3.043528 3.503535 0.4600062 0.03740694
[3,] 2.913456 3.372063 0.4586076 0.03960065
[4,] 2.824380 3.267316 0.4429355 0.03961456
[5,] 2.743791 3.175679 0.4318879 0.03962653
[6,] 2.667146 3.097106 0.4299601 0.03905863
[7,] 2.594961 3.024891 0.4299297 0.03916471
[8,] 2.504070 2.954818 0.4507475 0.03704226
Clustering Gap statistic ["clusGap"] from call:
clusGap(x = ruspini, FUNcluster = pam1, K.max = 8, B = 60)
B=60 simulated reference sets, k = 1..8; spaceH0="scaledPCA"
 --> Number of clusters (method 'firstSEmax', SE.factor=1): 4
         logW   E.logW         gap     SE.sim
[1,] 7.187997 7.127156 -0.06084119 0.04281827
[2,] 6.628498 6.766057  0.13755839 0.04852638
[3,] 6.261660 6.559480  0.29782002 0.04159768
[4,] 5.692736 6.380305  0.68756887 0.04975800
[5,] 5.580999 6.232089  0.65109047 0.04515626
[6,] 5.500583 6.111116  0.61053249 0.04582079
[7,] 5.394195 6.005277  0.61108186 0.04652737
[8,] 5.320052 5.914224  0.59417215 0.04661338

phyloseq documentation built on Nov. 8, 2020, 6:41 p.m.