2.1 Numeric data

For purely numeric inputs, pass the matrix or data frame directly to cluster(). Centering and scaling are available via center and scale:

X <- iris[, 1:4]  # 150 x 4 numeric matrix
head(X)
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width
#> 1          5.1         3.5          1.4         0.2
#> 2          4.9         3.0          1.4         0.2
#> 3          4.7         3.2          1.3         0.2
#> 4          4.6         3.1          1.5         0.2
#> 5          5.0         3.6          1.4         0.2
#> 6          5.4         3.9          1.7         0.4

2.2 Mixed-type data

When the data contain categorical variables, two strategies are available depending on the clustering method:

One-hot encoding (for k-means, PAM, DBSCAN, GMM):

iris_mixed <- iris          # 4 numeric + 1 factor column
iris_mixed$Group <- as.character(
  ifelse(iris$Sepal.Length > 5.8, "large", "small")
)
Xenc <- prepare_mixed_data(iris_mixed[, c(1:4, 6)], center = TRUE, scale = TRUE)
dim(Xenc)
#> [1] 150   6

Gower distance (for hierarchical methods):

d_gower <- mixed_distance(iris_mixed[, c(1:4, 6)])
class(d_gower)
#> [1] "dissimilarity" "dist"

The Gower dissimilarity between observations \(i\) and \(j\) is \[ d_G(i,j) = \frac{\sum_{f=1}^{p} w_{ijf}\,\delta_{ijf}\,s_{ijf}} {\sum_{f=1}^{p} w_{ijf}\,\delta_{ijf}}, \] where \(\delta_{ijf}\) flags comparable pairs, \(w_{ijf}\) is an optional feature weight, and \(s_{ijf}\) is the feature-specific partial similarity (scaled absolute difference for numeric, mismatch indicator for categorical).

3.1 K-means

Lloyd’s algorithm (Lloyd, 1982) minimizes the total within-cluster sum of squares: \[ \min_{\mathcal{C}} \sum_{j=1}^{k} \sum_{x_i \in C_j} \|x_i - \mu_j\|^2, \qquad \mu_j = \frac{1}{|C_j|} \sum_{x_i \in C_j} x_i. \] The algorithm alternates between nearest-centroid assignment and centroid recomputation until convergence. Multiple restarts (nstart) reduce sensitivity to random initialization.

fit_km <- cluster(X, method = "kmeans", k = 3, seed = 1)
fit_km
#> <cluster_fit>
#>   Method: kmeans
#>   Observations: 150
#>   Clusters: 3

3.2 PAM

Partitioning Around Medoids (Kaufman and Rousseeuw, 1990) selects \(k\) actual observations (medoids) \(m_1, \ldots, m_k\) to minimize \[ \min_{m_1,\ldots,m_k \in \mathcal{X}} \sum_{i=1}^{n} \min_{g} d(x_i, m_g). \] Unlike k-means, medoids are always observed data points, making PAM more robust to outliers.

if (requireNamespace("cluster", quietly = TRUE)) {
  fit_pam <- cluster(X, method = "pam", k = 3)
  fit_pam
}
#> <cluster_fit>
#>   Method: pam
#>   Observations: 150
#>   Clusters: 3

3.3 Hierarchical clustering

"hclust" and "agnes" build a dendrogram by successively merging the two closest clusters. The inter-cluster distance depends on the linkage:

The final partition is obtained by cutting the tree at height \(k\).

# From a distance object (supports mixed data via mixed_distance())
d <- mixed_distance(X)
fit_hc <- cluster(d, method = "hclust", k = 3)
fit_hc
#> <cluster_fit>
#>   Method: hclust
#>   Observations: 150
#>   Clusters: 3
plot_dendrogram(fit_hc)

AGNES (Kaufman and Rousseeuw, 1990) is the cluster package’s implementation and additionally reports the agglomerative coefficient:

if (requireNamespace("cluster", quietly = TRUE)) {
  fit_ag <- cluster(d, method = "agnes", k = 3)
  fit_ag
}
#> <cluster_fit>
#>   Method: agnes
#>   Observations: 150
#>   Clusters: 3

3.4 DBSCAN

DBSCAN (Ester et al., 1996) groups observations reachable through dense neighborhoods without requiring \(k\). A point \(p\) is a core point when \[ |N_\varepsilon(p)| = |\{q : d(p, q) \le \varepsilon\}| \ge \mathrm{MinPts}. \] Points not reachable from any core point are labeled noise (cluster 0).

if (requireNamespace("dbscan", quietly = TRUE)) {
  fit_db <- cluster(X, method = "dbscan", eps = 0.5, minPts = 5)
  fit_db
  # Noise points
  sum(clusters(fit_db) == 0)
}
#> [1] 17

3.5 Gaussian mixture models

A GMM models the data density as a weighted sum of Gaussians: \[ p(x) = \sum_{g=1}^{k} \pi_g \,\mathcal{N}(x \mid \mu_g, \Sigma_g). \] Parameters \((\pi_g, \mu_g, \Sigma_g)\) are estimated by the EM algorithm (Fraley and Raftery, 2002). Hard labels come from the MAP rule \(\hat{y}_i = \arg\max_g \pi_g \mathcal{N}(x_i \mid \mu_g, \Sigma_g)\); soft memberships are stored in membership(fit).

if (requireNamespace("mclust", quietly = TRUE)) {
  fit_gmm <- cluster(X, method = "gmm", k = 3, seed = 1)
  fit_gmm
  # Soft membership probabilities
  head(membership(fit_gmm))
}
#>      [,1]         [,2]         [,3]
#> [1,]    1 4.916819e-40 3.345948e-29
#> [2,]    1 5.846760e-29 1.599452e-23
#> [3,]    1 2.486168e-33 2.724243e-25
#> [4,]    1 2.050996e-28 2.180764e-22
#> [5,]    1 8.283066e-42 8.710458e-30
#> [6,]    1 2.118466e-41 1.791499e-29

3.6 K-prototypes

K-prototypes (Huang, 1998) extends k-means to mixed data by combining squared Euclidean distance on numeric features with a Hamming mismatch penalty on categorical features: \[ D(x_i, p_g) = \underbrace{\sum_{j \in \mathcal{N}} (x_{ij} - \mu_{gj})^2}_{\text{numeric}} +\, \lambda \underbrace{\sum_{j \in \mathcal{C}} I(x_{ij} \neq \nu_{gj})}_{\text{categorical}}. \] When lambda = NULL, \(\lambda\) is estimated automatically from the data.

if (requireNamespace("clustMixType", quietly = TRUE)) {
  fit_kp <- cluster(iris, method = "kproto", k = 3, seed = 1)
  fit_kp
  # Mixed prototypes (centroids + modes)
  prototypes(fit_kp)
}
#>           Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
#> cluster_1      0.07235    -0.29260       0.5067      0.1340 versicolor
#> cluster_2      0.77910    -0.07366       1.8150      0.8333  virginica
#> cluster_3     -0.83730     0.37070      -2.2960     -0.9533     setosa

3.7 KMM (K-Mixed-Modes)

KMM is phynotype’s native mixed-data algorithm. It minimizes the same weighted prototype distance as k-prototypes: \[ D(x_i, p_g) = \sum_{j \in \mathcal{N}} (x_{ij} - \mu_{gj})^2 + \lambda \sum_{j \in \mathcal{C}} I(x_{ij} \neq \nu_{gj}), \] with \(\lambda \ge 0\) estimated automatically (as the median numeric variance) when not supplied. KMM uses multiple random restarts and is designed for biostatistical applications where the number of categories is small relative to the number of observations.

iris_chr <- iris
iris_chr$Species <- as.character(iris_chr$Species)

