聚类分析应用实例

为了研究世界各国森林、草原资源的分布规律,共抽取了21个国家的数据,每份数据包括4项指标,原始数据见表5-2。试用该数据对国别进行聚类分析。[大谦MATLAB,dqmatlab点com]

表5-2 原始数据表

国 别 森林面积(万公顷) 森林覆盖率(%) 林木蓄积量(亿立方米) 草原面积(万公顷)
中 国 11978 12.5 93.5 31908
美 国 28446 30.4 202.0 23754
日 本 2501 67.2 24.8 58
德 国 1028 28.4 14.0 599
英 国 210 8.6 1.5 1147
法 国 1458 26.7 16.0 1288
意 大 利 635 21.1 3.6 514
加 拿 大 32613 32.7 192.8 2385
澳大利亚 10700 13.9 10.5 45190
俄 罗 斯 92000 41.1 841.5 37370
捷 克 458 35.8 8.9 168
波 兰 868 27.8 11.4 405
匈 牙 利 161 17.4 2.5 129
克罗地亚 929 36.3 11.4 640
罗马尼亚 634 26.7 11.3 447
保加利亚 385 34.7 2.5 200
印 度 6748 20.5 29.0 1200
印度尼西亚 2180 84.0 33.7 1200
尼日利亚 1490 16.1 0.8 2090
墨 西 哥 4850 24.6 32.6 7450
巴 西 57500 67.6 238.0 15900

MATLAB中提供了两种方法可以进行聚类分析。

(1)一次聚类。它的优点是可利用clusterdata函数对样本数据进行一次聚类。它的缺点是可供用户选择的范围比较窄,不能更改距离的计算方法。

(2)分步聚类。此方法可以分以下步骤进行分步聚类。

① 找到数据集合中变量两两之间的相似性和非相似性,并用pdist函数计算变量之间的距离;

② 用linkage函数定义变量之间的连接;

③ 用cophenetic函数评价聚类信息;

④ 用cluster函数创建聚类。

下面分别介绍这两种聚类方法。

1.一次聚类

在命令窗口中输入下面的命令行:

code.matlab
>> X=[11978	12.5	93.5	31908
28446	30.4	202.0	23754
2501	67.2	24.8	58
1028	28.4	14.0	599
210	8.6	1.5	1147
1458	26.7	16.0	1288
635	21.1	3.6	514
32613	32.7	192.8	2385
10700	13.9	10.5	45190
92000	41.1	841.5	37370
458	35.8	8.9	168
868	27.8	11.4	405
161	17.4	2.5	129
929	36.3	11.4	640
634	26.7	11.3	447
385	34.7	2.5	200
6748	20.5	29.0	1200
2180	84.0	33.7	1200
1490	16.1	0.8	2090
4850	24.6	32.6	7450
57500	67.6	238.0	15900];
>> T = clusterdata(X,0.9)
T =
    10
    10
     3
     8
     2
     9
     8
    10
    10
     7
     1
     8
     1
     8
     8
     1
     4
     9
     9
     5
     6

可见,MATLAB将数据集合分为10类。调整cutoff值,将有不同的分类。

2.分步聚类

(1)寻找变量之间的相似性

用pdist函数可以计算相似性矩阵或非相似性矩阵。在MATLAB中有多种方法可以计算距离。默认时用pdist函数计算欧氏距离,可以指定一种或多种选项(参见函数部分的内容)。在进行计算之前最好先将数据(用zscore函数)进行标准化。下面的代码返回变量之间的距离信息。

