聚类-R语言

R实战代码

1. 层次聚类-R语言

install.packages("flexclust",destdir = "D:\\Softwares\\R\\R-3.6.1\\packages")
install.packages("NbClust",destdir = "D:\\Softwares\\R\\R-3.6.1\\packages")
library(NbClust)
par(mfrow=c(1,1))

# 加载指定包中的指定数据集，以原名字存储
data(nutrient,package = "flexclust")   
# 将行名变为小写，= rownames() ;       列名 names(nutrient)  # 列名=colnames()
row.names(nutrient) <- tolower(row.names(nutrient))  

s1:预处理：标准化+距离矩阵

# 对变量标准化，scale(x, center = TRUE, scale = TRUE) 默认按列进行标准化w为均值0，方差为1；x:二维对象，
nutrient_scaled <- scale(nutrient)

d <- dist(nutrient_scaled)         # 计算二维表中所有行之间的距离，默认欧氏距离,
#返回一维向量，需要转换为矩阵方便查看，但是下面的层次聚类中使用的还得是转换之前的

s2: 层次聚类并可视化
对距离矩阵进行层次聚类，hclust(d,method):d:距离矩阵，method是方法，平均联动比较好

fit_average <- hclust(d,method = "average")      # 得到一个tree对象
plot(fit_average,hang = -1,            # hang:标签的悬挂方式，对齐悬挂在0的下面
     cex=0.8,main = "Average Linkage Clustering")  # cex:相对于默认字号的字号

s3:选择聚类个数（NbCluster）

NbClust：对标准化后的样本数据阵确定最佳聚类数，提供了30个用于确定聚类数量的指标，向用户提出最佳聚类方案。
NbClust(data , distance = “euclidean”, min.nc = 2, max.nc = 15, method = NULL, index = “all”, alphaBeale = 0.1)
允许的最少聚类个数，最大聚类的个数
method：使用的聚类分析方法“ ward.D”，“ ward.D2”，“ single”，“ complete”，“ average”，“ mcquitty”，“ median”，“ centroid”，“ kmeans”。

devAskNewPage(ask = T)  # 在新页面之前提示,控制（对于当前设备）在开始新的输出页面之前是否提示用户。
nc <- NbClust(nutrient_scaled,distance = "euclidean", min.nc = 2,max.nc = 15,method = "average")

table(nc$Best.nc[1,]) # 列出里面存储的第一行（26个指标的推荐类数）的频数表
# 画出频数的柱状图，直观显示推荐的类数，可以看出投票最多的聚类个数是2，3，5，15，可以挨个试试效果
barplot(table(nc$Best.nc[1,]),xlab = "Number of Clusters",ylab = "Number of Criteria",main = "Number of Clusters chosen by 26 Criteria")

s4: 获取最终聚类方案

cutree() : Cut a Tree into Groups of Data , tree可以是hclust中的树结果，通过制定划分个数k

clusters <- cutree(fit_average,k=5)     # 以5类显示结果,返回值是带标签的向量，标签是每一个样本，值为该样本对应的类别
table(clusters)                         # 显示类别和类内对应个数的频数表

# 对原始数据或标准化数据按照每行对应的类别进行分类聚合，描述不同类的中位数
aggregate(nutrient,by=list(cluster=clusters),median)             # 此处为原始度量形式
aggregate(nutrient_scaled,by=list(cluster=clusters),median)      # 此处为标准度量形式

# 结果重新绘图
plot(fit_average,hang = -1, cex=0.8,main = "Average Linkage Clustering")  # 原不加修饰的树状图
rect.hclust(fit_average,k=5)     # 在上幅图上，对原树状图叠加一个5类的分类框

rect.hclust(tree,k,border):在层次聚类的树状图的分支周围绘制矩形，以突出显示相应的簇。
tree:hclust得到的tree对象;k是类数
border是不同矩形的颜色向量

2. k-means聚类-R语言

rm(list = ls())
install.packages("rattle",destdir = "D:\\Softwares\\R\\R-3.6.1\\packages")
library(rattle)