fit_kmm <- cluster(iris_chr, method = "kmm", k = 3, seed = 1)
fit_kmm
#> <cluster_fit>
#>   Method: kmm
#>   Observations: 150
#>   Clusters: 3
summary(fit_kmm)
#> Cluster fit summary
#>   Method: kmm
#>   Observations: 150
#>   Clusters: 3
#>   Sizes: 1=50, 2=50, 3=50
prototypes(fit_kmm)
#>           Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
#> cluster_1        5.006       3.428        1.462       0.246     setosa
#> cluster_2        6.668       3.016        5.554       2.014  virginica
#> cluster_3        5.856       2.728        4.258       1.338 versicolor

5.1 Internal metrics

val <- validate(fit_km)
val$metrics_table
#>              metric    value               scale        direction
#> 1        silhouette   0.5528             -1 to 1 higher is better
#> 2 calinski_harabasz 561.6000 positive, unbounded higher is better
#> 3    davies_bouldin   0.6620 positive, unbounded  lower is better
#> 4      total_within  78.8500 positive, unbounded  lower is better
#> 5     bootstrap_ari   0.9710                <NA>             <NA>

The five internal metrics each capture a different aspect of cluster quality.

Silhouette width (Rousseeuw, 1987) measures how similar each observation is to its own cluster compared to the nearest alternative cluster. For observation \(i\), let \(a(i)\) be its mean distance to all other members of its cluster and \(b(i)\) be its mean distance to the closest other cluster. The silhouette width is: \[ s(i) = \frac{b(i) - a(i)}{\max\{a(i),\, b(i)\}}. \] This quantity lies in \([-1, 1]\). A value near \(+1\) means \(i\) is well inside its cluster; near \(0\) means it sits on a cluster boundary; negative means it is closer on average to a neighboring cluster than to its own. validate() reports the mean over all observations, a global summary of cluster cohesion and separation simultaneously.

Calinski-Harabasz index (Caliński and Harabasz, 1974) treats clustering as an ANOVA-like problem: it asks how large the between-cluster dispersion is relative to the within-cluster dispersion, normalizing by the corresponding degrees of freedom: \[ \mathrm{CH} = \frac{\mathrm{BSS} / (k-1)}{\mathrm{WSS} / (n-k)}, \] where \(\mathrm{BSS}\) is the between-cluster sum of squares and \(\mathrm{WSS}\) is the total within-cluster sum of squares. The ratio is unbounded above; higher values indicate tighter, more separated clusters. This index is most useful when comparing solutions with different values of \(k\).

Davies-Bouldin index (Davies and Bouldin, 1979) penalizes clusters that are internally diffuse and close to their neighbors. For each cluster \(j\), it finds the worst-case neighbor (the one that is both the most scattered and the closest) and averages these worst cases: \[ \mathrm{DB} = \frac{1}{k} \sum_{j=1}^{k} \max_{l \ne j} \frac{s_j + s_l}{d(\mu_j, \mu_l)}, \] where \(s_j\) is the mean distance from the cluster-\(j\) members to centroid \(\mu_j\), and \(d(\mu_j, \mu_l)\) is the Euclidean distance between centroids. Lower values are better; a perfect score of 0 would require clusters with zero internal scatter.

Total within-cluster sum of squares is the raw objective minimized by k-means, and is included as a direct measure of compactness: \[ \mathrm{WSS} = \sum_{j=1}^{k} \sum_{i \in C_j} \|x_i - \mu_j\|^2. \] It is most interpretable when compared across solutions with the same \(k\): the “elbow” in a WSS-vs-\(k\) plot signals a natural cluster count.

Bootstrap ARI (Fang and Wang, 2012) addresses stability rather than geometric quality. It refits the same clustering model on \(B\) bootstrap resamples of the data and asks: how similar are the resulting partitions to the original? Formally, \[ \mathrm{BootARI} = \frac{1}{B} \sum_{b=1}^{B} \mathrm{ARI}\!\left(\hat{y}_{I^{(b)}},\, \hat{y}^{(b)}\right), \] where \(I^{(b)}\) is the set of indices drawn in bootstrap \(b\), \(\hat{y}^{(b)}\) is the partition refitted on those indices, and ARI is evaluated only on the sampled observations. Values near 1 mean the cluster structure is robust to resampling; values below 0.6 suggest the solution is sensitive to the particular data points included.