code.matlab
>> Y=pdist(X)
Y =
   1.0e+04 *
  1 ~ 131.8376    3.3230    3.3169    3.2935    3.2377    3.3380    3.6020    1.3344    8.0212    3.3766    3.3405    3.3905    3.3163
  14 ~ 263.3444    3.3761    3.1150    3.2233    3.1609    2.5476    4.8255    3.5138    3.5888    3.6172    3.5116    3.6243    2.1771
  27 ~ 392.7829    6.4999    3.6601    3.6135    3.6854    3.5937    3.6287    3.6637    3.1297    3.4621    3.4583    2.8681    3.0097
  40 ~ 520.1570    0.2537    0.1613    0.1922    3.0202    4.5871    9.6969    0.2046    0.1670    0.2342    0.1677    0.1908    0.2121
  53 ~ 650.4398    0.1186    0.2270    0.7756    5.7236    0.0985    0.0812    0.0402    3.1636    4.5628    9.8126    0.0715    0.0251
  66 ~ 780.0986    0.0107    0.0422    0.0757    0.5752    0.1301    0.1561    0.7845    5.8509    0.1256    0.0763    3.2427    4.5275
  79 ~ 919.8682    0.1010    0.0992    0.1019    0.0880    0.0819    0.0963    0.6538    0.1972    0.1590    0.7827    5.9160    0.1130
  92 ~ 1043.1175    4.4864    9.7470    0.1502    0.1062    0.1739    0.0837    0.1177    0.1528    0.5291    0.0730    0.0803    0.7034
  105 ~ 1175.7916    3.2033    4.5796    9.8522    0.0389    0.0257    0.0611    0.0320    0.0068    0.0402    0.6151    0.1692    0.1793
  118 ~ 1300.8116    5.8910    4.8088    6.8929    3.2232    3.1807    3.2531    3.1733    3.2038    3.2303    2.5893    3.0457    3.1125
  131 ~ 1432.8222    2.8320    8.1679    4.6172    4.5852    4.6277    4.5609    4.5861    4.6157    4.4167    4.4808    4.4073    3.8191
  144 ~ 1565.5210    9.8816    9.8347    9.9106    9.8202    9.8548    9.8872    9.2611    9.6833    9.7147    9.2147    4.0640    0.0474
  157 ~ 1690.0300    0.0667    0.0330    0.0080    0.6374    0.2008    0.2182    0.8504    5.9172    0.0759    0.0243    0.0238    0.0525
  170 ~ 1820.5934    0.1535    0.1796    0.8093    5.8714    0.0923    0.0570    0.0236    0.6674    0.2287    0.2369    0.8694    5.9469
  183 ~ 1950.0353    0.0700    0.5846    0.1372    0.1555    0.7858    5.8593    0.0351    0.6160    0.1721    0.1853    0.8174    5.8929
  196 ~ 2080.6441    0.2056    0.2189    0.8515    5.9234    0.4568    0.5333    0.6532    5.2838    0.1129    0.6797    5.7240    0.6326
  209 ~ 2105.7688    5.3324

为了便于阅读,可以用squareform函数将距离向量转换为矩阵:

code.matlab
>> squareform(Y)
ans =
   1.0e+04 *
  1 ~ 130    1.8376    3.3230    3.3169    3.2935    3.2377    3.3380    3.6020    1.3344    8.0212    3.3766    3.3405    3.3905
    1.8376         0    3.5138    3.5888    3.6172    3.5116    3.6243    2.1771    2.7829    6.4999    3.6601    3.6135    3.6854
    3.3230    3.5138         0    0.1570    0.2537    0.1613    0.1922    3.0202    4.5871    9.6969    0.2046    0.1670    0.2342
    3.3169    3.5888    0.1570         0    0.0985    0.0812    0.0402    3.1636    4.5628    9.8126    0.0715    0.0251    0.0986
    3.2935    3.6172    0.2537    0.0985         0    0.1256    0.0763    3.2427    4.5275    9.8682    0.1010    0.0992    0.1019
    3.2377    3.5116    0.1613    0.0812    0.1256         0    0.1130    3.1175    4.4864    9.7470    0.1502    0.1062    0.1739
    3.3380    3.6243    0.1922    0.0402    0.0763    0.1130         0    3.2033    4.5796    9.8522    0.0389    0.0257    0.0611
    3.6020    2.1771    3.0202    3.1636    3.2427    3.1175    3.2033         0    4.8088    6.8929    3.2232    3.1807    3.2531
    1.3344    2.7829    4.5871    4.5628    4.5275    4.4864    4.5796    4.8088         0    8.1679    4.6172    4.5852    4.6277
    8.0212    6.4999    9.6969    9.8126    9.8682    9.7470    9.8522    6.8929    8.1679         0    9.8816    9.8347    9.9106
    3.3766    3.6601    0.2046    0.0715    0.1010    0.1502    0.0389    3.2232    4.6172    9.8816         0    0.0474    0.0300
    3.3405    3.6135    0.1670    0.0251    0.0992    0.1062    0.0257    3.1807    4.5852    9.8347    0.0474         0    0.0759
    3.3905    3.6854    0.2342    0.0986    0.1019    0.1739    0.0611    3.2531    4.6277    9.9106    0.0300    0.0759         0
    3.3163    3.5937    0.1677    0.0107    0.0880    0.0837    0.0320    3.1733    4.5609    9.8202    0.0667    0.0243    0.0923
    3.3444    3.6287    0.1908    0.0422    0.0819    0.1177    0.0068    3.2038    4.5861    9.8548    0.0330    0.0238    0.0570
    3.3761    3.6637    0.2121    0.0757    0.0963    0.1528    0.0402    3.2303    4.6157    9.8872    0.0080    0.0525    0.0236
    3.1150    3.1297    0.4398    0.5752    0.6538    0.5291    0.6151    2.5893    4.4167    9.2611    0.6374    0.5934    0.6674
    3.2233    3.4621    0.1186    0.1301    0.1972    0.0730    0.1692    3.0457    4.4808    9.6833    0.2008    0.1535    0.2287
    3.1609    3.4583    0.2270    0.1561    0.1590    0.0803    0.1793    3.1125    4.4073    9.7147    0.2182    0.1796    0.2369
    2.5476    2.8681    0.7756    0.7845    0.7827    0.7034    0.8116    2.8222    3.8191    9.2147    0.8504    0.8093    0.8694
    4.8255    3.0097    5.7236    5.8509    5.9160    5.7916    5.8910    2.8320    5.5210    4.0640    5.9172    5.8714    5.9469
  14 ~ 213.3163    3.3444    3.3761    3.1150    3.2233    3.1609    2.5476    4.8255
    3.5937    3.6287    3.6637    3.1297    3.4621    3.4583    2.8681    3.0097
    0.1677    0.1908    0.2121    0.4398    0.1186    0.2270    0.7756    5.7236
    0.0107    0.0422    0.0757    0.5752    0.1301    0.1561    0.7845    5.8509
    0.0880    0.0819    0.0963    0.6538    0.1972    0.1590    0.7827    5.9160
    0.0837    0.1177    0.1528    0.5291    0.0730    0.0803    0.7034    5.7916
    0.0320    0.0068    0.0402    0.6151    0.1692    0.1793    0.8116    5.8910
    3.1733    3.2038    3.2303    2.5893    3.0457    3.1125    2.8222    2.8320
    4.5609    4.5861    4.6157    4.4167    4.4808    4.4073    3.8191    5.5210
    9.8202    9.8548    9.8872    9.2611    9.6833    9.7147    9.2147    4.0640
    0.0667    0.0330    0.0080    0.6374    0.2008    0.2182    0.8504    5.9172
    0.0243    0.0238    0.0525    0.5934    0.1535    0.1796    0.8093    5.8714
    0.0923    0.0570    0.0236    0.6674    0.2287    0.2369    0.8694    5.9469
         0    0.0353    0.0700    0.5846    0.1372    0.1555    0.7858    5.8593
    0.0353         0    0.0351    0.6160    0.1721    0.1853    0.8174    5.8929
    0.0700    0.0351         0    0.6441    0.2056    0.2189    0.8515    5.9234
    0.5846    0.6160    0.6441         0    0.4568    0.5333    0.6532    5.2838
    0.1372    0.1721    0.2056    0.4568         0    0.1129    0.6797    5.7240
    0.1555    0.1853    0.2189    0.5333    0.1129         0    0.6326    5.7688
    0.7858    0.8174    0.8515    0.6532    0.6797    0.6326         0    5.3324
    5.8593    5.8929    5.9234    5.2838    5.7240    5.7688    5.3324         0

