Gap statistics calculation as a robustness check for testing cluster validity using Gower dissimilarity matrix and Wards linkage

[Anna Novaresio]

I would like to calculate and plot the gap statistic introduced by Tibshirani et al (2001) as a form of robustness check to test the cluster validity of a clustering exercise based on Gower method and Wards linkage.

I have successfully calculated the "within-cluster dispersion" on my real data (number of clusters k =10) implementing these lines of code (step 1):

 use "mydb", clear

local vars "cat_var_1 cat_var_2 cat_var_3 cat_var_4 num_var_w_1 num_var_w_2 num_var_w_3 num_var_w_4 num_var_w_5 "
local K_max = 10
local B = 10

matrix dissimilarity G_real = `vars', gower

mata:
D = st_matrix("G_real")
n = rows(D)
K_max = st_numscalar("K_max")  // pass local to Mata
results = J(10, 2, .)

for (k=1; k<=10; k++) {
    // generate clusters using Ward's linkage in Mata (vectorized)
    g = J(n,1,1)  // placeholder: assign all to cluster 1 (replace with clustering in Mata)
    // compute Wk safely using D and g
    Wk = 0
    for (i=1; i<=n; i++) {
        for (j=i+1; j<=n; j++) {
            if (g[i] == g[j]) Wk = Wk + D[i,j]
        }
    }
    if (Wk > 0) results[k,2] = ln(Wk)
    results[k,1] = k
}
st_matrix("results", results)
end

matrix list results

I am struggling to calculate the within cluster dispersion on simulated data (number of clusters k=10) using these lines of codes (step 2):

clear all
macro drop _all

local B = 10
local K_max = 10
local vars "cat_var_1 cat_var_2 cat_var_3 cat_var_4 num_var_w_1 num_var_w_2 num_var_w_3 num_var_w_4 num_var_w_5"

* Load your dataset
use "mydb", clear

* Number of observations
summarize `vars'
local n = r(N)

