Sample Usage
We follow the step-by-step demonstration from the Appendix A of Ghashti, Hare, and Thompson (2026). The example is a bivariate mixture with two well-separated clusters and a protected attribute that is imbalanced within each cluster.
- First we generate the demonstration data. There are two clusters
centreed at
(-3, -3)and(3, 3), and two protected groups offset in opposite directions on the second coordinate. The group counts are asymmetric within each cluster, one hundred against two hundred, so the overall group marginals are balanced while the baseline clustering is not fair.
set.seed(42)
demoData <- rbind(
data.frame(x1 = rnorm(100, -3, 1), x2 = rnorm(100, -3 - 0.25, 1), cluster = 1L, group = 0L),
data.frame(x1 = rnorm(200, -3, 1), x2 = rnorm(200, -3 + 0.25, 1), cluster = 1L, group = 1L),
data.frame(x1 = rnorm(200, 3, 1), x2 = rnorm(200, 3 - 0.25, 1), cluster = 2L, group = 0L),
data.frame(x1 = rnorm(100, 3, 1), x2 = rnorm(100, 3 + 0.25, 1), cluster = 2L, group = 1L)
)
dataMat <- as.matrix(demoData[, c("x1", "x2")])
groupVec <- demoData$group
trueCluster <- demoData$cluster
cat("n =", nrow(dataMat),
" group counts =", paste(table(groupVec), collapse = "/"),
" cluster counts =", paste(table(trueCluster), collapse = "/"), "\n")
#> n = 600 group counts = 300/300 cluster counts = 300/300We plot the data by protected group. The two groups overlap heavily in each cluster, so the imbalance is not visually obvious, yet it is enough to make an unconstrained clustering unfair.
ggplot(demoData, aes(x1, x2, shape = factor(group), fill = factor(group))) +
geom_point(size = 2, colour = "black", stroke = 0.3, alpha = 0.7) +
scale_shape_manual("Group", values = c("0" = 21, "1" = 24)) +
scale_fill_manual("Group", values = c("0" = "#8ABF69", "1" = "#D08890")) +
labs(x = expression(x[1]), y = expression(x[2])) +
coord_fixed() + theme_bw() +
theme(panel.grid = element_blank(), legend.position = "bottom")
- We construct the transport interpolation.
buildTransport()computes the Wasserstein-2 barycenter on the full sample and the group transport maps, and returns a closurefnthat produces the transported data at any intensity.
alphaVec <- resolveAlpha("uniform", groupVec, sort(unique(groupVec)))
transport <- buildTransport(dataMat, groupVec, alphaVec, verbose = FALSE)
cat("barycenter atoms =", nrow(transport$barycenter),
" converged =", transport$baryConverged,
" iterations =", transport$baryIter, "\n")
#> barycenter atoms = 600 converged = TRUE iterations = 4We can compare the soft-fairness violation of a baseline fuzzy c-means fit on the original data against a fit on the fully transported data at intensity one.
baseFit <- fcm(dataMat, numClusters = 2, numStart = 5)
fullFit <- fcm(transport$fn(1), numClusters = 2, numStart = 5)
cat("Delta_soft at t = 0:", round(softViolation(baseFit$membership, groupVec), 3), "\n")
#> Delta_soft at t = 0: 0.154
cat("Delta_soft at t = 1:", round(softViolation(fullFit$membership, groupVec), 3), "\n")
#> Delta_soft at t = 1: 0.001- We now run the full centroid-family procedure with
fitSKM(), sweeping the transport intensity on a grid at the tolerancedeltaFair = 0.05. The returned object reports the selected intensitytOptimal, the baseline and final violations, and the swept diagnostics inhistory.
set.seed(1)
fitCentroid <- fitSKM(dataMat, groupVec, numClusters = 2, deltaFair = 0.05,
tSeq = seq(0, 1, by = 0.02), verbose = FALSE)
cat("t* =", fitCentroid$tOptimal,
" Delta_soft:", round(fitCentroid$violationBaseline, 3),
"->", round(fitCentroid$violation, 3), "\n")
#> t* = 0.82 Delta_soft: 0.154 -> 0.046- The
historydata frame records the soft and hard violations at every grid point. Plotting the soft-violation trace against the intensity shows where it first crosses the tolerance, which is the selectedtOptimal.
hist <- fitCentroid$history
ggplot(hist, aes(t, violationSoft)) +
geom_line() + geom_point(size = 1) +
geom_hline(yintercept = 0.05, linetype = "dashed", colour = "#E31A1C") +
geom_vline(xintercept = fitCentroid$tOptimal, linetype = "dotted", colour = "grey30") +
labs(x = expression(t), y = expression(Delta[soft](t))) +
theme_bw() + theme(panel.grid = element_blank())
- We can count how many observations are reassigned between the baseline and the fair solution, after aligning the two labelings so that cluster identities match. This is the practical footprint of the fairness intervention.
alignLabels <- function(current, reference) {
overlap <- table(current, reference)
mapping <- apply(overlap, 1, which.max)
as.integer(mapping[as.character(current)])
}
baseLabels <- fitCentroid$clustersBaseline
fairLabels <- alignLabels(fitCentroid$clusters, baseLabels)
cat("reassigned:", sum(fairLabels != baseLabels),
"of", length(baseLabels),
sprintf("(%.1f%%)", 100 * mean(fairLabels != baseLabels)), "\n")
#> reassigned: 37 of 600 (6.2%)- Finally we visualize the fair clustering in the original coordinates, with observations colored by their fair cluster label and shaped by protected group.
plotDF <- data.frame(x1 = dataMat[, 1], x2 = dataMat[, 2],
cluster = factor(fairLabels), group = factor(groupVec))
ggplot(plotDF, aes(x1, x2, shape = group, fill = cluster)) +
geom_point(size = 2, colour = "black", stroke = 0.3) +
scale_shape_manual("Group", values = c("0" = 21, "1" = 24)) +
scale_fill_manual("Cluster", values = c("1" = "#4E9BC7", "2" = "#F4A460")) +
labs(x = expression(x[1]), y = expression(x[2])) +
coord_fixed() + theme_bw() +
theme(panel.grid = element_blank(), legend.position = "bottom") +
guides(fill = guide_legend(override.aes = list(shape = 22)),
shape = guide_legend(override.aes = list(fill = "grey60")))