Calculating distances in SAS PROC CLUSTER method=centroid

Question

I'm trying to understand the calculations behind the centroid method of cluster analysis in SAS's PROC CLUSTER. I am able to calculate the reported distances at each join when I set the options to method=centroid nonorm, but I am unable to figure out how the reported distances at each join are calculated when I set the options to method=centroid nonorm nosquare. In the documentation, it states that "Distance data are squared unless you specify the NOSQUARE option."

Below is some code using a very simple example with 5 points in a plane (from Afifi and Clark, 1996).

data ac;
    input x y;
    cards;
1 1
1 2
3 6
5 7
8 5
;

* this gives centroid distances of 1, 2.2361, 4.272, and 6.2472 ;
proc cluster data=ac method=centroid nonorm;
    run;

* this gives centroid distances of 1, 2.2361, 3.7933, and 5.0593 ;
proc cluster data=ac method=centroid nonorm nosquare;
    run;

Here are my calculations of the distances "by hand" for the first PROC CLUSTER statement.

CL4: (1, 1) and (1, 2)
    distance = sqrt((1-1)^2 + ((1-2)^2) = 1
    centroid = ((1+1)/2, (1+2)/2) = (1, 1.5)

CL3: (3, 6) and (5, 7)
    distance = sqrt((3-5)^2 + ((6-7)^2) = 2.236
    centroid = ((3+5)/2, (6+7)/2) = (4, 6.5)

CL2: CL3 and (8, 5)
    distance = sqrt((4-8)^2 + ((6.5-5)^2) = 4.272
    centroid = ((3+5+8)/3, (6+7+5)/3) = (5.3, 6)

CL1: CL4 and CL2
    distance = sqrt((1-5.3)^2 + ((1.5-6)^2) = 6.247
    centroid = ((1+1+3+5+8)/5, (1+2+6+7+5)/5) = (3.6, 4.2)

Can anyone help me calculate the distances by hand for the second PROC CLUSTER statement?

Answer 1

This is likely related to the implementation details of the clustering algorithm, specifically the Lance-Williams dissimilarity update (Original paper; Descriptions: 1, 2, wiki; SAS states this is what it uses), which is a generalized formula to recursively calculate the distances/dissimilarity between any two clusters without needing the original coordinates of any of the points - only the distance/dissimilarity matrix between the points.

For any method that uses euclidean distances, this update algorithm needs the distances to be squared to match the real distances that you get when calculating by using the original point coordinates (someone more well-versed than will have to explain why).

To see this explicitly, this is what your example looks like in R:

x = structure(list(V1 = c(1L, 1L, 3L, 5L, 8L), V2 = c(1L, 2L, 6L, 7L, 5L)), .Names = c("X", "Y"), class = "data.frame", row.names = c(NA, -5L))
d = dist(x)

h1 = hclust(d, method='centroid')
h1$height

Which gives you the same results as NOSQUARE

[1] 1.000000 2.236068 3.793268 5.059300

And then using squared distances:

h2 = hclust(d^2, method='centroid')
sqrt(h2$height)

Will give you the the expected

[1] 1.000000 2.236068 4.272002 6.247222