(2)定义变量之间的连接

将数据中变量之间的近似性计算出来以后,使用linkage函数就可以决定数据集合中哪些变量可以被归入某类中。该函数利用pdist函数生成的距离信息,连接相近的成对变量,形成二分类。然后,linkage函数将新生成的类与其他变量相连接,生成更大的类,直到原始数据中所有的变量被连接到系统聚类树中。

code.matlab
>> Z=linkage(Y)
Z =
   1.0e+04 *
    0.0007    0.0015    0.0068
    0.0011    0.0016    0.0080
    0.0004    0.0014    0.0107
    0.0013    0.0023    0.0236
    0.0012    0.0022    0.0238
    0.0024    0.0026    0.0243
    0.0025    0.0027    0.0330
    0.0006    0.0018    0.0730
    0.0005    0.0028    0.0763
    0.0019    0.0029    0.0803
    0.0030    0.0031    0.0812
    0.0003    0.0032    0.1186
    0.0017    0.0033    0.4398
    0.0020    0.0034    0.6326
    0.0001    0.0009    1.3344
    0.0002    0.0036    1.8376
    0.0008    0.0037    2.1771
    0.0035    0.0038    2.5476
    0.0021    0.0039    2.8320
    0.0010    0.0040    4.0640

在输出信息中,每一行代表一次连接,前两列代表连接的变量,第3列代表连接对象之间的距离。例如,首先将国家7和国家15相连接,形成新类22,它们之间的距离是68;然后将国家16和国家11相连,形成新类23,它们之间的距离为80,比前一组略大;接着,将国家4和国家14相连,形成新类24;再将新类23与国家13相连,形成新类25;依此类推,最后将新类40和国家10相连,形成整个系统。整个输出信息展示了聚类的过程。

(3)评价聚类信息

连接变量生成聚类树以后,用户可能会想修改聚类树或希望了解变量连接的更多信息。下面分别介绍具体实现方法。

① 修改聚类树:衡量通过linkege函数生成的聚类信息的有效性的一种方法是将它与pdist函数生成的原始相似数据对比。如果有效,则二者强相关。可以使用cophenet函数比较这两组值,并且计算它们的相关性,返回cophenetic相关系数。该值越接近于1,表示聚类效果越好。

