自己整理編寫的R語言常用資料分析模型的模板，原檔案為Rmd格式，直接複製貼上過來，作為個人學習筆記儲存和分享。部分參考薛毅的《統計建模與R軟體》和《R語言實戰》

聚類分析是一類將資料所研究物件進行分類的統計方法，這一類方法的共同特點是：事先不知道類別的個數和結構，據以進行分析的資料是物件之間的相似性或相異性的資料。將這些相似（相異）性資料看成是物件之間的“距離”遠近的一種度量，將距離近的變數歸為一類，不同類之間的物件距離較遠。這就是聚類分析方法的共同思路。

setwd("C://Users//DELL//Documents//WD")

需要安裝的包

install.packages("NbClust"
 
)
install.packages("flexclust")

---

1. 資料的預處理

1.1 R裡的距離計算公式

在度量空間中，可自行定義距離公式。常用的定量變數間的距離公式有：絕對值距離、Euclide距離、Minkowski距離、Chebyshev距離、Mahalanobis距離、Lance和Williams距離。除此之外還有對定型變數的距離公式。

在R軟體中，dist()函式給出了各種距離的計算結果，其使用格式為：

  dist(x, method = " ",diag = FALSE, upper = FALSE, p = 2)

其中x是樣本構成的資料矩陣（樣本按行輸入）或資料框，method表示計算距離選用的方法，預設值為Euclide距離，其他可選用的方法有：

• “euclidean” – Euclide距離
• “maximum” – Chebyshev距離
• “manhattan” – 絕對值距離
• “canberra” – Lance距離
• “minkowski” – Minkowski距離，其中引數p是Minkowski距離的階數，即公式中的q
• “binary” – 定型變數的距離
  （以上距離的數學公式可參見薛毅書P467）

diag是邏輯變數，當diag = TRUE時,給出對角線上的距離；upper是邏輯變數,當upper = TRUE時，給出上三角矩陣的值（預設值僅給出下三角矩陣的值）

1.2 資料中心化與標準化變換

在作聚類分析過程中，大多數資料往往是不能直接參與運算的，需要先將資料作中心化或標準化處理。

* 1.2.1 中心化變換*

Xij* = Xij - Xj均， i=1,2,…,p 為中心化變換，
變換後資料的均值為0，方差陣不變

* 1.2.2 標準化變換*

Xij* = (Xij - Xj均)/Sj, i=1,2,…,p 為中心化變換，
變換後資料，每個變數的樣本均值均為0，標準差為1，而且標準化後的資料與變數的量綱無關

在R軟體中，用scale()函式作資料的中心化或標準化，使用格式為：

 scale(x, center = TRUE, scale = TRUE)

其中x是樣本構成的資料矩陣，center是邏輯變數，TRUE（預設值）表示對資料作中心化變換，FALSE表示不作變換；scale是邏輯變數，TRUE（預設值）表示表示對資料作標準化變換，FALSE表示不作變換。
1.2.1中公式的計算函式為：X* = scale(x,scale = FALSE);
1.2.2中公式的計算函式為：X* = scale(x)

* 1.2.3 極差標準化變換*

Xij* = (Xij - Xj均)/Rj, i=1,2,…,n，j=1,2,…,p 為極差標準化變換，
變換後的資料，每個變數的樣本均值都為0，極差為1，且|Xij*|<1，在以後的分析計算中可以減少誤差的產生，同時變換後的資料也是無量綱的量

在R軟體中，可用sweep()函式作極差標準化變換，變換過程如下：

 center <- sweep( x, 2, apply(x, 2, mean) )
 R <- apply(x, 2, max) - apply(x, 2, min)
 x_star <- sweep(center, 2, R, "/")

其中x是樣本構成的資料矩陣，第一行是將資料中心化，第二行是計算極差Rj，j=1,2,…,p，第三行是將中心化後的資料除以極差，得到資料的極差標準化資料。

