一、判別分析

? ? ? ? 判別分析是一種分類技術，其通過一個已知類別的“訓練樣本”來建立判別準則，并通過預測變量來為未知類別的數據進行分類。根據判別的模型分為線性判別和非線性判別，線性判別中根據判別準則又分為Fisher判別，Bayes判別和距離判別。本文介紹最基礎的Fisher判別，又稱線性判別，R中可用MASS包內的lda()函數進行。

? ? ? ? 注：線性判別的基礎假設是數據服從正態分布

1.1 協差陣齊性檢驗

? ? ? ? 進行線性判別分析之前需要進行協差陣齊性檢驗，通過協差陣齊性檢驗的數據才可以進行線性判別。若未通過協差陣齊性檢驗，可以進行二次判別，R中可以進行二次判別的函數是MASS包內的qda()函數。

? ? ? ? 進行協差陣檢驗使用heplots包內的boxM()函數，其使用方法在我的上一個文章也有提到，用法為

• boxM()用法：boxM(formula,data,...)或boxM(Y,group,...)

????????formula類參數表示為Y~X1+X2...類型，即指定Y為因變量，X1、X2等為自變量。

????????boxM()中的formula參數通常是? “data~group”? 的形式，data是變量，group是分組變量，分組變量通常要先轉為因子

????????另外一種用法中，Y是數據矩陣或數據框，group是指定的分組變量

---

例1：

? ? ? ??一個城市的居民家庭，按其有無割草機可分為兩組，有割草機的一組記為，沒有的一組記為，割草機工廠欲判斷家庭(26,10)是否將買割草機。

：有割草機家庭		：無割草機家庭?
（$1000）	(1000)	（$1000）?	(1000)??
20.0	9.2	25.0	9.8
28.5	8.4	17.6	10.4
21.6	10.8	21.6	8.6
20.5	10.4	14.4	10.2
29.0	11.8	28.0	8.8
36.7	9.6	16.4	8.8
36	8.8	19.8	8.0
27.6	11.2	22.0	9.2
23.0	10.0	15.8	8.2
31.0	10.4	11.0	9.4
17.0	11.0	17.0	7.0
27.0	10.0	21.0	7.4

?注：輸入數據時通常將相同變量的數據放在一列，并設置分組變量標識每一個觀測的類別

> data1=read.csv('Table1.csv') > head(data1)x1 x2 y 1 20.0 9.2 1 2 28.5 8.4 1 3 21.6 10.8 1 4 20.5 10.4 1 5 29.0 11.8 1 6 36.7 9.6 1 > library(heplots) > boxM(cbind(x1,x2)~y,data=data1) #formula參數用法Box's M-test for Homogeneity of Covariance Matricesdata: Y Chi-Sq (approx.) = 0.99346, df = 3, p-value = 0.8028

????????結果顯示不同組的協差陣無顯著性差異，即通過協差陣齊性檢驗

---

1.2 建立判別模型

? ? ? ? 對通過協差陣齊性檢驗的訓練集數據可以用lda()建立線性判別模型，未通過協差陣檢驗的訓練集數據可以用qda()建立二次線性判別模型。qda()的用法與lda()相同，此處介紹lda()函數

? ? ? ? lda()和qda()都是MASS包內的函數

• lda()用法：lda(formula,data,...)或lda(x,grouping,prior=proportions,...)

????????lda()函數中的formula參數與boxM()函數不同，左側是分組變量，即groups~x1+x2+...的形式

????????另一個用法中，x是各個指標變量，grouping指定分組變量，prior指定先驗概率。當未指定prior先驗概率時，默認以邊緣概率作為先驗概率，即各水平的觀測數占總觀測數的比例為先驗概率。若在進行判別分析前已經對研究對象有一定了解，知曉各個水平出現的概率，則可以指定proior先驗概率，相當于做的是Bayes判別。

? ? ? ? 另外再提一提CV參數，CV是邏輯值，當CV=T時返回留一法交叉驗證的各項列表。其中包括訓練集數據各個觀測的分類結果和各個觀測的后驗概率。但是此時lda()返回的是列表而不是模型，并不能用predict()函數進行新數據的預測。

---

續上例：

