Details
In this vignette, we want to show you how to use the new package to estimate cluster numbers. Of course, we start by loading the library.
Data
Now we simulate a simple dataset, with four predefined clusters.
set.seed(1)
x <- rbind(matrix(rnorm(30, sd = 0.1), ncol = 2),
matrix(rnorm(30, mean = 1, sd = 0.2), ncol = 2),
matrix(rnorm(30, mean = 5, sd = 0.1), ncol = 2),
matrix(rnorm(30, mean = 7, sd = 0.2), ncol=2))
Simulated Data Set
Available Indices
First, we want look at the set of available indices for estimating cluster numbers. We include some totally unnecessary manipulations simply to be able to order these method by the year of publication.
ai <- availableIndices()
desc <- ai$desc
nc <- nchar(desc)
year <- substring(desc, nc-3, nc)
ai[order(year),]## tag description
## ball ball Ball and Hall 1965
## rubin rubin rubin; Friedman and Rubin 1967
## friedman friedman friedman; Friedman and Rubin 1967
## beale beale Beale 1969
## marriot marriot Marriot 1971
## scott scott Scott and Symons 1971
## frey frey Frey and Van Groenewoud 1972
## duda duda Duda and Hart 1973
## pseudot2 pseudot2 pseudot2; Duda and Hart 1973
## gplus gplus G+; Rohlf 1974
## dunn dunn Dunn 1974
## tau tau Tau; Rohlf 1974
## ch ch Calinski and Harabasz 1974
## hartigan hartigan Hartigan 1975
## gamma gamma Baker and Hubert 1975
## mcclain mcclain McClain and Rao 1975
## hubert hubert Hubert and Levin 1976
## cindex cindex cindex; Hubert and Levin 1976
## ratkowsky ratkowsky Ratkowsky and Lance 1978
## db db Davies and Bouldin 1979
## ptbiserial ptbiserial Milligan 1980, 1981
## ccc ccc Sarle 1983
## tracew tracew TraceW; Milligan and Cooper 1985
## trcovw trcovw TrCovW; Milligan and Cooper 1985
## silhouette silhouette Rousseeuw 1987
## kl kl Krzanowski and Lai 1988
## sdindex sdindex Halkidi et al. 2000
## dindex dindex Lebart et al. 2000
## sdbw sdbw Halkidi and Vazirgiannis 2001
## gap gap Tibshirani et al. 2001Simple Test
In the original NbClust package, you could choose one test, or “all”. In the newer NumClust version, you can still select “all”, but you can also list any subset that you want to try. Here is an example, using two common indices:
NumClust(x, distance = "euclidean", method = "ward.D2",
min.nc = 2, max.nc = 6,
index = c("gap", "silhouette"))## $Stat
## gap silhouette
## 2 0.5983558 0.8384990
## 3 1.2651727 0.8882648
## 4 2.8787625 0.8789639
## 5 2.6350839 0.7564874
## 6 2.5135645 0.6717751
##
## $Critical
## gap
## 2 -0.6515477
## 3 -1.5931357
## 4 0.2687306
## 5 0.1509136
## 6 0.3239187
##
## $Best.nc
## gap silhouette
## N 4.0000000 3.0000000
## score.4 0.2687306 0.8882648Compare With Original
if(!requireNamespace("NbClust", quietly = TRUE)) {
stop("NumClust not available.\n")
}
library(NbClust)dmat<- dist(x, method = "euclidean", diag=FALSE)
dexes <- availableIndices()$tag
bb <- list()
NbBest <- list()
NbVal <- list()
numBest <- list()
numVal <- list()
for (dex in sort(dexes)) {
cat("\n\n-----------------\n", dex, "\n", file = stderr())
foo <- NumClust(x, diss = dmat, distance = NULL, method = "ward.D2",
min.nc = 2, max.nc = 6,
index = dex)
numBest[[dex]] <- foo$Best.nc
numVal[[dex]] <- foo$Best.nc[2, 1]
goo <- try(NbClust(x, diss = dmat, distance = NULL, method = "ward.D2",
min.nc = 2, max.nc = 6, index = dex))
if (inherits(goo, "try-error")) {
Nb <- rep(NA, length(foo$Stat))
best.nc <- NA
} else {
Nb <- goo$All.index
best.nc <- goo$Best.nc
if (is.null(best.nc)) {
best.nc <- c(NA, NA)
}
}
NbBest[[dex]] <- best.nc[[1]]
NbVal[[dex]] <- best.nc[[2]]
print(daft <- data.frame(Num=foo$Stat, Nb = Nb))
bb[[dex]] <- best.nc[[1]]
cat("Best:", best.nc, "\n")
cat("NumBest:", foo$Best.nc, "\n")
}## ball Nb
## 2 39.1973780 39.1974
## 3 5.8896703 5.8897
## 4 0.5282767 0.5283
## 5 0.3377401 0.3377
## 6 0.2260937 0.2261
## Best: 3 33.3077
## NumBest: 3 -33.30771
## beale Nb
## 2 51.7047317 51.7047
## 3 15.3395648 15.3396
## 4 0.8884561 0.8885
## 5 0.7721213 0.7721
## 6 0.9386913 0.9387
## Best: 4 0.8885
## NumBest: 4 0.4234108
## ccc Nb
## 2 14.80567 14.8057
## 3 16.02273 16.0227
## 4 23.84542 23.8454
## 5 21.67956 21.6796
## 6 20.29839 20.2984
## Best: 4 23.8454
## NumBest: 4 23.84542
## ch Nb
## 2 671.277 671.277
## 3 1561.453 1561.453
## 4 8688.900 8688.900
## 5 8012.296 8012.296
## 6 7836.778 7836.778
## Best: 4 8688.9
## NumBest: 4 8688.9
## cindex Nb
## 2 0.3571498 0.3571
## 3 0.3439237 0.3439
## 4 0.3201271 0.3201
## 5 0.3153900 0.3154
## 6 0.3621455 0.3621
## Best: 5 0.3154
## NumBest: 5 0.31539
## db Nb
## 2 0.2800896 0.2801
## 3 0.1356308 0.1356
## 4 0.1869418 0.1869
## 5 0.4510373 0.4510
## 6 0.5499630 0.5500
## Best: 3 0.1356
## NumBest: 3 0.1356308
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## dindex Nb
## 2 1.0799245 1.0799
## 3 0.4513364 0.4513
## 4 0.1610390 0.1610
## 5 0.1415145 0.1415
## 6 0.1274647 0.1275
## Best: NA NA
## NumBest: 2 -0.9945253
## duda Nb
## 2 0.01833136 0.0183
## 3 0.05921572 0.0592
## 4 0.51103872 0.5110
## 5 0.54599599 0.5460
## 6 0.49199030 0.4920
## Best: 4 0.511
## NumBest: 4 0.5110387
## dunn Nb
## 2 1.5189096 1.5189
## 3 1.2401601 1.2402
## 4 1.5072536 1.5073
## 5 0.1571841 0.1572
## 6 0.1996308 0.1996
## Best: 2 1.5189
## NumBest: 2 1.51891
## [1] "Frey index : No clustering structure in this data set"
## frey Nb
## 2 1.741933 1.7419
## 3 2.095140 2.0951
## 4 9.587710 9.5877
## 5 6.469208 6.4692
## 6 9.346128 9.3461
## Best: NA NA
## NumBest: NA NA
## friedman Nb
## 2 1193.002 1193.003
## 3 1281.868 1281.869
## 4 2182.392 2182.392
## 5 2704.176 2704.176
## 6 3405.087 3405.087
## Best: 4 900.5235
## NumBest: 4 900.5235
## gamma Nb
## 2 1.0000000 1.0000
## 3 1.0000000 1.0000
## 4 1.0000000 1.0000
## 5 0.9863777 0.9864
## 6 0.9846886 0.9847
## Best: 2 1
## NumBest: 2 1
## gap Nb
## 2 0.5983558 0.5984
## 3 1.2651727 1.2652
## 4 2.8787625 2.8788
## 5 2.6350839 2.6351
## 6 2.5135645 2.5136
## Best: 4 2.8788
## NumBest: 4 0.2687306
## gplus Nb
## 2 0.000000 0.0000
## 3 0.000000 0.0000
## 4 0.000000 0.0000
## 5 1.977401 1.9774
## 6 2.055367 2.0554
## Best: 2 0
## NumBest: 2 0
## hartigan Nb
## 2 199.337317 199.3373
## 3 419.612704 419.6127
## 4 14.073998 14.0740
## 5 13.466108 13.4661
## 6 9.140195 9.1402
## Best: 4 405.5387
## NumBest: 4 -405.5387
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
##
## hubert Nb
## 2 0.001084892 0.0011
## 3 0.001174784 0.0012
## 4 0.001233452 0.0012
## 5 0.001241350 0.0012
## 6 0.001242841 0.0012
## Best: NA NA
## NumBest: 2 -0.9999936
## kl Nb
## 2 7.987146e+00 7.9871
## 3 2.329332e+00 2.3293
## 4 4.991667e+03 4991.6666
## 5 2.934864e-02 0.0293
## 6 1.682832e+01 16.8283
## Best: 4 4991.667
## NumBest: 4 4991.667
## marriot Nb
## 2 287.54002 287.5400
## 3 139.48598 139.4860
## 4 17.54545 17.5454
## 5 17.68035 17.6803
## 6 16.24845 16.2484
## Best: 4 122.0754
## NumBest: 4 122.0754
## mcclain Nb
## 2 0.1631230 0.1631
## 3 0.1696218 0.1696
## 4 0.1284306 0.1284
## 5 0.1440709 0.1441
## 6 0.1517364 0.1517
## Best: 4 0.1284
## NumBest: 4 0.1284306
## pseudot2 Nb
## 2 1499.43722 1499.4372
## 3 444.84738 444.8474
## 4 12.43839 12.4384
## 5 10.80970 10.8097
## 6 10.32560 10.3256
## Best: 4 12.4384
## NumBest: 4 12.43839
## ptbiserial Nb
## 2 0.9214413 0.9214
## 3 0.8301886 0.8302
## 4 0.6781922 0.6782
## 5 0.6239792 0.6240
## 6 0.5877972 0.5878
## Best: 2 0.9214
## NumBest: 2 0.9214413
## ratkowsky Nb
## 2 0.6784105 0.6784
## 3 0.5721525 0.5722
## 4 0.4994636 0.4995
## 5 0.4468303 0.4468
## 6 0.4079672 0.4080
## Best: 2 0.6784
## NumBest: 2 0.6784105
## rubin Nb
## 2 28.91285 28.9128
## 3 128.28198 128.2820
## 4 1072.64599 1072.6460
## 5 1342.22488 1342.2249
## 6 1670.85297 1670.8530
## Best: 4 -674.7851
## NumBest: 4 -674.7851
## scott Nb
## 2 203.6820 203.6820
## 3 295.7417 295.7417
## 4 454.6537 454.6537
## 5 480.9714 480.9714
## 6 507.9174 507.9174
## Best: 4 158.912
## NumBest: 4 158.912
## sdbw Nb
## 2 0.079577346 0.0796
## 3 0.012767814 0.0128
## 4 0.002198366 0.0022
## 5 0.001682014 0.0017
## 6 0.001280817 0.0013
## Best: 6 0.0013
## NumBest: 6 0.001280817
## sdindex Nb
## 2 0.7795364 0.7795
## 3 0.9114994 0.9115
## 4 1.5979048 1.5979
## 5 6.5939910 6.5940
## 6 6.5728540 6.5729
## Best: 2 0.7795
## NumBest: 2 0.7795364
## silhouette Nb
## 2 0.8384990 0.8385
## 3 0.8882648 0.8883
## 4 0.8789639 0.8790
## 5 0.7564874 0.7565
## 6 0.6717751 0.6718
## Best: 3 0.8883
## NumBest: 3 0.8882648
## tau Nb
## 2 442.3729 442.3729
## 3 409.9576 409.9576
## 4 320.3390 320.3390
## 5 286.3638 286.3638
## 6 264.3638 264.3638
## Best: 2 442.3729
## NumBest: 2 442.3729
## tracew Nb
## 2 78.394756 78.3948
## 3 17.669011 17.6690
## 4 2.113107 2.1131
## 5 1.688700 1.6887
## 6 1.356562 1.3566
## Best: 3 45.1698
## NumBest: 3 45.16984
## trcovw Nb
## 2 2.7304255 2.7304
## 3 0.8701784 0.8702
## 4 0.8576112 0.8576
## 5 0.8461776 0.8462
## 6 0.3490341 0.3490
## Best: 3 1.8602
## NumBest: 3 -1.860247## ball beale ccc ch cindex db dindex
## 0.000 -0.465 0.000 0.000 0.000 0.000 NA
## duda dunn frey friedman gamma gap gplus
## 0.000 0.000 NA 0.000 0.000 -2.610 0.000
## hartigan hubert kl marriot mcclain pseudot2 ptbiserial
## 0.000 NA 0.000 0.000 0.000 0.000 0.000
## ratkowsky rubin scott sdbw sdindex silhouette tau
## 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## tracew trcovw
## 0.000 0.000dfr <- data.frame(Num = sapply(numBest, function(B) B[1]),
Nb = sapply(NbBest, function(B) B[[1]]))
apply(dfr, 1, diff)## ball beale ccc ch cindex db dindex
## 0 0 0 0 0 0 NA
## duda dunn frey friedman gamma gap gplus
## 0 0 NA 0 0 0 0
## hartigan hubert kl marriot mcclain pseudot2 ptbiserial
## 0 NA 0 0 0 0 0
## ratkowsky rubin scott sdbw sdindex silhouette tau
## 0 0 0 0 0 0 0
## tracew trcovw
## 0 0