其中用到的sweep()函式是對陣列或矩陣運算，格式為：

 sweep(x, MARGIN, STATS, FUN = "-" , ...)

其中x是陣列或矩陣，MARGIN是運算的區域，對於矩陣來講，1表示行，2表示列。STATS是統計量，如apply(x,2,mean)表示各列的均值。FUN表示函式的運算，預設值為減法運算。

從sweep()函式的用法可知，將上述命令第三行改為x_star <- sweep(center, 2, sd(x), “/”)，得到的就是(普通)標準化變換後的結果

* 1.2.4 極差正則化變換*

Xij* = [Xij - MIN 1<=k<=n (Xkj)]/Rj, i=1,2,…,n，j=1,2,…,p 為極差正則化變換
變換後資料 0 <= Xij* <= 1，極差為1，也是無量綱的量

利用sweep()函式，可以得到資料的極差正則化變換，其變換過程如下：

 center <- sweep( x, 2, apply(x, 2, min) )
 R <- apply(x, 2, max) - apply(x, 2, min)
 x_star <- sweep(center, 2, R, "/")

1.3 相似係數

聚類分析不僅用來對樣本進行分類，也可以對變數來進行分類。在對變數進行分類時，常用相似係數來度量變數之間的相似程度。

* 1.3.1 夾角餘弦*

計算公式與計算向量餘弦的演算法一樣，變數Xi的n次觀測值作為一個向量，Xj的n個觀測值做另一個向量，這裡不作展示

在R軟體中，可用scale()函式完成兩向量餘弦的計算，公式如下：

 y <- scale(x, center = F, scale = T)/sqrt(nrow(x)-1)
 C <-t(y) %*% y

其中x是樣本構成的資料矩陣，C是計算得出的相似係數構成的矩陣

* 1.3.2 相關係數*

相關係數就是對資料作標準化處理後的夾角餘弦，也就是變數Xi和變數Xj的相關係數rij，這裡記作cij，公式這裡不作展示

在R軟體中，cij的計算更加方便，即樣本的相關矩陣：C <- cor(x)

變數之間常藉助於相似係數來定義距離，如令 dij^2 = 1 - cij^2，有時也用相似係數來度量變數間的相似程度。 （上述兩公式見薛毅書P473）

---

2.層次聚類法(也叫系統聚類法)

既可處理分類變數，也可處理連續變數，但不能同時處理兩種變數型別，不需要指定類別數。聚類結果間存在著巢狀，或者說層次的關係。系統聚類方法是聚類分析諸方法中用得最多的一種，其基本思想是：開始將n個樣本各自作為一類，並規定樣本之間的距離和類和類之間的距離，然後將距離最近的兩類合併為一個新類，計算新類與其他類的距離；重複進行兩個類的合併，每次減少一類，直至所有的樣本合併為一類。可用於定義“距離”的統計量包括了歐氏距離(euclidean)、馬氏距離(manhattan)、兩項距離(binary)、明氏距離(minkowski)。還包括相關係數和夾角餘弦。在計算類間距離時則有六種不同的方法，分別是最短距離法、最長距離法、類平均法、重心法、中間距離法、離差平方和法。

在R軟體中，hclust()函式提供了系統聚類的計算，plot()函式可畫出系統聚類的樹形圖（或稱為譜系圖)

hclust()函式使用格式為：

 hclust(d, method = "complete", members = NULL)

其中d是由“dist”構成的結構。method是系統聚類方法（預設值是最長距離法），可選用的其他引數有：

• “single” – 最短距離法
• “complete” – 最長距離法
• “median” – 中間距離法
• “mcquitty” – Mcquitty相似法
• “average” – 類平均法
• “centroid” – 重心法
• “ward“ – 離差平方和法

members的預設值為NULL，或與d有相同變數長度的向量，具體使用方法見幫助