data(wine,package = "rattle") # 载入葡萄酒数据集，第一列是分为的3类，聚类时丢弃这一列
data <- scale(wine[,-1])  # 标准化

确定聚类的个数

# 作出类内平方和对聚类数量的曲线
wssplot <- function(data,nc=15,seed=1234)    # nc:考虑的最大聚类数目，默认值15，seed：随机数种子，因为该算法是随机初始化
{
  wss <- (nrow(data)-1)*sum(apply(data,2,var))  # 当类数是1时，总平方和
  for (i in 2:nc) {       # 对不同的类数
    set.seed(seed)
    wss[i] <- sum(kmeans(data,centers = i)$withinss)  # 得到该类数下的kmeans聚类后的组内平方和
  }
  plot(1:nc,wss,type = "b",xlab = "Number of Clusters",ylab = "Within groups Sum of Squares")
}

wssplot(data)

法二：

library(NbClust)
devAskNewPage(ask = T)   # 同时生成4张图时，是两个两个显示的，一般就显示最后的两张，设置这个询问后，每次都让你按键进入下两个的显示
nc <- NbClust(data,min.nc = 2,max.nc = 15,method = "kmeans")
par(mfrow=c(1,1))
table(nc$Best.nc[1,])
barplot(table(nc$Best.nc[1,]),
        xlab = "推荐的聚类个数",ylab="支持的指标个数",
        main="26个指标推荐的聚类个数")

k-means Cluster

set.seed(1234)
fit <- kmeans(data,centers = 3,nstart = 25)  # 随机初始化25次，选择效果最好的那一次初始化
fit$size        # 查看聚类个数
fit$centers     # 查看类中心（标准化后的）
aggregate(wine[,-1],by=list(fit$cluster),mean)  # 查看原始数据的类心（类均值），原始数据按聚类结果的cluster分组

评价聚类效果

t <- table(wine$Type,fit$cluster)  # 类似于混淆矩阵的交叉列表,显示分错类的个数
library(flexclust)
randIndex(t)   # 兰德指数=0.897
# 等价于如下
comPart(wine$Type,fit$cluster,type = "ARI")  # ARI：调整后的RI，其他类型："J":Jaccard Index,"FM":Fowlkes-Mallows
# 等价于如下，默认调整后的RI
randIndex(wine$Type,fit$cluster)  # 参数correct=T

3. 基于中心点的划分聚类 (k-means的稳健版本)

由于k-means基于均值，所以对异常值敏感，如果围绕中心点（最有代表性的观测值）划分则更稳健，还适用于混合数据类型，不仅限于连续型

library(cluster)
set.seed(1234)
fit_pam <- pam(wine[,-1],k=3,stand=T)  # stand:自动标准化，-均值/标准差
fit_pam$medoids    # 查看中心点，实际的观测值

评价聚类效果

t_pam <- table(wine$Type,fit_pam$clustering)
randIndex(t_pam)        # 0.699

NbCluster中的CCC指标可以验证数据是否可以聚类

nc <- NbClust(data,min.nc = 2,max.nc = 15,method = "kmeans")
plot(nc$All.index[,4])

案例：旅游用户评分聚类分析

一. 选题介绍

近年来，随着旅游市场的不断发展，选择旅游的人数大幅增长。但是在前往未知的旅行地时，游客往往并不明确自己喜欢去什么景点，只能在网上看别人的攻略和评论。而旅游网站上有大量的评价信息和用户信息，网站可以根据这些已有的信息对用户进行细分，判断不同的用户喜欢什么类型的景点或建筑，从而有针对性的向不同的用户推荐不同的符合他们口味的地点，起到个性化推荐的作用。
常用的个性化推荐算法包括协同过滤和矩阵分解，但是均比较简单，适合小型系统。而聚类作为无监督学习的代表，适合在大型系统中将用户区分开。因此本案例意图对旅行中给出评分的用户进行聚类，以期对不同的用户推荐不同的旅行内容。

