Details
As usual, we start by loading the library.
Data
Now we simulate a simple dataset, with four predefined clusters.
set.seed(2468)
x <- rbind(matrix(rnorm(60, sd = 0.1), ncol = 4),
matrix(rnorm(60, mean = 1, sd = 0.2), ncol = 4),
matrix(rnorm(60, mean = 5, sd = 0.1), ncol = 4),
matrix(rnorm(60, mean = 7, sd = 0.2), ncol=4))
Simulated Data Set
Hierarchical Node-Based Assessment
We are now going to describe the “NodeTickler” algorithm. We assume that we will define clusters using the usual hierarchical clustering algorithm (hclust) and usintg cutree to define and assign samples to clusters. The underlying motivation for the new algorithm is that we want to examine the splits at each node as we descend through the dendrogram. Since each node is binary, at such a node, we have three sets of samples:
top: the samples that belong to the parent node;left: the samples assigned to the left-hand node by the split, andright`: the samples assigned to the right-hand node by the split.
The assignments to “left” or “right” are neither arbitrary nor symmetric. The documentation for the hclust function tells us that: “[t]he algorithm used … is to order the subtree so that the tighter cluster is on the left (the last, i.e., most recent, merge of the left subtree is at a lower value than the last merge of the right subtree).”
We compute the centroids (mean vector of features) across the samples in each of these three sample sets. We expect all three of these to be fairly highly correlated (which is basically why the left and right branches were joined together during the agglomerative clustering step that built the tree). So, we can measure the correlation coefficients between these three vectors. The following code chunk devines a function that carries out these computations at the Kth node during the descent.
parseSplit <- function(dset, HC, K) {
high <- cutree(HC, k = K)
low <- cutree(HC, k = K+1)
tbl <- table(high, low)
v <- which(apply(tbl, 1, function(X) sum(X > 0)) == 2)
w <- which(tbl[v,] > 0)
top <- apply(dset[high == v, ], 2, mean)
left <- apply(dset[low == w[1], , drop = FALSE], 2, mean)
right <- apply(dset[low == w[2], , drop = FALSE], 2, mean)
cor(cbind(top, left, right))
}Now we explicitly apply this function across a range of values for the number of clusters.
dmat <- dist(x)
HC <- hclust(dist(x), "ward.D2")
N <- 2:12
ps <- lapply(N, function(K) {
parseSplit(x, HC, K)
})
RT <- sapply(ps, function(X) X[1,3])
LT <- sapply(ps, function(X) X[1,2])
LR <- sapply(ps, function(X) X[2,3])opar <- par(mfrow = c(3,2))
plot(RT, LT, pch = 16)
abline(0,1)
plot(LR, pch = 16, main = "Left-Right")
fit <- loess(LR ~ N)
lines(fit$x, fit$fitted, lwd =2, col = "red")
Q = 40
abline(v = Q)
plot(pmax(LT, RT), pch = 16, main = "Best-Top")
plot(pmin(LT, RT), pch = 16, main = "Worst-Top")
fit <- loess(pmin(LT, RT) ~ N)
lines(fit$x, fit$fitted, lwd =2, col = "red")
abline(v = Q)
plot(RT, pch = 16, main = "Right-Top")
plot(LT, pch = 16, main = "Left-Top")
fit <- loess(LT ~ N)
lines(fit$x, fit$fitted, lwd =2, col = "red")
abline(v = Q)
Compartiosons of corrleation values between node components.
We can now wrap a function around the parseSplit function so it will accept the general argl argument list documented in the “registries” man page.
nodeTickler <- function(argl) {
dset <- argl$jeu
HC = hclust(argl$md, "ward.D2")
K = length(unique(argl$cl1))
step <- parseSplit(dset, HC, K)
LT <- step[1, 2]
LR <- step[2, 3]
LT - LR
}We also write a function to select the best number of clusters. Here, we want to stop one level before the score begins to decrease.
bestTickler <- function(vals, ignore) {
dlev <- c(NA, diff(vals))
dlev[is.na(dlev)] <- max(abs(dlev), na.rm = TRUE) + 1
if (any(dlev < 0)) {
index <- max(1, which(dlev < 0)[1] - 1)
score <- dlev[index+1]
} else {
index <- NA
score <- NA
}
list(index = index, score = score)
}Then we register the method…
registerIndice("nodeTickler", "Node Based Evaluation 2026",
nodeTickler, bestTickler)
availableIndices()## tag description
## kl kl Krzanowski and Lai 1988
## gplus gplus G+; Rohlf 1974
## tracew tracew TraceW; Milligan and Cooper 1985
## hartigan hartigan Hartigan 1975
## dunn dunn Dunn 1974
## duda duda Duda and Hart 1973
## tau tau Tau; Rohlf 1974
## hubert hubert Hubert and Levin 1976
## nodeTickler nodeTickler Node Based Evaluation 2026
## ccc ccc Sarle 1983
## silhouette silhouette Rousseeuw 1987
## sdindex sdindex Halkidi et al. 2000
## ch ch Calinski and Harabasz 1974
## cindex cindex cindex; Hubert and Levin 1976
## frey frey Frey and Van Groenewoud 1972
## rubin rubin rubin; Friedman and Rubin 1967
## ratkowsky ratkowsky Ratkowsky and Lance 1978
## beale beale Beale 1969
## db db Davies and Bouldin 1979
## marriot marriot Marriot 1971
## ptbiserial ptbiserial Milligan 1980, 1981
## dindex dindex Lebart et al. 2000
## trcovw trcovw TrCovW; Milligan and Cooper 1985
## scott scott Scott and Symons 1971
## gamma gamma Baker and Hubert 1975
## pseudot2 pseudot2 pseudot2; Duda and Hart 1973
## ball ball Ball and Hall 1965
## mcclain mcclain McClain and Rao 1975
## friedman friedman friedman; Friedman and Rubin 1967
## sdbw sdbw Halkidi and Vazirgiannis 2001
## gap gap Tibshirani et al. 2001… and apply it to our simulated data set, along with several older emthods.
dexes <- c("nodeTickler", "silhouette", "gap", "ball")
nc <- NumClust(x, diss = dmat, distance = NULL, method = "ward.D2",
min.nc = 2, max.nc = 10,
index = dexes)
nc## $Stat
## nodeTickler silhouette gap ball
## 2 0.11074295 0.8366636 0.5944488 78.1847377
## 3 0.18010332 0.8802667 1.2352380 12.0575653
## 4 1.18386582 0.8559884 2.6221024 1.3567373
## 5 1.11584443 0.6904229 2.2906955 0.9469275
## 6 0.70193590 0.5675375 2.1054701 0.6756001
## 7 0.14207758 0.5751713 1.7731798 0.5186317
## 8 0.81713780 0.5778402 1.6322999 0.4034194
## 9 0.38965092 0.5787063 1.5327358 0.3258970
## 10 0.04556285 0.5706051 1.4673370 0.2654282
##
## $Critical
## gap
## 2 -0.6255201
## 3 -1.3664102
## 4 0.3564588
## 5 0.2146196
## 6 0.3700355
## 7 0.1843476
## 8 0.1493877
## 9 0.1201936
## 10 0.1568599
##
## $Best.nc
## nodeTickler silhouette gap ball
## N 4.00000000 3.0000000 4.0000000 3.00000
## score.5 -0.06802139 0.8802667 0.3564588 -66.12717