plot()函式畫出譜系圖的格式為：

 plot(x, labels = NULL, hang = 0.1, axes = TRUE, frame.plot = FALSE, ann = TRUE, 
 main = "Cluster Dendrogram", sub = NULL, xlab = NULL, ylab = "Height",...)

其中x是由hclust()函式生成的物件，hang是表明譜系圖中各類所在的位置，當hang取負值時，譜系圖中的類從底部畫起。其他引數見plot幫助檔案

下面用例項來說明系統聚類方法

2.1 例1. 城鎮居民消費水平（對變數進行分類） （cluster1）

城鎮居民消費水平通常用8項指標來描述：人均糧食之處（元/人）x1，人均副食支出（元/人）x2，人均煙、酒、茶支出（元/人）x3，人均其他副食支出（元/人）x4，人均衣著商品支出（元/人）x5，人均日用品支出（元/人）x6，人均燃料支出（元/人）x7，人均非商品支出（元/人）x8。這8項指標間存在著一定的相關性。為了研究城鎮居民的消費結構，需將相關性強的指標歸併到一起，這實際上是對指標聚類。

* 2.1.1 繪製譜系圖*

#用sas中相同的資料來源，輸入資料，列名稱是x1-x8，行名稱是1-30的序號
x <- read.csv("cluster1.csv")
#相關係數矩陣
r <- cor(x)
#生成距離結構
d<-as.dist(1-r)  #method預設值預設是Euclide距離
#用四種不同的距離方法生成系統聚類
hc1<-hclust(d, "single"); hc2<-hclust(d, "complete") #最短距離法和最長距離法
hc3<-hclust(d, "median"); hc4<-hclust(d, "mcquitty") #中間距離法和Mcquitty相似法
#繪製出以上四種方法的樹形結構圖，並以2x2的形式畫在一張圖上.影象距離邊界的距離
opar <- par(mfrow = c(2, 2))
#hang是表明譜系圖中各類所在的位置 當hang取負值時，譜系圖中的類從底部畫起  生成譜系圖
plot(hc1,hang=-1); plot(hc2,hang=-1)
plot(hc3,hang=-1); plot(hc4,hang=-1)
par(opar)

R中與繪製譜系圖有關的函式還有as.dendrogram(object, hang = -1, …)

