應用統計學與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的交互式小地圖（老規矩）。