> library(MASS) > model=lda(y~x1+x2,data=data1) > model Call: lda(y ~ x1 + x2, data = data1)Prior probabilities of groups: #先驗概率0 1 0.5 0.5 Group means: #各水平均值x1 x2 0 19.13333 8.816667 1 26.49167 10.133333Coefficients of linear discriminants: #線性判別函數的系數LD1 x1 0.1453404 x2 0.7590457

---

1.3 新數據的判別

? ? ? ? 建立判別模型后就可以對新數據進行判別了，對新數據的判別用predict()函數，但是在進行判別之前可以用table()函數對訓練集數據回代判別模型的分類結果和訓練集數據的實際類別做列聯表，相當于查看模型的殘差

• predict()用法：predict(object,...)

????????predict()函數的用處很廣泛，其用于根據各種模型擬合函數的結果進行預測，參數object就是進行預測時的要用的模型

????????第二個參數通常是newdata，即進行預測的新數據。newdata需是列表或數據框，而且這個列表或數據框的變量名需要和訓練集的變量名相同。

????????對lda模型，若未指定newdata或newdata指定不規范，則返回的是訓練集數據回代模型的結果，包括分類結果、后驗概率和線性判別函數值

• table()用法：table(x,y,...)

????????table()函數用來每個因子水平組合處構建計數的列聯表。

????????在查看lda模型的殘差時，通常指定x為訓練集數據回代判別模型的分類結果，指定y為訓練集數據各個觀測的實際類別。這樣生成的列聯表行列數相同，對角線上的是預測正確的觀測數，非對角線上的是預測錯誤的觀測數，可以由此確定模型的判對率或錯誤率。

---

續上例：

> table(data1$y,predict(model)$class)0 10 10 21 1 11 #可見24個觀測值中，有3個預測錯誤 #占訓練集總觀測數的12.5% > predict(model,newdata=list(x1=26,x2=10)) $class #判斷的分類結果 [1] 1 Levels: 0 1$posterior #新數據的后驗概率0 1 1 0.1439459 0.8560541$x #新數據的線性判別函數值LD1 1 0.8617716

????????家庭(26,10)的判別結果為1類，即有割草機家庭。那么工廠可以認為這個家庭是欲購買割草機的。

---

二、聚類分析

? ? ? ? 聚類分析是將對象的集合分組為由類似的對象組成的多個類的分析過程。本文介紹R中的系統聚類和k均值聚類

2.1 系統聚類

? ? ? ? R中的系統聚類方法用hclust()函數實現，但是在此之前需要對數據進行標準化，并且生成距離陣，以距離陣出發進行系統聚類。數據標準化可用scale()函數，生成距離陣用dist()函數

2.1.1 生成距離陣

? ? ? ? 用scale標準化數據后，生成距離陣需要定義各個樣品之間的距離，在dist()函數的method參數可以選擇樣品距離的定義，下面介紹兩個函數的用法

• scale()用法：scale(x,center=TRUE,scale=TRUE)

? ? ? ? 其中x為數據矩陣或數據框，返回的是一個對每列變量標準化后的矩陣

• dist()用法：dist(x,method="euclidean",diag=FALSE,upper=FALSE,p=2)

? ? ? ? x為數據矩陣或數據框，函數返回的是距離陣，在R中的數據類型是dist而不是matrix

? ? ? ? 參數method指定樣品間距離的定義，默認為"euclidean"——歐幾里得距離，可選的距離還有

? ? ? ? ? ? ? ? ? ? ? ? "maximum"——切比雪夫距離

? ? ? ? ? ? ? ? ? ? ? ? "manhattan"——絕對距離

? ? ? ? ? ? ? ? ? ? ? ? "canberra"——蘭氏距離

? ? ? ? ? ? ? ? ? ? ? ? "minkowski"——閔可夫斯基距離

? ? ? ? ? ? ? ? ? ? ? ? "binary"——定性變量的距離

? ? ? ? 參數diag和upper是邏輯值，當diag=TRUE時給出對角線上的距離，也就是樣品到自身的距離；當upper=TRUE時給出上三角的值，距離陣是對稱矩陣，對角線上的距離也一定為0，進行系統聚類時這兩個參數可以不管。

? ? ? ? 參數p是當method="minkowski"時，閔可夫斯基距離的階數，p=1時為絕對距離，p=2時為歐幾里得距離，p=＋∞時為切比雪夫距離。