5.2 External metrics

When reference labels are available:

val_ext <- validate(fit_km, truth = iris$Species)
val_ext$metrics_table
#>              metric    value               scale        direction
#> 1        silhouette   0.5528             -1 to 1 higher is better
#> 2 calinski_harabasz 561.6000 positive, unbounded higher is better
#> 3    davies_bouldin   0.6620 positive, unbounded  lower is better
#> 4      total_within  78.8500 positive, unbounded  lower is better
#> 5               ari   0.7302                <NA>             <NA>
#> 6               nmi   0.7582                <NA>             <NA>
#> 7     bootstrap_ari   0.9710                <NA>             <NA>

Adjusted Rand index (Hubert and Arabie, 1985) counts the proportion of observation pairs on which two labelings agree (both in the same cluster, or both in different clusters) and corrects that count for the level of agreement expected by chance alone. Given a contingency table \((n_{ij})\) with row sums \(a_i\) and column sums \(b_j\): \[ \mathrm{ARI}(U, V) = \frac{ \sum_{ij}\binom{n_{ij}}{2} - E }{ \tfrac{1}{2}\!\left[\sum_i\binom{a_i}{2}+\sum_j\binom{b_j}{2}\right] - E }, \quad E = \frac{\sum_i\binom{a_i}{2}\cdot\sum_j\binom{b_j}{2}}{\binom{n}{2}}. \] A value of 1 indicates perfect agreement; the expected value under random labelings is 0; negative values are possible when agreement is worse than chance.

Normalized mutual information (Strehl and Ghosh, 2002) measures how much knowing one partition reduces uncertainty about the other, normalized so that the scale is comparable across different cluster counts: \[ \mathrm{NMI}(U, V) = \frac{I(U; V)}{\sqrt{H(U)\,H(V)}}, \] where \(I(U; V) = \sum_{u,v} p_{uv} \log(p_{uv} / p_u p_v)\) is mutual information and \(H(U) = -\sum_u p_u \log p_u\) is Shannon entropy. Unlike ARI, NMI is not corrected for chance, but it is symmetric and bounded in \([0, 1]\), making it useful for reporting alongside ARI.

5.3 Grid search

To compare cluster counts, pass raw data directly:

grid_val <- validate(X, method = "kmeans", k = 2:6)
grid_val$metrics_table
#>               metric    value               scale        direction k
#> 1         silhouette   0.6810             -1 to 1 higher is better 2
#> 2  calinski_harabasz 513.9000 positive, unbounded higher is better 2
#> 3     davies_bouldin   0.4043 positive, unbounded  lower is better 2
#> 4       total_within 152.3000 positive, unbounded  lower is better 2
#> 5      bootstrap_ari   0.9973                <NA>             <NA> 2
#> 6         silhouette   0.5528             -1 to 1 higher is better 3
#> 7  calinski_harabasz 561.6000 positive, unbounded higher is better 3
#> 8     davies_bouldin   0.6620 positive, unbounded  lower is better 3
#> 9       total_within  78.8500 positive, unbounded  lower is better 3
#> 10     bootstrap_ari   0.9494                <NA>             <NA> 3
#> 11        silhouette   0.4981             -1 to 1 higher is better 4
#> 12 calinski_harabasz 530.8000 positive, unbounded higher is better 4
#> 13    davies_bouldin   0.7803 positive, unbounded  lower is better 4
#> 14      total_within  57.2300 positive, unbounded  lower is better 4
#> 15     bootstrap_ari   0.8195                <NA>             <NA> 4
#> 16        silhouette   0.4887             -1 to 1 higher is better 5
#> 17 calinski_harabasz 495.5000 positive, unbounded higher is better 5
#> 18    davies_bouldin   0.8060 positive, unbounded  lower is better 5
#> 19      total_within  46.4500 positive, unbounded  lower is better 5
#> 20     bootstrap_ari   0.8458                <NA>             <NA> 5
#> 21        silhouette   0.3648             -1 to 1 higher is better 6
#> 22 calinski_harabasz 473.9000 positive, unbounded higher is better 6
#> 23    davies_bouldin   0.9142 positive, unbounded  lower is better 6
#> 24      total_within  39.0400 positive, unbounded  lower is better 6
#> 25     bootstrap_ari   0.8523                <NA>             <NA> 6