其中object是是由hclust得到的物件，此時plot()函式用法為：

 plot(x, type = c("rectangle","triangle"), center = FALSE, 
 edge.root = is.leaf(x)||!is.null(attr(x,"edgetest"), 
 nodePar = NULL, edgePar = list(),
 leaflab = c("perpendicular","textlike","none"),
 dLeaf = NULL, xlab = "", ylab = "", xaxt = "n", yaxt = "s",
 horiz = FALSE, frame.plot = FALSE,...)

其中x是由dendrogram得到的物件，type表示畫譜系圖的型別，“rectangle”表示矩形（預設值），“triangle”為三角形；horiz是邏輯變數，當horiz = TRUE時，表示譜系圖水平放置。其他引數見幫助檔案

* 2.1.2 設定譜系圖形狀*

#以hc1為例繪製不同引數下的譜系圖
dend1<-as.dendrogram(hc1)
#生成兩行一列  影象距離邊界的距離
opar <- par(mfrow = c(2, 2),mar = c(4,3,1,2))
#通過設定引數分別畫出四種不同的圖
plot(dend1)  #原始影象
plot(dend1, nodePar=list(pch = c(1,NA), cex=0.8, lab.cex = 0.8),
     type = "t", center=TRUE)  #斜線
plot(dend1, edgePar=list(col = 1:2, lty = 2:3), dLeaf=1, edge.root = TRUE)  #虛線
plot(dend1, nodePar=list(pch = 2:1,cex=.4*2:1, col = 2:3), horiz=TRUE)      #縱向
par(opar)

---

2.2 例2. 女中學生體型指標分類（對變數進行分類）（一個小函式讓圖形好看一點）

對305名女中學生測量8個體型指標，給出相應的相關矩陣。將相關係數看成相似係數，定義距離為：dij = 1 - rij，用最長距離法作系統聚類分析

# 輸入相關矩陣. 
x<- c(1.000, 0.846, 0.805, 0.859, 0.473, 0.398, 0.301, 0.382,
      0.846, 1.000, 0.881, 0.826, 0.376, 0.326, 0.277, 0.277, 
      0.805, 0.881, 1.000, 0.801, 0.380, 0.319, 0.237, 0.345, 
      0.859, 0.826, 0.801, 1.000, 0.436, 0.329, 0.327, 0.365, 
      0.473, 0.376, 0.380, 0.436, 1.000, 0.762, 0.730, 0.629, 
      0.398, 0.326, 0.319, 0.329, 0.762, 1.000, 0.583, 0.577, 
      0.301, 0.277, 0.237, 0.327, 0.730, 0.583, 1.000, 0.539, 
      0.382, 0.415, 0.345, 0.365, 0.629, 0.577, 0.539, 1.000)
names<-c("身高 x1", "手臂長 x2", "上肢長 x3", "下肢長 x4", "體重 x5", 
         "頸圍 x6", "胸圍 x7", "胸寬 x8")
r<-matrix(x, nrow=8, dimnames=list(names, names))
# 作系統聚類分析, 
# as.dist()的作用是將普通矩陣轉化為聚類分析用的距離結構. 
d<-as.dist(1-r); hc<-hclust(d); dend<-as.dendrogram(hc)
# 寫一段小程式, 其目的是在繪圖命令中呼叫它, 使譜系圖更好看.
nP<-list(col=3:2, cex=c(2.0, 0.75), pch= 21:22, 
         bg= c("light blue", "pink"),
         lab.cex = 1.0, lab.col = "tomato")
addE <- function(n){
   if(!is.leaf(n)) {
       attr(n, "edgePar") <- list(p.col="plum")
       attr(n, "edgetext") <- paste(attr(n,"members"),"members")
   }
   n
}
# 畫出譜系圖.
op<-par(mfrow=c(1,1), mar=c(4,3,0.5,0))
de <- dendrapply(dend, addE); plot(de, nodePar= nP)
par(op)

從圖中容易看出，變數x2（手臂長）與x3（上肢長）最先合併成一類，解析來是變數x1（身高）與x4（下肢長）合併成一類.再合併就是將新得到的兩類合併成一類（可稱為“長”類）.後面要合併的是x5（體重）與x3（頸圍）.再合併就是將x7（胸圍）加到新類中.再加就是x8（胸寬）.最後合併為一類.

---

2.3 類的個數的確定

在聚類問題中如何確定合適的分類個數，是一個非常困難的問題，至今沒有特別令人滿意的答案。

* 2.3.1 rect.hclust()函式確定分類個數*

在R軟體中，與確定類的個數有關的函式是rect.hclust()函式，它的本質是由給定類的個數或給定閾值（類與類間的距離）來確定分類個數

 rect.hclust(tree, k = NULL, which = NULL, x = NULL, h = NULL, border = 2, cluster = NULL)

其中tree是由hclust生成的結構，k是類的個數，h是譜系圖中的閾值，要求分成的各類的距離大於h.border是數或向量，表示矩形框的顏色

對2.2中例2的八個體型指標的聚類分析中，將變數分為三類，即k=3，其程式和計算結果如下：

plclust(hc,hang = -1)   #與plot類似用法
re <- rect.hclust(hc,k=3)   #加入k=3的類

根據譜系圖確定分類個數的準則：
A.各類重心的距離必須很大
B.確定的類中，各類所包含的元素都不要太多
C.類的個數必須符合實用目的
D.若採用集中不同的聚類方法處理，則在各自的聚類途中因發現相同的類。

得到身高（x1），手臂長（x2），上肢長（x3），下肢長（x4）分為第一類；胸寬（x8）為第二類；體重（x5），頸圍（x6），胸圍（x7）分為第三類。

上邊程式中，plclust()函式是另一種繪譜系圖的函式，與plot()函式所畫圖形略有差別，使用格式如下：

 plclust(tree, hang = 0.1, unit = FALSE, level = FALSE, hmin = 0,
         square = TRUE, labels = NULL, plot. = TRUE,
         axes = TRUE, frame.plot = FALSE, ann = TRUE,
         main = "", sub = NULL, xlab = NULL, ylab = "Height")

其中tree是由hclust()函式生成的物件，其他引數設定與plot()函式中的相同

* 2.3.2 例3.NbClust包確定在一個聚類分析裡類的最佳數目*（用flexclust包裡自帶資料作為例子,clust1是對變數聚類,資料是相關係數矩陣，無法使用）

par(ask=TRUE)
opar <- par(no.readonly=FALSE)
#輸入flexclust包自帶資料nutrient，是27個不同種類的肉的5個成分含量
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")
#用NbClust()函式判斷生成類的個數
library(NbClust)
nc <- NbClust(nutrient.scaled, distance="euclidean", 
              min.nc=2, max.nc=15, method="average")
#設定圖形
par(opar)
table(nc$Best.n[1,])
#標準圖
barplot(table(nc$Best.n[1,]), xlab="Numer of Clusters", ylab="Number of Criteria",main="Number of Clusters Chosen by 26 Criteria") 

函式給出最佳類別數量為2，選擇2類

#使用前先清除歷史圖片，或者將第二個程式碼塊執行兩次
#每個種類的個數
clusters <- cutree(fit.average, k=2) 
table(clusters)
#每個種類的中位數
aggregate(nutrient, by=list(cluster=clusters), median) 
#標準化後的中位數
aggregate(as.data.frame(nutrient.scaled), by=list(cluster=clusters),median)
#畫圖，k選擇2
plot(fit.average, hang=-1, cex=.8,main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=2)

根據實際情況，選擇5類

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)
plot(fit.average, hang=-1, cex=.8,  main="Average Linkage Clustering\n5 Cluster Solution")
rect.hclust(fit.average, k=5)   #k選擇5

