《精通机器学习：基于R 第二版》学习笔记
1、无监督学习
前面的案例中都有一个响应变量Y，Y是X的函数，或表示为y = f(x)。我们的数据中有Y的实际值，以Y为依据训练X。这种方式称为监督式学习。但是，从数据中学习知识的时候通常没有响应变量Y，或者故意忽略了Y，这便是无监督学习。无监督学习建立和选择算法的基准是算法能够明确满足业务需求的程度，而不是算法的正确率。
为什么要研究无监督学习？首先，无监督学习可以帮助你理解并识别出数据中的模式，这可能非常有价值。其次，你可以通过无监督学习转换数据，提高监督式学习技术的效果。
2、数据理解与数据准备
> library(pacman)
> p_load(dplyr, HDclassif)
> data(wine)
> str(wine)
## 'data.frame': 178 obs. of 14 variables:
## $ class: int 1 1 1 1 1 1 1 1 1 1 ...
## $ V1 : num 14.2 13.2 13.2 14.4 13.2 ...
## $ V2 : num 1.71 1.78 2.36 1.95 2.59 1.76 1.87 2.15 1.64 1.35 ...
## $ V3 : num 2.43 2.14 2.67 2.5 2.87 2.45 2.45 2.61 2.17 2.27 ...
## $ V4 : num 15.6 11.2 18.6 16.8 21 15.2 14.6 17.6 14 16 ...
## $ V5 : int 127 100 101 113 118 112 96 121 97 98 ...
## $ V6 : num 2.8 2.65 2.8 3.85 2.8 3.27 2.5 2.6 2.8 2.98 ...
## $ V7 : num 3.06 2.76 3.24 3.49 2.69 3.39 2.52 2.51 2.98 3.15 ...
## $ V8 : num 0.28 0.26 0.3 0.24 0.39 0.34 0.3 0.31 0.29 0.22 ...
## $ V9 : num 2.29 1.28 2.81 2.18 1.82 1.97 1.98 1.25 1.98 1.85 ...
## $ V10 : num 5.64 4.38 5.68 7.8 4.32 6.75 5.25 5.05 5.2 7.22 ...
## $ V11 : num 1.04 1.05 1.03 0.86 1.04 1.05 1.02 1.06 1.08 1.01 ...
## $ V12 : num 3.92 3.4 3.17 3.45 2.93 2.85 3.58 3.58 2.85 3.55 ...
## $ V13 : int 1065 1050 1185 1480 735 1450 1290 1295 1045 1045 ...
数据包括178种葡萄酒，有13个变量表示酒中的化学成分，还有一个标号变量Class，表示品种等级或葡萄种植品种。
V1 ：酒精
V2 ：苹果酸
V3 ：灰分
V4 ：灰分碱性
V5 ：镁
V6 ：总酚含量
V7 ：黄酮类化合物
V8 ：非黄酮类酚类
V9 ：原花青素
V10 ：色彩强度
V11 ：色调
V12 ：OD280/OD315
V13 ：脯氨酸
给变量添加有意义的名称：
> names(wine) <- c("class", "alcohol", "malic_acid", "ash", "alk_ash", "magnesium",
+ "t_phenols", "flavanoids", "non_flav", "proantho", "c_intensity", "hue", "od280_315",
+ "proline")
> names(wine)
## [1] "class" "alcohol" "malic_acid" "ash" "alk_ash"
## [6] "magnesium" "t_phenols" "flavanoids" "non_flav" "proantho"
## [11] "c_intensity" "hue" "od280_315" "proline"
数据标准化：
> df <- scale(wine[, -1]) %>% as_tibble()
> head(df)
## # A tibble: 6 x 13
## alcohol malic_acid ash alk_ash magnesium t_phenols flavanoids non_flav
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1.51 -0.561 0.231 -1.17 1.91 0.807 1.03 -0.658
## 2 0.246 -0.498 -0.826 -2.48 0.0181 0.567 0.732 -0.818
## 3 0.196 0.0212 1.11 -0.268 0.0881 0.807 1.21 -0.497
## 4 1.69 -0.346 0.487 -0.807 0.928 2.48 1.46 -0.979
## 5 0.295 0.227 1.84 0.451 1.28 0.807 0.661 0.226
## 6 1.48 -0.516 0.304 -1.29 0.858 1.56 1.36 -0.176
## # ... with 5 more variables: proantho <dbl>, c_intensity <dbl>, hue <dbl>,
## # od280_315 <dbl>, proline <dbl>
查看品种等级class分布情况：
> table(wine$class)
##
## 1 2 3
## 59 71 48
3、层次聚类
层次聚类算法的基础是观测之间的相异度测量。我们使用的是通用的测量方式——欧氏距离，当然还有其他方式。
Ward距离 （Ward's-Linkage）使总的簇内方差最小，使用簇中的点到质心的误差平方和作为测量方式
最大距离（Complete linkage） 两个簇之间的距离就是两个簇中的观测之间的最大距离
最小距离（Single linkage） 两个簇之间的距离就是两个簇中的观测之间的最小距离
平均距离（Average linkage） 两个簇之间的距离就是两个簇中的观测之间的平均距离
质心距离（Centroid linkage） 两个簇之间的距离就是两个簇的质心之间的距离
最后需要注意的是，要对你的数据进行标准化的缩放操作，使数据的均值为0，标准差为1，这样在计算距离时才能进行比较。否则，变量的测量值越大，对距离的影响就越大。
> p_load(stats,NbClust)
>
> num.complete <- NbClust(df,
+ # 测距方式为欧氏距离
+ distance = "euclidean",
+ # 簇的最小数量
+ min.nc = 2,
+ # 簇的最大数量
+ max.nc = 6,
+ # 使用所有有效性指标
+ method = "complete",index = "all")
Hubert index统计图
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
##
Hubert 指数图，在左侧的图中，需要找出一个明显的拐点；在右侧的图中，要找到峰值。
Dindex统计图
Dindex 指数图也提供了同样的信息。
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 1 proposed 2 as the best number of clusters
## * 11 proposed 3 as the best number of clusters
## * 6 proposed 5 as the best number of clusters
## * 5 proposed 6 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************
根据少数服从多数原则，我们应该选择3个簇作为最优解。另外一种确定簇数目的著名方法是gap统计量。
查看每种有效性指标的最优簇数量和与之对应的指标值：
> num.complete$Best.nc
## KL CH Hartigan CCC Scott Marriot TrCovW
## Number_clusters 5.0000 3.0000 3.0000 5.000 3.0000 3.000000e+00 3.00
## Value_Index 14.2227 48.9898 27.8971 1.148 340.9634 6.872632e+25 22389.83
## TraceW Friedman Rubin Cindex DB Silhouette Duda
## Number_clusters 3.0000 3.0000 5.0000 3.0000 6.0000 3.0000 5.0000
## Value_Index 256.4861 10.6941 -0.1489 0.3551 1.6018 0.2038 0.8856
## PseudoT2 Beale Ratkowsky Ball PtBiserial Frey McClain
## Number_clusters 5.0000 5.0000 3.0000 3.0000 6.0000 1 2.0000
## Value_Index 6.3314 1.1253 0.3318 462.0304 0.5877 NA 0.7687
## Dunn Hubert SDindex Dindex SDbw
## Number_clusters 6.0000 0 6.0000 0 6.0000
## Value_Index 0.1892 0 0.9951 0 0.5031
聚类有效性指标：表现最好的前5种指标是CH指数、Duda指数、Cindex、Gamma和Beale指数。在R中，我们可以使用NbClust包中的NbClust()函数，求出23种聚类有效性指标的结果。
可以看到第一个指标（KL）的最优簇数量是5，第二个指标（CH）的最优簇数量是3。选择3个簇进行聚类：
> # 计算距离矩阵
> dis <- dist(df, method = "euclidean")
>
> # 建立层次聚类模型 hclust() 函数需要两个基本输入：距离矩阵和聚类方法
> hc <- hclust(dis, method = "complete")
画出聚类结果树状图：
> # 参数 hang=-1 表示将观测排列在图的底部
> plot(hc, hang = -1, labels = F, main = "Complete-Linkage")
最大距离聚类结果树状图
对树状图剪枝后，使用sparcl包生成彩色树状图：
> p_load(sparcl)
> com3 <- cutree(hc, 3)
> # 参数 branchlength=50要根据实际数据确定
> ColorDendrogram(hc, y = com3, main = "Complete", branchlength = 50)
彩色树状图
比较聚类结果和品种等级的标号：
> table(com3, wine$class)
##
## com3 1 2 3
## 1 51 18 0
## 2 8 50 0
## 3 0 3 48
下面试试Ward距离法，hclust() 函数使用的默认方法是最大距离法。在最大距离法得出的结果中，两个簇之间的距离等于两个簇中的观测之间的最大距离。Ward距离法对观测进行聚集的目标是，使簇内观测的误差平方和最小。在R中，ward.D2方法使用欧氏距离的平方来测量距离，这是真正的Ward距离法。
> num.ward <- NbClust(df, diss = NULL, distance = "euclidean", min.nc = 2, max.nc = 6,
+ method = "ward.D2", index = "all")
Hubert 统计图
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
##
Dindex统计图
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 2 proposed 2 as the best number of clusters
## * 18 proposed 3 as the best number of clusters
## * 2 proposed 6 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************
跟前面的结果一样，簇的数目还是3。
聚类并生成树状图：
> hc.ward <- hclust(dis, method = "ward.D2")
> ward3 <- cutree(hc.ward, 3)
> ColorDendrogram(hc.ward, y = ward3, main = "Ward's-Linkage", branchlength = 50)
ward距离法聚类树状图
与品种等级标号进行比较：
> table(ward3, wine$class)
##
## ward3 1 2 3
## 1 59 5 0
## 2 0 58 0
## 3 0 8 48
比较两种距离计算方法的观测匹配情况：
> table(com3, ward3)
## ward3
## com3 1 2 3
## 1 53 11 5
## 2 11 47 0
## 3 0 0 51
查看簇中的观测，计算出数据中13个变量在各个簇中的均值：
> # 计算均值统计量
> aggregate(wine[, -1], list(com3), mean)
## Group.1 alcohol malic_acid ash alk_ash magnesium t_phenols flavanoids
## 1 1 13.40609 1.898986 2.305797 16.77246 105.00000 2.643913 2.6689855
## 2 2 12.41517 1.989828 2.381379 21.11724 93.84483 2.424828 2.3398276
## 3 3 13.11784 3.322157 2.431765 21.33333 99.33333 1.675686 0.8105882
## non_flav proantho c_intensity hue od280_315 proline
## 1 0.2966667 1.832899 4.990725 1.0696522 2.970000 984.6957
## 2 0.3668966 1.678103 3.280345 1.0579310 2.978448 573.3793
## 3 0.4443137 1.164314 7.170980 0.6913725 1.709804 622.4902
> aggregate(wine[, -1], list(ward3), mean)
## Group.1 alcohol malic_acid ash alk_ash magnesium t_phenols flavanoids
## 1 1 13.66922 1.970000 2.463125 17.52812 106.15625 2.850000 3.0096875
## 2 2 12.20397 1.938966 2.215172 20.20862 92.55172 2.262931 2.0881034
## 3 3 13.06161 3.166607 2.412857 21.00357 99.85714 1.694286 0.8478571
## non_flav proantho c_intensity hue od280_315 proline
## 1 0.2910937 1.908125 5.450000 1.071406 3.158437 1076.0469
## 2 0.3553448 1.686552 2.895345 1.060000 2.862241 501.4310
## 3 0.4494643 1.129286 6.850179 0.721000 1.727321 624.9464
各个数值都非常接近，使用统计图有助于理解和解释：
> par(mfrow = c(1, 2))
> boxplot(wine$proantho ~ com3, data = wine, main = "Proline by Complete Linkage")
> boxplot(wine$proantho ~ ward3, data = wine, main = "Proline by Ward's Linkage")
两种聚类方法对比
最大距离法的疑似离群点更少一些。查看每个变量的箱线图有助于和领域专家确定最好的层次聚类方法。
4、K均值聚类
使用K均值聚类时，需要明确指定所需的簇的数目，然后算法开始迭代，直到每个观测都属于某个簇。算法的目标是使簇内的差异最小，簇内差异由欧氏距离的平方定义。所以，第k个簇的簇内差异等于簇内所有两个观测之间的欧氏距离的平方和，再除以簇内观测的数量。
因为第1步中的初始分配是随机的，所以会造成每次聚类结果不一致。因此，重要的一点是，要进行多次初始分配，让软件找出最优的解。
> num.kmeans <- NbClust(df, min.nc = 2, max.nc = 15, method = "kmeans")
Hubert统计图
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
##
Dindex统计图
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 2 proposed 2 as the best number of clusters
## * 19 proposed 3 as the best number of clusters
## * 1 proposed 14 as the best number of clusters
## * 1 proposed 15 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************
从结果可以看出，3个簇再次成为最优解。
> set.seed(123)
> # 随机分配的参数nstart表示初始分配
> km <- kmeans(df, centers = 3, nstart = 25)
查看观测在簇之间的分布情况：
> table(km$cluster)
##
## 1 2 3
## 51 62 65
查看簇中心点矩阵，它保存了每个簇中每个变量的中心点值：
> km$centers
## alcohol malic_acid ash alk_ash magnesium t_phenols
## 1 0.1644436 0.8690954 0.1863726 0.5228924 -0.07526047 -0.97657548
## 2 0.8328826 -0.3029551 0.3636801 -0.6084749 0.57596208 0.88274724
## 3 -0.9234669 -0.3929331 -0.4931257 0.1701220 -0.49032869 -0.07576891
## flavanoids non_flav proantho c_intensity hue od280_315
## 1 -1.21182921 0.72402116 -0.77751312 0.9388902 -1.1615122 -1.2887761
## 2 0.97506900 -0.56050853 0.57865427 0.1705823 0.4726504 0.7770551
## 3 0.02075402 -0.03343924 0.05810161 -0.8993770 0.4605046 0.2700025
## proline
## 1 -0.4059428
## 2 1.1220202
## 3 -0.7517257
第二个簇具有相对较高的酒精成分。我们做一个箱线图看看酒精成分的分布，用Ward距离法进行对比：
> par(mfrow = c(1, 2))
> boxplot(wine$alcohol ~ km$cluster, data = wine, main = "Alcohol Content,K-Means")
> boxplot(wine$alcohol ~ ward3, data = wine, main = "Alcohol Content,Ward's")
Alcohol含量对比图
对比K均值聚类结果和品种等级：
> table(km$cluster, wine$class)
##
## 1 2 3
## 1 0 3 48
## 2 59 3 0
## 3 0 65 0
与层次聚类法的结果区别还是很大的。
5、果瓦系数与围绕中心的划分
无论层次聚类还是K均值聚类，都不是为分析混合数据而专门设计的。混合数据是既包括定量数据又包括定性数据的数据，说得更具体一些，就是包括名义数据、定序数据、定距数据和定比数据。
这时可以使用果瓦相异度系数将混合数据转换为适当的特征空间。在这种方法中，你甚至可以使用因子作为输入变量。在聚类算法方面，推荐使用PAM聚类算法。
> p_load(cluster)
>
> # 将酒精变量转换为因子类型
> wine$alcohol <- ifelse(df$alcohol > 0, "High", "Low") %>% as.factor()
>
> # 建立相异度矩阵
> dis.matrix <- daisy(wine[, -1], metric = "gower")
>
> # 建立聚类对象
> set.seed(123)
> pam.fit <- pam(dis.matrix, k = 3)
>
> # 与品种等级对比
> table(pam.fit$clustering, wine$class)
##
## 1 2 3
## 1 57 5 0
## 2 2 64 5
## 3 0 2 43
使用compareGroups包建立用于商业展示的描述性统计表格：
> p_load(compareGroups)
>
> wine$cluster <- pam.fit$clustering
>
> # 建立compareGroups对象，保存聚类结果中的描述性统计
> group <- compareGroups(cluster ~ ., data = wine)
>
> # 转换为表格
> clustab <- createTable(group)
> clustab
商业统计表格
这张表显示了因子水平在各个簇中的比例，对于数值型变量，显示了均值以及括号中的标准差。要将表格导出为.csv文件，使用 export2csv() 函数即可：
> export2csv(clustab, file = "wine_clusters.csv")
6、随机森林
使用随机森林方法处理混合数据有以下优点：
对异常值和高度扭曲的变量的鲁棒性更好；
不需要对数据进行转换和比例缩放；
可以处理混合数据（数值型和因子）；
可以兼容缺失值；
可以在有大量变量的数据上使用，实际上，通过检查变量重要性，可以消除那些无用的特征；
生成的相异度矩阵可以作为我们前面讨论过的一些技术的输入（层次聚类、K均值聚类和PAM）。
> p_load(randomForest)
> set.seed(123)
> rf <- randomForest(x = wine[, -1], ntree = 2000, proximity = T)
> rf
##
## Call:
## randomForest(x = wine[, -1], ntree = 2000, proximity = T)
## Type of random forest: unsupervised
## Number of trees: 2000
## No. of variables tried at each split: 3
查看矩阵维数：
> dim(rf$proximity)
## [1] 178 178
检查矩阵的前5行和前5列：
> rf$proximity[1:5, 1:5]
## 1 2 3 4 5
## 1 1.0000000 0.3527273 0.3840580 0.3442029 0.1881533
## 2 0.3527273 1.0000000 0.1854545 0.1724138 0.0701107
## 3 0.3840580 0.1854545 1.0000000 0.2934783 0.2326389
## 4 0.3442029 0.1724138 0.2934783 1.0000000 0.1048689
## 5 0.1881533 0.0701107 0.2326389 0.1048689 1.0000000
对这些值的一种理解方式是，将它们看作两个观测出现在同一个终端节点的次数百分比。
查看变量重要性：
> importance(rf)
## MeanDecreaseGini
## alcohol 3.869152
## malic_acid 11.482949
## ash 9.655096
## alk_ash 10.598505
## magnesium 9.744547
## t_phenols 16.810291
## flavanoids 20.779212
## non_flav 10.866666
## proantho 13.239989
## c_intensity 13.889777
## hue 13.795760
## od280_315 16.919988
## proline 14.734792
## cluster 11.135772
创建相异度矩阵（聚类算法的输入特征）：
> diss.mat <- sqrt(1 - rf$proximity)
>
> diss.mat[1:2, 1:2]
## 1 2
## 1 0.0000000 0.8045326
## 2 0.8045326 0.0000000
PAM聚类：
> set.seed(123)
> pam.rf <- pam(diss.mat, k = 3)
>
> # 与分类特征值对比
> table(pam.rf$clustering, wine$class)
##
## 1 2 3
## 1 55 0 0
## 2 4 70 5
## 3 0 1 43
如果聚类数据是一团乱麻，那么可以考虑随机森林方法。上面的结果和使用其他技术得到的结果很相似。还能通过对随机森林参数进行调优来改善。