---

例2：

????????試對下列數據分別用系統聚類法和均值法進行聚類，并對結果進行比較分析

序號	省份	工資性收入	家庭性收入	財產性收入	轉移性收入
1	北京	4524.25	1778.33	588.04	455.64
2	天津	2720.85	2626.46	152.88	79.64
3	河北	1293.50	1988.58	93.74	105.81
4	山西	1177.94	1563.52	62.70	86.49
5	內蒙古	504.46	2223.26	73.05	188.10
6	遼寧	1212.20	2163.49	113.24	201.28?

注：僅展示部份表格

> data2=read.csv('Table_0.csv',encoding='UTF-8') > rownames(data2)=data2[,1] > data2=data2[,-1] #處理數據框 > head(data2)工資性收入 家庭性收入 財產性收入 轉移性收入 北京 4524.25 1778.33 588.04 455.64 天津 2720.85 2626.46 152.88 79.64 河北 1293.50 1988.58 93.74 105.81 山西 1177.94 1563.52 62.70 86.49 內蒙古 504.46 2223.26 73.05 188.10 遼寧 1212.20 2163.49 113.24 201.28 > st.data2=scale(data2) #數據標準化 > distance2=dist(st.data2) #生成距離陣，存入distance2中

---

2.1.2? 進行系統聚類

????????用hclust()函數進行系統聚類，函數返回聚類結果，我們可以通過查看聚類結果的各個組件得到聚類過程，也會用到plot()或者plclust()函數畫出聚類的系譜圖、用rect.hclust()函數根據聚類結果將樣品劃分類別

• hclust()用法：hclust(d,method="complete",...)

? ? ? ? d是距離陣，method指定類間距離的定義，默認為"complete"——最長距離法，可選的其他類間距離的定義有

? ? ? ? ? ? ? ? ? ? ? ? "single"——最短距離法

? ? ? ? ? ? ? ? ? ? ? ? "median"——中間距離法

? ? ? ? ? ? ? ? ? ? ? ? "centroid"——重心法

? ? ? ? ? ? ? ? ? ? ? ? "average"——類平均法

? ? ? ? ? ? ? ? ? ? ? ? "ward.D"——離差平方和法

等等

? ? ? ? 聚類完成后，聚類結果的merge組件記錄了聚類過程，height組件記錄了聚類過程中每次合并的兩類間的距離，order組件記錄了繪制聚類系譜圖的節點次序

? ? ? ? merge組件返回的聚類過程，每一行代表一次合并的兩類編號，編號為負值表示在原始數據的編號，編號為正值表示在此表中的對應行合并后的類

• rect.hclust()用法：rect.hclust(tree,k=NULL,...)

? ? ? ? rect.hclust()用于在系譜圖中繪制分類矩形，在同一個矩形內的為一類

? ? ? ? 通常通過系譜圖大致確定分類的個數，rect.hclust()函數中的k是分類個數，tree是hlust()返回的聚類結果。若將rect.hlust()函數的結果賦給某個變量，那么這個變量會存儲每個類的樣品名(原始數據的行名)

---

續上例：

> hclu=hclust(distance2,method='centroid') #重心法 > hclu$merge[1:6,] #展示聚類過程的前6步[,1] [,2][1,] -18 -23 #第一次合并的是第18號和第23號樣品，形成類1[2,] -4 -12 #第二次合并的是第4號和第12號樣品，形成類2[3,] -24 -27[4,] -28 3 #第四次合并的是第28號樣品和類3，類3是第3步合并形成的類[5,] -22 1[6,] 2 5 #第六次合并的是類2和類5，形成類6 > hclu2=hclust(distance2) #作為比較，用最遠距離法再聚類 > par(mfrow=c(1,2)) #一頁多圖 > plot(hclu) > a=rect.hclust(hclu,4) #看系譜圖，大致可分為3、4類 > plot(hclu2) > b=rect.hclust(hclu2,4)

????????左圖為重心法的聚類結果，右圖為最遠距離法的聚類結果。根據每個類內的樣品不宜過多的原則，通常最遠距離法的聚類效果更好，最遠距離法通俗的理解是，以樣品間最遠的距離作為類間距離進行合并，合并的兩個類中，最遠的兩個樣品的距離都比其他類樣品的距離都小，說明這兩個類距離足夠近，通常更有說服力。

