R語言數據分析與挖掘(第八章):判別分析(1)——距離判別法

  • 2019 年 12 月 13 日
  • 筆記

判別分析概述

判別分析是判斷個體所屬類別的一種多元統計分析方法。它在醫學領域有着廣泛的應用,主要有疾病診斷、疾病預測和病因學分析。例如,根據病人的癥狀、生化指標判斷病人得的是什麼疾病,根據病人癥狀的嚴重程度或者指標的高低預測病人的預後等等。比如,高血壓、高血糖、動脈硬化程度這些都是腦血管疾病的患病危險因素;那麼如果知道了人體的這些指標,並對這些數據進行分析,就可以對尚未明確診斷的人是否發生腦血管疾病進行預測;對於很可能是腦血管疾病的人就可以事先給予預防,或者在入院後儘快得到救治,提高診療有效率。

判別分析也屬於對事物現象進行分類的統計分析方法,它和聚類分析不同的地方在於:聚類分析(後面會講)事先並不知道分型情況,而判別分析需要事先知道分型情況,已知的分型數據又叫訓練數據。判別分析需要事先得到一些已經明確知道診斷結果的訓練數據,利用這些數據建立判別準則,然後依據準則對未知類別的預測值進行判別。如果是對於分類不明的數據,可以先用聚類分析對這組數據進行分類,然後再用判別分析對新建立的類別進行判斷。

 在判別分析中,因為判別準則的不同,可分為多種判別分析法。常用的有費歇爾(Fisher)判別分析、貝葉斯(Bayes)判別分析和距離判別分析。我們這裡先介紹距離判別法。

距離判別的基本思想是樣品X離哪個總體的距離最近,就判斷X屬於哪個總體。在判別法中根據不同的功能需求,會經常用到dist()、mahalanobis()和wmd()這三個函數。

函數介紹

1

dist()函數

進行距離判別的第一步是計算"距離",最常用的函數為dist(),該函數按照指定方法計算數據矩陣之間的距離,書寫格式為:

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

參數介紹:

X:指定用於計算距離的數據對象,可以為一個矩陣、數據框或者「dist" 類的對象;

Method:指定測度距離的方法,"euclidean" 表示歐氏距離,"maximum" 表示切比雪夫距離,"manhattan"表示絕對值距離,"canberra"表示 Lance距離,"binary"表示定性變量距離,"minkowski"表示明科夫斯基距離,使用時要指定p值,默認值為"euclidean";

Diag:邏輯值,為TRUE時表示輸出距離矩陣的對角線,默認值為FALSE;

Upper:邏輯值,為TRUE時表示輸出距離矩陣的上三角部分,默認值為FALSE;

p:指定明科夫斯基距離的權數,當參數method設置為"minkowski"時, 此參數需要進行設定,默認值為2。

2

mahalanobis()函數

不難發現,函數dist()不能用於計算馬氏距離,下面介紹一個專門用於計算馬氏距離的函數: mahalanobis(), 其基本書寫格式為:

mahalanobis(x, center, cov, inverted = FALSE, ...)

參數介紹:

x:指定用於計算距離的數據對象,p維的數據向量或矩陣;

Center:指定分佈的均值,即總體均值;

Cov:指定分佈的協力差,即總體協方差,一般用樣本的協方差進行評估;

inverted:邏輯值,若為TRUE,則表示參數cov應該包括協方差的逆。

3

wmd()函數

上述介紹的兩個函數均返回距離值,而不能直接判別,下面介紹一個可直接用於判別的函數: wmd(), 該函數存在於WMDR包中,可用於實現加權馬氏距離的判別,它利用函數mahalanobis()計算出馬氏距離,然後進行判別分析,最終返回包含結果和準確度的表單,其基本書寫格式為:

wmd(TrnX,TrnG,Tweight = NUL, TstX = NULL, var.equal = F)

參數介紹:

TmX:指定訓練集的數據對象,可以為矩陣或數據框;

TrnG:一個因子類的向量,用於指定已知的訓練樣本的分類;

