计算CH指数出现负值的原因排查及代码优化咨询

Hey there! Let's break down why your Calinski-Harabasz (CH) Index is showing negative values and how to troubleshoot the issue. First off, a quick clarification: the CH Index can indeed be negative—this happens when your clustering performs worse than randomly assigning samples to groups (i.e., the between-cluster variance is smaller than the within-cluster variance). So this is a red flag that something's off with your data prep, distance metric, or clustering workflow.

Let's dig into potential issues and fixes #

1. Critical: Your data might not be formatted correctly for JSD distance

Jensen-Shannon Divergence (JSD) is designed to compare probability distributions—meaning each sample (column in your matrix) should sum to 1 (relative abundances), not raw counts. Looking at your code, you didn't include a step to convert raw counts to relative abundances, which is likely distorting your distance calculations. This is probably the biggest culprit behind bad clustering and negative CH values.

2. JSD Distance Function Improvements

Your current JSD function uses nested loops which are slow and error-prone. Here's an optimized version that includes relative abundance normalization and vectorized operations:

dist.JSD <- function(inMatrix, pseudocount = 1e-6) {
  # Convert raw counts to relative abundances (each column sums to 1)
  inMatrix <- sweep(inMatrix, 2, colSums(inMatrix), "/")
  # Add pseudocount to avoid log(0) errors
  inMatrix[inMatrix == 0] <- pseudocount
  # Re-normalize to ensure columns still sum to 1 after adding pseudocount
  inMatrix <- sweep(inMatrix, 2, colSums(inMatrix), "/")
  
  # Define KLD and JSD helper functions
  KLD <- function(x, y) sum(x * log(x/y))
  JSD <- function(x, y) sqrt(0.5 * KLD(x, (x+y)/2) + 0.5 * KLD(y, (x+y)/2))
  
  # Compute pairwise JSD using vectorized outer() instead of loops
  jsd_matrix <- outer(seq_len(ncol(inMatrix)), seq_len(ncol(inMatrix)), 
                      function(i, j) JSD(inMatrix[,i], inMatrix[,j]))
  colnames(jsd_matrix) <- colnames(inMatrix)
  rownames(jsd_matrix) <- colnames(inMatrix)
  # Convert to dist object for clustering
  as.dist(jsd_matrix, method = "dist")
}

3. Clustering Stability & CH Index Calculation

PAM (Partitioning Around Medoids) can get stuck in local minima, especially for larger k values. Running it multiple times and picking the best result (highest CH index) will give you more reliable values:

library(fpc)
set.seed(123) # Ensure reproducibility

nclusters <- rep(NA, 20)
for (k in 2:20) {
  # Run PAM 5 times to avoid local minima
  pam_runs <- lapply(1:5, function(x) pam.clustering(data.dist, k))
  # Calculate CH index for each run
  ch_scores <- sapply(pam_runs, function(cluster_labels) {
    index.G1(x = t(data), cl = cluster_labels, d = data.dist, centrotypes = "medoids")
  })
  # Keep the highest CH score for this k (best clustering)
  nclusters[k] <- max(ch_scores)
}

# Plot the results with cleaner labels
plot(nclusters, type="h", xlab="Number of clusters (k)", ylab="Calinski-Harabasz Index", 
     lwd=2, col="steelblue", main="Optimal Cluster Number via CH Index")

Key Takeaways #

• Always normalize your microbial data to relative abundances before calculating JSD—this is non-negotiable for valid distance comparisons.
• Negative CH values mean that clustering with that k is worse than random, so you can safely ignore those k values and focus on the peaks in your plot.
• Running PAM multiple times improves stability, especially for higher k values.

备注：内容来源于stack exchange，提问作者esti