R語言中的線性判别分析

      在R語言中，線性判别分析（Liner Discriminant Analysis，簡稱LDA），依靠軟體包MASS中有線性判别函數lqa()來實作。該函數有三種調用格式：

1）當對象為資料框data.frame時

   lda(x,grouping,prior = propotions,tol = 1.0e-4,method,CV = FALSE,nu,...)

2) 當對象為公式Formula時

   lda(formula,data,...,subnet,na.action)

3) 當對象為矩陣Matrix時

   lda(x,group,...,subnet,na.action)

   對于第一種情況，grouping表示每個觀測樣本的所屬類别；

   prior表示各類别的先驗機率，預設取訓練集中各樣本的比例；

                   tol表示篩選變量，預設取0.0001

   對于第二種情況，formula表示判别公式，比如，y~x1+x2+x3,或者y~x1*x1

                   data表示資料集

                   subnet表示樣本

                   na.action表示處理缺失值的方法，預設為“如果樣本中有缺失值，則lda()函數無法運作”；如果設定為na.omit，則表示“自動删除樣本中的缺失值，然                    後，進行計算”

   對于第三種情況，x表示矩陣                

   data表示資料集

                   subnet表示樣本

                   na.action表示處理缺失值的方法，預設為“如果樣本中有缺失值，則lda()函數無法運作”；如果設定為na.omit，則表示“自動删除樣本中的缺失值，然                    後，進行計算”

    下面，舉一個例子，來說明線性判别分析。我們選用kknn軟體包中的miete資料集進行算法示範。miete資料集記錄了1994年慕尼黑的住房傭金标準中一些有趣變量，比如房子的面積、是否有浴室、是否有中央供暖、是否供應熱水等等，這些因素都影響着傭金的高低。

    1.資料概況

    首先，簡單了解一下，miete資料集。

> library(kknn)
> data(miete)
> head(miete)
           

 nm wfl     bj bad0 zh ww0 badkach fenster kueche mvdauer bjkat wflkat

1  693.29  50 1971.5    0  1   0       0       0      0       2     4      1

2  736.60  70 1971.5    0  1   0       0       0      0      26     4      2

3  732.23  50 1971.5    0  1   0       0       0      0       1     4      1

4 1295.14  55 1893.0    0  1   0       0       0      0       0     1      2

5  394.97  46 1957.0    0  0   1       0       0      0      27     3      1

6 1285.64  94 1971.5    0  1   0       1       0      0       2     4      3

       nmqm rooms nmkat adr wohn

1 13.865800     1     3   2    2

2 10.522857     3     3   2    2

3 14.644600     1     3   2    2

4 23.548000     3     5   2    2

5  8.586304     3     1   2    2

6 13.677021     4     5   2    2

> dim(miete)
           

[1] 1082   17

  我們看到，該資料集一共有1082條樣本，和17個變量。下面，我們利用summary（）來檢視變量的資訊。

> summary(miete)
           

       nm              wfl               bj       bad0     zh      ww0     

 Min.   : 127.1   Min.   : 20.00   Min.   :1800   0:1051   0:202   0:1022  

 1st Qu.: 543.6   1st Qu.: 50.25   1st Qu.:1934   1:  31   1:880   1:  60  

 Median : 746.0   Median : 67.00   Median :1957                            

 Mean   : 830.3   Mean   : 69.13   Mean   :1947                            

 3rd Qu.:1030.0   3rd Qu.: 84.00   3rd Qu.:1972                            

 Max.   :3130.0   Max.   :250.00   Max.   :1992                            

 badkach fenster  kueche     mvdauer      bjkat   wflkat       nmqm       

 0:446   0:1024   0:980   Min.   : 0.00   1:218   1:271   Min.   : 1.573  

 1:636   1:  58   1:102   1st Qu.: 2.00   2:154   2:513   1st Qu.: 8.864  

                          Median : 6.00   3:341   3:298   Median :12.041  

                          Mean   :10.63   4:226           Mean   :12.647  

                          3rd Qu.:17.00   5: 79           3rd Qu.:16.135  

                          Max.   :82.00   6: 64           Max.   :35.245  

     rooms       nmkat   adr      wohn   

 Min.   :1.000   1:219   1:  25   1: 90  

 1st Qu.:2.000   2:230   2:1035   2:673  

 Median :3.000   3:210   3:  22   3:319  

 Mean   :2.635   4:208                   

 3rd Qu.:3.000   5:215                   

 Max.   :9.000                           

      我們可以選擇nmkat(淨租金)作為 待判别變量---一是，由于該變量在含義上容易受其他變量影響，為被解釋變量；二是，nmkat自身就有5個等級類别，其相應的樣本量依次為 219、230、210、208、215，即每一類的樣本量都為200個左右，分布比較均勻。

      2.資料預處理

      下面，我們将miete資料集，分為訓練集和測試集。

      為了提高判别效果，我們考慮采用分層抽樣的方式，由于miete資料集中，待判别變量nmkat的5個等級分布比較均勻，是以采用5個等級按等量抽取樣本。（如果分布不均勻，則采用按比例抽取樣本）。具體如下：