5.4 Per-cluster silhouette widths

Individual silhouette widths reveal which observations are well-assigned (positive width, comfortably inside their cluster) and which are borderline or potentially misclassified (width near zero or negative):

if (requireNamespace("cluster", quietly = TRUE)) {
  plot_silhouette(fit_km)
}

8.1 How it works

For \(B\) candidate partitions, the co-association matrix \(C \in [0,1]^{n\times n}\) is \[ C_{ij} = \frac{1}{B} \sum_{b=1}^{B} I\!\left\{u^{(b)}_i = u^{(b)}_j\right\}, \] the proportion of candidates that co-cluster observations \(i\) and \(j\). The consensus dissimilarity \(D = 1 - C\) is then hierarchically clustered (average linkage) and the optimal \(k\) is chosen by maximizing the mean silhouette width: \[ k^* = \operatorname*{arg\,max}_{k \in \mathcal{K}} \bar{s}_k(D). \]

Ensemble stability is reported as the mean pairwise partition agreement (PPA): \[ \mathrm{PPA}(U, V) = \frac{1}{\binom{n}{2}} \sum_{i < j} I\!\left\{(u_i = u_j) = (v_i = v_j)\right\}. \]

8.2 Usage

mfit <- metacluster(
  X,
  methods = c("kmeans", "pam", "hclust"),
  k       = 2:5,
  seed    = 1
)
mfit
#> <metacluster_fit>
#>   Methods: kmeans, pam, hclust
#>   Candidate fits: 12
#>   Final clusters: 5

Candidate partition table:

head(mfit$candidate_table)
#>   candidate method k n_clusters
#> 1         1 kmeans 2          2
#> 2         2 kmeans 3          3
#> 3         3 kmeans 4          4
#> 4         4 kmeans 5          5
#> 5         5    pam 2          2
#> 6         6    pam 3          3

Ensemble stability:

mfit$stability_summary
#>                         metric mean_agreement min_agreement max_agreement
#> 1 pairwise_partition_agreement         0.8186        0.6761             1

Consensus \(k\) selection scores:

mfit$selection_summary
#>   k silhouette
#> 1 2     0.6567
#> 2 3     0.7659
#> 3 4     0.8603
#> 4 5     0.9076

8.3 Validation and exploration of the consensus

validate(mfit)$metrics_table
#>                         metric    value               scale        direction
#> 1                   silhouette   0.4926             -1 to 1 higher is better
#> 2            calinski_harabasz 494.1000 positive, unbounded higher is better
#> 3               davies_bouldin   0.8168 positive, unbounded  lower is better
#> 4 pairwise_partition_agreement   0.8186                <NA>             <NA>
plot_coassoc(mfit)

plot_consensus(mfit)

Consensus dendrogram:

plot_dendrogram(mfit)

9.1 The unsupervised interpretation problem

Classical interpretability methods (permutation importance, LIME, individual conditional profiles) were designed for supervised models where a response variable \(y\) is available. Clustering produces no such variable. The challenge is that we want to understand the clustering rule \(f: \mathcal{X} \to C\), not the cluster structure in the data.

phynotype resolves this by treating the fitted model as a black-box predictor. Every supported method with a predict() implementation maps new observations to cluster labels (and optionally to distance or membership scores). This map is well-defined, deterministic, and can be probed exactly like any classifier. The three interpretability functions below operate exclusively through this map: they never look at the raw training data or the clustering criterion directly. As a result, the explanations describe the decision boundary of the fitted model, not the geometry of the original point cloud. This distinction matters: a feature can be highly informative for separating the data yet contribute little to the fitted boundary if the algorithm did not use it, and vice versa.

