应用统计学与R语言实现学习笔记（十）——聚类分析

Chapter 10 Cluster Analysis

本篇是第十章，内容是聚类分析。因为以后的几章是典型的分析方法。并且在14章的案例里面可能不会体现，因此内容里会渗透较多的R语言操做。因为开源中国不支持Latex公式，这篇彻头彻尾又是各类数学公式，为防整篇博客都变成截图，我就把内容进行了删减。 这篇博客的完整内容包含各种数学表达。能够见我CSDN和hexo搭的我的博客。

CSDN博客

hexo的我的博客

1 多元分布基本概念

在研究实际问题的时候，咱们常常遇到的是多变量的问题，因为指标间相互不独立，单独割裂开来分别研究分析，不能从总体上把握所研究问题的实质。因此咱们必须对多元变量及其分布进行统计和分析，在地学领域这种问题比比皆是，这里就不展开阐述了，接下来是一堆纯数学概念，数学恐惧者慎入，这部分的重点应该是关于协方差矩阵。通常来讲，假设所研究的问题有p个指标，进行了n次独立观测，获得了np个数据。 分布函数的性质

• 非降的右连续函数
• 分布函数的取值范围为[0，1]
• 分布函数当变量取值为无穷大时，函数值收敛到1

多元几率密度函数 协方差矩阵 相关系数矩阵

2 数据的变换处理

数据变换是将原始数据矩阵中的每一个元素按照某种特定的运算把它变成为一个新值，并且数值的变化不依赖于原始数据集合中其它数据的新值。事实上多元数据的变换处理一般是为了消除不一样量纲的差别。 较经常使用的数据变换以下：

• 中心化变换 中心化变换是一种坐标轴平移处理方法，它是先求出每一个变量的样本平均值，再从原始数据中减去该变量的均值，就获得中心化变换后的数据。 中心化变换的结果是使每列数据之和均为0，即每一个变量的均值为0，并且每列数据的平方和是该列变量样本方差的(n-1)倍，任何不一样两列数据之交叉乘积是这两列变量样本协方差的(n-1)倍，因此这是一种很方便地计算方差与协方差的变换。
• 极差规格化变换 极差规格化变换是从数据矩阵的每个变量中找出其最大值和最小值，这二者之差称为极差，而后从每一个变量的每一个原始数据中减去该变量中的最小值，再除以极差，就获得规格化数据。通过极差规格化变换后，数据矩阵中每列即每一个变量的最大数值为1，最小数值为0，其他数据取值均在0和1之间；而且变换后的数据都再也不具备量纲，便于不一样的变量之间的比较。
• 标准化变换 标准化变换首先对每一个变量进行中心化变换，而后用该变量的标准差进行标准化。通过标准化变换处理后，数据矩阵中每列数据即每一个变量的平均值为0，方差为1，且再也不具备量纲，便于不一样变量之间的比较。变换后，数据矩阵中任何两列数据乘积之和是所对应的两个变量相关系数的（ n-1）倍，因此这是一种很方便地计算相关矩阵的变换。
• 对数变换 对数变换是将各个原始数据取对数，将原始数据的对数值做为变换后的新值。

3 聚类分析

聚类分析是一种分类技术。与多元分析的其余方法相比，该方法较为粗糙，理论上还不完善，但应用方面取得了很大成功。与回归分析、判别分析一块儿被称为多元分析的三大方法。 聚类的目的——根据已知数据（ 一批观察个体的许多观测指标） ， 按照必定的数学公式计算各观察个体或变量（指标）之间亲疏关系的统计量（距离或相关系数等）。 根据某种准则（ 最短距离法、最长距离法、中间距离法、重心法等），使同一类内的差异较小，而类与类之间的差异较大，最终将观察个体或变量分为若干类。 聚类的种类—— 根据分类的方法可将聚类分析分为：系统聚类、快速聚类、有序聚类。 根据分类的对象可将聚类分析分为：Q型——样品聚类clustering for individuals；R型——指标聚类clustering for variables。 数据结构

4 样品间亲疏程度的测度

样品间亲疏程度的测度 研究样品或变量的亲疏程度的数量指标有两种，一种叫类似系数，性质越接近的变量或样品，它们的类似系数越接近于1，而彼此无关的变量或样品它们的类似系数则越接近于0，类似的为一类，不类似的为不一样类；另外一种叫距离，它是将每个样品看做p维空间的一个点，并用某种度量测量点与点之间的距离，距离较近的归为一类，距离较远的点属于不一样的类。 变量之间的聚类即R型聚类分析，经常使用类似系数来测度变量之间的亲疏程度。而样品之间的聚类即Q型聚类分析，则经常使用距离来测度样品之间的亲疏程度。 距离 假使每一个样品有p个变量，则每一个样品均可以当作p维空间中的一个点， n个样品就是p维空间中的n个点，则第i样品与第j样品之间的距离能够进行计算。 定义距离的准则 定义距离要求知足第i个和第j个样品之间的距离的四个条件 经常使用距离

