本章代碼
聚類分析是一種數(shù)據(jù)歸約技術迫筑，旨在揭露一個數(shù)據(jù)集中觀測值的子集。它可以把大量的觀測 值歸約為若干個類叮喳。這里的類被定義為若干個觀測值組成的群組睁宰，群組內觀測值的相似度比群間 相似度高。這不是一個精確的定義嘱么，從而導致了各種聚類方法的出現(xiàn)狮含。
聚類分析被廣泛用于生物和行為科學、市場以及醫(yī)學研究中曼振。例如几迄，一名心理學研究員可能 基于抑郁癥病人的癥狀和人口統(tǒng)計學數(shù)據(jù)對病人進行聚類，試圖得出抑郁癥的亞型冰评，以期通過亞型來找到更加有針對性和有效的治療方法映胁，同時更好地了解這種疾病。營銷研究人員根據(jù)消費者 的人口統(tǒng)計特征與購買行為的相似性制定客戶細分戰(zhàn)略集索，并基于此對其中的一個或多個子組制定 相應的營銷戰(zhàn)略屿愚。醫(yī)學研究人員通過對DNA微陣列數(shù)據(jù)進行聚類分析來獲得基因表達模式汇跨，從而 幫助他們理解人類的正常發(fā)育以及導致許多疾病的根本原因。
最常用的兩種聚類方法是層次聚類（hierarchical agglomerative clustering）和劃分聚類 （partitioning clustering）妆距。在層次聚類中穷遂，每一個觀測值自成一類，這些類每次兩兩合并娱据，直到所有的類被聚成一類為止蚪黑。在劃分聚類中，首先指定類的個數(shù)K中剩，然后觀測值被隨機分成K類忌穿，再重新形成聚合的類。
這兩種方法都對應許多可供選擇的聚類算法结啼。對于層次聚類來說掠剑，最常用的算法是單聯(lián)動 （single linkage）、全聯(lián)動（complete linkage ）郊愧、平均聯(lián)動（average linkage）朴译、質心（centroid）和 Ward方法。對于劃分聚類來說属铁，最常用的算法是K均值（K-means）和圍繞中心點的劃分（PAM）眠寿。 每個聚類方法都有它的優(yōu)點和缺點，我們將在本章討論焦蘑。
這一章的例子圍繞食物和酒（也是我的愛好）盯拱。我們對 flexclust 包中的營養(yǎng)數(shù)據(jù)集 nutrient 作層次聚類，以期回答以下問題例嘱。
- 基于五種營養(yǎng)標準的27類魚狡逢、禽、肉的相同點和不同點是什么拼卵？
- 是否有一種方法能把這些食物分成若干個有意義的類甚侣？
我們再用劃分聚類來分析178種意大利葡萄酒樣品的13種化學成分。數(shù)據(jù)在 rattle 包的 wine 數(shù)據(jù)集中间学。這里要解決的問題如下：
- q 這些意大利葡萄酒樣品能繼續(xù)分成更細的組嗎殷费？
- 如果能，有多少子組低葫？它們的特征是什么详羡？
事實上，樣品中共有三個品種的酒（記為類）嘿悬。這可以幫助我們評估聚類分析能否辨別這一 結構实柠。
盡管聚類方法種類各異，但是它們通常遵循相似的步驟善涨。我們在16.1節(jié)描述了這些步驟窒盐。16.2 節(jié)主要探討層次聚類草则，16.3節(jié)則探討劃分聚類。最后蟹漓，16.4節(jié)提出了一些相關建議炕横。為了保證本 章的代碼能正常運行，你必須確保安裝了以下軟件包： cluster 葡粒、 NbClust 份殿、 flexclust 、 fMultivar 嗽交、 ggplot2 和 rattle 卿嘲。第17章也將用到 rattle 包。
1 聚類分析的一般步驟
像因子分析（第14章）一樣夫壁，有效的聚類分析是一個多步驟的過程拾枣，這其中每一次決策都可能影響聚類結果的質量和有效性。本節(jié)介紹一個全面的聚類分析中的11個典型步驟盒让。
(1) 選擇合適的變量
第一（并且可能是最重要的）步是選擇你感覺可能對識別和理解數(shù)據(jù)中不同觀測值分組有重要影響的變量放前。例如，在一項抑郁癥研究中糯彬，你可能會評估以下一個或多個方面：心理學癥狀，身體癥狀葱她，發(fā)病年齡撩扒，發(fā)病次數(shù)、持續(xù)時間和發(fā)作時間吨些，住院次數(shù)搓谆，自理能力，社會和工作經歷豪墅，當前的年齡泉手，性別，種族偶器，社會經濟地位斩萌，婚姻狀況，家族病史以及對以前治療的反應屏轰。高級的聚類方法也不能彌補聚類變量選不好的問題颊郎。
(2) 縮放數(shù)據(jù)
如果我們在分析中選擇的變量變化范圍很大，那么該變量對結果的影響也是最大的霎苗。這往往是不可取的姆吭，分析師往往在分析之前縮放數(shù)據(jù)。最常用的方法是將每個變量標準化為均值為0和標準差為1的變量唁盏。其他的替代方法包括每個變量被其最大值相除或該變量減去它的平均值并除以變量的平均絕對偏差内狸。這三種方法能用下面的代碼來解釋：
df1 <- apply(mydata, 2, function(x){(x-mean(x))/sd(x)})
df2 <- apply(mydata, 2, function(x){x/max(x)})
df3 <- apply(mydata, 2, function(x){(x – mean(x))/mad(x)})
在本章中检眯，你可以使用 scale() 函數(shù)來將變量標準化到均值為0和標準差為1的變量。這和第一個代碼片段（ df1 ）等價昆淡。
(3) 尋找異常點
許多聚類方法對于異常值是十分敏感的锰瘸，它能扭曲我們得到的聚類方案。 你可以通過 outliers包中的函數(shù)來篩選（和刪除）異常單變量離群點瘪撇。mvoutlier 包中包含了能識別多元變量的離群點的函數(shù)获茬。一個替代的方法是使用對異常值穩(wěn)健的聚類方法，圍繞中心點的劃分（16.3.2節(jié)）可以很好地解釋這種方法倔既。
(4) 計算距離
盡管不同的聚類算法差異很大恕曲，但是它們通常需要計算被聚類的實體之間的距離。兩個觀測值之間最常用的距離量度是歐幾里得距離渤涌，其他可選的量度包括曼哈頓距離佩谣、蘭氏距離、非對稱二元距離实蓬、最大距離和閔可夫斯基距離（可使用 ?dist 查看詳細信息）茸俭。在這一 章中，計算距離時默認使用歐幾里得距離安皱。計算歐幾里得距離的方法見16.1.1節(jié)调鬓。
(5) 選擇聚類算法
接下來選擇對數(shù)據(jù)聚類的方法，層次聚類對于小樣本來說很實用（如150 個觀測值或更少）酌伊，而且這種情況下嵌套聚類更實用腾窝。劃分的方法能處理更大的數(shù)據(jù)量，但是需要事先確定聚類的個數(shù)居砖。一旦選定了層次方法或劃分方法虹脯，就必須選擇一個特定的聚類算法。這里再次強調每個算法都有優(yōu)點和缺點奏候。16.2節(jié)和16.3節(jié)中介紹了最常用的幾種方法循集。你可以嘗試多種算法來看看相應結果的穩(wěn)健性。
(6) 獲得一種或多種聚類方法
這一步可以使用步驟(5)選擇的方法蔗草。
(7) 確定類的數(shù)目
為了得到最終的聚類方案咒彤，你必須確定類的數(shù)目。對此研究者們也提出了很多相應的解決方法咒精。常用方法是嘗試不同的類數(shù)（比如2～K）并比較解的質量蔼紧。在 NbClust包中的NbClust() 函數(shù)提供了30個不同的指標來幫助你進行選擇（也表明這個問題有多么難解）。本章將多次使用這個包狠轻。
(8) 獲得最終的聚類解決方案
一旦類的個數(shù)確定下來奸例，就可以提取出子群，形成最終的聚類方案了。
(9) 結果可視化
可視化可以幫助你判定聚類方案的意義和用處查吊。層次聚類的結果通常表示 為一個樹狀圖谐区。劃分的結果通常利用可視化雙變量聚類圖來表示。
(10) 解讀類
一旦聚類方案確定逻卖，你必須解釋（或許命名）這個類宋列。一個類中的觀測值有何相似之處？不同的類之間的觀測值有何不同评也？這一步通常通過獲得類中每個變量的匯總統(tǒng)計來完成炼杖。對于連續(xù)數(shù)據(jù)，每一類中變量的均值和中位數(shù)會被計算出來盗迟。對于混合數(shù)據(jù)（數(shù)據(jù)中包含 分類變量）坤邪，結果中將返回各類的眾數(shù)或類別分布。
(11) 驗證結果
驗證聚類方案相當于問：“這種劃分并不是因為數(shù)據(jù)集或聚類方法的某種特性罚缕，而是確實給出了一個某種程度上有實際意義的結果嗎艇纺？”如果采用不同的聚類方法或不同的樣本，是否會產生相同的類邮弹？ fpc 黔衡、clv 和 clValid 包包含了評估聚類解的穩(wěn)定性的函數(shù)。 因為觀測值之間距離的計算是聚類分析的一部分腌乡，所以我們將在下一節(jié)詳細討論盟劫。
2 計算距離
聚類分析的第一步都是度量樣本單元間的距離、相異性或相似性与纽。兩個觀測值之間的歐幾里得距離定義為：
這里i和j代表第i和第j個觀測值侣签，p是變量的個數(shù)。
查看在 flexclust 包中的營養(yǎng)數(shù)據(jù)集渣锦，它包括對27種肉、魚和禽的營養(yǎng)物質的測量氢哮。最初的幾個觀測值由下面的代碼給出：
data(nutrient, package="flexclust")
head(nutrient, 4)
energy protein fat calcium iron
BEEF BRAISED 340 20 28 9 2.6
HAMBURGER 245 21 17 9 2.7
BEEF ROAST 420 15 39 7 2.0
BEEF STEAK 375 19 32 9 2.6
前兩個觀測值（ BEEF BRAISED 和 HAMBURGER ）之間的歐幾里得距離為：
R軟件中自帶的 dist() 函數(shù)能用來計算矩陣或數(shù)據(jù)框中所有行（觀測值）之間的距離袋毙。格式是 dist(x, method=)，這里的 x 表示輸入數(shù)據(jù)冗尤，并且默認為歐幾里得距離听盖。函數(shù)默認返回一個下三角矩陣，但是 as.matrix()函數(shù)可使用標準括號符號得到距離裂七。
對于營養(yǎng)數(shù)據(jù)集的數(shù)據(jù)框來說皆看，前四行的距離為：
d <- dist(nutrient)
as.matrix(d)[1:4,1:4]
BEEF BRAISED HAMBURGER BEEF ROAST BEEF STEAK
BEEF BRAISED 0.00000 95.6400 80.93429 35.24202
HAMBURGER 95.64000 0.0000 176.49218 130.87784
BEEF ROAST 80.93429 176.4922 0.00000 45.76418
BEEF STEAK 35.24202 130.8778 45.76418 0.00000
觀測值之間的距離越大，異質性越大背零。觀測值和它自己之間的距離是0腰吟。
混合數(shù)據(jù)類型的聚類分析
歐幾里得距離通常作為連續(xù)型數(shù)據(jù)的距離度量。但是如果存在其他類型的數(shù)據(jù)，則需要相異的替代措施毛雇，你可以使用 cluster 包中的 daisy() 函數(shù)來獲得包含任意二元（binary）嫉称、名義（nominal）、有序（ordinal）灵疮、連續(xù)（continuous）屬性組合的相異矩陣织阅。 cluster 包中的其他函數(shù)可以使用這些異質性來進行聚類分析。例如 agnes() 函數(shù)提供了層次聚類震捣， pam() 函數(shù)提供了圍繞中心點的劃分的方法荔棉。
需要注意的是，在營養(yǎng)數(shù)據(jù)集中蒿赢，距離很大程度上由能量（ energy ）這個變量控制润樱，這是因為該變量變化范圍更大∷咧玻縮放數(shù)據(jù)有利于均衡各變量的影響祥国。在下一節(jié)中，你可以對該數(shù)據(jù)集應用層次聚類分析晾腔。
3 層次聚類分析
如前所述舌稀，在層次聚類中，起初每一個實例或觀測值屬于一類灼擂。聚類就是每一次把兩類聚成新的一類壁查，直到所有的類聚成單個類為止，算法如下：
(1) 定義每個觀測值（行或單元）為一類剔应；
(2) 計算每類和其他各類的距離睡腿；
(3) 把距離最短的兩類合并成一類，這樣類的個數(shù)就減少一個峻贮；
(4) 重復步驟(2)和步驟(3)席怪，直到包含所有觀測值的類合并成單個的類為止。
在層次聚類算法中纤控，主要的區(qū)別是它們對類的定義不同（步驟(2)）挂捻。表16-1給出了五種最常見聚類方法的定義和其中兩類之間距離的定義。
單聯(lián)動聚類方法傾向于發(fā)現(xiàn)細長的船万、雪茄型的類刻撒。它也通常展示一種鏈式的現(xiàn)象，即不相似的觀測值分到一類中耿导，因為它們和它們的中間值很相像声怔。全聯(lián)動聚類傾向于發(fā)現(xiàn)大致相等的直徑緊湊類。它對異常值很敏感舱呻。平均聯(lián)動提供了以上兩種方法的折中醋火。相對來說，它不像鏈式，而 且對異常值沒有那么敏感胎撇。它傾向于把方差小的類聚合介粘。
Ward法傾向于把有少量觀測值的類聚合到一起，并且傾向于產生與觀測值個數(shù)大致相等的類晚树。它對異常值也是敏感的姻采。質心法是一種很受歡迎的方法，因為其中類距離的定義比較簡單爵憎、 易于理解慨亲。相比其他方法，它對異常值不是很敏感宝鼓。但是它可能不如平均聯(lián)動法或Ward方法表現(xiàn) 得好刑棵。
層次聚類方法可以用 hclust()函數(shù)來實現(xiàn)，格式是hclust(d, method=)愚铡，其中 d是通過 dist() 函數(shù)產生的距離矩陣蛉签，并且方法包括"single" 、 "complete" 沥寥、 "average" 碍舍、 "centroid"和 "ward"。 在本節(jié)中邑雅，你可以使用平均聯(lián)動聚類方法處理16.1.1節(jié)介紹的營養(yǎng)數(shù)據(jù)片橡。我們的目的是基于27種食物的營養(yǎng)信息辨別其相似性、相異性并分組淮野。下面的代碼清單提供了實施聚類的代碼捧书。
營養(yǎng)數(shù)據(jù)的平均聯(lián)動聚類
data(nutrient, package="flexclust")
row.names(nutrient) <- tolower(row.names(nutrient))
nutrient.scaled <- scale(nutrient)
d <- dist(nutrient.scaled) fit.average <- hclust(d, method="average") plot(fit.average, hang=-1, cex=.8, main="Average Linkage Clustering")
首先載入數(shù)據(jù)，同時將行名改為小寫（因為我討厭大寫的標簽）骤星。由于變量值的變化范圍很大经瓷，我們將其標準化為均值為0、方差為1洞难。27種食物之間的距離采用歐幾里得距離舆吮，應用的方法是平均聯(lián)動。最后廊营，結果用樹狀圖來展示（參見圖16-1）歪泳。 plot() 函數(shù)中的 hang 命令展示觀測值的標簽（讓它們在掛在0下面）萝勤。
樹狀圖應該從下往上讀露筒，它展示了這些條目如何被結合成類。每個觀測值起初自成一類敌卓，然 后相距最近的兩類（beef braised和smoked ham）合并慎式。其次，pork roast和pork simmered合并， chicken canned和tuna canned合并瘪吏。再次癣防，beef braised/smoked ham這一類和pork roast/pork simmered 這一類合并（這個類目前包含四種食品）。合并繼續(xù)進行下去掌眠，直到所有的觀測值合并成一類蕾盯。高度刻度代表了該高度類之間合并的判定值。對于平均聯(lián)動來說蓝丙，標準是一類中的點和其他類中 的點的距離平均值级遭。
如果你的目的是理解基于食物營養(yǎng)成分的相似性和相異性，圖16-1就足夠了渺尘。它提供了27種 食物之間的相似性/異質性的層次分析視圖挫鸽。tuna canned和chicken canned是相似的，但是都和clams canned有很大的不同鸥跟。但是丢郊，如果最終目標是這些食品分配到的類（希望有意義的）較少，我們 需要額外的分析來選擇聚類的適當個數(shù)医咨。
NbClust包提供了眾多的指數(shù)來確定在一個聚類分析里類的最佳數(shù)目枫匾。不能保證這些指標得出的結果都一致。事實上腋逆，它們可能不一樣婿牍。但是結果可用來作為選擇聚類個數(shù)K值的一個參考。NbClust()函數(shù)的輸入包括需要做聚類的矩陣或是數(shù)據(jù)框惩歉，使用的距離測度和聚類方法等脂，并考慮 最小和最大聚類的個數(shù)來進行聚類。它返回每一個聚類指數(shù)撑蚌，同時輸出建議聚類的最佳數(shù)目上遥。
下面的代碼清單使用該方法處理營養(yǎng)數(shù)據(jù)的平均聯(lián)動聚類。
nc <- NbClust(nutrient.scaled, distance="euclidean",
min.nc=2, max.nc=15, method="average")
table(nc$Best.n[1,])
0 1 2 3 4 5 9 10 13 14 15
2 1 4 4 2 4 1 1 2 1 4
barplot(table(nc$Best.n[1,]),
xlab="Numer of Clusters", ylab="Number of Criteria",
main="Number of Clusters Chosen by 26 Criteria")
這里争涌，四個評判準則贊同聚類個數(shù)為2粉楚，四個判定準則贊同聚類個數(shù)為3，等等亮垫。結果在圖16-2中模软。
你可以試著用“投票”個數(shù)最多的聚類個數(shù)（2、3饮潦、5和15）并選擇其中一個使得解釋最有意義燃异。下面的代碼清單展示了五類聚類的方案。
clusters <- cutree(fit.average, k=5)
table(clusters)
aggregate(nutrient, by=list(cluster=clusters), median)
aggregate(as.data.frame(nutrient.scaled), by=list(cluster=clusters),median)
pdf(file = "test.pdf")
plot(fit.average, hang=-1, cex=.8,
main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=5)
dev.off()
cutree() 函數(shù)用來把樹狀圖分成五類?继蜡。第一類有7個觀測值回俐，第二類有16個觀測值逛腿，等等。aggregate() 函數(shù)用來獲取每類的中位數(shù)?仅颇，結果有原始度量和標準度量兩種形式单默。最后，樹狀圖被重新繪制忘瓦， rect.hclust() 函數(shù)用來疊加五類的解決方案?搁廓。結果展示在圖16-3中。
4 劃分聚類分析
在劃分方法中耕皮，觀測值被分為K組并根據(jù)給定的規(guī)則改組成最有粘性的類枚抵。這一節(jié)討論兩種。方法：K均值和基于中心點的劃分（PAM）明场。
4.1 K 均值聚類
最常見的劃分方法是K均值聚類分析汽摹。從概念上講，K均值算法如下：
(1) 選擇K個中心點（隨機選擇K行）苦锨；
(2) 把每個數(shù)據(jù)點分配到離它最近的中心點逼泣；
(3) 重新計算每類中的點到該類中心點距離的平均值（也就說，得到長度為p的均值向量舟舒，這里的p是變量的個數(shù)）拉庶；
(4) 分配每個數(shù)據(jù)到它最近的中心點；
(5) 重復步驟(3)和步驟(4)直到所有的觀測值不再被分配或是達到最大的迭代次數(shù)（R把10次 作為默認迭代次數(shù)）秃励。
這種方法的實施細節(jié)可以變化氏仗。R軟件使用Hartigan & Wong（1979）提出的有效算法，這種 算法是把觀測值分成k組并使得觀測值到其指定的聚類中心的平方的總和為最小夺鲜。也就是說皆尔，在 步驟(2)和步驟(4)中，每個觀測值被分配到使下式得到最小值的那一類中：
表示第i個觀測值中第j個變量的值币励。
表示第
個類中第
個變量的均值慷蠕，其中
是變量的個數(shù)。
K均值聚類能處理比層次聚類更大的數(shù)據(jù)集食呻。另外流炕，觀測值不會永遠被分到一類中。當我們 提高整體解決方案時仅胞，聚類方案也會改動每辟。但是均值的使用意味著所有的變量必須是連續(xù)的，并 且這個方法很有可能被異常值影響干旧。它在非凸聚類（例如U型）情況下也會變得很差渠欺。
在R中K均值的函數(shù)格式是 kmeans(x,centers) ，這里 x 表示數(shù)值數(shù)據(jù)集（矩陣或數(shù)據(jù)框）莱革， centers 是要提取的聚類數(shù)目峻堰。函數(shù)返回類的成員、類中心盅视、平方和（類內平方和悬赏、類間平方和堵泽、 總平方和）和類大小。
由于K均值聚類在開始要隨機選擇k個中心點，在每次調用函數(shù)時可能獲得不同的方案诀诊。使用 set.seed() 函數(shù)可以保證結果是可復制的。此外义起，聚類方法對初始中心值的選擇也很敏感蔽挠。 kmeans() 函數(shù)有一個 nstart 選項嘗試多種初始配置并輸出最好的一個。例如断箫，加上 nstart=25 會生成25個初始配置拂酣。通常推薦使用這種方法。
不像層次聚類方法仲义，K均值聚類要求你事先確定要提取的聚類個數(shù)婶熬。同樣，NbClust包可以用來作為參考埃撵。另外赵颅，在K均值聚類中，類中總的平方值對聚類數(shù)量的曲線可能是有幫助的暂刘〗让可根據(jù)圖中的彎曲（類似于14.2.1節(jié)描述的卵石試驗彎曲）選擇適當?shù)念惖臄?shù)量。
圖像可以用下面的代碼產生：
wssplot <- function(data, nc=15, seed=1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))
for (i in 2:nc){
set.seed(seed)
wss[i] <- sum(kmeans(data, centers=i)$withinss)}
plot(1:nc, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares")
}
data 參數(shù)是用來分析的數(shù)值數(shù)據(jù)谣拣， nc 是要考慮的最大聚類個數(shù)募寨，而 seed 是一個隨機數(shù)種子。
讓我們用K均值聚類來處理包含178種意大利葡萄酒中13種化學成分的數(shù)據(jù)集森缠。該數(shù)據(jù)最初來自于UCI機器學習庫（http://www.ics.uci.edu/~mlearn/MLRepository.html）绪商，但是可以通過 rattle包獲得。在這個數(shù)據(jù)集里辅鲸，觀測值代表三種葡萄酒的品種格郁，由第一個變量（類型）表示。你可以放棄這一變量独悴，進行聚類分析例书，看看是否可以恢復已知的結構。
> data(wine, package="rattle")
> head(wine)
Type Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color Hue Dilution Proline
1 1 14.23 1.71 2.43 15.6 127 2.80 3.06 0.28 2.29 5.64 1.04 3.92 1065
2 1 13.20 1.78 2.14 11.2 100 2.65 2.76 0.26 1.28 4.38 1.05 3.40 1050
3 1 13.16 2.36 2.67 18.6 101 2.80 3.24 0.30 2.81 5.68 1.03 3.17 1185
4 1 14.37 1.95 2.50 16.8 113 3.85 3.49 0.24 2.18 7.80 0.86 3.45 1480
5 1 13.24 2.59 2.87 21.0 118 2.80 2.69 0.39 1.82 4.32 1.04 2.93 735
6 1 14.20 1.76 2.45 15.2 112 3.27 3.39 0.34 1.97 6.75 1.05 2.85 1450
標準化數(shù)據(jù)：
> df <- scale(wine[-1])
> df[1:4,1:4]
Alcohol Malic Ash Alcalinity
[1,] 1.5143408 -0.56066822 0.2313998 -1.1663032
[2,] 0.2455968 -0.49800856 -0.8256672 -2.4838405
[3,] 0.1963252 0.02117152 1.1062139 -0.2679823
[4,] 1.6867914 -0.34583508 0.4865539 -0.8069748
決定聚類的個數(shù)：
> wssplot(df)
Hit <Return> to see next plot: library(NbClust)
> set.seed(1234)
> devAskNewPage(ask=TRUE)
> nc <- NbClust(df, min.nc=2, max.nc=15, method="kmeans")
Hit <Return> to see next plot: table(nc$Best.n[1,])
*** : 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.
Hit <Return> to see next plot: barplot(table(nc$Best.n[1,]),xlab="Number of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria")
*** : 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
*******************************************************************
進行K均值聚類分析：
> set.seed(1234)
> fit.km <- kmeans(df, 3, nstart=25)
> fit.km$size
[1] 62 65 51
> fit.km$centers
Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color
1 0.8328826 -0.3029551 0.3636801 -0.6084749 0.57596208 0.88274724 0.97506900 -0.56050853 0.57865427 0.1705823
2 -0.9234669 -0.3929331 -0.4931257 0.1701220 -0.49032869 -0.07576891 0.02075402 -0.03343924 0.05810161 -0.8993770
3 0.1644436 0.8690954 0.1863726 0.5228924 -0.07526047 -0.97657548 -1.21182921 0.72402116 -0.77751312 0.9388902
Hue Dilution Proline
1 0.4726504 0.7770551 1.1220202
2 0.4605046 0.2700025 -0.7517257
3 -1.1615122 -1.2887761 -0.4059428
> aggregate(wine[-1], by=list(cluster=fit.km$cluster), mean)
cluster Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids Nonflavanoids Proanthocyanins Color Hue
1 1 13.67677 1.997903 2.466290 17.46290 107.96774 2.847581 3.0032258 0.2920968 1.922097 5.453548 1.0654839
2 2 12.25092 1.897385 2.231231 20.06308 92.73846 2.247692 2.0500000 0.3576923 1.624154 2.973077 1.0627077
3 3 13.13412 3.307255 2.417647 21.24118 98.66667 1.683922 0.8188235 0.4519608 1.145882 7.234706 0.6919608
Dilution Proline
1 3.163387 1100.2258
2 2.803385 510.1692
3 1.696667 619.0588
因為變量值變化很大刻炒，所以在聚類前要將其標準化?决采。下一步，使用wssplot()和 Nbclust()函數(shù)確定聚類的個數(shù)?坟奥。圖16-4表示從一類到三類變化時树瞭，組內的平方總和有一個明顯的下降趨勢拇厢。三類之后，下降的速度減弱晒喷，暗示著聚成三類可能對數(shù)據(jù)來說是一個很好的擬合孝偎。在圖16-5中， NbClust 包中的24種指標中有14種建議使用類別數(shù)為三的聚類方案凉敲。需要注意的是并非30個指標都可以計算每個數(shù)據(jù)集衣盾。
使用 kmeans()函數(shù)得到的最終聚類中势决，聚類中心也被輸出了?。因為輸出的聚類中心是基于標準化的數(shù)據(jù)蓝撇，所以可以使用 aggregate() 函數(shù)和類的成員來得到原始矩陣中每一類的變量均值果复。
K均值可以很好地揭示類型變量中真正的數(shù)據(jù)結構嗎？交叉列表類型（葡萄酒品種）和類成員由下面的代碼表示：
> ct.km <- table(wine$Type, fit.km$cluster)
> ct.km
1 2 3
1 59 0 0
2 3 65 3
3 0 0 48
你可以用 flexclust 包中的蘭德指數(shù)（Rand index）來量化類型變量和類之間的協(xié)議：
> library(flexclust)
> randIndex(ct.km)
[1] 0.897
調整的蘭德指數(shù)為兩種劃分提供了一種衡量兩個分區(qū)之間的協(xié)定渤昌，即調整后機會的量度据悔。它的變化范圍是從–1（不同意）到1 （完全同意）。葡萄酒品種類型和類的解決方案之間的協(xié)定是0.9耘沼。結果不壞极颓，那我們應該喝點酒？
4.2 圍繞中心點的劃分
因為K均值聚類方法是基于均值的群嗤，所以它對異常值是敏感的菠隆。一個更穩(wěn)健的方法是圍繞中心點的劃分（PAM）。與其用質心（變量均值向量）表示類狂秘，不如用一個最有代表性的觀測值來 表示（稱為中心點）骇径。K均值聚類一般使用歐幾里得距離，而PAM可以使用任意的距離來計算者春。 因此破衔，PAM可以容納混合數(shù)據(jù)類型，并且不僅限于連續(xù)變量钱烟。
PAM算法如下：
(1) 隨機選擇K個觀測值（每個都稱為中心點）晰筛；
(2) 計算觀測值到各個中心的距離/相異性；
(3) 把每個觀測值分配到最近的中心點拴袭；
(4) 計算每個中心點到每個觀測值的距離的總和（總成本）读第；
(5) 選擇一個該類中不是中心的點，并和中心點互換拥刻； (6) 重新把每個點分配到距它最近的中心點怜瞒；
(7) 再次計算總成本；
(8) 如果總成本比步驟(4)計算的總成本少般哼，把新的點作為中心點吴汪；
(9) 重復步驟(5)～(8)直到中心點不再改變惠窄。
在PAM方法中應用基礎數(shù)學的一個例子可以在這里找到：http://en.wikipedia.org/wiki/k- medoids（我不經常引用維基百科，但這確實是一個很好的例子）漾橙。
你可以使用 cluster 包中的 pam() 函數(shù)使用基于中心點的劃分方法杆融。格式是pam(x, k, metric="euclidean", stand=FALSE)，這里的 x 表示數(shù)據(jù)矩陣或數(shù)據(jù)框近刘，k 表示聚類的個數(shù)， metric 表示使用的相似性/相異性的度量臀晃，而 stand 是一個邏輯值觉渴，表示是否有變量應該在計算 該指標之前被標準化。圖16-6中列出了使用PAM方法處理葡萄酒的數(shù)據(jù)徽惋。
pdf(file = "test.pdf")
set.seed(1234)
fit.pam <- pam(wine[-1], k=3, stand=TRUE)
fit.pam$medoids
clusplot(fit.pam, main="Bivariate Cluster Plot")
dev.off()
注意案淋，這里得到的中心點是葡萄酒數(shù)據(jù)集中實際的觀測值。在這種情況下险绘，分別選擇36踢京、107和175個觀測值來代表三類。通過從13個測定變量上得到的前兩個主成分繪制每個觀測的坐標來創(chuàng)建二元圖（參見第14章）宦棺。每個類用包含其所有點的最小面積的橢圓表示瓣距。還需要注意的是，PAM在如下例子中的表現(xiàn)不如K均值：
> ct.pam <- table(wine$Type, fit.pam$clustering)
> 1 2 3
> 1 59 0 0
> 2 16 53 2
> 3 0 1 47
> randIndex(ct.pam)
> [1] 0.699
調整的蘭德指數(shù)從（K均值的）0.9下降到了0.7
5 避免不存在的類
在結束討論前代咸，還要提出一點注意事項蹈丸。聚類分析是一種旨在識別數(shù)據(jù)集子組的方法，并且在此方面十分擅長呐芥。事實上逻杖，它甚至能發(fā)現(xiàn)不存在的類。
請看以下的代碼:
library(fMultivar)
set.seed(1234)
df <- rnorm2d(1000, rho=.5)
df <- as.data.frame(df)
plot(df, main="Bivariate Normal Distribution with rho=0.5")
隨后思瘟，使用 wssplot()和Nbclust()函數(shù)來確定當前聚類的個數(shù)：
wssplot <- function(data, nc=15, seed=1234){
wss <- (nrow(data)-1)*sum(apply(data,2,var))
for (i in 2:nc){
set.seed(seed)
wss[i] <- sum(kmeans(data, centers=i)$withinss)}
plot(1:nc, wss, type="b", xlab="Number of Clusters",
ylab="Within groups sum of squares")
}
pdf(file = "test1.pdf")
wssplot(df)
wssplot() 函數(shù)建議聚類的個數(shù)是3，然而 NbClust 函數(shù)返回的準則多數(shù)支持2類或3類滨攻。如
果利用PAM法進行雙聚類分析：
pdf(file = "test.pdf")
library(ggplot2)
library(cluster)
fit <- pam(df, k=2)
df$clustering <- factor(fit$clustering)
p <- ggplot(data=df, aes(x=V1, y=V2, color=clustering, shape=clustering))
p <- p + geom_point() + ggtitle("Clustering of Bivariate Normal Data")
print(p)
dev.off()
你將得到如圖16-10所示的雙聚類圖像（ ggplot() 函數(shù)是綜合圖像 ggplot2 包的一部分够话，第19章包含 ggplot2 的詳細內容）。
很明顯劃分是人為的。實際上在這里沒有真實的類奇钞。那么你怎樣避免這種錯誤呢澡为？雖然并非萬無一失，但是我發(fā)現(xiàn) NbClust 包中的立方聚類規(guī)則（Cubic Cluster Criteria景埃，CCC）往往可以幫助我們揭示不存在的結構媒至。代碼是：
plot(nc$All.index[,4], type="o", ylab="CCC",
xlab="Number of clusters", col="blue")
結果展示在圖16-11中顶别。當CCC的值為負并且對于兩類或是更多的類遞減時，就是典型的單峰分布拒啰。
聚類分析（或你對它的解讀）找到錯誤聚類的能力使得聚類分析的驗證步驟很重要驯绎。如果你試圖找出在某種意義上“真實的”類（而不是一個方便的劃分），就要確保結果是穩(wěn)健的并且是可重復的谋旦。你可以嘗試不同的聚類方法剩失，并用新的樣本復制結果。如果同一類持續(xù)復原册着，你就可以對得出的結果更加自信拴孤。
6 小結
在這一章里，我們回顧了把觀測值凝聚成子組的常見聚類方法甲捏。首先演熟，我們回顧了常見聚類分析的一般步驟。接下來司顿，描述了層次聚類和劃分聚類的常見方法芒粹。
最后，如果尋求的不只是更方便的劃分方法大溜，我強調了需要驗證所產生的類化漆。聚類分析是一個寬泛的話題，而R有一些最全面的方案來實施現(xiàn)有的方法钦奋。想要了解更多获三，可以參考CRAN任務視圖的聚類分析和有限混合模型部分（http://cran.r-project.org/web/views/Cluster.html）。
除此之外锨苏，Tan疙教、Steinbach & Kumar（2006）寫了一本關于數(shù)據(jù)挖掘技術的好書。它有一章清晰地講解了聚類分析伞租，可以自由下載（www-users.cs.umn.edu/~kumar/dmbook/ch8.pdf）贞谓。最后，Everitt葵诈、Landau裸弦、Leese & Stahl（2011）就這個問題寫了一本實用性的課本，評價很高作喘。聚類分析方法用于發(fā)現(xiàn)數(shù)據(jù)集中有凝聚力的觀測值分組理疙。在下一章中，我們將考慮已經定義好分組的情況泞坦，你的目標是用一個準確的方法把觀測值分入這些組窖贤。