All three functions require a fit trained on row-by-feature data with prediction support ("kmeans", "pam", "hclust", "agnes", "gmm", "kproto", "kmm").

9.2 Global feature importance

What it answers: which features does the fitted clustering rule rely on most, globally across all observations?

The approach adapts permutation importance (Breiman, 2001; Fisher et al., 2019) to the unsupervised setting. The idea is simple: if a feature is important to the model, randomly shuffling its values across observations should disrupt the cluster assignments. If it is irrelevant, shuffling should leave assignments unchanged.

Concretely, let \(\hat{c} = f(X)\) be the baseline partition produced by the fitted model on the evaluation data \(X\). For each feature \(j\), we independently permute its column \(M\) times and re-predict: \[ \mathrm{FI}_j = \frac{1}{M} \sum_{m=1}^{M} \frac{1}{n} \sum_{i=1}^{n} I\!\left\{\hat{c}_i \neq f\!\left(X^{(m)}_{\pi(j)}\right)_i\right\}. \] Here \(X^{(m)}_{\pi(j)}\) is the data matrix with column \(j\) replaced by a random permutation of its values. The importance \(\mathrm{FI}_j\) is the mean fraction of observations whose predicted cluster changes when column \(j\) is scrambled. A value of 0 means the model ignores feature \(j\) entirely; a value of 1 means every assignment changes when \(j\) is removed.

Two alternative metrics measure the same idea through internal scores rather than label changes. With metric = "silhouette", importance is the mean drop in silhouette width after permutation: \[ \mathrm{FI}_j^{\mathrm{sil}} = \frac{1}{M} \sum_{m=1}^{M} \left[S\!\left(f, X\right) - S\!\left(f, X^{(m)}_{\pi(j)}\right)\right], \] so a positive value means permuting \(j\) degrades within-cluster cohesion. With metric = "total_within", it is the mean increase in within-cluster dispersion, capturing the same idea from a compactness perspective.

imp <- feature_importance(fit_km, n_repeats = 5, seed = 1)
imp$summary
#>        feature importance std_error n_repeats
#> 1 Petal.Length   0.526700  0.005963         5
#> 2 Sepal.Length   0.077330  0.004522         5
#> 3  Petal.Width   0.052000  0.009286         5
#> 4  Sepal.Width   0.009333  0.003399         5
plot(imp)

The bar chart ranks features by their mean importance. In this k-means fit on iris, Petal.Length and Petal.Width dominate: permuting either of them strongly disrupts which cluster each observation is sent to, confirming that the fitted boundary is primarily organised along the petal dimensions.

9.3 Local explanations (LIME)

What it answers: for a specific observation, which features pushed it into its assigned cluster, and in which direction?

Global importance tells us which features matter on average across all observations, but it cannot explain an individual assignment. LIME (Local Interpretable Model-agnostic Explanations; Ribeiro et al., 2016) addresses this by fitting a simple linear model locally around each observation being explained.

In the supervised setting LIME explains a classifier’s predicted class. Here the “classifier” is the clustering rule \(f\). For each explained observation \(x\), we define a target \(g\{f(\cdot)\}\) that converts model output into a scalar:

• With target = "cluster" (default): \(g\{f(z)\} = I\{f(z) = c\}\), a binary indicator of whether the perturbed sample \(z\) is predicted to belong to the same cluster \(c\) as the original observation. The local linear model then identifies which features explain membership in that cluster.
• With target = "score": \(g\{f(z)\}\) is the cluster-specific membership probability (for GMM) or a normalized radial similarity score derived from distances (for centroid/prototype methods). This gives a continuous response and is appropriate when soft assignment information is meaningful.

