I am interested in determining the optimal number of clusters calculated by the PAM clustering algorithm using the Calinski-Harabasz (CH) index. To that end, I found 2 different R functions calculating CH values for a given clustering, but which returned different results: ?cluster.stats (in the fpc package), and ?index.G1 (in the clusterSim package).

First one is called via:

pam.res <- pam(dist.matrix, 2, diss=TRUE)
ch1     <- cluster.stats(dist.matrix, pam.res$clustering, silhouette=TRUE)$ch

Second one is called via:

ch2 <- index.G1(t(dataframe), pam.res$clustering, d=dist.matrix)

Data may be found here: dataframe.RData, or here: dist.matrix.RData [dead links].

  • Can anybody explain the difference between these two CH index calculations to me?

    Using cluster.stats(), the highest CH index is obtained for 2 clusters ($\approx32$); while using index.G1(), the highest CH index is obtained for 3 clusters ($\approx60$, and the value for 2 clusters is totally different from the previous, $\approx54$).

  • Which function is normally used to calculate the CH index?

  • 1
    $\begingroup$ Not everyone here uses R to read in your file format. Please give your data in more convenient format (plain text, Excel?) or publish it in your question if it's not too lengthy. $\endgroup$
    – ttnphns
    Commented Sep 21, 2011 at 11:36
  • $\begingroup$ You are right. Here you may find the files in txt format: dataframe.txt and dist.matrix.txt Loading dist.matrix.txt must be followed by converting it to dist object as follows: otu <- read.table("dataframe.txt") dist.matrix <- as.dist(read.table("dist.matrix.txt")) $\endgroup$
    – anat
    Commented Sep 21, 2011 at 16:10
  • $\begingroup$ This question might be better suited for the r-help listserve, where (presumably) the package maintainers would be notified of the discrepancy. $\endgroup$ Commented Sep 23, 2013 at 15:04
  • $\begingroup$ This question appears to be off-topic because it is about potentially conflicting R functions. $\endgroup$ Commented Jan 26, 2014 at 4:08
  • $\begingroup$ This question requires statistical knowledge to answer, so I think it's on topic here (although contacting the package maintainers directly might be best). $\endgroup$
    – Peter Flom
    Commented Jan 26, 2014 at 10:13

2 Answers 2


There is one method of calculating Caliński & Harabasz (1974) index for the same distance matrix, so if two R functions show different results one of them is wrong. Hence your question is off-topic.

  • Look how Caliński & Harabasz index is calculated, in their original paper [1] or e.g. here.

  • Then check the source code of both R functions, find a bug and report it to package creators.

Here are fpc and clusterSim example sites on GitHub where you can view the source code: https://github.com/cran/fpc/tree/master/R, https://github.com/cran/clusterSim.

[1] Caliński, T., and J. Harabasz. "A dendrite method for cluster analysis." Communications in Statistics. Vol. 3, No. 1, 1974, pp. 1–27.


Using a synthetic, two dimensional dataset of 200 points, euclidean distance and complete linkage I am not able to reproduce the discrepancies which you encountered. Also the clusterCrit package and another implementation return the same values

> # fpc
> ch1 <- calinhara(X, pc, cn=max(pc))
> # clusterSim
> ch2 <- index.G1 (X,pc,d=NULL,centrotypes="centroids")
> # clusterCrit
> ch3 <- as.numeric(intCriteria(X,pc,"Calinski_Harabasz"))
> cat('fpc: ', ch1, '\nclusterSim: ', ch2, '\nclusterCrit: ', ch3)
fpc:  369.0315 
clusterSim:  369.0315 
clusterCrit:  369.0315


>>> itn.calinski_harabasz(X, pc)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.