# Cubic Clustering Criterion using R update

I have seen other users ask about recreating SAS's CCC output in other programs. This question, Cubic clustering criterion in R, has an answer that says to use NbClust to calculate, but that function does not handle large datasets well. It makes a call to dist that must allocate a 50 gig object. I have tried replacing the function with cluster::daisy, and proxy::dist from this SO question with the same memory problems.

Avoiding the dist call altogether may be the best option. I am looking to other options to recreate it. In this question How to define number of clusters in K-means clustering?, a user goes through the math provided by SAS. But I do not have the stats chops to translate that into R code.

Keeping it simple, I have kmeans output that provides total sum of squares (tot.ss), within.ss, between.ss, and I also calculated the $R^2$.

kmeans(x = mydata, centers = 23, iter.max = ITER)
Within cluster sum of squares by cluster:
[1]  91248.77  72122.06  78680.32  90402.25  86341.35 153533.51  73988.63  64903.32
[9]  38334.98  84125.14  92366.93  74721.24 110313.76  96859.55  84516.37  56068.08
[17]  76201.69  86194.35  59526.00  53709.75  72503.21  50767.36  80531.94
(between_SS / total_SS =  36.5 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"
[7] "size"         "iter"         "ifault"

Can I calculate the CCC using these measures?

The second question has a long description from the SAS pdf. But I saw a simplified equation here.

where $E(R^2)$ is the expected $R^2$, and $R^2$ is the observed $R^2$, and $K$ is the variance-stabilizing transformation.

*Can this equation be completed by R's kmeans output and a calculated $R^2$

Edit

One reason why I am focusing on kmeans is that SAS users utilize PROC FASTCLUS when running large datasets. It is equivalent to R's kmeans function. The package NbClust calculates the CCC that I'm looking for, but it does it on the full data with euclidean distance, which is impossible for most computers. That is equivalent to SAS's PROC CLUSTER.

• The CCC is a comparison of how much better the clustering is than a uniform random. If you had a decent first estimator of kernel bandwidth, then you could use a kernel method to tell you "where your data is". Also called a "nonparametric density". You could restrict your uniform sampling to those regions if you use a threshold. (0.001 aka 0.1%?) It might require a stronger signal to show, but might also work against things like "the curse of dimensionality". – EngrStudent - Reinstate Monica May 12 '16 at 21:19
• Can you say that in English? : ) Is the kernel bandwidth and uniform sampling thresholds in the R output or input parameters? – Pierre Lafortune May 12 '16 at 21:21

I figured it out. The solution I chose is to rewrite and extract what NbClust was doing but to exclude the dist matrix call and everything else that I did not need. I check my custom CCC function against the actual output to be sure that the output is the same:

> NbClust(iris[-5], min.nc=3, max.nc=3, index="ccc", method="kmeans")
$All.index [1] 37.6701$Best.nc
Number_clusters     Value_Index
3.0000         37.6701

> CCC(iris[-5], 3)
[1] 37.67012

The measures are the same.

Here's the full function for anyone else interested in CCC in R.

As some may notice, this function also calculates other disgnostics like c("scott", "rubin", "marriot", "friedman"). I only needed the CCC for my purposes but the others can also be extracted:

CCC <- function(data, nc) {
Indices.WBT <- function(x,cl,P,s,vv)
{
n <- dim(x)[1]
pp <- dim(x)[2]
qq <- max(cl)
z <- matrix(0,ncol=qq,nrow=n)
clX <- as.matrix(cl)

for (i in 1:n)
for (j in 1:qq)
{
z[i,j]==0
if (clX[i,1]==j)
{z[i,j]=1}
}

xbar <- solve(t(z)%*%z)%*%t(z)%*%x
B <- t(xbar)%*%t(z)%*%z%*%xbar
W <- P-B
marriot <- (qq^2)*det(W)
trcovw <- sum(diag(cov(W)))
tracew <- sum(diag(W))
if(det(W)!=0)
scott <- n*log(det(P)/det(W))
else {cat("Error: division by zero!")}
friedman <- sum(diag(solve(W)*B))
rubin <- sum(diag(P))/sum(diag(W))

R2 <- 1-sum(diag(W))/sum(diag(P))
v1 <- 1
u <- rep(0,pp)
c <- (vv/(qq))^(1/pp)
u <- s/c
k1 <- sum((u>=1)==TRUE)
p1 <- min(k1,qq-1)
if (all(p1>0,p1<pp))
{
for (i in 1:p1)
v1 <- v1*s[i]
c <- (v1/(qq))^(1/p1)
u <- s/c
b1 <- sum(1/(n+u[1:p1]))
b2 <- sum(u[p1+1:pp]^2/(n+u[p1+1:pp]),na.rm=TRUE)
E_R2 <- 1-((b1+b2)/sum(u^2))*((n-qq)^2/n)*(1+4/n)
ccc <- log((1-E_R2)/(1-R2))*(sqrt(n*p1/2)/((0.001+E_R2)^1.2))
}else
{
b1 <- sum(1/(n+u))
E_R2 <- 1-(b1/sum(u^2))*((n-qq)^2/n)*(1+4/n)
ccc <- log((1-E_R2)/(1-R2))*(sqrt(n*pp/2)/((0.001+E_R2)^1.2))
}
results <- list(ccc=ccc,scott=scott,marriot=marriot,trcovw=trcovw,tracew=tracew,friedman=friedman,rubin=rubin)
return(results)
}

nc <- nc
jeu1 <- as.matrix(data)
numberObsBefore <- dim(jeu1)[1]
jeu <- na.omit(jeu1) # returns the object with incomplete cases removed
nn <- numberObsAfter <- dim(jeu)[1]
pp <- dim(jeu)[2]
TT <- t(jeu)%*%jeu
sizeEigenTT <- length(eigen(TT)$value) eigenValues <- eigen(TT/(nn-1))$value

for (i in 1:sizeEigenTT)
{
if (eigenValues[i] < 0) {
#cat(paste("There are only", numberObsAfter,"nonmissing observations out of a possible", numberObsBefore ,"observations."))
stop("The TSS matrix is indefinite. There must be too many missing values. The index cannot be calculated.")
}
}
s1 <- sqrt(eigenValues)
ss <- rep(1,sizeEigenTT)
for (i in 1:sizeEigenTT)
{
if (s1[i]!=0)
ss[i]=s1[i]
}
vv <- prod(ss)

set.seed(1)
cl1 <- kmeans(jeu,nc)$cluster TT <- t(jeu)%*%jeu Indices.WBT(x=jeu, cl=cl1, P=TT,s=ss,vv=vv)$ccc
}