---

2.4 例4. 一個具體的例項來總結前面介紹的聚類分析的方法（cluster2、對樣本分類）

列出1999年全國31個省市自治區，的城鎮居民家庭平均每人全年消費性支出的八個主要指標（變數）資料，這八個變數是x1食品，x2衣著，x3家庭裝置用品及服務，x4醫療保健，x5交通與通訊，x6娛樂教育文化服務，x7居住，x8雜項商品和服務

例. 某研究者收集了24種菌株，其中17~22號為已知的標準菌株，它們分別取自牛、羊、犬、豬、鼠、綿羊，其他的為xx菌株。測得各菌株的16種脂肪酸百分含量，試作樣品聚類分析，以便了解哪些為止菌株與已知的標準菌株在全部指標上最為接近。

#用sas中相同的資料來源，輸入資料，列名稱是x1-x8，行名稱是1-30的序號
y <- read.csv("cluster2.csv")
#生成距離結構
bacterium <- dist(scale(y))  #method預設值預設是Euclide距離,scale函式先將y作標準化變換
#用四種不同的距離方法生成系統聚類
hc1<-hclust(bacterium, "complete")  #最長距離法
hc2<-hclust(bacterium, "average")   #類平均法
hc3<-hclust(bacterium, "centroid")  #重心法
hc4<-hclust(bacterium, "ward.D2")      #離差平方和法
#以上聚類模型畫圖

* 2.4.1 rect.hclust()函式，根據譜系圖，結合實際情況選定類別個數*