二. 数据获取与描述

本案例使用的数据来自于UCI机器学习数据库，是谷歌用户对欧洲不同类别的景点的评分数据。用户评分范围在0-5之间，并且每个类别的景点已计算了用户的平均评分。由于原始数据太大，共有5454条样本，所以对其进行无放回随机抽样，得到500条样本。
经过数据预处理后的样本前六行如图1所示，共有18个变量，分别是用户ID，教堂、度假村、沙滩、公园、剧院、博物馆、购物商场、动物园、饭店、酒吧、当地服务、旅店住宿、艺术画廊、舞蹈俱乐部、观光景点、纪念碑和花园。

图 1 评分数据前六行

三.模型建立-K均值聚类

本案例的目的是想将评分相似的用户聚成一类，一类的用户可能有相似的旅行场所的兴趣，以便在后期可以针对用户进行个性化推荐，即对不同类的用户推荐不同的内容。聚类属于无监督学习，是基于距离来判断样本间的相似程度。常用聚类方法有层次聚类和划分聚类，包括K均值聚类和围绕中心点的划分聚类等。由于层次聚类和Ｋ中心点聚类均适用于小数聚集的情况，而K均值聚类由于简洁与效率应用最为广泛，所以本案例采用了K均值聚类方法。
K均值聚类方法的基本思想是：
(1). 首先在所有数据中随机选择K个数据点，作为初始的聚类中心
(2). 然后计算其余各点到各聚类中心的距离，将这些点划分到距离它最近的聚类中心
(3). 重新计算新的类中心，并将数据重新划分到最近的类中
(4). 经过不断迭代，直至损失函数不再减少。
R软件中的K均值聚类中的第k类采用的目标函数是：

由于算法中要进行均值的计算，所以要求变量的取值必须是连续的，而且算法容易受到异常值的影响。本案例的数据经过预处理后基本符合要求，使用K均值聚类方法如下：

1. 确定聚类的个数K

由于算法需要在开始时随机选择K个点作为初始聚类中心，所以要提前确定聚类的个数。图2展示了数据被分为不同的聚类个数时总的类内的平方和，可以看出当聚类数为3时，图中有一个小的拐点，之后曲线的下降变缓，所以适当的类数可以是3。

图 2 聚类个数与组内组内平方和的对比图
R还提供了一种更准确的判断方法，使用NbClust包提供的更多指标来选择类数。如图3所示。可以发现图中显示了共有23个指标给出了自己建议的类数，其中有7个指标支持聚成3类，因此根据“投票”原则，可以选择推荐个数最多的聚类方案。本案例综合两幅图的结果，将聚类个数设为3。

图 3 NbClust包确定聚类个数
由于初始中心点是随机选择的，并且对结果是否收敛影响较大，所以聚类时设定了尝试25次不同的初始配置，最终采用最好的一次作为初始配置。

四.结果分析

聚类后的总平方和为8483，每一类的组内平方和分别为1080、2870、2379。每一类的个数分别为105、196、199。根据前两个主成分画出所有数据点的聚类图，如图4所示，可以发现前两个成分解释了36.56%的可能性。数据被很清晰的分为了三类。

图 4 聚类结果图
将原始数据的聚类中心输出，如表1所示。