The LIME objective is to find a sparse linear surrogate \(\hat{h}\) that approximates \(g\{f(\cdot)\}\) in a weighted neighborhood around \(x\): \[ \hat{h} = \operatorname*{arg\,min}_{h \in \mathcal{H}} \sum_{i=1}^{N} \pi_x(z_i) \left[g\!\left\{f(z_i)\right\} - h(z_i)\right]^2 + \Omega(h), \] where \(z_1, \ldots, z_N\) are perturbed samples drawn around \(x\) (numeric features are perturbed by Gaussian noise calibrated to their training standard deviation; categorical features are randomly resampled from observed levels), and \(\pi_x(z_i) = \exp(-\|x - z_i\|^2 / \sigma^2)\) is an exponential kernel that down-weights samples far from \(x\). Sparsity \(\Omega(h)\) is enforced operationally by returning only the top n_features absolute effects from the fitted weighted linear model.

The output is a signed coefficient for each feature: a positive coefficient means that increasing the feature value makes membership in cluster \(c\) more likely under the local linear approximation; a negative coefficient means the opposite. These are local effects: they describe the boundary in the immediate neighborhood of the explained observation and should not be extrapolated globally.

lx <- lime_explain(fit_km, iris[c(1, 51, 101), 1:4],
                   n_permutations = 100, n_features = 4, seed = 1)
lx$explanations[, c("observation", "feature", "estimate", "direction")]
#>    observation      feature  estimate direction
#> 1            1 Petal.Length -0.244300  negative
#> 2            1 Sepal.Length -0.025140  negative
#> 3            1  Petal.Width -0.022190  negative
#> 4            1  Sepal.Width  0.006828  positive
#> 5            2 Petal.Length -0.427600  negative
#> 6            2  Petal.Width -0.149000  negative
#> 7            2 Sepal.Length -0.058490  negative
#> 8            2  Sepal.Width -0.042810  negative
#> 9            3 Petal.Length  0.435000  positive
#> 10           3 Sepal.Length  0.160600  positive
#> 11           3  Sepal.Width  0.081570  positive
#> 12           3  Petal.Width  0.076890  positive
plot(lx)

The plot shows one panel per explained observation. Reading the bars: observation 1 (setosa, cluster 1) is explained primarily by small Petal.Length pushing it away from the other clusters; observation 51 (versicolor) and observation 101 (virginica) are separated mainly along Petal.Width and Petal.Length, consistent with the global importance result.

9.4 Ceteris paribus profiles

What it answers: how does the predicted cluster assignment (or cluster score) change as one feature varies, holding everything else fixed at its observed value for a specific observation?

LIME gives a linear summary of the local boundary. Ceteris paribus (CP) profiles give the full nonlinear picture for each feature of interest. The name comes from the Latin for “all other things being equal”: for observation \(x_i\) and feature \(j\), we slide \(j\) across a grid of values \(z_1, \ldots, z_G\) while clamping all other features at their observed values, then record the model’s prediction at each grid point: \[ \mathrm{CP}_{i,j}(z) = g\!\left\{f\!\left(x_{i,-j},\, z\right)\right\}, \] where \(x_{i,-j}\) denotes the feature vector of observation \(i\) with feature \(j\) replaced by \(z\), \(f\) is the fitted clustering rule, and \(g\) extracts either the cluster label (producing a step function, target = "cluster") or a continuous cluster-specific score (target = "score").

CP profiles answer the question: “if this particular patient had a different value for Petal.Length, would their phenotype assignment change, and at what threshold?” They are individual diagnostic curves, useful for understanding borderline cases and for communicating individual-level model behaviour to domain experts.

cp <- ceteris_paribus(
  fit_km,
  iris[c(1, 51, 101), 1:4],
  features  = c("Petal.Length", "Petal.Width"),
  grid_size = 20
)
plot(cp)

Each panel shows one feature; each line corresponds to one explained observation. The dashed vertical line marks the observation’s actual feature value, and the dot marks its baseline prediction. Step changes in the profile (for target = "cluster") reveal the exact feature thresholds at which the model switches cluster assignment, directly interpretable as phenotype decision boundaries.