SAS学习系列35 聚类分析.docx
《SAS学习系列35 聚类分析.docx》由会员分享，可在线阅读，更多相关《SAS学习系列35 聚类分析.docx（26页珍藏版）》请在冰点文库上搜索。
SAS学习系列35聚类分析
35.聚类分析
（一）概述
聚类分析，相当于“物以类聚”，用于对事物的类别面貌尚不清楚，甚至在事前连总共有几类都不能确定的情况下对数据进行分类。
而判别分析，必须事先知道各种判别的类型和数目，并且要有一批来自各判别类型的样本，才能建立判别函数来对未知属性的样本进行判别和归类。
聚类分析是把分类对象按一定规则分成组或类，这些组或类不是事先给定的而是根据数据特征而定的。
在同类的对象在某种意义上倾向于彼此相似，而在不同类里的这些对象倾向于不相似。
根据这种相似性的不同定义，聚类分析也有不同的方法。
聚类分析分为：
对样品的聚类，对变量的聚类。
样品聚类：
其统计指标是类与类之间距离，把每一个样品看成空间中的一个点，用某种原则规定类与类之间的距离，将距离近的点聚合成一类，距离远的点聚合成另一类。
变量聚类：
其统计指标是相似系数，将比较相似的变量归为一类，而把不怎么相似的变量归为另一类，用它可以把变量的亲疏关系直观地表示出来。
（二）原理
一、距离和相似系数
1.距离
设有n组样品，每组样品有p个变量的数据如下：
变量
样品
V1
V2
…
Vp
X1
x11
x12
…
x1p
X2
x21
x22
…
x2p
…
…
…
…
…
Xn
xn1
xn2
…
xnp
例如，Xi到Xj的闵科夫斯基距离定义为：
q=2时为欧几里得距离；还有马氏距离：
dij=（Xi-Xj）TS-1（Xi-Xj）
其中，Xi=（xi1,…,xip），S-1为n个样品的p×p的协方差矩阵的逆矩阵。
注：
马氏距离考虑了观测变量之间的相关性和变异性（不再受各指标量纲的影响）。
距离选择的基本原则：
（1）要考虑所选择的距离公式在实际应用中有明确的意义。
如欧氏距离就有非常明确的空间距离概念。
马氏距离有消除量纲影响的作用。
（2）要综合考虑对样本观测数据的预处理和将要采用的聚类分析方法。
如在进行聚类分析之前已经对变量作了标准化处理，则通常就可采用欧氏距离。
（3）应根据研究对象的特点不同做出具体分折。
实际中，聚类分析前不妨试探性地多选择几个距离公式分别进行聚类，然后对聚类分析的结果进行对比分析，以确定最合适的距离测度方法。
2.相似系数
变量间的相似性，可以从它们的方向趋同性或“相关性”进行考察，从而得到“夹角余弦法”和“相关系数”两种度量方法。
（1）夹角余弦
Xi=（xi1,…,xip）和Xj作为p维向量，按通常内积中夹角的概念来定义
（2）相关系数
显然|cij|≤1，借助于相似系数，可以定义变量之间的距离。
例如，采用非相似测度距离为dij=1-|cij|,或dij2=1-rij2.
二、类的特征、类间的距离
1.类的定义
给定某阈值T，①若dij≤T，对任意i,j∈G；或者
②若（∑j∈Gdij）/（k-1）≤T，对任意i∈G；或者
③若对任意i∈G，都存在j∈G使得dij≤T；
则称G为一个类。
条件从强到弱的顺序为：
①②③
2.类的特征
设类G的元素为x1,…,xm,m为G中样品数。
描述类的特征有：
（1）均值或重心
（2）样品协方差阵
（3）直径
3.类间的距离
设Gp和Gq中分别有p和q个样品，它们的重心分别记为
和
.
（1）最短距离
它不对类的形状加以限制，保证了对拉长和不规则类的检测，例如（b）（c）（e）（d）形式的类；但它却牺牲了恢复压缩类的性能，另外它也趋向于在分开主要类之前去掉分布的尾部。
（2）最长距离
严重地倾向于产生直径粗略相等的类，而且可能被异常值严重地扭曲，例如（a）.
（3）重心法距离
即两个重心之间欧氏距离的平方。
重心法在处理异常值上比其他谱系方法更稳健，但是在其他方便不如Ward或类平均距离法的效果好。
（4）类平均距离
Gp和Gq中所有两个样品对之间距离的平均。
类平均距离法趋向于合并具有较小偏差的类，而且稍微有点倾向与产生相同方差的类。
例如（a）.
（5）Ward最小方差法或Ward离差平方和距离
D表示类的直径，也有
若样品间的距离采用欧氏距离，上式可表示为
该方法在每次合并类Gp和Gq为Gp+q时，总是选择使得合并后的Dw（p,q）值最小的类Gp和Gq（最小方差法）.合并后增加的最小方差Dw（p,q）除以合并后总的离差平方和TSS的比值（即半偏R2）的统计意义是容易解释的。
Ward方法趋向于合并具有少量观察的类，并倾向于形成具有大约相同数目观察的类，一般是在多元正态混合型、等球形协方差、等抽样概率假设下合并类。
例如（a）.Ward方法对异常值也很敏感。
（6）密度估计法
例如（b）.使用一种基于密度估计的新的非相似测度d*来计算样品xi和xj的近邻关系；
然后根据基于d*方法计算的距离，采用最小距离法进行聚类。
有三种不同的密度估计法：
k最近邻估计法、均匀核估计法、Wong混合法（适用于大数据集而不适用于小数据集）。
（7）两阶段密度估计法
在密度估计法中，众数类经常在尾部所有点聚类之前就已经被合并掉了。
两阶段密度估计法对密度估计法这一缺点进行了修正，以确保众数类被合并之前，所有点能被分配到众数类中去。
在第一阶段中，互不相交的众数类被生成，只有在两个类中至少有一个类中的样品个数少于proccluster过程的选项“mode=值”时才能合并。
在第二阶段中，把前面形成的众数类再按最短距离法谱系地合并（某一类能够完全地包含在别的类里，但在类之间不允许其他形式的重迭）。
（8）最大似然法
最大似然估计法是在三个基本假设：
多元正态混合分布、等球面协方差和不等抽样概率的前提下，使得在谱系的每个水平上似然值最大。
该方法除了偏向于生成不等大小的类外，与Ward最小方差法很相似。
proccluster过程的偏度选项penalty=p值可用于调整偏向的程度。
修正公式可将公式中的2改为p值。
三、类的统计量
（1）类Gp的均方根标准差
即类内的离差平方和（定义为直径），点越分散直径越大。
（2）R2
R2用于评价每次合并成m个类时的聚类效果，R2越大说明聚类效果（能被分开）越好。
R2=1，说明类被完全分开；当所有样品最后合并成一类时，R2=0，说明类被完全混合在一起了，分不开。
R2的值总是随着分类个数的减少而变小。
那么如何根据R2的值来确定所有样品应分为几类最合适呢？
首先，最合适分类的R2的值不能太小，最好能达到0.7以上；其次，不能孤立地看合并后R2值的大小，应该看R2值的变化。
类的合并总是使R2值变小，如果这种变小是均匀的，合并是应该的，但当某次合并使R2值减少很多，这次合并就不应该，即我们找到了最合适的分类数。
例如，从10类合并到5类时，R2值每次减少大约在0.01到0.02之间，R2值减小到0.8，从5类合并到4类时，R2值减少了0.3，即R2=0.8-0.3=0.5，这时，我们就可以得出分5类是最合适的。
（3）半偏R2统计量
合并类Gp和类Gq为类Gm时，可以用半偏R2统计量评价这次合并的效果：
表示合并后类内离差平方和的增量（=上次R2-本次R2）。
（4）伪F统计量
其中，G为聚类的个数，n为观察总数，v为观察样品的变量维数。
该统计量用于评价分为G个聚类的效果。
伪F值越大越表示这些观察可显著地分为G个类。
（5）伪t2统计量
该统计量用以评价合并类Gp和类Gq为类Gm的效果。
该值很大说明合并后离差平方和的增加量相对于原来两类的类内离差平方和很大。
这表明上一次合并的两个类Gp和Gq是很分开的，也就是说上一次聚类的效果是很好的。
否则，这两个类应该在再合并时的离差平方和的增加量很小。
四、类数的确定
至今无太满意的结果。
（1）一般情况下，是看每个变量的R2值及所有变量的累积R2值，观察这些R2值与聚类个数所画的折线图形，以此来判断；
（2）通过对数据的典型变量（可用proccandisc典型判别分析选出）做散点图来判断；
（3）立方聚类准则CCC；
通过数据模拟研究，发现三个最好的准则是：
①伪F统计量；②伪t2统计量；③立方聚类准则CCC。
在cluster过程中这三个准则都被应用和输出，而在fastclus过程中仅伪F统计量和CCC统计量被应用和输出。
建议寻找这三种统计量之间的一致性，即：
CCC统计量和F统计量的局部峰值所对应的聚类数，与这个聚类数伪t2统计量的一个较小值和下一个聚类数的一个较大伪t2统计量相吻合。
还必须强调这些准则仅仅适用于紧凑的或略微拉长的类，也较好地适合于略微正态分布的类。
（三）聚类的方法
一、系统聚类法
目前使用最多的一种方法。
基本步骤是：
（1）首先各样品自成一类，这样对
组样品就相当于有
类（即一类包括一个样品）；
（2）计算各类间的距离（可选8种距离见前文），将其中最近的两类进行合并；
（3）计算新类与其余各类的距离，再将距离最近的两类合并；
（4）重复上述的步骤，直到所有的样品都聚为一类时为止。
注意：
①在整个聚类的过程中，如果在某一步的距离矩阵中最小元素不止一个时，则可以将其同时合并；
②系统聚类法由于每次都计算距离，计算量（内存占用）较大；
③研究表明：
最好综合特性的聚类方法为类平均法或Ward最小方差法，而最差的则为最短距离法。
Ward最小方差法倾向于寻找观察数相同的类。
类平均法偏向寻找等方差的类。
具有最小偏差的聚类方法为最短距离法和密度估计法。
拉长的或无规则的类使用最短距离法比其他方法好。
最没有偏见的聚类方法为密度估计法。
二、K均值聚类法（动态聚类法）
基本步骤：
（1）选择K个观察作为“凝聚点”或称类的中心点，作为这些类均值的第一次猜测值，将所有的样品分成K个初始类；
（2）通过欧氏距离将某个样品划入离中心最近的类中，并对获得样品与失去样品的类，重新计算中心坐标；
（3）重复步骤2，直到所有的样品都不能再分配时为止（反复迭代，直到再也无法降低样本与重心离差之和为止）。
优点：
速度快省内存；缺点：
需要事先确定类数。
注意：
该聚类法对异常点很敏感，通常会把异常点分配为单独的一个类（可用于检测异常点procfastclus）。
（三）SAS实现
SAS提供了5个聚类过程，即cluster，fastclus，modeclus、varclus和tree过程。
cluster为系统聚类过程，可使用十一种聚类方法进行谱系聚类，包括重心法、Ward离差平方和法和欧氏平均距离法等。
fastclus为动态聚类过程，使用K-均值算法寻找不相交的聚类，适宜于大样本分析，观察值可多达10万个。
modeclus为动态聚类过程，使用非参数密度估计法寻找不相交的聚类。
varclus过程可用于系统或动态聚类，通过斜交多组分量分析对变量作“谱系的”和“不相交的”两种聚类。
cluster过程、fastclus过程和modeclus过程用于对样品聚类，varclus过程用于对变量聚类。
tree过程将cluster或varclus过程的聚类结果画出树形结构图或谱系图。
一、PROCCLUSTER过程步
基本语法：
procclusterdata=数据集 <可选项>;
var变量列表;
id变量;
freq变量;
copy变量列表;
rmsstd变量;
by变量列表;
说明：
（1）可选项
outtree=输出数据集
供tree过程调用，用来输出聚类结果的树状图；
method=算法
ward（离差平方和法），average（类平均法），centroid（重心法），complete（最长距离法），single（最短距离法），median（中间距离法），density（密度法），flexible（可变类平均法），twostage（两阶段密度法），eml（最大似然法），mcquitty（相似分析法）；
standard/std——对变量实施标准化；
nonorm——阻止距离被正态化成均数为1或均方根为1；
nosquare——阻止过程在method=average/centroid/median/ward方法中距离数据被平方；
mode=n——当合并两个类时，规定对被指定的众数类中的每个类至少有n个成员，该选项只能在method=density/twostage时使用；
penalty=p——指定用于method=eml中的惩罚系数（p>0,默认p=2）；
trim=p——要求从分析中删去那些概率密度估计较小的点（0≤p<100,被当作百分比），在method=ward/complete时，因为类可能被异常值严重地歪曲，最好使用这个选项（也可用于method=single）；
dim=n——用于method=density/twostage时指定使用的维数（n≥1），若是坐标数据，缺省值为变量个数；若是距离数据，缺省值为1；
hybrid——要求用Wong混合聚类方法，其中密度用k均值法的初始聚类分析中的均值计算得到。
这个选项只能在method=density/twostage时使用；
k=n——指定k最近邻估计法中近邻的个数（2≤n<观察数）；
r=n——指定均匀核密度估计法的支撑球半径（n>0）；
notie——阻止cluster过程在聚类历史过程中检查每次产生的类间最小距离连结（ties）的情况，此选项可以减少过程执行的时间和空间；
rsquare——输出R2和半偏R2；
rmsstd——输出每一类的均方根标准差；
ccc——输出在均匀的原假设下判断聚类分成几类合适的立方聚类准则统计量ccc和近似期望值R2；同时输出选项rsquare有关的R2和半偏R2；此选项不适合于method=single（容易删掉分布的结尾部分）；
pseudo——输出伪F统计量（PSF）和伪t2统计量（PST2），当分类数目不同时，它们有不同的取值；
simple——输出简单统计数；
在输出报表中，可以根据输出的ccc、psf和pst2统计量确定多少分类数较合适，当ccc和psf值出现峰值所对应的分类数较合适，而pst2值是在出现峰值所对应的分类数减1较合适。
copy语句——指定输入数据集中的一些变量拷贝到outtree=的输出数据集中；
rmsstd语句——当输入数据集中的坐标数据代表类的均值时，定义表示均方根标准差变量，通常与freq语句中的变量配合使用。
例1用cluster和tree过程1976年74个国家每100000人口的出身率和死亡率数据（C:
\MyRawData\populations.txt）进行聚类分析。
代码：
datavital;
infile'C:
\MyRawData\populations.txt';
title'CRUDEBIRTHANDDEATHRATESIN1976';
inputcountry$19.birthdeath;
run;
procprintdata=vital;
run;
procplotformchar='|-----|--'hpct=100vpct=100;
plotdeath*birth/hpos=86vpos=26;
run;
procclusterdata=vitalout=treemethod=averageprint=15cccpseudo;
varbirthdeath;
title'ClusterAnalysisofBirthandDeathRatesin74Countries';
run;
proctreedata=treenoprintout=outncl=8;
copybirthdeath;
run;
procplotformchar='|-----|--';
plotdeath*birth=cluster/hpos=86vpos=26;
title2'PlotofClustersFormMethod=Average';
run;
运行结果及说明：
A表示此处有一个散点，而B表示此处有二个重叠散点，C表示此处有三个重叠散点；直观上看散点图可以区分成三类：
低出身率和低死亡率的国家（富裕发达的国家），高出身率和高死亡率的国家（贫穷落后的国家），中出身率和低死亡率的国家（中等发展中的国家）。
输出包括协方差阵的特征值，两相邻特征值之差，各特征值占总方差的百分比，和累计百分比。
聚类是从74国家作为74类，先聚类成73=74-1类开始，直至将所有国家聚为１类的聚类过程和统计数。
“print=15”，聚类数列只输出显示聚类成15类后各合并过程；
“连接聚类列”下的两列为每一次聚类成1个新类的2个样品，标有obn表示是原始样品中的第n号样品，标有cln表示是在上面聚类过程中已经聚成的第n类；
例如，第14新类（聚类数=14）由前面形成的第20旧类（cl20）和原始样品的第8号（ob41）合并产生，共有4个（freq=4）原始样品；同时也可求出前面形成的cl20旧类有4-1=3个原始样品；最后由旧的cl2类和cl3类合并形成1类，其中cl2类由cl4类和cl5类合并，cl3类由cl6类和cl14类合并。
频数列为新类中所含的样品数；
半偏R2，表示每一次合并对信息的损失程度；
R2，表示累计聚类结果；
近似期望R2为在均匀零假设下近似期望的R2；
ccc、psf、psf2列为3个判断分成几类较合适的统计量；
NormRMSDist列为距离被正态化后两样品或类间的平均距离；
Tie（连接）列里用T来指出最小距离的一个连结，而空白表明没有连结；
“out=tree”将聚类过程和聚类结果及相应的统计量输出到tree数据集中。
关于如何选择合适的聚类数。
主要方法是查看“三次聚类准则值CCC”在ncl聚类数为多少时达到峰值，以及“伪F值PSF”和“伪t2值PST2”在ncl聚类数为多少时达到峰值。
同时可以通过再合并成新类时R2值的减少最多（即每一次合并时信息的损失程度“半偏R2值”最大）来验证，还可以直接通过树状图决定分类数。
通常CCC、PSF、PST2值在不同的ncl处有多处峰值，为便于比较和观测，将CCC、PSF和PST2统计量按聚类数大小顺序连接成曲线（则必须先对tree数据集中表示聚类数的字段_ncl_进行排序），在同一张图中绘制相应的散布图。
代码：
datavital;
infile'C:
\MyRawData\populations.txt';
inputcountry$19.birthdeath;
run;
procclusterdata=vitalout=treemethod=averageprint=15cccpseudo;
varbirthdeath;
title'ClusterAnalysisofBirthandDeathRatesin74Countries';
run;
procprintdata=tree;
title'tree';
run;
procsortdata=tree;
by_ncl_;
goptionsreset=globalgunit=pctcback=whiteborderhtitle=6htext=3ftext=swissbcolors=（back）;
procgplotdata=tree;
plot_ccc_*_ncl_/haxis=axis1vaxis=axis2;
plot2_psf_*_ncl__pst2_*_ncl_/vaxis=axis3overlay;
symbol1i=joinv='C'l=3h=2.5w=1.5;
symbol2i=joinv='F'l=1h=2.5w=1.5;
symbol3i=joinv='T'l=1h=2.5w=1.5;
axis1label=（'NumberofClusters'）order=（1to15by1）;
axis2label=（'CCC'）order=（0to6by0.5）;
axis3label=（'PSFandPST2'）order=（0to300by25）;
title2"Plotof_CCC_*_NCL_.Symbolusedis'C'";
title3"Plotof_PSF_*_NCL_.Symbolusedis'F'";
title4"Plotof_PST2_*_NCL_.Symbolusedis'T'";
run;
运行结果：
CCC值在ncl=2,4,8,12处有峰值，这些是立方聚类准则认为可考虑取的分类数，但ncl=8时更陡些；
psf在ncl=3,8,12处有峰值，这些是根据伪F统计量可考虑取的分类数，但ncl=8时更陡些；
pst2在ncl=5,9,13处有峰值，通常伪t2统计量认为取再减1的分类较合适，或者直接取伪t2统计量最小值处的分类数，因为再次合并将会有比较大伪t2统计量，所以根据伪t2统计量可考虑的分类数为4,8,12;
从半偏R2值分析，从12类合并成11类时，有较大的信息损失值0.003198，从8类合并成7类时，有较大的信息损失值0.009841，这项指标特别指出从3类合并成2类时，有仅次于分成1类的信息损失值0.096203，半偏R2指标值支持分成2,3,8,12类；
最后，从带有聚类过程和合并距离的树状图中，也能较直观地分成2,3,4,8类；
综合结论，可以把74国家分成2类、3类、4类、8类、12类，尤其要注意的是分成8类，因为这种分类已从聚类运算结果分析中显示了良好的指标性能。
注：
完全分开的球状类，适合用average（类平均法）；
没有完全分开的球状类，适合用ward（离差平方和法）；
大小和离散程度不等的多元正态类，适合用eml（最大似然法）；
拉长的多元正态类，适合用twostage,k=10（两阶段密度法）；
中心相同大小不等的圆环类，适合用single（最短距离法）。
二、PROCFASTCLUS过程步
选择“凝聚点”，根据距离接近的原则归类，适合对大样品数据做快速的“粗”分类。
初始“凝聚点”可以指定或由过程自动选择；可以指定最大的“凝聚点”个数，及“凝聚点”之间的必须分隔开的最小距离。
注意：
fastclus过程对异常点很敏感，通常会把异常点分配为单独的一个类，因此也是检测异常点的有效过程。
基本语法：
procfastclusdata=数据集<可选项>;
var变量列表;
id变量;
freq变量;
weight变量;
by变量列表;
可选项：
seed=输入数据集——指定该数据集作为初始“凝聚点”；
mean=输出数据集——包含每个类的均值和一些统计量；
out=输出数据集——包含原始数据和分类标志的聚类结果；
maxclusters/maxc=n——指定允许的最大聚类数（默认为100）；
radius=t——为选择新的“凝聚点”指定最小距离准则；当一个观察点与已有“凝聚点”的最小距离都大于t值时，该观察可
升级会员 