R语言入门--第十四节（聚类分析）

（1）定义每一个观测值为一类；

（2）计算每一类和其它各类的距离；

（3）把“距离”最短的两类合并成一类，这样类的个数就减少一个；

（4）重复步骤1和步骤2，直到包含所有观测值的类合并成单个类为止。

基于5种营养标准含量（变量）的27种食物（观测）进行层次聚类分析，探索不同食物的相同点与不同点，并分成有意义的类。此处层次聚类算法以平均联动(average)为例。

（1）数据预处理--归一化

（2）计算欧几里得距离

（3）平均联动层次聚类分析

（1）确定聚类个数

NbClust包提供了众多的指数来确定在一个聚类分析里类的最佳数目。

（2）获取最终的聚类方案

由上图，尝试解释每类变量的含义：

K均值聚类为最常见的划分方法。

（1）选择K个中心点（随机选择K个观测），K数值就是我们预期的聚类数。

（2）把每个数据点分配给离它最近的中心点；第一次中心点是随机选择的，但也可以设置参数，选择最优的初始值。

（3）重新计算每类中的点到该类中心点距离的平均值；此时的中心点应该为每一类的均值中心点，对异常值敏感（之后都是如此）

（4）分配每个数据到它最近的中心点；

（5）重复步骤3、4，直到所有的观测值不在被分配或是达到最大的迭代次数（默认10次）

（1）数据预处理：去除第一列干扰数据，并归一化数据。

（2）确定待提取的聚类个数，同样可用NbClust包判断（顺序与层次聚类分析不同，如前所述，层次聚类分析在最后才确定聚类个数）

（3）K均值聚类分析

（4）最后将聚类结果与原始数据标准结果（第一列数据）进行比对，看看分析质量如何。

兰德指数接近0.9，看来K均值聚类算法还不错~

K均值法对均值异常敏感，相比来说，PAM为更稳健的方法。

（1）随机选择K个观测（每个都称为中心点）；

（2）计算观测值到各个中心的距离；

（3）把每个观测值分配到最近的中心点；

（4）计算每个中心点到每个观测值的距离的总和（总成本）；

（5）选择一个该类中不是中心的点，并和中心点互换；

（6）重新把每个点分配到距它最近的中心点；

（7）再次计算总成本；

（8）若新的总成本比步骤4计算的总成本少，就把新的点作为中心点；

（9）重复步骤5-8，直到中心点不变。

群落按照物种相似形组成进行聚类分析，可以用树状图较好的表现物种的组成关系。受到很多植被学家的重视。这里以R软件实现聚类分析为例。如果按照物种组成的相似性做聚类分析，那么可以用Jaccard指数（经过转换的）。Jaccard指数只考虑物种在两个样方间是否重复出现，盖度在分析的过程中并不起什么作用。但是如果对乔木和灌木进行分析，就可以考虑个体的数量，计算样方物种组成的相似性的时候用Bray-Curtis指数。Jaccard指数和Bray-Curtis指数在众多生态学相关的程序包中都是可以计算的。下面说一下在R软件中，结合vegan程序包，对草本样方的物种组成进行聚类分析。下面是在R中的具体 *** 作过程：#第一步#是矩阵的整理，建议先整理一下各样地的名录，成如下格式，再用R整理成物种矩阵。

R语言学习笔记之聚类分析

使用k-means聚类所需的包：

factoextra

cluster #加载包

library(factoextra)

library(cluster)l

#数据准备

使用内置的R数据集USArrests

#load the dataset

data("USArrests")

#remove any missing value (i.e, NA values for not available)

#That might be present in the data

USArrests <- na.omit(USArrests)#view the first 6 rows of the data

head(USArrests, n=6)

在此数据集中，列是变量，行是观测值

在聚类之前我们可以先进行一些必要的数据检查即数据描述性统计，如平均值、标准差等

desc_stats <- data.frame( Min=apply(USArrests, 2, min),#minimum

Med=apply(USArrests, 2, median),#median

Mean=apply(USArrests, 2, mean),#mean

SD=apply(USArrests, 2, sd),#Standard deviation

Max=apply(USArrests, 2, max)#maximum

)

desc_stats <- round(desc_stats, 1)#保留小数点后一位head(desc_stats)

变量有很大的方差及均值时需进行标准化

df <- scale(USArrests)

#数据集群性评估

使用get_clust_tendency()计算Hopkins统计量

res <- get_clust_tendency(df, 40, graph = TRUE)

res$hopkins_stat

## [1] 0.3440875

#Visualize the dissimilarity matrix

res$plot

Hopkins统计量的值<0.5，表明数据是高度可聚合的。另外，从图中也可以看出数据可聚合。

#估计聚合簇数

由于k均值聚类需要指定要生成的聚类数量，因此我们将使用函数clusGap()来计算用于估计最优聚类数。函数fviz_gap_stat()用于可视化。

set.seed(123)

## Compute the gap statistic

gap_stat <- clusGap(df, FUN = kmeans, nstart = 25, K.max = 10, B = 500)

# Plot the result

fviz_gap_stat(gap_stat)

图中显示最佳为聚成四类（k=4）

#进行聚类

set.seed(123)

km.res <- kmeans(df, 4, nstart = 25)

head(km.res$cluster, 20)

# Visualize clusters using factoextra

fviz_cluster(km.res, USArrests)

#检查cluster silhouette图

Recall that the silhouette measures (SiSi) how similar an object ii is to the the other objects in its own cluster versus those in the neighbor cluster. SiSi values range from 1 to - 1:

A value of SiSi close to 1 indicates that the object is well clustered. In the other words, the object ii is similar to the other objects in its group.

A value of SiSi close to -1 indicates that the object is poorly clustered, and that assignment to some other cluster would probably improve the overall results.

sil <- silhouette(km.res$cluster, dist(df))

rownames(sil) <- rownames(USArrests)

head(sil[, 1:3])

#Visualize

fviz_silhouette(sil)

图中可以看出有负值，可以通过函数silhouette()确定是哪个观测值

neg_sil_index <- which(sil[, "sil_width"] <0)

sil[neg_sil_index, , drop = FALSE]

##          cluster    neighbor     sil_width

## Missouri    3          2        -0.07318144

#eclust():增强的聚类分析

与其他聚类分析包相比，eclust()有以下优点：

简化了聚类分析的工作流程

可以用于计算层次聚类和分区聚类

eclust()自动计算最佳聚类簇数。

自动提供Silhouette plot

可以结合ggplot2绘制优美的图形

#使用eclust()的K均值聚类

# Compute k-means

res.km <- eclust(df, "kmeans")

# Gap statistic plot

fviz_gap_stat(res.km$gap_stat)

# Silhouette plotfviz_silhouette(res.km)

##    cluster size ave.sil.width

## 1     1     13      0.31

## 2     2     29      0.38

## 3     3      8      0.39

#使用eclust（）的层次聚类

# Enhanced hierarchical clustering

res.hc <- eclust(df, "hclust") # compute hclust

fviz_dend(res.hc, rect = TRUE) # dendrogam

#下面的R代码生成Silhouette plot和分层聚类散点图。

fviz_silhouette(res.hc) # silhouette plot

##   cluster size ave.sil.width

## 1    1     19      0.26

## 2    2     19      0.28

## 3    3     12      0.43

fviz_cluster(res.hc) # scatter plot

#Infos

This analysis has been performed using R software (R version 3.3.2)

---