1

I posted this question over on Esri's GeoNet, but have not yet gotten an answer.

Using ArcGIS v10.2.1, I have been using the Grouping Analysis tool to generate a range of clusters from my data. I wrote a Python script to run Grouping Analysis a large number of times (with different clustering algorithms and spatial constraints and with the number of clusters (k) ranging from 3 to 250) and report the 'grouping effectiveness' statistic (the Calinski-Harabasz pseudo-F statistic: see ArcGIS Help (http://resources.arcgis.com/en/help/main/10.2/#/How_Grouping_Analysis_works/005p0000004w000000/), or CH Index for short ) so I can try to identify the 'best' number(s) of clusters for my data. This has all been working fine.

To compare these results with some from additional clustering methods, I have also used another application (REDCAP: http://www.spatialdatamining.org/software/redcap) on my data as well. Fantastic though the software is, REDCAP does not report the CH Index to help evaluate the 'best' number of clusters in the data. I would like to be able to compare/evaluate my results from ArcGIS's Grouping Analysis with those from REDCAP using the same 'grouping effectiveness' statistic, so I have written another Python script to try to calculate the CH Index 'by hand', as it were, and this is where I've run into a problem.

I wrote my script to try to replicate the approach set out in the Help file for 'How Grouping Analysis works' (see the link to Esri's online Help above) and in the Python script the tool calls (Partition.py). I tried out my script on some of my output feature classes from Grouping Analysis, to check that my script was producing the same results. Here (finally) is the problem: when the Grouping Analysis was based on a single variable, I get exactly the same CH Index result from my script as from Esri's own tool; when the Grouping Analysis was based on more than one variable, I get a different (higher) CH Index result from my script. I am stumped as to why this is happening.

Here are some examples of the different CH Index values I get: For the data referenced in my code below, where the grouping is based on 3 variables and 9 groups:

CH Index from Esri tool: 1911.75679228486

CH Index from my script: 2448.47996865

For a similar dataset, using the same 3 variables but with 31 groups:

CH Index from Esri tool: 1069.1115254454

CH Index from my script: 1114.27024452

Compare this with the results from a dataset where the grouping is based on 1 variable, with 4 groups:

CH Index from Esri tool: 6349.09734505356

CH Index from my script: 6349.09734505

This is my Python 2.7 code:

# Get and print the starting time for running the script
import datetime
start_time = datetime.datetime.now()
print("Script started at " +str(start_time))

import arcpy
import numpy as np

arcpy.env.workspace = r"E:\Pr6653_GIS_Data_Maps\SKATER_ClusterOutputs_1021FGDB.gdb"
# The input feature class
inFC = "Clstrs_DstNcl_All_DSHC_RC_K9"
# The fields used to generate the clusters using Grouping Analysis
inFields = ["DstNcl_All_SCL01INV", "DSInterp_SCL01INV", "HCInterp_SCL01INV"]
# The field holding the values identifying the clusters
clstField = "SS_GROUP"
# Create a new list from inFields and add the cluster field, to be used when
#     creating the structured Numpy array from the input FC table
arrFields = list(inFields)
arrFields.append(clstField)

# Convert the attribute table of the input FC to a Numpy array
baseArr = arcpy.da.TableToNumPyArray(inFC, arrFields)
n = len(baseArr) * 1.0 # This is the number of features in the input FC

# Get values based on the clustering field
clstArr = np.unique(baseArr[clstField])
nc = len(clstArr) * 1.0 # This is the number of groups
nc_min = clstArr.min() # This is the minimum value for cluster ID
nc_max = clstArr.max() # This is the maximum value for cluster ID

SST = 0.0 # This is the total sum of squared differences
for fld in inFields: # For each field/variable...
    varMean = baseArr[fld].mean() # Get the overall mean for the variable
    varSS = ((baseArr[fld] - varMean)**2.0).sum() # Calculate the sum of squared differences from the variable mean for each value
    SST += varSS # Add that sum to the total sum of squared differences

SSE = 0.0 # This is the sum of squared differences within each group
for group in range(nc_min, nc_max + 1): # For each group...
    grpArr = baseArr[baseArr[clstField] == group] # Create a new array just for that group
    for fld in inFields: # For each field/variable...
        grpMean = grpArr[fld].mean() # Get the mean for the variable in that group
        grpSS = ((grpArr[fld] - grpMean)**2.0).sum() # Calculate the sum of squared differences from the group mean for each value in the group
        SSE += grpSS # Add that sum to the within-group sum of squared differences

# Calculate the pseudo-F statistic
R2 = (SST - SSE) / SST
fStat = float((R2 / (nc - 1.0)) / ((1.0 - R2) / (n - nc)))

# Print the results
print("For input " + inFC + ":")
print("Number of features = " + str(n))
print("Number of clusters = " + str(nc))
print("SST = " + str(SST))
print("SSE = " + str(SSE))
print("R2 = " + str(R2))
print("CH Index = " + str(fStat))
# Clean up
del baseArr, grpArr, clstArr
del n, nc, nc_max, nc_min
del SST, SSE, R2, fStat
del inFields, arrFields

print("Script complete. Elapsed time: " + str(datetime.datetime.now() - start_time))
# End of script

If it would help, I'd be happy to post some example data as well.

This code seems to work fine in a univariate situation, but not in a multivariate situation. I should also point out that I don't have a strong background in either statistics or programming.

Improve this question

1 Answer 1

Reset to default
1

I figured out my problem. I needed to standardise the values in each of the input fields to z-scores (subtract the mean and divide by the standard deviation).

My new code is as follows:

# Get and print the starting time for running the script
import datetime
start_time = datetime.datetime.now()
print("Script started at " +str(start_time))

import arcpy
import numpy as np

def zScore(inTable, inFld, outFld):
    """Calculate z-scores for values in one field and write them to another field"""
    fldArr = arcpy.da.TableToNumPyArray(inTable, inFld)
    fldMean = fldArr[inFld].mean()
    fldStdDev = fldArr[inFld].std()

    with arcpy.da.UpdateCursor(inTable, [inFld, outFld]) as cur:
            for x, y in cur:
                y = (float(x) - fldMean)/fldStdDev
                cur.updateRow([x,y])
    del fldArr
    del fldMean, fldStdDev

arcpy.env.workspace = r"E:\Pr6653_GIS_Data_Maps\SKATER_ClusterOutputs_1021FGDB.gdb"

inFC = "Clstrs_DstNcl_All_DSHC_RC_K9"
inFields = ["DstNcl_All_SCL01INV", "DSInterp_SCL01INV", "HCInterp_SCL01INV"]

arrFields = []

# Copy input dataset into memory
print("Copying source features into memory")
arcpy.CopyFeatures_management(inFC, "in_memory\\inMemSrcCopy")

# Add and populate fields for z-Scores
print("Adding and populating z-score fields for in-memory FC")
for fld in inFields:
    arcpy.AddField_management("in_memory\\inMemSrcCopy", fld + "_Z", "DOUBLE")
    zScore("in_memory\\inMemSrcCopy", fld, fld + "_Z")
    arrFields.append(fld + "_Z")

clstField = "SS_GROUP"
arrFields.append(clstField)

# Convert the attribute table of the in-memory input FC to a Numpy array
baseArr = arcpy.da.TableToNumPyArray("in_memory\\inMemSrcCopy", arrFields)
n = len(baseArr) # This is the number of features in the input FC

# Get values based on the clustering field
clstArr = np.unique(baseArr[clstField])
nc = len(clstArr) # This is the number of groups
nc_min = clstArr.min() # This is the minimum value for cluster ID
nc_max = clstArr.max() # This is the maximum value for cluster ID


SST = 0.0 # This is the total sum of squared differences
print("Calculating the total sum of squared differences")

for group in range(nc_min, nc_max + 1): # For each group...
    grpArr = baseArr[baseArr[clstField] == group] # Create a new array just for that group
    for fld in inFields: # For each field/variable...
        fld = fld + "_Z"
        varMean = baseArr[fld].mean() # Get the overall mean for the variable
        varSS = ((grpArr[fld] - varMean)**2.0).sum() # Calculate the sum of squared differences from the variable mean for each value in the group
        SST += varSS # Add that sum to the total sum of squared differences

SSE = 0.0 # This is the sum of squared differences within each group
print("Calculating the sum of squared differences within each group")
##lstSSE = []
for group in range(nc_min, nc_max + 1): # For each group...
# for group in range(nc): # For each group...
    grpArr = baseArr[baseArr[clstField] == group] # Create a new array just for that group
    for fld in inFields: # For each field/variable...
        fld = fld + "_Z"
        grpMean = (grpArr[fld]).mean() # Get the mean for the variable in that group
        grpSS = ((grpArr[fld] - grpMean)**2.0).sum() # Calculate the sum of squared differences from the group mean for each value in the group
        SSE += grpSS # Add that sum to the within-group sum of squared differences

R2 = (SST - SSE) / SST
fStat = (R2 / (nc - 1)) / ((1 - R2) / (n - nc))

print("For input " + inFC + ":")
print("\tNumber of features = " + str(n))
print("\tNumber of clusters = " + str(nc))
print("\tSST = " + str(SST))
print("\tSSE = " + str(SSE))
print("\tR2 = " + str(R2))
print("\tCH Index = " + str(fStat))

del baseArr, grpArr, clstArr
del n, nc, nc_max, nc_min
del SST, SSE, R2, fStat
del inFields, arrFields
arcpy.Delete_management("in_memory")

print("Script complete. Elapsed time: " + str(datetime.datetime.now() - start_time))
# End of script

This revised code produces results identical to those reported by Esri's tool.

Problem solved.

Improve this answer

Your Answer

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

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