当前位置：首页 > article >正文

层次聚类R复现

article 2025/3/4 13:07:56

在模型建立的过程中，我们一般都要把许多X变量进行聚类，进而分析模型的实际意义。

其中大概可分为3类

1. K-Means聚类

2. 层次聚类

3. DBSCAN聚类

在这里，我们以层次聚类为例子进行分享

一安装所需包

#================================================================
# R内置数据集变量聚类分析示例
# 目的：通过层次聚类识别冗余变量并从每个聚类选择代表性变量
#================================================================

# 安装并加载所需包
if(!require(cluster)) install.packages("cluster")
if(!require(corrplot)) install.packages("corrplot")
if(!require(dendextend)) install.packages("dendextend")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(dplyr)) install.packages("dplyr")
if(!require(factoextra)) install.packages("factoextra")
if(!require(MASS)) install.packages("MASS")  # 包含Boston数据集

library(cluster)# 聚类分析
library(corrplot)# 可视化相关性矩阵
library(dendextend)# 树状图可视化
library(ggplot2)# 数据可视化
library(dplyr)# 数据处理
library(factoextra)# 最佳聚类数量选择
library(MASS)# 包含Boston数据集

二载入和介绍数据包

#================================================================
# 1. 数据准备
#================================================================

# 加载数据 (Boston数据集包含波士顿郊区房价及相关特征)
data(Boston)

# 查看数据结构
str(Boston)

# 清理数据 (移除缺失值)
boston_clean <- na.omit(Boston)

# 数据简单描述
summary(boston_clean)

# Boston数据集变量解释
# crim: 犯罪率 per capita by town
# zn: 住宅用地比例 for lots over 25,000 sq. ft.
# indus: 非零售业务用地比例 per town
# chas: Charles River（贫民区） dummy variable (= 1 if tract bounds river; 0 otherwise)
# nox: 氮氧化物浓度 (parts per 10 million)
# rm: 每个住宅的平均房间数
# age: 1940年之前建成的自住单位比例
# dis: 距离五个波士顿就业中心的加权距离
# rad: 距离高速公路的便利指数
# tax: 每 10,000 美元的全值财产税率
# ptratio: 每个镇的师生比例
# black: 1000(Bk - 0.63)^2，其中 Bk 是每个镇的黑人比例
# lstat: 低收入人口比例
# medv: 自住房的中位数价值 (以千美元计)

使用较为经典的Boston数据包，包含波士顿房价等因素

如图所示

三变量相关性分析

#================================================================
# 2. 相关性分析
#================================================================

# 计算变量间的相关性矩阵
corr_matrix <- cor(boston_clean, use = "pairwise.complete.obs")

# 可视化相关性矩阵
corrplot(corr_matrix, 
         method = "color", 
         type = "upper", 
         order = "hclust", 
         tl.col = "black", 
         tl.cex = 0.7,
         title = "波士顿房价数据变量相关性矩阵",
         mar = c(0,0,2,0))

如图，得到了不同变量之间的相关性，我们的目的就是基于这些相关性，将变量分为不同的组

四层次聚类分析

#================================================================
# 3. 层次聚类分析
#================================================================

# 计算变量间的距离矩阵 (基于1-|相关系数|)
dist_matrix <- as.dist(1 - abs(corr_matrix))

# 应用层次聚类
hc <- hclust(dist_matrix, method = "complete")

# 转换为树状图对象
dend <- as.dendrogram(hc)

# 可视化树状图
par(mar = c(5, 4, 4, 10)) #控制绘图边距
plot(dend, 
     main = "波士顿房价数据变量层次聚类",
     ylab = "高度",
     horiz = TRUE)

可以看到，我们依据变量之间的相关性，绘制出树状图。但是我们如何确定最佳的聚类数呢？这就引出了不同的筛选标准

五最佳聚类数筛选标准手肘法则&AIC判别

手肘法则

手肘法则：WSS（簇内方差平方和）关于聚类数k的曲线，选择梯度发生变化的那个点，如图所示

#================================================================
# 4. 确定最佳聚类数量
#================================================================

# 使用肘部法则确定最佳聚类数量
wss <- fviz_nbclust(t(as.matrix(boston_clean)), 
                   FUNcluster = hcut,
                   method = "wss", 
                   k.max = 10)
print(wss)

在这里可以把k选作2、也可以选为3

AIC判别法则

AIC判别：AIC（Akaike Information Criterion）即赤池信息准则，其计算公式为： $AIC = 2k - 2ln(L)$