Tweight:指定權重,若沒有進行指定,則軟件默認使用主成分分析中的相應貢獻率作為權重,默認值為NUll,表示不進行加權,採用傳統的馬氏距離判別法;

TstX: 指定測試集的數據對象,可以為向量、矩陣或數據框,若為向量,則將被識別為單個案例的行向量,默認值為NULL,表示直接對訓練集進行判別;

Var.equal: 指定不同類別之間是否具有相同的協方差,默認值為F。

需要注意的是,函數wmd()中訓練集的樣本量與測試集的樣本量相等,否則R語言會報錯。

綜合案例:基於距離判別的iris數據集分類

下面利用iris 數據集進行操作演練,由於該數據集中鳶尾花的種類有三種,下面將原始數據分為訓練集和測試集,分別包含隨機抽取的100個樣本,具體操作如下:

> set.seed(1234)  > sa<-sample(1:150,100)  > sa    [1]  18  93  91  92 126 149   2  34  95  73  98  76  40 127 138 114  39   [18]  36  25  31  42 136  21   6  28 102  66 113 125 137  55  32 133  60   [35]  22  88  23  30 112  90  61  71 143  67  35  53 109  50 132  78   8   [52] 131 104  49  15  48  47  70  17 101 148   4 146 144 128 110  26  43   [69]   5  46  10 140  87  85   7 134  29 121  24 142  65  33  80  37  13   [86]  59 122  20  68 120  62  54 100  58 141  74 147 111 145  38  > dtrain<-iris[sa,1:5]  > head(dtrain)      Sepal.Length Sepal.Width Petal.Length Petal.Width    Species  18           5.1         3.5          1.4         0.3     setosa  93           5.8         2.6          4.0         1.2 versicolor  91           5.5         2.6          4.4         1.2 versicolor  92           6.1         3.0          4.6         1.4 versicolor  126          7.2         3.2          6.0         1.8  virginica  149          6.2         3.4          5.4         2.3  virginica  > dtest<-iris[-sa,1:4]  > head(dtest)     Sepal.Length Sepal.Width Petal.Length Petal.Width  1           5.1         3.5          1.4         0.2  3           4.7         3.2          1.3         0.2  9           4.4         2.9          1.4         0.2  11          5.4         3.7          1.5         0.2  12          4.8         3.4          1.6         0.2  14          4.3         3.0          1.1         0.1

上面代碼表示:首先利用函數set.seed()釋放隨機種子,然後利用函數sample()從1到150中隨機抽取100個數據,將抽取的數據對應到數據中觀測樣本的編號,得到訓練數據集,剩下的數據集即為測試集。下面對訓練集中的觀測樣本安裝變量Species分類:

> d1<-subset(dtrain,Species=="setosa");dim(d1)  [1] 38  5  > d2<-subset(dtrain,Species=="versicolor");dim(d2)  [1] 29  5  > d3<-subset(dtrain,Species=="virginica");dim(d3)  [1] 33  5

上述輸出結果顯示:在訓練集中種類為setosa的鳶尾花共有38條記錄,種類為versicolor()和virginica的鳶尾花分別有29和33條記錄。

下面利用函數mahalanobis()計算馬氏距離:

> ma1<-mahalanobis(dtest,colMeans(d1[,1:4]),cov(d1[,1:4]))  > ma2<-mahalanobis(dtest,colMeans(d2[,1:4]),cov(d2[,1:4]))  > ma3<-mahalanobis(dtest,colMeans(d3[,1:4]),cov(d3[,1:4]))  > (distance<-cbind(ma1,ma2,ma3,iris[-sa,5]))               ma1         ma2        ma3  1      0.3590061 135.3154437 271.444570 1  3      1.5544105 112.3495585 242.943354 1  9      3.4278559  86.6318539 203.172929 1  11     1.8993357 151.9381072 290.189235 1  12     1.9068175 117.1927197 235.675433 1  14     9.0397389 110.4697730 241.640736 1  16     9.8088108 211.1672458 359.672511 1  19     6.3669561 147.5440882 282.272697 1  27     4.0876272 104.1682597 226.561182 1  41     1.8477792 132.0256229 270.980300 1  44    17.8963186 102.9131186 222.375954 1  45    12.0006760 125.8813043 238.336324 1  51   519.2519291   5.8057277  26.778995 2  52   483.9369049   2.7504623  21.480586 2  56   440.2499733   3.2297851  14.484601 2  57   550.1907724   3.7024849  16.423962 2  63   314.6573331   6.6213187  26.256011 2  64   505.8581500   1.9220704  11.529449 2  69   509.6105640  14.5675407  11.236568 2  72   347.7106464   2.1919434  28.310069 2  75   406.4106931   1.7786523  25.078831 2  77   548.0360755   3.6213575  14.421208 2  79   485.3833731   0.7923886  12.021558 2  81   278.4175024   1.7677099  25.346488 2  82   246.8996580   2.8295498  30.772019 2  83   308.7761180   1.1340186  28.096945 2  84   653.7060829   8.5350447   3.495719 2  86   506.4900210   6.8295042  22.639212 2  89   358.5380439   3.0143070  26.291627 2  94   189.8000871   4.8847665  37.271789 2  96   357.0214367   5.0103361  27.493277 2  97   378.1639201   1.3599246  21.536322 2  99   169.1113735  11.7826462  52.221092 2  103 1041.4600884  22.9217779   1.827084 3  105 1041.9268906  28.4050506   2.251143 3  106 1286.1243510  36.4572358  10.796066 3  107  560.0828547  18.2689904  10.068365 3  108 1082.7402091  29.1658678   8.484063 3  115  951.7079875  57.3638916   6.100539 3  116  943.8604134  31.8419254   1.768798 3  117  814.6300092  10.5481905   1.706876 3  118 1336.8526085  31.3734151  11.353519 3  119 1489.6049927  69.4397271  29.160649 3  123 1302.1372419  44.5459977  16.435043 3  124  666.3834809  11.6169643   2.831568 3  129  948.9251799  26.6719664   1.741693 3  130  864.8827533  15.9814974   6.508066 3  135  747.0655076  29.1338728  11.601745 3  139  635.8783895   7.2968399   4.865241 3  150  708.1090088   9.5614468   3.787757 3

上述代碼表示:分別對訓練集計算三種類別的馬氏距離,其中函數colMeans()表示按列計算均值;訓練集中每一個觀測樣本分別對應三個馬氏距離,然後利用函數cbind()將三個馬氏距離值與原始數據集中測試樣本對應的分類合併在一起,輸出結果如上所示。對於測試集中的每一個觀測樣本而言,三個馬氏距離中最小的那一個所對應的類別即為測試樣本屬於的類別,如第一條記錄中,第一個馬氏距離的值明顯小於另外兩個,故第一條記錄應歸為第一類,即該鳶尾花屬於setosa類,與原始數據集中的分類一致,說明分類正確。

上述分類過程沒有直接得到結果,需要比較每一個測試樣本對應的三個馬氏距離進行分類,如需了解該分類的效率,還要進行後續操作,分類過程較為煩瑣,下面利用函數wmd()直接得到分類結果:

install.packages("WMDR")# 對高版本的R已經不適用  library(WMDR)  dta<-iris[,1:4]  species<-gl(3,50)  wmd(dta,species)  wmd(dta,species,diag(rep(0.25,4)))#設定權重,增加準確率。

函數wmd()的輸出結果中,第一部分表示對150個觀測值進行分類的結果(由於函數中沒有指定測試集和訓練集,故軟件默認訓練集和測試集均為同一個); "num of wrong judgement"表示判別錯誤的樣本編號,可以看到有8個樣品錯判; "samples divided to"和"samples actually belongs to" 分別表示上述樣品的錯誤判別歸類和真實類別;"percent of right judgement" 表示判別結果的正確率,由(150-8) /150得到正確率為94.7%。