• 明氏距离(Minkowski) 明氏距离、欧氏距离、绝对距离、Chebychev距离。 明氏距离主要有如下两个缺点 明氏距离的值与各指标的量纲有关，而各指标计量单位的选择有必定的人为性和随意性，各变量计量单位的不一样不只使此距离的实际意义难以说清，并且，任何一个变量计量单位的改变都会使此距离的数值改变从而使该距离的数值依赖于各变量计量单位的选择。 明氏距离的定义没有考虑各个变量之间的相关性和重要性。实际上，明氏距离是把各个变量都同等看待，将两个样品在各个变量上的离差简单地进行了综合。
• 兰氏距离(Lance & Williams) 这是兰思和维廉姆斯(Lance & Williams)所给定的一种距离，这是一个自身标准化的量，因为它对大的奇异值不敏感，这样使得它特别适合于高度偏倚的数据。虽然这个距离有助于克服明氏距离的第一个缺点，但它也没有考虑指标之间的相关性。
• 马氏距离(Mahalanobis) 这是印度著名统计学家马哈拉诺比斯(P.C． Mahalanobis)所定义的一种距离，里面包含了协方差矩阵，Σ表示观测变量之间的协方差短阵。在实践应用中，若整体协方差矩阵Σ未知，则可用样本协方差矩阵做为估计代替计算。马氏距离又称为广义欧氏距离。显然，马氏距离与上述各类距离的主要不一样就是马氏距离考虑了观测变量之间的相关性。若是假定各变量之间相互独立，即观测变量的协方差矩阵是对角矩阵，则马氏距离就退化为用各个观测指标的标准差的倒数做为权数进行加权的欧氏距离。所以，马氏距离不只考虑了观测变量之间的相关性，并且也考虑到了各个观测指标取值的差别程度。
• 斜交空间距离 因为各变量之间每每存在着不一样的相关关系，用正交空间的距离来计算样本间的距离易变形，因此能够采用斜交空间距离。当各变量之间不相关时，斜交空间退化为欧氏距离。
• 配合距离 适用于分类变量，尤为是名义尺度变量。 类似系数 研究样品间的关系经常使用距离，研究指标间的关系经常使用类似系数。类似系数经常使用的有夹角余弦和相关系数。
• 夹角余弦(Cosine) 夹角余弦是从向量集合的角度所定义的一种测度变量之间亲疏程度的类似系数。
• 相关系数是将数据标准化后的夹角余弦 距离和类似系数选择的原则 通常说来，同一批数据采用不一样的亲疏测度指标，会获得不一样的分类结果。产生不一样结果的缘由，主要是因为不一样的亲疏测度指标所衡量的亲疏程度的实际意义不一样，也就是说，不一样的亲疏测度指标表明了不一样意义上的亲疏程度。所以咱们在进行聚类分析时，应注意亲疏测度指标的选择。一般，选择亲疏测度指标时，应注意遵循的基本原则主要有：所选择的亲疏测度指标在实际应用中应有明确的意义。如在经济变量分析中，经常使用相关系数表示经济变量之间的亲疏程度。适当地考虑计算工做量的大小。如对大样本的聚类问题，不适宜选择斜交空间距离，因采用该距离处理时，计算工做量太大。 亲疏测度指标的选择要综合考虑已对样本观测数据实施了的变换方法和将要采用的聚类分析方法。如在标准化变换之下，夹角余弦实际上就是相关系数；又如若在进行聚类分析以前已经对变量的相关性做了处理，则一般就可采用欧氏距离，而没必要选用斜交空间距离。此外，所选择的亲疏测度指标，还须和所选用的聚类分析方法一致。如聚类方法若选用离差平方和法，则距离只能选用欧氏距离。 样品间或变量间亲疏测度指标的选择是一个比较复杂且带主观性的问题，咱们应根据研究对象的特色做具体分析，以选择出合适的亲疏测度指标。实践中，在开始进行聚类分析时，不妨试探性地多选择几个亲疏测度指标，分别进行聚类，而后对聚类分析的结果进行对比分析，以肯定出合适的亲疏测度指标。

5 类与类之间的距离