k <- 3   #根據譜系圖，結合實際情況，自行選定類別個數
#生成兩行一列  影象距離邊界的距離
opar<-par(mfrow=c(2,1), mar=c(2,4,2,2))
#hang是表明譜系圖中各類所在的位置 當hang取負值時，譜系圖中的類從底部畫起  生成譜系圖
plot(hc1,hang=-1);re1<-rect.hclust(hc1,k,border="red")  #將分類結果分成5類 用紅色矩形筆跡標記
plot(hc2,hang=-1);re2<-rect.hclust(hc2,k,border="red")
#在活動裝置中返回所有圖形引數和他們的值
par(opar)
opar<-par(mfrow=c(2,1), mar=c(2,4,2,2))
plot(hc3,hang=-1);re3<-rect.hclust(hc3,k,border="red")
plot(hc4,hang=-1);re4<-rect.hclust(hc4,k,border="red")
par(opar)

由結果的譜系圖看，k=3時的最長距離法分類較為合理，分類結果是：
第一類：3、4、5、7、8、12、14、15、16、17、18、22、23、24、26、27、28、29、30； 第二類：1、2、6、10、11、25； 第三類：9、13、19、20、21

* 2.4.2 NbClust包確定類別個數*

#library(NbClust)
#預設index引數為all，會報錯，出現極小值需要設定tolerance。但是這個函式過於複雜，不作手動修改
#nc <- NbClust(scale(y),min.nc=2,max.nc=8,method = "average")
#index=“hubert”，通過圖形判斷類別個數
nc <- NbClust(scale(y),min.nc=2,max.nc=8,method = "average",index = "hubert")

沒有查詢到Hubert統計值的具體含義，猜想左邊圖k=3時極小值，右圖k=3到k=3+1有明顯下降趨勢，所以設定k=3.

---

3.動態聚類法（這裡介紹的是用K-均值方法的kmeans()函式）例5. K-均值法作cluster2的聚類分析

 系統聚類法以此形成類之後就不能改變，這就要求一次分類分得比較準確，對分類的方法要求較高，相應的計算量也會增大。尤其是對於大樣本問題，需要很長的計算時間，還會佔據足夠大的計算機記憶體。基於這種情況，產生了動態聚類，即動態聚類法。

 動態聚類法又稱逐步聚類法，其基本思想是，開始先粗略的分一下類，然後按照某種最優原則修改不合理的分類，直至類分得比較合理為止，形成最終的分類結果。這種方法具有計算量較小，佔計算機記憶體較少和方法簡單的優點，適用於大樣本的Q型聚類分析。

 這裡介紹用於動態聚類的R函式kmeans()函式，它採用的是K-均值方法。針對連續變數，也可處理有序分類變數，運算很快，但需要指定類別數。K-均值聚類法不會自動對資料進行標準化處理，需要先自己手動進行標準化分析

 kmeans(x, centers, iter.max = 10, nstart = 1, algorithm = c("Hartigan-Wong","Lloyd","Forgy","MacQueen"))

其中x是由資料構成的矩陣或資料框，centers是聚類的個數或者是初始類的中心，iter.max為最大迭代次數（預設值為10），nstart隨機集合的個數（當centers為聚類的個數時），algorithm為動態聚類的演算法（預設值為Hartigan-Wong方法）

3.1 令類的個數k=5作聚類分析

與例4.一樣，為消除資料數量級的影響，先對資料作標準化處理，然後再用kmeans()函式作動態聚類，為與前面的方法作比較，類的個數選為5

y <- read.csv("cluster2.csv")
#演算法選擇Hartigan-Wong方法，即預設狀態
km <- kmeans(scale(y),5,nstart = 20);km
#便於看出聚類後的分類情況，sort()函式對分類情況排序
sort(km$cluster)
library(fpc)
plotcluster(scale(y), km$cluster)    #生成聚類圖

“cluster”是一個整數向量，用於表示記錄所屬的聚類
“centers”是一個矩陣，表示每聚類中各個變數的中心點
“totss”表示所生成聚類的總體距離平方和
“withinss”表示各個聚類組內的距離平方和
“tot.withinss”表示聚類組內的距離平方和總量
“betweenss”表示聚類組間的聚類平方和總量
“size”表示每個聚類組中成員的數量

