在模型建立的过程中，我们一般都要把许多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）即赤池信息准则，其计算公式为：
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 <- 0for(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]])
}
具体两种法则如何选择，可以分别进行，然后通过最后展示结果进行人工选择