基于K均值聚类的葡萄酒品种判别
特别注意：主要思路、程序和分析过程来源于：https://www.kaggle.com/xvivancos/tutorial-clustering-wines-with-k-means。本文在此基础上做了部分改动。
一. K均值聚类原理简介
1.1 K均值聚类的主要步骤如下[1]：
1）选择n个数值型变量参与聚类分析，并选择聚类数为k
2）由系统选择k个样品作为聚类的种子
3）按照到这些类重心的距离最小的原则把所有的样品分类到各类重心所在的类中，完成初始分类
4）重新计算各类的重心（即现分的各子类中所包含的样品的均值）
5）重复第三和第四步，直到达到迭代要求为止
1.2 迭代要求可以是：
1） 没有（或者小于某个数值的）对象被重新分配给不同的类
2） 没有（或者小于某个数值的）类中心发生变化
3） 误差平方和（SSE）达到局部最小
4） 达到指定的迭代次数
此外，K均值算法必须在平均值有意义的情况下才能使用，因此不适用于分类变量。需要给定聚类数目，并且对异常数据和数据噪声比较敏感。
二. 数据简介及预处理
2.1 数据简介
本文所选取的数据wines 下载于https://archive.ics.uci.edu/ml/datasets/wine（简称 UCI），其提供者是Stefan Aeberhard。wines数据集提供了178条关于葡萄酒信息，这些信息是对意大利同一地区三种不同品种葡萄酒的化学分析结果，原数据有30个变量，由于某些原因，UCI上只提供了14列数据。变量的具体信息如表2.1所示。
2.2 数据可视化
首先，简单了解不同葡萄酒种类的比例，如图2.1所示，不同种类分别占比33.15%、26.97%、39.89%，可见种类分布的差异比较均衡。对于较均衡数据的分类效果，可以用误判率指标来比较。
由于聚类分析是一种无监督学习，因此在分析前现将class 这一列去掉。对于处理后的数据，我们分别探索各变量的分布形状（图2.2）及变量间的关系（图2.3）。
由图2.2，我们大致能够对各变量的分布有一个直观的印象。
由图2.3，我们看到，除了变量Total_Phenols 和变量 Flavanoids 有很强的相关性之外，其他变量间的相关性不是很强。
2.3 数据预处理
由于异常值和特殊的变量对聚类有较大影响，因此当分类变量的测量尺度不一致时，需要事先做标准化处理。这里使用scale()函数对数据集进行标准化处理。现以变量Alcolhol 和 变量Malic_acid 为例来可视化标准化的结果，如图2.4所示。可以看出，标准化后的横纵坐标不一样，但两变量的关系没有改变。
三. K均值聚类应用
3.1 K的选择
类似于方差分析，一个数据集分为K个聚类后，存在组间平方和、组内平方和。一个理想的K应该在使得组间平方和尽可能大的同时，组内平方和尽可能少。
为找出理想的K，当K取1:10时，分别求出其组间和组内平方和，如图3.1所示。
综合考虑组间和组内平方和，以及K太多失去聚类意义这三个因素，选择两幅图中变化幅度较大的点（K=3时）。
3.2 利用K均值聚类将数据分成三类
使用Kmeans()可达到聚类的目的，这里选取数据集的前六列可视化聚类后各变量的关系，如图3.2所示。可以看到，数据比较好的被分为了三类。不同种类的分布也各有特点，具体将在3.3小节部分讨论。
为更好的量化聚类结果，可计算聚类后与聚类前的混淆矩阵如下。由混淆矩阵可知，聚类错误率仅为3.37%，说明聚类效果还不错。
3.3 聚类结果的特点
由图3.2 可知，不同种类葡萄酒的Alcohol 均值不一样，按颜色分类的话，均值大小依次是：红、蓝、绿。对变量 Total_Phenols 而言，均值差异也比较明显，其大小依次是：红、绿、蓝。此外，对变量 Malic_add 来说，红色种类葡萄酒的分布非常集中。
四．基于单变量的 K均值聚类
为了分析不同类别间各含量的差异，考虑绘制各类别的不同变量间的差异。由于13个变量的比较较复杂，这里尝试用基于原始数据的主成分方法筛选出不同种类间差异较大的因素，得到各主成分的方差解释度及载荷矩阵如图3.3。由图3.3 可知，wines 各变量由“第一主成分”解释的比例高达99.81%，且“第一个主成分”完全可以由变量 Proline 代表，因此，可选择变量 Proline进行不同类别间的葡萄酒差异比较。
通过计算得到变量 Proline 在不同类别间的均值差异如表3.1 所示。由表3.1可知，种类1的 Proline 值较大，种类3次之，种类2最小，且不同种类之间的Proline 均值含量相差较大。陈勇（2004）等[2]指出，脯氨酸（Proline）是所有葡萄酒中氨基酸的主体，也是各样品含量差异最大的氨基酸，这与葡萄酒品种、气候、发酵工艺及酒的陈酿时间有关。结合我们的分析结果，可尝试以脯氨酸作为特征变量对葡萄酒进行分类。
由于在上一部分已经选出合适的聚类簇K，因此这里直接基于 Proline 将数据分为三类。得到混淆矩阵如下。由混淆矩阵可知，对类别1，2, 3的误判率分别为22.03%，71.83%，60.41%。
基于变量 Proline 进行K均值聚类的方法只对类别1有较好的判别效果，原因在于变量Proline 在类别1和非类别1间差别较大，但在类别2和类别3中不太明显。当数据信息不完整（包括变量 Proline 的信息，但缺少其他信息），可用基于Proline 的K均值聚类对类别1与非类别1做初步判别。
四. 总结
第一，本文通过使用K均值聚类判别葡萄酒种类，发现脯氨酸含量在类别1与非类别中差异较大。
第二，K均值聚类在 wines 数据集上都非常好的甄别出葡萄酒的种类，这和 wines 数据集本身的结构特点关系（种类分布较均衡，变量间相关性不是特别高）。对于不均衡分类的数据K均值聚类的算法如何有待实践。
五. 参考文献
[1] 吴密霞, 刘春玲.多元统计分析(M).北京：科学出版社, 2014:181-182.
[2] 陈勇,曾新安,董新平,肖利民.中国主产干红葡萄酒中氨基酸含量对照与探讨[J].食品与发酵工业,2004(01):107-110.
[3] 薛毅，陈立萍.统计建模与R软件(M).北京：清华出版社，2007.
[4] https://www.kaggle.com/xvivancos/tutorial-clustering-wines-with-k-means.（主要程序和分析部分参考这篇文章）
代码：
# 加载包
library(tidyverse)
library(corrplot)
library(gridExtra)
library(GGally)
library(caret)
library(ggmap)
library(dplyr)
library(ggplot2)
library(knitr)
#读取数据
wines <- read.csv("Wine.csv")
head(wines)
#数据可视化
#绘制葡萄酒种类的饼图
a1 <- nrow(wines[wines$Class == 1,])
a2 <- nrow(wines[wines$Class == 2,])
a3 <- nrow(wines[wines$Class == 3,])
value <-c(a1, a2, a3)
per.value <- paste(round(100 * value / sum(value),2),"%")
Group = c("种类1", "种类2", "种类3")
slice.col <- rainbow(10)
pie(value, labels = per.value, col= slice.col, main = "葡萄酒种类饼图")
legend("topright",Group, cex=0.85, fill = slice.col)
#为进行聚类分析，去掉类别列
#去掉类别类
wine <- wines[,-1]
#对每个属性画直方图
# Histogram for each Attribute
wine %>%
gather(Attributes, value, 1:13) %>%
ggplot(aes(x=value, fill=Attributes)) +
geom_histogram(colour="black", show.legend=FALSE) +
facet_wrap(~Attributes, scales="free_x") +
labs(x="Values", y="Frequency",
title="Wine Attributes - Histograms") +
theme_bw()
#各属性的相关系数矩阵
corrplot(cor(wine), type="upper", method="ellipse", tl.cex=0.9)
cor(wine)
#可以看出属性 Total_Phenols 和属性 Flavanoids 有很强的相关性
#数据预处理
# 标准化
wineNorm <- as.data.frame(scale(wine))
# 原始数据
p1 <- ggplot(wine, aes(x=Alcohol, y=Malic_acid)) +
geom_point() +
labs(title="Original data") +
theme_bw()
# 标准化数据
p2 <- ggplot(wineNorm, aes(x=Alcohol, y=Malic_acid)) +
geom_point() +
labs(title="Normalized data") +
theme_bw()
# Subplot
grid.arrange(p1, p2, ncol=2)
#K-means
#确定k
bss <- numeric()
wss <- numeric()
set.seed(1234)
for(i in 1:10){
# 对每个k，计算组内平方和 和 组间平方和
bss[i] <- kmeans(wineNorm, centers=i)$betweenss
wss[i] <- kmeans(wineNorm, centers=i)$tot.withinss
}
# 根据组间平方和的大小选择K
p3 <- qplot(1:10, bss, geom=c("point", "line"),
xlab="Number of clusters", ylab="Between-cluster sum of squares") +
scale_x_continuous(breaks=seq(0, 10, 1)) +
theme_bw()
# 根据组内平方和的大小选择k
p4 <- qplot(1:10, wss, geom=c("point", "line"),
xlab="Number of clusters", ylab="Total within-cluster sum of squares") +
scale_x_continuous(breaks=seq(0, 10, 1)) +
theme_bw()
grid.arrange(p3, p4, ncol=2)
#选择k=3
set.seed(1234)
wine_k3 <- kmeans(wineNorm, centers=3)
aggregate(wines, by=list(wine_k3$cluster), mean)
ggpairs(cbind(wine, Cluster=as.factor(wine_k3$cluster)),
columns=1:6, aes(colour=Cluster, alpha=0.5),
lower=list(continuous="points"),
upper=list(continuous="blank"),
axisLabels="none", switch="both") +
theme_bw()
wine_k3$cluster
confusionMatrix(as.factor(wines$Class), as.factor(wine_k3$cluster))
#利用主成分分析降维
wine_pr <- princomp(wine)
summary(wine_pr, loading = T)
screeplot(wine_pr, type = "lines")
wine$class <- wine_k3$cluster
wine_data <- wine[,c("Proline", "class")]
#分组统计各类别的均值
mean_pr <- tapply(wine_data[,1], INDEX = wine_data$class, FUN = mean)
#利用变量 proline 进行K-means聚类
#选择数据
wine_p <- wine[,"Proline"]
set.seed(1234)
wine_3 <- kmeans(wine_p, centers=3)
wine_3$cluster
library(e1071)
wines$Class
confusionMatrix(as.factor(wines$Class), as.factor(wine_3$cluster))