應用案例
1 線性判別分析
執行線性判別分析可使用lda()函數,且該函數有三種執行形式,依次嘗試使用。
(1)公式formula格式
我們使用nmkat變量作為待判別變量,其他剩余的變量作為特征變量,根據公式nmkat~使用訓練集數據來運行lda()函數:
library(MASS)
library("MASS")
fitlda1<-lda(nmkat~.,datatrain) #以公式格式執行判別分析
names(fitlda1) #查看lda輸出項名稱
結果分析:我們看到,可以根據lda()函數得到10項輸出結果,分別為執行過程中所使用的先驗概率prior、數據集中各類別的樣本量counts、各變量在每一類別中的均值 means等。
fitlda1$prior #查看本次執行過程中所使用的先驗概率
fitlda1$counts #查看數據集datatrain中各類別的樣本量
結果分析:由於我們在之前的抽樣過程中采用的是nmkat各等級的等概率分層抽樣方式,因此如上各類別的先驗概率和樣本量在5個等級中都是相等的。具體的,5類的先驗概率都為0.2,之和為1,且訓練集中每一類都抽出了144個樣本。
fitlda1$means
結果分析:在如上的均值輸出結果中,我們可以看到一些很能反映現實情況的數據特征。比如,對於占地面積wfl變量,它明顯隨着租金nmkat的升高而逐步提高,我們看到在租金為等級1(少於500馬克)時,占地面積的均值僅為55.53平方米,而對於租金等級5(租金不低於1150 馬克),平均占地面積則達到了92.24平方米。面積越大的房屋租金越貴,這是十分符合常識的。
執行fitlda1可直接將判別結果輸出。
(2)數據框data.frame及矩陣matrix格式
由於這兩種函數格式的主體參數都為x與grouping,我們放在一起實現,程序代碼如下:
fitlda2<-lda(datatrain[,-12],datatrain[,12])
#設置屬性變量(除第12個變量nmkat外)與待判別變量(第12個變量nmkat)的取值
fitlda2
2.判別規則可視化
我們首先使用plot()直接以判別規則fit_ldal為對象輸出圖形,如下圖所示:
plot(fitlda1)
結果分析:從圖可以看到,在所有4個線性判別式(Linear Discriminants,即 LD)下1至5這5個類別的分布情況,不同類別樣本已用相應數字標出。
我們可以通過dimen參數的設定來控制輸出圖形中所使用的判別式個數當參數取值大於總共的判別式個數式,則默認取所有樣式,比如下面分別以 dimen=1和2為例生成圖形,分別如圖所示:
plot(fitlda1,dimen=1) #變量太多,畫布顯示不出來
plot(fitlda1,dimen=2)
結果分析:1類和5類比較分明,2,3,4類重疊很多。
3.對測試集待判別變量取值進行預測
下面我們來使用之前得到的判別規則fitldal來對測試集 datatest中的待判別變量nmkat 的類別進行預測。
prelda1<-predict(fitlda1,datatest) #使用類別規則fitlda1預測datatest中nmkat變量的類別
prelda1$class #輸出datatest中個樣本的預測結果
prelda1$posterior #輸出個樣本屬於每一類別的后驗概率
結果分析:將class與 posterior的輸出結果相結合來看,我們知道每一樣本屬於各類別的后驗概率最高者為該樣本被判定的類別。比如,posterior輸出項中序號為“3”的樣本屬於第2類的概率最高,約為0.501,因此該樣本在此次判別中被歸為類別2。
為了進一步評價本次判別的效果,我們可以生成測試集中nmkat變量的預測結果與其實際類別的混淆矩陣,代碼如下:
table(datatest$nmkat,prelda1$class)
table(datatest$nmkat)
結果分析:由上面的混淆矩陣,行表示實際的類別,列表示預測判定的類別。在362個測試樣本中,實際屬於第1類的有75個,而由判定結果,75個樣本中,有57個判定正確,14個錯判成第2類,3個錯判成第3類。且該矩陣的非對角線上的元素之和為120,也就是說120個樣本被判錯了,錯誤率則可以通過計算得到,為0.33=120/362。
其中,屬於第3類的樣本被錯分的個數最多,共有31個(約占總量的一半)被錯誤分類。這之中有17個被錯分入類別2,13個分入類別4。這很可能是由於類別3與類別2、4的樣本之間相似度太高,表現在圖形中即為有較大的重疊區域所導致的分類困難,正如我們在上圖所看到的,2、3、4這三個中等租金額的樣本點聚集在一起難以分割。
可以計算錯誤率來評價分類效果:
errorlda1<-sum(as.numeric(as.numeric(prelda1$class)!=as.numeric(datatest$nmkat)))/nrow(datatest)
errorlda1