(between_SS / total_SS = 80.9 %)，組間的距離平方和佔了整體距離平方和的的80.9 %，也就是說各個聚類間的距離做到了最大

結合實際情況17-22為標準類，發現k=5並不符合實際情況。

3.2 類別個數k的選擇

在SAS中PROC CLUSTER會在輸出聚類結果的同時還是輸出一大坨統計量或者偽統計量，主要是這四個:
RSQ、SPRSQ、PSF以及PST2，分別指的是R^2、半偏R^2（偽R^2）、偽F和偽t^2 。

先簡單介紹下這四個量，如果把樣本分為k個類，R2統計量的思想其實和迴歸分析中的可決係數係數，即類間偏差平方和的總和與總離差平方總和的比值，如果值越大就說明k個類越能夠區分開，即聚類效果越好。但需要注意的是，跟迴歸中類似，R2的值總是隨著k的減少而減小，也就是說單純的根據R2來判斷是沒有意義，因此我們更關注的是他的變化情況，如果從k到k-1有非常顯著的下降，那麼可以認為分為k類是比較合適。基於此有人就又提出了半偏R2統計量，其實就是相鄰兩個R2的差值。F的值越大分類效果越顯著，但唯一的缺點是他只是看起來像F統計量，但實際上卻不服從F分佈。當然偽t2同樣如此。
數學公式參考網頁 https://www.52ml.net/8663.html/comment-page-1#comment-185

下面提供上述統計量的計算和影象

* 3.2.1 先看within-cluster sum of squares組內平方和的趨勢圖*

y <- read.csv("cluster2.csv")
mydata <- y
wss <- (nrow(mydata)-1)*sum(apply(mydata,2,var))
  for (i in 2:13) wss[i] <- sum(kmeans(mydata,centers=i)$withinss)