可以通过cophenetic相关系数来比较用不同距离计算方法或聚类算法对同一套数据得到的聚类结果。对于本例:

code.matlab
>> C=cophenet(Z,Y)
C =
    0.9393

将pdist函数中距离的计算方法分别指定为“Mahal”“SEuclid”和“CityBlock”,重新使用pdist函数计算以后,用cophenet函数计算,得到的C值分别等于0.7486、0.8018和0.9291。均小于距离为欧氏距离时的计算结果,所以使用默认设置时的效果最佳。

② 了解与聚类连接相关的更多信息:决定数据集合中聚类的方法之一是比较聚类树中每一个连接的长度与相邻次一级连接的长度。如果二者相近,则表示在此水平上变量之间是相似的,这些连接被认为具有较高水平的连续性;如果二者不相近,则表示聚类树中变量之间不相似,则称该连接与周围是不连续的。cluster函数用不连续性的度量决定在哪里对数据进行分类。

下面的命令行生成聚类树:

code.matlab
>> dendrogram(Z)
ans =
  151.0055
  153.0055
  154.0055
  155.0055
  156.0055
  157.0055
  158.0055
  159.0055
  160.0055
  161.0055
  162.0055
  163.0055
  164.0055
  165.0055
  166.0055
  167.0055
  168.0055
  169.0055
  170.0055
  171.0055

图5-1中的纵轴为变量之间的距离,横轴为聚类变量。系统聚类树从左到右反映了聚类的先后次序。该图可以与linkage函数的计算结果相对应。

Document Image
\[\]

图5-1 本问题的冰柱图

在系统聚类树中,每一个连接的相对连续性可用不连续性系数定量表达,该函数用于比较某个连接的长度与相邻连接长度的均值。若该变量与周围变量连续,则不连续性系数较低,反之则反。

用inconsistent函数可以生成聚类树中每一个连接的不连续系数表。该函数可以比较每一个连接与以下两个水平的连接,这里称为比较的深度。用该函数可以指定其他深度。聚类树底部的对象被称为叶节点,不连续系数为0。

code.matlab
>> I=inconsistent(Z)
I =
  1.0e+004*
    0.0068         0   0.0001         0
    0.0080         0   0.0001         0
    0.0107         0   0.0001         0
    0.0158   0.0110   0.0002   0.0001
    0.0153   0.0120   0.0002   0.0001
    0.0196   0.0077   0.0003   0.0001
    0.0270   0.0053   0.0003   0.0001
    0.0730        0   0.0001      0
    0.0546   0.0306   0.0002   0.0001
    0.0766   0.0052   0.0002   0.0001
    0.0793   0.0026   0.0003   0.0001
    0.0999   0.0265   0.0002   0.0001
    0.2792   0.2271   0.0002   0.0001
    0.5362   0.1363   0.0002   0.0001
    1.3344        0   0.0001       0
    1.5860   0.3559   0.0002   0.0001
    2.0074   0.2401   0.0002   0.0001
    1.7858   1.0157   0.0003   0.0001
    2.6898   0.2011   0.0002   0.0001
    3.4480   0.8711   0.0002   0.0001

在矩阵中,第1列为所有连接长度的均值,第2列为所有连接长度的标准差,第3列为计算中所包含的连接数,第4列为不连续系数。该输出信息可以与linkage函数的输出信息进行对比,第1行代表国家7和国家15相连,因为二者均为叶节点,所以不连续系数为0。

(4)创建聚类

创建二分聚类树以后,可以用cluster函数将系统分成更大的类。有以下两种方式可以创建聚类:

① 找到原始数据集合的自然分界:如果用cluster函数对数据进行分类,指定一个不连续系数0.9作为cutoff参数的值,则cluster函数将所有变量分为10类,如下所示。

code.matlab
>> T=cluster(Z, 0.9)
T =
    10
    10
     3
     8
     2
     9
     8
    10
    10
     7
     1
     8
     1
     8
     8
     1
     4
     9
     9
     5
     6

如果选择不连续系数0.7作为cutoff参数的值,则所有数据分为16类。

code.matlab
>> T=cluster(Z, 0.7)
T =
    16
     9
     6
     3
     4
    15
    14
    10
    16
    12
    13
     2
     1
     3
    14
    13
     7
    15
     5
     8
    11

② 通过指定类数进行聚类:将所有变量归为3类,将cutoff参数的值设置为3,用cluster函数进行分类。

code.matlab
>> T=cluster(Z,3)
T =
     2
     2
     2
     2
     2
     2
     2
     2
     2
     3
     2
     2
     2
     2
     2
     2
     2
     2
     2
     2
     1