K：模型内自由参数的数量

ln（L）：似然函数的自然对数

AIC越小，说明效果越好

# 测试不同聚类数量的AIC
k_max <- 10
aic_values <- numeric(k_max)

for(k in 1:k_max) {
  clusters <- cutree(hc, k = k)
  
  # 计算每个聚类的组内方差
  wss_cluster <- 0
  for(i in 1:k) {
    if(sum(clusters == i) > 0) {
      cluster_vars <- names(clusters[clusters == i])
      if(length(cluster_vars) > 1) {
        cluster_data <- boston_clean[, cluster_vars]
        wss_cluster <- wss_cluster + sum(scale(cluster_data, center = TRUE, scale = FALSE)^2)
      }
    }
  }
  
  # AIC = WSS + 2 * k * p (p是参数数量，这里简化为变量数)
  aic_values[k] <- wss_cluster + 2 * k * ncol(boston_clean)
}

# 绘制AIC图
# 使用ggplot2创建更美观的AIC图
aic_df <- data.frame(k = 1:k_max, aic = aic_values)
min_aic_k <- which.min(aic_values)

ggplot(aic_df, aes(x = k, y = aic)) +
    geom_line(color = "#3366CC", linewidth = 1) +
    geom_point(color = "#3366CC", size = 3) +
    geom_point(data = aic_df[min_aic_k, ], aes(x = k, y = aic), 
                         color = "#CC3333", size = 4, shape = 18) +
    geom_text(data = aic_df[min_aic_k, ], 
                        aes(label = paste("最小AIC:", round(aic, 1), "(k=", k, ")")),
                        hjust = 0.5, vjust = -1, color = "#CC3333") +
    labs(title = "基于AIC确定最佳聚类数量",
             subtitle = "较低的AIC值表示更好的模型平衡度",
             x = "聚类数量 (k)",
             y = "AIC 值") +
    scale_x_continuous(breaks = 1:k_max) +
    theme_minimal() +
    theme(
        plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 11, color = "darkgrey"),
        axis.title = element_text(face = "bold"),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(color = "gray90"),
        panel.border = element_rect(fill = NA, color = "gray80")
    )

# 基于肘部法则和AIC分析确定最佳聚类数
# 观察肘部法则图形和AIC图后手动设置
# 使用factoextra多种方法确定最佳聚类数量

在这里，可将k选为7

optimal_k_elbow <- 2  # 这里选择3个聚类作为示例
optimal_k_AIC <- 7  # 这里选择4个聚类作为示例

对手肘法则的最佳聚类结果进行展示

#================================================================
# 5A. 基于肘部法则的聚类划分
#================================================================

# 基于肘部法则切割树状图
clusters_elbow <- cutree(hc, k = optimal_k_elbow)

# 为每个聚类着色
dend_colored_elbow <- dend %>%
  set("branches_k_color", k = optimal_k_elbow) %>%
  set("branches_lwd", 2)

# 可视化聚类结果
par(mar = c(5, 4, 4, 10))
plot(dend_colored_elbow, 
     main = paste0("肘部法则聚类结果 (k = ", optimal_k_elbow, ")"),
     horiz = TRUE)

# 查看每个聚类的变量
cluster_vars_elbow <- list()
for(i in 1:optimal_k_elbow) {
  cluster_vars_elbow[[i]] <- names(clusters_elbow[clusters_elbow == i])
  cat("\n肘部法则聚类", i, "包含的变量:\n")
  print(cluster_vars_elbow[[i]])
}

对AIC判别法则的最佳聚类结果进行展示

#================================================================
# 5B. 基于AIC的聚类划分
#================================================================

# 基于AIC切割树状图
clusters_aic <- cutree(hc, k = optimal_k_AIC)

# 为每个聚类着色
dend_colored_aic <- dend %>%
  set("branches_k_color", k = optimal_k_AIC) %>%
  set("branches_lwd", 2)

# 可视化聚类结果
par(mar = c(5, 4, 4, 10))
plot(dend_colored_aic, 
     main = paste0("AIC聚类结果 (k = ", optimal_k_AIC, ")"),
     horiz = TRUE)

# 查看每个聚类的变量
cluster_vars_aic <- list()
for(i in 1:optimal_k_AIC) {
  cluster_vars_aic[[i]] <- names(clusters_aic[clusters_aic == i])
  cat("\nAIC聚类", i, "包含的变量:\n")
  print(cluster_vars_aic[[i]])
}