层次聚类R复现

在模型建立的过程中，我们一般都要把许多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 plete.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）即赤池信息准则，其计算公式为：

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]]) }

具体两种法则如何选择，可以分别进行，然后通过最后展示结果进行人工选择