可以发现，第一类的用户对观光景点、纪念碑、花园、沙滩这些旅游地点评分较高，均值均超过了2.5分，而对旅店、饭店、博物馆、舞蹈俱乐部评分最低。说明这类人群可能爱好优美的自然风光和室外景点，而相对不喜欢去旅途中室内的场所。对这类人群可以推荐他们室外的风景好的观光景点。
对第二类用户可以发现他们对剧院的平均评分最高，达到了4.14分，对博物馆、购物商场、公园、饭店的评分也很高，超过了3分，而对教堂、画廊、舞蹈俱乐部、纪念碑这些地点评分最低，低于1.6分。说明这类人群可能爱好休闲娱乐较强的可以放松减压的场所，而相对不喜欢去一些气氛肃穆的场所。对于这类人群就可以推荐各种休闲减压的旅行场所。
在第三类用户的评分中，可以看出他们对饭店的评分最高为4.17分，其次还有购物商场、酒吧和当地服务。而其他场所均比较低，舞蹈俱乐部和纪念碑评分最低，在1分以下。说明这类用户最喜欢吃吃吃和买买买，也注重当地服务的体验，很可能是不差钱的外地游客，喜欢当地的新鲜的风土人情。对于他们就可以推荐各种当地特色美食和特产。
综合三类用户来看，公园、剧院和商场这三种场所评分均较高，均超过2分，说明人们都比较喜欢去这三类地方。而舞蹈俱乐部的评分均最低，均在2分以下，说明一般的普通游客在旅游时并不是很喜欢这种场所。因此就可以避免向普通用户推荐这类场所，或者再进一步细化，只向爱好跳舞的游客推荐。

五.总结

本案例使用了K均值聚类根据对旅游场所的评分将用户聚为了三类，在以后可以对不同类的用户进行个性化推荐。但是聚类是一种比较弱的个性化推荐，因为这种方法的本质是识别一类的用户，对该类的所有用户推荐相同的内容。在实际中可以考虑将用户聚类作为第一步，缩小不同用户的范围，然后再针对一类的用户使用其他推荐算法进行更细致精确的推荐。

六. 代码

library(NbClust)
library(ggplot2)
library(cluster)
library(ggfortify)
library(dplyr)
rm(list = ls())

data <- read.csv("E:\\Bing\\研究生\\多元\\旅行评分\\谷歌旅行评分\\google_review_ratings.csv",header=T,stringsAsFactors=F)
data$X <- NULL
data <- na.omit(data)
data$Category.11 <- as.numeric(data$Category.11)
names(data) <- c("id","church","resort","beach","park","theatre",
                 "museum","mall","zoo","restaurant","pub_bar","localService",
                 "burger_pizza_shop","hotel","juiceBar","artGallery","danceClub",
                 "swimmingPool","gym","bakery","beauty_spas","cafe","viewPoint","monument","garden")
set.seed(1234)
da <- data[sample(1:nrow(data),500),]       # 随机抽样，抽一半样本，490
da <- da[order(as.numeric(row.names(da))),]  # 索引排序
summary(da)

## 筛选变量
da <- da[,-c(13,15,18:22)]
head(da)
da <- da[,-1]
df <- scale(da)

## 确定聚类的个数
# 作出类内平方和对聚类数量的曲线
# # nc:考虑的最大聚类数目，seed：随机数种子
wssplot <- function(data,nc=15,seed=1234) 
{
  wss <- (nrow(data)-1)*sum(apply(data,2,var))  # 当类数是1时，总平方和
  for (i in 2:nc) {       # 对不同的类数
    set.seed(seed)
    wss[i] <- kmeans(data,centers = i)$tot.withinss  # 得到该类数下的kmeans聚类后的组内平方和
  }
  plot(1:nc,wss,type = "b",xlab = "Number of Clusters",ylab = "Within groups Sum of Squares")
}
wssplot(df)

## 法二：
set.seed(123)
devAskNewPage(ask = T)   
nc <- NbClust(df,min.nc = 2,max.nc = 15,method = "kmeans")
par(mfrow=c(1,1))
table(nc$Best.nc[1,])
barplot(table(nc$Best.nc[1,]),
        xlab = "推荐的聚类个数",ylab="支持的指标个数",
        main="23个指标推荐的聚类个数",col="lightblue")

## kmeans聚类
set.seed(1234)
fit <- kmeans(df,3,nstart = 25)
fit$centers               # 类中心，标准化后的
options(digits = 3)
group <- aggregate(da,by=list(fit$cluster),mean) %>% t()
autoplot(fit,df,frame=T)    # 聚类结果可视化
fit$cluster               # 显示类别
fit$totss   　　　　　　　# 总平方和
fit$withinss
fit$size