• 最短距离法(Nearest Neighbor)
• 最长距离法(Furthest Neighbor)
• 重心法(Centroid method)
• 类平均法(average linkage)——组间链接法(Between-groups Linkage)和组内链接法(Within-groups Linkage)
• Ward离差平方和法(Ward's minimumvariance method)

6 系统聚类(hierarchical clustering method)

系统聚类的步骤

（1） 开始将n个样品各做为一类。 （2） 根据样品的特征，选择合适的距离公式，计算n个样品两两之间的距离，构成距离矩阵。 （3） 选择距离矩阵中最小的非对角线元素$d_{pq}$，将相应的两类$G_p$和$G_q$合并为一新类$G_r$={$G_p, G_q$}。 （4） 利用递推公式计算新类与当前各种的距离。 分别删除原矩阵的第p，q行和第p，q列，并新增一行和一列添上的结果，产生新的距离矩阵。 （5） 再合并、计算，直至只有一类为止。 （6） 画聚类图，解释。

• 最短距离法
• 最长距离法
• 中间距离法
• 可变方法
• 重心法
• 类平均方法
• 可变类平均法
• 离差平方和法 相似于方差分析的想法，若是类分得恰当，同类内的样品之间的离差平方和应较小，而类间的离差平方和应当较大。当k固定时，选择使SST达到最小的分类。分类可能指数级增加，寻找最优难以完成。离差平方和法的思路：先让n个样品各自成一类，而后缩小一类，每缩小一类离差平方和就要增大，选择使SST增长最小的两类合并，直到全部的样品归为一类为止（局部最优）。 以上聚类方法的计算步骤彻底相同，仅类与类之间距离的定义不一样。 Lance和Williams于1967年将其统一成对应的形式。

肯定类的个数 从系统聚类的计算机结果能够获得任何可能数量的类。可是，聚类的目的是要使各种之间的距离尽量地远，而类中点的距离尽量的近， 而且分类结果还要有使人信服的解释。每每作系统聚类的时候，大部分状况下咱们都是依靠人的主观判断肯定最后分类的个数。这里给出了一些统计方法来肯定类的个数。

• 给定阈值 经过观测聚类图， 给出一个合适的阈值T。要求类与类之间的距离要超过T值。 例如咱们给定T=0.35， 当聚类时， 类间的距离已经超过了0.35， 则聚类结束。
• 统计量R² 总离差平方和的分解，能够证实：总离差平方和＝组内离差平方和＋组间离差平方和。令T为总离差平方和，令$P_G$为分为G类的组内离差平方和。统计量为R²=1-PG/T其中T是数据的总离差平方和， 是组内离差平方和。R²比较大，说明分G个类时组内离差平方和比较小，也就是说分G类是合适的。可是，分类越多，每一个类的组内离差平方和就越小，R²也就越大；因此咱们只能取合适的G，使得R²足够大，而G自己很小，随着G的增长，R²的增幅不大。好比，假定分4类时，R²=0.8；下一次合并分三类时，降低了许多，R²=0.32，则分4类是合适的。
• 伪F统计量 伪F统计量用于评价聚为G类的效果。若是聚类的效果好，类间的离差平方和相对于类内的离差平方和大，因此应该取伪F统计量较大而类数较小的聚类水平。
• 伪t²统计量 伪t²统计量用来评价合并第K和L类的效果，伪t²统计量大说明不该该合并这两类，应该取合并前的水平。

7 快速聚类(k-means clustering method)

系统聚类法的缺陷——系统聚类法是一种比较经常使用的聚类方法。然而当样本点数量十分庞大时，则是一件很是繁重的工做，聚类的计算速度也比较慢。好比在市场抽样调查中，有4万人就其对衣着的偏好做了回答，但愿能迅速将他们分为几类。这时， 用系统聚类法计算的工做量极大，做出的树状图也十分复杂，不便于分析。 quick cluster method， k-means method也叫动态聚类、逐步聚类、迭代聚类、k-均值聚类，快速聚类适用于大型数据。 用一个简单的例子来讲明快速聚类法的工做过程。 例如咱们要把图中的点分红两类。 快速聚类的步骤： 一、随机选取两个点做为聚核。 二、对于任何点，分别计算到两个点的距离。 三、若到第一个点的距离小于到第二个点的距离，则将该点划为第一类，不然划给第二类。因而得图（b）的两个类。 四、分别计算两个类的重心，则得新的两个点，以其为新的聚核，对空间中的点进行从新分类，获得新分类。

• 选择凝聚点 初始凝聚点（聚类种子、initial cluster seeds/clustercenters）就是一批有表明性的点，是欲造成类的中心。初始凝聚点的选择直接决定初始分类，对分类结果也有很大的影响，因为凝聚点的不一样选择，其最终分类结果也将出现不一样，故选择时要慎重。一般选择初始凝聚点的方法有： 人为选择，当人们对所欲分类的问题有必定了解时，根据经验，预先肯定分类个数和初始分类，并从每一类中选择一个有表明性的样品做为凝聚点。 将数据人为地分为A类，计算每一类的重心，就将这些重心做为凝聚点。 用密度法选择凝聚点。以某个正数d为半径，以每一个样品为球心，落在这个球内的样品数（不包括做为球心的样品）就叫作这个样品的密度。计算全部样品点的密度后，首先选择密度最大的样品做为第一凝聚点，而且人为地肯定一个正数D（通常D＞d，常取D＝2d）。而后选出次大密度的样品点，若它与第一个凝聚点的距离大于D，则将其做为第二个凝聚点；不然舍去这点，再选密度次于它的样品。这样，按密度大小依次考查，直至所有样品考查完毕为止．此方法中，d要给的合适，太大了使凝聚点个数太少，过小了使凝聚点个数太多。 人为地选择一正数d，首先以全部样品的均值做为第一凝聚点。而后依次考察每一个样品，若某样品与已选定的凝聚点的距离均大于d，该样品做为新的凝聚点，不然考察下一个样品。 随机地选择，若是对样品的性质毫无所知，可采用随机数表来选择，打算分几类就选几个凝聚点。或者就用前A个样品做为凝聚点（假设分A类）。这方法通常不提倡使用。
• 衡量聚类结果的合理性指标 该指标为全部K个类中全部元素与其重心的距离的平方和。若分类不合理时，指标会很大，随着分类的过程，逐渐降低并趋于稳定。 算法终止的标准 定义算法终止的标准是：相对偏差小于ε，ε是事前给定的一个充分小量。 快速聚类步骤 第一，选择若干个观测值点为“凝聚点”； 第二，经过分配每一个“凝聚点”最近的类来造成临时分类。每一次对一个观测值点进行归类，“凝聚点”更新 为这一类目前的均值；全部的观测值点分配完后，这些类的“凝聚点”用临时类的均值代替；该步骤能够一直进行直到“凝聚点”的改变很小或为零时止； 第三，最终的分类由分配每个观测到最近的“凝聚点”而造成。

8 有序聚类

有序样本聚类法

有序样本聚类法又称为最优分段法。该方法是由费歇在1958年提出的。它主要适用于样本由一个变量描述，或者将多变量综合成为一个变量来分析的状况。对于有序样本聚类，其实是须要找出一些分点，将它们划分为几个分段，每一个分段看做一类，这样的分类又称分割。分点位置不一样获得的分割不一样，有序样本聚类是要找到一个分割使得各段内部样本差别很小，而各段之间样本的差别很大。有序样本聚类法经常被用于系统的评估问题，被用来对样本点进行分类划级。这种行政上的规定每每是不客观、不合理的。合理的分类应该把发展状况最近似的地区划入同一类。这就是有序样本聚类的工做思路。系统聚类开始n个样品各自自成一类，而后逐步并类，直至全部的样品被聚为一类为止。而有序聚类则相反，开始全部的样品为一类，而后分为二类、三类等，直到分红n类。每次分类都要求产生的离差平方和最小。

有序样本聚类算法步骤

1. 定义类的直径；
2. 定义分类的损失函数L[p(n,k)]； 由损失函数的构造能够看出，损失函数是各种的直径之和。若是分类很差，则各种的直径之和大，不然比较小。当n和k固定时， L[b(n,k)]越小表示各种的离差平方和越小，分类是合理的。所以要寻找一种分法b(n,k)，使分类损失函数L[b(n,k)]达到最小。记该分法为p[n,k]。
3. L[p(n,k)]的递推公式； 若是要找到n个样品分为k个类的最优分割，应创建在将j-1（j＝2,3,…,n)个样品分为k-1类的最优分割的基础上。
4. 寻找最优解。