> library(sampling)
> n = round(2/3*nrow(miete)/5)
> n
           

[1] 144

    可以看到，訓練集占總樣本的2/3，測試集占總樣本的1/3，從訓練集中nmkat變量的每一個等級抽取的樣本數是144個。

> #以nmkat變量的5個等級劃分層次，進行分層抽樣
> sub_train = strata(miete,stratanames = "nmkat",size=rep(n,5),method="srswor")
> head(sub_train)
           

 #顯示訓練集抽取的情況，包括nmkat變量取值、該樣本在資料集中的序号、被抽取的機率、以及被抽取的層次。

   nmkat ID_unit      Prob Stratum

1      3       1 0.6857143       1

2      3       2 0.6857143       1

3      3       3 0.6857143       1

16     3      16 0.6857143       1

20     3      20 0.6857143       1

22     3      22 0.6857143       1

>

> #擷取如上ID_unit所對應的樣本構成訓練集，并删除變量1、3、12
> data_train = getdata(miete[,c(-1,-3,-12)],sub_train$ID_unit)
> data_test = getdata(miete[,c(-1,-3,-12)],-sub_train$ID_unit)
> dim(data_train); dim(data_test) #分别顯示訓練集、測試集的次元
           

[1] 720  14

[1] 362  14

> head(data_test)
           

   wfl bad0 zh ww0 badkach fenster kueche mvdauer bjkat      nmqm rooms nmkat

7   28    0  1   0       0       1      1       9     4 17.011071     1     1

8   36    0  1   0       0       0      1       3     4 19.710278     1     3

9   33    0  1   0       0       0      0       1     4 25.840606     1     4

10  57    0  1   0       1       0      1       9     6 11.534035     2     2

11  75    0  1   0       1       0      1       3     6 16.504533     3     5

17  79    0  1   0       0       0      0      20     4  7.507215     3     2

   adr wohn

7    2    2

8    2    2

9    2    2

10   2    2

11   2    2

17   2    2

   至此，資料了解和資料預處理過程結束，得到可以直接使用的訓練集data_train和測試集data_test。

      3.線性判别

      這裡使用公式formula格式，進行判别。首先要加載軟體包MASS,接着使用nmkat作為待判别變量，其他變量作為特征變量，根據公式nmkat~. （如果變量為y，則公式為y~. ） ，使用訓練集data_train來運作 lda()函數。

> library(MASS)
> fit_lda1 = lda(nmkat~., data_train) #以公式格式進行線性判别
> names(fit_lda1)
> fit_lda1$prior 
           

  1   2   3   4   5 

0.2 0.2 0.2 0.2 0.2 

  我們可以看到，各類别的先驗機率在5個等級中都為0.2，之和為1，即它們都相等，這與它們分别均勻對應。

> fit_lda1$means
           

      wfl       bad01       zh1        ww01  badkach1   fenster1    kueche1

1 54.87500 0.055555556 0.6041667 0.138888889 0.3888889 0.06944444 0.04166667

2 60.59722 0.013888889 0.8125000 0.027777778 0.5486111 0.08333333 0.04166667

3 66.76389 0.013888889 0.8194444 0.041666667 0.5208333 0.06944444 0.07638889

4 74.70833 0.013888889 0.8750000 0.041666667 0.6458333 0.03472222 0.07638889

5 90.10417 0.006944444 0.9375000 0.006944444 0.7708333 0.04166667 0.20833333

    mvdauer     bjkat.L      bjkat.Q      bjkat.C    bjkat^4     bjkat^5

1 14.444444 -0.21580517 -0.104562739  0.031056500 0.17323372 -0.17585847

2 11.923611 -0.12450298 -0.211398581 -0.002588042 0.20604313 -0.12642562

3 11.847222 -0.12782306 -0.145478593  0.049690399 0.16273470 -0.07349309

4 10.347222 -0.08964215 -0.127293769 -0.035197366 0.12861291  0.01137393

5  5.333333 -0.04482107 -0.009092412  0.018633900 0.02624753 -0.01574852

       nmqm    rooms        adr.L      adr.Q     wohn.L     wohn.Q

1  8.231574 2.173611 -0.019641855 -0.7654655 0.03437325 -0.3997431

2 10.988627 2.416667 -0.034373246 -0.7569604 0.08838835 -0.4252586

3 12.495436 2.597222 -0.019641855 -0.7654655 0.11294067 -0.3487121

4 14.107342 2.861111  0.004910464 -0.7909811 0.16695577 -0.3912379

5 17.108865 3.250000  0.019641855 -0.7484552 0.27498597 -0.2041241

     從上面的結果中，可以看到一些很能反映現實情況的資料特征。比如，住房面積wfl變量，它明顯随着租金nmkat的升高而逐漸提高。這與我們的常識“房子的面積越大，租金就越貴”是十分吻合的。