> a #重心法系統聚類的分類結果 [[1]] 浙江 11 [[2]]天津 河北 山西 內蒙古 遼寧 吉林 黑龍江 江蘇 安徽 福建 江西 2 3 4 5 6 7 8 10 12 13 14 山東 河南 湖北 湖南 廣東 廣西 海南 重慶 四川 貴州 云南 15 16 17 18 19 20 21 22 23 24 25 西藏 陜西 甘肅 青海 寧夏 新疆 26 27 28 29 30 31 [[3]] 北京 1 [[4]] 上海 9 > b #最遠距離法系統聚類的分類結果 [[1]] 山西 安徽 湖南 廣西 重慶 四川 貴州 云南 西藏 陜西 甘肅 青海 寧夏 4 12 18 20 22 23 24 25 26 27 28 29 30 [[2]]天津 河北 內蒙古 遼寧 吉林 黑龍江 江蘇 福建 江西 山東 河南 2 3 5 6 7 8 10 13 14 15 16 湖北 廣東 海南 新疆 17 19 21 31 [[3]] 上海 9 [[4]] 北京 浙江 1 11

? ? ? ? ?根據最遠距離法系統聚類的分類結果，將我國各省、直轄市、自治區的收入水平劃分為4個水平，山西、安徽等水平相似，天津、河北等水平相似，上海獨成一類，北京、浙江為一類。說明我國西部地區收入水平普遍低于東部地區，其中沿海的北京、浙江的收入水平更高，上海的收入水平最高。

---

2.2? k均值聚類

? ? ? ? 系統聚類法的計算過程繁瑣，計算量大，當樣品量大時需要占用大量計算機內存空間。k均值聚類是一種快速聚類方法，其基本思想是將每一個樣品分配給最近中心(均值)的類中。每分配一個樣本，聚類的聚類中心會根據聚類中現有的對象被重新計算。這個過程將不斷重復直到滿足某個終止條件。

? ? ? ? k均值聚類需要將樣品隨意分成k個初始類，終止條件通常為所有的樣品到最近距離的類都是自身所在的類。在R中進行k均值聚類的函數是kmeans()，這個函數需要指定最終分類個數

• kmeans()用法：kmeans(x,centers,iter.max=10,nstart=1,...)

? ? ? ? x為數據矩陣或數據框，在進行k均值聚類前也需要對數據進行標準化

? ? ? ? centers用于指定最終聚類個數或者初始類的中心，當centers是最終聚類個數時，nstart用于指定初始類的個數。iter.max是最大迭代次數

? ? ? ? kmeans()返回的對象具有以下組件：

cluster——分類結果

centers——各個類的聚類中心

totss——總距離平方和

withinss——各組的組內距離平方和

tot.withinss——組內距離平方和的合計

betweenss——組間距離平方和

2.2.1? 聚類數目的確定

? ? ? ? 隨著聚類數目的增加，每個組內的距離平方和將會越來越小并趨于0，當聚類數目合適時，組內的距離平方和越小表示分類的效果越好。但是聚類數目不是越大越好，當通常聚類數目增加到一定大小時，每增加一個聚類數目的組內距離平方和減小的程度也會越來越小。因此我們可以畫出聚類數目與組內距離平方和合計值的碎石圖確定聚類數目。

? ? ? ? k均值聚類碎石圖的實現如下：

kmeans.screen=function(x,k=15,nstart=1,iter.max=10){twss=numeric()for(i in 1:k){twss[i]=kmeans(x,i,iter.max=iter.max,nstart=nstart)$tot.withinss}plot(1:k,twss,type='b',pch=19,xlab='聚類數目',ylab='組內平方和總計',main='k均值聚類碎石圖') }

? ? ? ? 聚類數目不宜過多，通常在2~5個之間。當組內平方和總計下降到一個較小的水平，或聚類數目加一時組內平方和總計下降幅度很小時，可以確定聚類數目。

---

續上例，進行k均值聚類：

> source("kmeans.screen.R") #加載函數 > kmeans.screen(st.data2,nstart=3,iter.max=20)

? ? ? ? ?在2~5個間選擇聚類數目，結合碎石圖可以確定聚類數目為4個，因為4個到5個的組內平方和總計下降幅度很小。

---