* Store min/max for each variable
matrix mins = J(1, `: word count `vars'', .)
matrix maxs = J(1, `: word count `vars'', .)
local i = 1
foreach v of local vars {
    quietly summarize `v'
    matrix mins[1,`i'] = r(min)
    matrix maxs[1,`i'] = r(max)
    local ++i
}

* Initialize cumulative logWk matrix
matrix logWk_sim_total = J(10,1,0)

*-----------------------------------------------
* Run simulations fully in Mata
*-----------------------------------------------
mata:
B = st_numscalar("B")
K_max = st_numscalar("K_max")
n = st_numscalar("n")
mins = st_matrix("mins")
maxs = st_matrix("maxs")

// Import logWk_sim_total safely as column vector
logWk_sim_total = st_matrix("logWk_sim_total")
logWk_sim_total = vec(logWk_sim_total)

vars_count = cols(mins)

for (b=1; b<=B; b++) {

    // Generate uniform simulated data
    X = J(n, vars_count, 0)
    for (j=1; j<=vars_count; j++) {
        X[,j] = runiform(n,1)*(maxs[1,j]-mins[1,j]) + mins[1,j]
    }

    // Compute Gower distance
    D = J(n,n,0)
    for (i=1; i<=n; i++) {
        for (j=i+1; j<=n; j++) {
            D[i,j] = mean(abs(X[i,.]-X[j,.]) / (maxs - mins))
            D[j,i] = D[i,j]
        }
    }

    // Loop over cluster counts
    for (k=1; k<=10; k++) {

        // Simple cluster assignment: split into k equal-sized clusters
        cluster = ceil((1::n)'*k/n)

        // Compute Wk
        Wk = 0
        for (i=1; i<=n; i++) {
            for (j=i+1; j<=n; j++) {
                if (cluster[i,1]==cluster[j,1]) Wk = Wk + D[i,j]
            }
        }

        // Add safely using single index
        if (Wk>0) logWk_sim_total[k] = logWk_sim_total[k] + ln(Wk)
    }
}

// Return cumulative logWk to Stata
st_matrix("logWk_sim_total", logWk_sim_total)
end
*-----------------------------------------------
* View cumulative simulated logWk
*-----------------------------------------------
matrix list logWk_sim_total

With the above mentioned lines of code I get this error:

  <istmt>:  3200  conformability error (3 lines skipped)

right after this:

 // Add safely using single index
>         if (Wk>0) logWk_sim_total[k] = logWk_sim_total[k] + ln(Wk)
>     }
> }

I have noticed that the first step of the analysis works only if the max number of clusters to be tested is inputed in numerical form and I am applying this very same rule also to the second step of the analysis.
I wish I could perform then the comparison between real data and simulated data using these final lines of code (step 3):

matrix list logWk_sim_total

* Convert cumulative logWk_sim_total to average over B
matrix logWk_sim_avg = logWk_sim_total / `B'

* Retrieve logWk_real from results matrix
matrix logWk_real = J(`K_max',1,.)
forvalues k = 1/`K_max' {
    logWk_real[`k',1] = results[`k',2]
}

* Compute Gap = average logWk_sim - logWk_real
matrix Gap = logWk_sim_avg - logWk_real

* Optional: store results together
matrix final = logWk_real, logWk_sim_avg, Gap
matrix colnames final = "logWk_real" "logWk_sim_avg" "Gap"

* View Gap statistic table
matrix list final

* Plot Gap vs k
twoway (line Gap 1/`K_max'), xlabel(1/`K_max') ///
    title("Gap Statistic by Number of Clusters") ///
    ytitle("Gap") xtitle("Number of clusters (k)")

I am very new to MATA and I would really appreciate if anybody could help me identifying which commands should be changed/improved to perform a complete and successful gap statistic using Stata 19.
Many thanks in advance.
Anna

[Choonjoo Lee]

Dear Anna Novaresio

The conformability error is only one of several issues. Your B, K_max, and n are Stata local macros, whereas st_numscalar() reads named Stata scalars. You could use strtoreal(st_local("B")) inside Mata, or define a Stata scalar first. Also, mean() of a row vector returns a row vector, which cannot fill one D[i,j], and (1::n)' gives a row rather than a column of cluster labels.

More fundamentally, the first step assigns everyone to one cluster for every k, and the second assigns clusters by observation order. Neither computes Ward groups. Clearing all between steps also destroys the matrices needed for the comparison. Finally, the paper's within-cluster W_k includes a cluster-size divisor; the raw pairwise sum in your code omits it. The graph command also needs data variables, rather than a Stata matrix.

The code below creates 40 synthetic observations, so you can run it unchanged first. For your own analysis, replace the demo-data block with a use command pointing to your dataset and edit the variable lists. It uses the same Gower distance and Ward tree on each real and reference dataset. Each 0/1 variable is sampled independently from its observed levels; each continuous variable is sampled uniformly over its observed range. That defines a product-marginal reference distribution. B=100 gives a less noisy Monte Carlo standard error than B=10. I ran the displayed settings (K=10, B=100) in Stata 16.1; your Stata 19 dataset has not been available to test.

version 16.1
set more off
clear all
* Self-contained demo data: four 0/1 variables and five continuous ones.
* For your analysis, replace this demo block with: use "full/path/mydb.dta", clear
set seed 84327
set obs 40
generate byte demo_group = _n>20
forvalues j=1/4 {
    generate byte cat_var_`j' = runiform() < (0.15 + 0.7*demo_group)
}
forvalues j=1/5 {
    generate double num_var_w_`j' = demo_group + rnormal(0,.3)
}
drop demo_group

* For real data, also edit these lists. catvars must be 0/1.
local catvars "cat_var_1 cat_var_2 cat_var_3 cat_var_4"
local contvars "num_var_w_1 num_var_w_2 num_var_w_3 num_var_w_4 num_var_w_5"
local vars "`catvars' `contvars'"
local K 10
local B 100
set seed 20261025

* Use the same complete observations for every real/reference calculation.
egen byte nmiss = rowmiss(`vars')
keep if nmiss==0
drop nmiss
assert _N > `K'
foreach v of local catvars {
    assert inlist(`v', 0, 1)
}

* W_k = sum over within-cluster pairs of Gower distance / cluster size.
* This adapts Tibshirani et al.'s pairwise form to Gower dissimilarity.
mata:
real scalar gap_within(string scalar dname, string scalar gname)
{
    real matrix D
    real colvector g, nr
    real scalar n, i, j, w
    D = st_matrix(dname)
    g = st_data(., gname)
    n = rows(D)
    nr = J(n, 1, 0)
    for (i=1; i<=n; i++) nr[i] = sum(g:==g[i])
    w = 0
    for (i=1; i<=n; i++) {
        for (j=i+1; j<=n; j++) {
            if (g[i]==g[j]) w = w + D[i,j]/nr[i]
        }
    }
    return(w)
}

void gap_resample_binary(string scalar names)
{
    string rowvector v
    real matrix C
    real scalar n, j
    v = tokens(names)
    C = st_data(., v)
    n = rows(C)
    for (j=1; j<=cols(C); j++) {
        st_store(., v[j], C[ceil(n*runiform(n,1)), j])
    }
}
end

capture program drop gapward
program define gapward
    syntax varlist, K(integer) OUT(name)
    quietly matrix dissimilarity G_gap = `varlist', Gower
    quietly clustermat wardslinkage G_gap, add name(gap_tree)
    quietly cluster generate gapg = groups(1/`k'), name(gap_tree)
    matrix `out' = J(`k', 1, .)
    forvalues j=1/`k' {
        mata: st_numscalar("gap_w", gap_within("G_gap", "gapg`j'"))
        assert scalar(gap_w)>0
        matrix `out'[`j',1] = ln(scalar(gap_w))
        drop gapg`j'
    }
    cluster drop gap_tree
end

* Observed log W_k: one Ward tree, cut at every k.
gapward `vars', k(`K') out(logW_real)
matrix logW_ref = J(`K', `B', .)

* Reference null: sample each binary column independently from its observed
* levels, and each continuous column uniformly over its observed range.
* Then repeat the SAME Gower + Ward + W_k procedure.
forvalues b=1/`B' {
    preserve
    mata: gap_resample_binary("`catvars'")
    foreach v of local contvars {
        quietly summarize `v', meanonly
        local lo = r(min)
        local hi = r(max)
        quietly replace `v' = `lo' + (`hi'-`lo')*runiform()
    }
    quietly gapward `vars', k(`K') out(logW_b)
    forvalues k=1/`K' {
        matrix logW_ref[`k',`b'] = logW_b[`k',1]
    }
    restore
}

* Gap(k) = mean_b log W*_kb - log W_k; s_k includes sqrt(1+1/B).
mata:
L = st_matrix("logW_real")
S = st_matrix("logW_ref")
R = J(rows(L),5,.)
for (k=1; k<=rows(L); k++) {
    v = S[k,.]'
    m = mean(v)
    se = sqrt((1+1/rows(v))*sum((v:-m):^2)/(rows(v)-1))
    R[k,.] = (k,L[k,1],m,m-L[k,1],se)
}
st_matrix("gap_result",R)
end
matrix colnames gap_result = k logW_real mean_logW_ref Gap se
matrix list gap_result

* Tibshirani's 1-SE rule: first k satisfying Gap(k) >= Gap(k+1)-se(k+1).
local chosen `K'
forvalues k=1/`=`K'-1' {
    local next = `k'+1
    if `chosen'==`K' & gap_result[`k',4] >= ///
        gap_result[`next',4] - gap_result[`next',5] local chosen `k'
}
display as result "First k satisfying 1-SE rule: `chosen'"

preserve
clear
svmat double gap_result, names(col)
twoway connected Gap k, xlabel(1(1)`K') ///
    xtitle("Number of clusters") ytitle("Gap statistic") ///
    title("Gower + Ward: gap statistic")
graph export "gap_gower_ward.png", replace width(1600)
restore
display as result "GAP_DEMO_COMPLETE"

One methodological caution: Stata permits Ward linkage with Gower dissimilarity, but its manual warns that Ward's usual minimum-variance interpretation is tied to squared Euclidean distance. Thus I would describe the result as a gap calculation for the specified Gower-plus-Ward procedure, with the reference distribution stated explicitly. It estimates k relative to that null; it is not a p-value for cluster validity. Stata's Gower treats 0/1 columns as binary and other numeric columns as continuous. For a nominal variable with three or more levels, encode the factor appropriately and sample the original factor jointly before creating its indicators. Independently sampling its indicator columns could create impossible combinations.

Sources: Tibshirani, Walther, and Hastie (2001), https://doi.org/10.1111/1467-9868.00293 ; Stata's Gower definition, https://www.stata.com/manuals/mvmeasure_option.pdf ; Stata's Ward warning and linkage syntax, https://www.stata.com/manuals/mvclusterlinkage.pdf .