###這裡的wss(within-cluster sum of squares)是組內平方和
plot(1:13, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")

* 3.2.2 R方和偽R方*

#先對k=2的時候計算，再構造迴圈，計算k=1-7的情況，並畫圖
#分類個數i為2時，R方的公式為：
R2 <- kmeans(mydata,centers=2)$betweenss / kmeans(mydata,centers=2)$tot.withinss
   for (i in 2:7)R2[i] <- kmeans(mydata,centers=i)$betweenss / kmeans(mydata,centers=i)$tot.withinss
fakeR2 <- R2[2]-R2[1]
   for (i in 2:7)fakeR2[i] <- R2[i]-R2[i-1]
plot(R2,type="b",xlab="Number of Clusters", ylab="R Square")
plot(fakeR2,type="b",xlab="Number of Clusters", ylab="Fake R Square")

* 3.2.3 偽F和偽t^2*

#先對k=2的時候計算，再構造迴圈，計算k=1-7的情況，並畫圖
fakeF <- 1 -  kmeans(mydata,centers=2)$betweenss*(nrow(mydata)-2) /
             (kmeans(mydata,centers=2)$tot.withinss-kmeans(mydata,centers=2)$betweenss)*(2-1)
   for (i in 2:7)fakeF[i] <- 1 -  kmeans(mydata,centers=i)$betweenss*(nrow(mydata)-i) /
                              (kmeans(mydata,centers=i)$tot.withinss-kmeans(mydata,centers=i)$betweenss)*(i-1)
#偽t方沒有找到計算公式，略
plot(fakeF,type="b",xlab="Number of Clusters", ylab="fake F")

* 3.2.4 結合以上統計量一起畫圖比較*

#四行一列的圖，調整圖形左右上下的間距大小
opar<-par(mfrow=c(4,1), mar=c(4,5,3,2))
plot(1:13, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")
plot(R2,type="b",xlab="Number of Clusters", ylab="R Square")
plot(fakeR2,type="b",xlab="Number of Clusters", ylab="Fake R Square")
plot(fakeF,type="b",xlab="Number of Clusters", ylab="fake F")
par(opar)

圖一wws看出，分類個數為2-4個。從3類變成2類時，偽R方急劇減小，聚類結果的解釋度迅速下降，說明2類不合適，3類較合適。此時wss趨於平穩，R方也在3類變成2類時明顯波動，偽F值也是波動明顯，最後選定分類個數k=3

* 3.2.5 NbClust包來確定類別數量*

library(NbClust)
#預設index引數為all，會報錯，出現極小值需要設定tolerance。但是這個函式過於複雜，不作手動修改
#nc <- NbClust(scale(y),min.nc=2,max.nc=15,method = "kmeans")
#index=“hubert”，通過圖形判斷類別個數
nc <- NbClust(scale(y),min.nc=2,max.nc=15,method = "kmeans",index = "hubert")

沒有查詢到Hubert統計值的具體含義，猜想左邊圖k=3時極小值，右圖k=3到k=3+1有明顯上升趨勢，所以設定k=3.

---

3.3 選擇3.2裡的類別個數，類的個數k=3作聚類分析**

y <- read.csv("cluster2.csv")
#演算法選擇Hartigan-Wong方法，即預設狀態
km <- kmeans(scale(y),centers = 3, nstart = 20);km
#便於看出聚類後的分類情況，sort()函式對分類情況排序
sort(km$cluster)
library(fpc)
plotcluster(scale(y), km$cluster)   #生成聚類圖

24個菌株分為三類，第一類：1、2、3、4、5、7、9、10、11、16、18、22
第二類：6、8、12、13、14、15、17、19、21、23、24
第三類：20
和張佳玉sas的結果完全一致。

---

4. 結果匯出（對應sas三個過程）

* 4.1 例1.(2.1) cluster1資料的對變數分類* (k=3)

x <- read.csv("cluster1.csv")
r <- cor(x)
d<-as.dist(1-r) 
hc1<-hclust(d, "single"); hc2<-hclust(d, "complete")
hc3<-hclust(d, "median"); hc4<-hclust(d, "mcquitty")
n1 <- cutree(hc1,k=3);n2 <- cutree(hc2,k=3);n3 <- cutree(hc3,k=3);n4 <- cutree(hc4,k=3)
n1=as.data.frame(n1);n2=as.data.frame(n2);n3=as.data.frame(n3);n4=as.data.frame(n4)
n <- rbind(t(n1),t(n2),t(n3),t(n4))
x[31,] <- n[1,];x[32,] <- n[2,];x[33,] <- n[3,];x[34,] <- n[4,]
write.csv(x,"cluster1_result.csv")

* 4.2 例4.(2.4) cluster2資料對樣本分類* (k=3)

y <- read.csv("cluster2.csv")
bacterium <- dist(scale(y))
hc1<-hclust(bacterium, "complete");hc2<-hclust(bacterium, "average")
hc3<-hclust(bacterium, "centroid");hc4<-hclust(bacterium, "ward.D2")
n1 <- cutree(hc1,k=3);n2 <- cutree(hc2,k=3);n3 <- cutree(hc3,k=3);n4 <- cutree(hc4,k=3)
n1=as.data.frame(n1);n2=as.data.frame(n2);n3=as.data.frame(n3);n4=as.data.frame(n4)
n <- cbind(n1,n2,n3,n4)
y$n1 <- n$n1;y$n2 <- n$n2;y$n3 <- n$n3;y$n4 <- n$n4
write.csv(y,"cluster2_result1.csv")

* 4.3 例5.(3.3) K-均值法作cluster2的聚類分析* (k=3)

y <- read.csv("cluster2.csv")
km <- kmeans(scale(y),centers = 3, nstart = 20)
y$cluster <- km$cluster
write.csv(y,"cluster2_result2.csv")

n1,n2,n3,n4(31,32,33,34)分別為最短距離法，最長距離法，中間距離法，Mcquitty相似法的分類結果