2.2.2? 進行k均值聚類

? ? ? ? 確定好聚類數目就可以進行k均值聚類了，建議畫碎石圖時的參數和進行k均值聚類時的參數一致，否則碎石圖的結果可能會不太可信。通過查看kmeans()返回的cluster模塊可以查看分類結果。

? ? ? ? kmeans()并沒有類似rect.hclust()的將分類結果歸類的函數，但是可以編寫一個函數實現這一簡單的功能：

rect.kmeans=function(kmean){ #kmean是k均值聚類的返回對象clu=kmean$cluster;class=levels(as.factor(clu))numeric.class=as.numeric(class);results=list()for(i in numeric.class){results[[class[i]]]=names(clu)[clu==i]}results }

---

續上例：

> k.clu=kmeans(st.data2,4,nstart=3,iter.max=20) > k.clu K-means clustering with 4 clusters of sizes 14, 5, 10, 2Cluster means:工資性收入 家庭性收入 財產性收入 轉移性收入 1 -0.4255246 -0.7442161 -0.4629740 -0.4096606 2 0.9116437 1.0425867 0.4206866 0.4127901 3 -0.4612330 0.7463783 -0.1949991 -0.2512104 4 3.0057285 -1.1288454 3.1640973 3.0917010Clustering vector: #這個是分類結果，也可以通過k.clu$cluster獲取北京 天津 河北 山西 內蒙古 遼寧 吉林 黑龍江 上海 江蘇 浙江 4 2 3 1 3 3 3 3 4 2 2 安徽 福建 江西 山東 河南 湖北 湖南 廣東 廣西 海南 重慶 1 2 1 3 3 3 1 2 1 3 1 四川 貴州 云南 西藏 陜西 甘肅 青海 寧夏 新疆 1 1 1 1 1 1 1 1 3 Within cluster sum of squares by cluster: [1] 5.823392 6.871055 5.067501 6.804320(between_SS / total_SS = 79.5 %)Available components:[1] "cluster" "centers" "totss" "withinss" "tot.withinss" [6] "betweenss" "size" "iter" "ifault" > source("rect.kmeans.R") #加載將分類結果歸類的函數 > rect.kmeans(k.clu) $`1`[1] "山西" "安徽" "江西" "湖南" "廣西" "重慶" "四川" "貴州" "云南" "西藏" "陜西" [12] "甘肅" "青海" "寧夏"$`2` [1] "天津" "江蘇" "浙江" "福建" "廣東"$`3`[1] "河北" "內蒙古" "遼寧" "吉林" "黑龍江" "山東" "河南" "湖北" [9] "海南" "新疆" $`4` [1] "北京" "上海"

? ? ? ? 可看到k均值聚類的結果，將我國各省、直轄市、自治區的收入水平分為4類，山西、安徽等地為一類；天津、江蘇等地區為一類；河北、內蒙古等地區為一類；北京和上海為一類。

? ? ? ? k均值聚類的結果也有說明我國中西部地區與東部地區收入水平低的趨向，相比于系統聚類的結果，k均值聚類對南北方的省份也進行了一定的劃分。根據現實各省份總體的發展情況，中西部偏南地區的省份如山西、安徽（第1類）收入水平較低，大部分其余中西部和北方以及海南（第3類）的收入水平次低，而沿海的天津、江蘇等省份（第2類）的收入水平較高，作為我國政治中心和經濟中心的北京和上海（第4類）的收入水平最高。

---

三、小結

對線性判別和聚類分析所用到的函數以及應用場景進行總結：

判別分析
函數	應用場景
boxM()	協差陣齊性檢驗
lda()	建立線性判別模型
qda()	建立二次判別模型
predict()	模型預測
table()	做列聯表

系統聚類		k均值聚類
函數	應用場景	函數	應用場景
scale()	標準化數據	scale()	標準化數據
dist()	生成距離陣	kmeans.screen()	自編，畫碎石圖
hclust()	進行系統聚類	kmeans()	進行k均值聚類
plot()	畫圖，系統聚類中用于畫系譜圖	rect.kmeans()	自編，將k均值聚類的分類結果歸類
rect.hclust()	返回分類結果并在系譜圖中圈出

總結

以上是生活随笔為你收集整理的可用于 线性判别、聚类分析 的R语言函数总结的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。