9 聚类分析的主要步骤

（1）选择变量 和聚类分析的目的密切相关；在不一样研究对象上的值有明显的差别；变量之间不能高度相关。 （2）计算类似性 类似性是聚类分析中的基本概念，它反映了研究对象之间的亲疏程度，聚类分析就是根据对象之间的类似性来分类的。 （3）聚类 选定了聚类的变量，计算出样品或指标之间的类似程度后，构成了一个类似程度的矩阵。这时主要涉及两个问题：选择聚类的方法和肯定造成的类数。 （4）聚类结果的解释和证明 对聚类结果进行解释是但愿对各个类的特征进行准确的描述，给每类起一个合适的名称。这一步能够借助各类描述性统计量进行分析，一般的作法是计算各种在各聚类变量上的均值，对均值进行比较，还能够解释各种差异的缘由。 （5）有关问题 几种聚类方法得到的结果不必定相同，指标聚类采用类似系数，类似系数大或距离小则表示类间关系密切。 （6）变量聚类分析 对于变量聚类分析，聚类分析作完以后，各种中有较多的指标。 为了达到降维的目的， 须要在每类中选出一个表明指标。 具体作法是：假设某类中有k个指标， 首先分别计算类内指标之间的相关指数， 而后计算某个指标与类内其它指标之间相关指数的平均数。

