基于分组变量的多变量配对聚类（1kNN）实现需求

现有一个包含id列、g分组变量及300个数值变量的tibble数据集，需完成以下操作：

• 在每个分组内对数据行进行两两配对聚类（1kNN），若分组内行数为奇数，剩余1行不参与聚类；
• 优先采用**马氏距离（Mahalanobis Distance）**进行聚类，也可接受其他合理距离方案；
• 生成簇内马氏距离作为诊断变量；
• 无需像MatchIt工具那样强制对数据行进行二元分类。

示例数据集 #

library(tibble)
set.seed(123) # 设置随机种子保证结果可复现
tibble(
  id = c(1:8),
  g = rep(c("A","B"),4),
  v1 = rnorm(8),
  v2 = rnorm(8),
  v3 = rnorm(8)
) -> obs

解决方案代码 #

library(dplyr)
library(stringr)

# 定义马氏距离矩阵计算函数（基于分组内协方差矩阵）
mahalanobis_dist_matrix <- function(data, num_vars) {
  num_mat <- as.matrix(data[num_vars])
  cov_mat <- cov(num_mat)
  
  # 处理协方差矩阵奇异的情况（添加微小扰动保证可逆）
  if(any(eigen(cov_mat)$values < 1e-6)) {
    cov_mat <- cov_mat + diag(ncol(cov_mat)) * 1e-6
  }
  
  # 生成两两之间的马氏距离矩阵
  as.matrix(dist(num_mat, method = "mahalanobis", cov = cov_mat))
}

# 分组执行配对聚类逻辑
processed_data <- obs %>%
  group_by(g) %>%
  group_modify(function(group_data, group_key) {
    row_count <- nrow(group_data)
    # 奇数行时剔除最后一行，只保留偶数行参与配对
    if(row_count %% 2 != 0) {
      group_data <- group_data %>% slice_head(n = row_count - 1)
      row_count <- row_count - 1
    }
    
    # 提取所有数值变量列
    numeric_vars <- setdiff(names(group_data), c("id", "g"))
    dist_mat <- mahalanobis_dist_matrix(group_data, numeric_vars)
    
    # 贪心配对：每次选取距离最近的未配对行组成簇
    is_paired <- logical(row_count)
    cluster_ids <- integer(row_count)
    cluster_dist <- numeric(row_count)
    current_cluster <- 1
    
    while(any(!is_paired)) {
      # 找到第一个未配对的行索引
      idx1 <- which(!is_paired)[1]
      is_paired[idx1] <- TRUE
      
      # 找到与该行距离最近的未配对行
      available_dists <- dist_mat[idx1, !is_paired]
      idx2 <- as.integer(names(available_dists)[which.min(available_dists)])
      is_paired[idx2] <- TRUE
      
      # 分配簇ID并记录簇内距离
      cluster_ids[c(idx1, idx2)] <- current_cluster
      cluster_dist[c(idx1, idx2)] <- dist_mat[idx1, idx2]
      current_cluster <- current_cluster + 1
    }
    
    # 添加聚类相关列
    group_data %>%
      mutate(
        cluster = cluster_ids,
        mahalanobis_distance = cluster_dist,
        pair = str_c(group_key$g, cluster)
      )
  }) %>%
  ungroup() %>%
  arrange(pair)

# 查看处理后的数据
print(processed_data)

代码说明 #

1. 马氏距离计算：自定义函数处理分组内的协方差矩阵，针对奇异矩阵添加微小扰动确保计算有效；
2. 分组预处理：自动识别分组内的奇数行并剔除最后一行，保证仅偶数行参与配对；
3. 贪心配对逻辑：通过循环选取距离最近的未配对行组成簇，避免强制二元分类；
4. 输出字段：生成cluster（簇ID）、mahalanobis_distance（簇内马氏距离）、pair（分组-簇唯一标识）三个新增字段，并按pair排序。

内容的提问来源于stack exchange，提问作者GiulioGCantone