10 R语言中聚类分析实现

R语言自带的聚类分析函数包括了hclust和k-means。因此本篇主要介绍这两个函数的使用。 而首先hclust是基于距离进行的聚类分析，因此事实上在作层次聚类的时候，第一步是先计算距离。 固然前期说明下，这里的样例数据是北京市12个大气污染监测站点在2017年6月7日和6月8日全天的PM2.5数据（数据来自笔者本身写的代码获取而得，调用了环境云的API），样例数据连同完整的代码会在笔记写完后统一给出。

环境云官网：http://www.envicloud.cn/

数据：

dist.pm25<-dist(airnew[,-1],method='euclidean')
heatmap(as.matrix(dist.pm25),labRow=stationname,labcol=F)

咱们作的分析是对一天内24小时下12个站点的PM2.5聚类分析。因此这个问题的多元变量，是不一样时间段的PM2.5值，前期已经把数据结构成功作成矩阵形式，接下来就须要计算距离了。 距离矩阵在R里面是比较好求取的。dist函数。 dist函数的参数事实上有很多，可是其实通常重点用的就是输入矩阵的参数（代码中的airnew[,-1]，-1表明去掉第一列数据（站点名称）），还有计算距离的方式——method。这里选的是欧氏距离。这个参数的可选取值还包括maximum（最大距离）、manhattan（曼哈顿距离）、canberra（兰氏威廉姆斯距离）、binary（定性距离，其实就是配合距离）、minkowski（闵可夫斯基距离——明氏距离）。还有用得多些的参数——diag和upper。diag为TRUE的时候给出对角线上的距离。upper为TURE的时候给出上三角矩阵上的值。默认都是FLASE。 函数计算完以后获得的是一个距离矩阵。咱们用热力图的方式进行可视化，这就是上面的第二句代码。 heatmap函数是个热图可视化函数，要求输入一个矩阵。labRow实际上是输入列名，labcol是与labRow相关，用来映射输入的值的。结果以下图。

计算完矩阵，便可进行聚类分析了。hclust函数的必要参数与前面距离的参数相似——输入矩阵参数，方法参数（这里聚类的方法前面也有提到，这里就不赘述了，有兴趣的能够本身看官方帮助文档）。而聚类完的结果存储在model1里面，用plot便可画出聚类谱系图。事实上，plclust也是相同的做用，参数基本是统一的，labels填写咱们聚类的变量。而聚类完的结果则能够用cutree来得到，输入的model1——聚类结果，k是要求的类数。

model1=hclust(dist.pm25,method="ward")
plot(model1,labels=stationname,hang=-1,las=1)
plclust(model1,labels=stationname,hang=-1)

对聚类结果作个简单可视化。以0点和1点的PM2.5值分别为x和y轴，以聚类结果作划分。

result=cutree(model1,k=3)
plot(airnew[,2],airnew[,3],col=result,pch=as.integer(result))

接下来是K-means的方法。函数并不复杂，输入数据框或者矩阵（作聚类的数据），center就是聚类数，nstart是迭代次数。迭代次数高，聚类可信度高些。后面的这个函数是聚类可视化的函数，是fpc包下面的，使用前请先确认是否安装。

kres<-kmeans(airnew[,-1],centers=3,nstart=10)
plotcluster(airnew[,-1],kres$cluster)

对比了两者的分类结果，是一致的。

聚类结束后，咱们就这个数据和结果作些简单的分析。事实上做为地学人员，咱们就简单地画个站点分布图来对应看看具体状况。从这张图来看，PM2.5的聚类结果显示了它具备很好的空间分异性。固然下面的图有点简陋，给出一个对比的，基于leaflet和R Notebook的交互式小地图（老规矩）。