【原】【R分享|實戰】LDA 線性判別分析

科白君 2021-12-19

展開全文

“ No one konws everything, and you don't have to.” --科白君

"R數據分析"專題·第22篇

編輯 |Linnaea borealis

4502字 | 14分鐘閱讀

本期推送內容

之前我們分享了不少降維相關的分析方法，例如PCoA，NMDS，PCA分析，它們都是無監督學習。無監督學習指我們事先沒有任何訓練樣本，直接對數據進行建模。無監督學習的主要算法是聚類，聚類目的在于把相似的東西聚在一起，主要通過計算樣本間和群體間距離得到。與之相對的便是有監督學習，它通過已有的訓練樣本得到一個最優模型，再利用這個模型將所有的輸入映射為相應的輸出，對輸出進行簡單的判斷從而實現預測和分類。在這一期我們將與大家分享有監督學習中LDA分析的基本知識，以及如何在R語言中實現LDA分析與預測。

線性判別分析(LDA)的基本介紹

本期內容提到的LDA分析全稱是Linear discriminant Analysis，即線性判別分析。最早由Fisher在1936年提出，多用于數據降維以及分類預測，例如：①根據給出的性狀指標，區分物種類別。②判斷一個人信用標準的好壞，③判斷學生是否能被高校錄取。LDA與回歸分析類似，但它的解釋變量是分類的而不是連續的。LDA的中心思想可以用一句話概括："投影后類內方差最小，類間方差最大"，換句話說就是我們將不同種類的高維數據投影到低維度上，希望投影結果中相同種類數據的投影點盡可能接近，而不同種類數據的中心點盡可能遠離。

如下圖所示的藍紅兩類數據，我們試圖將他們投影到一維上，保證同類相近，不同類分離。仔細觀察兩種投影方案，我們可以發現第二種方案的投影效果要比第一種好，因為它不僅將兩類數據完全分離開，且二者在自己的位置更為集中。以上就是LDA的主要思想了，在實際應用中，我們的數據是多個類別的，我們的原始數據一般也是超過二維的，投影后的也一般不是直線，而是一個低維的超平面。

相比于DA(判別分析)，LDA突出的是“Linear(線性)”，它試圖按預先分類找到能夠分離總體樣本的最佳線性組合(函數)。

Z便是上文中提到最佳線性函數。

LDA與PCA的異同點

作為常用的線性降維方法，LDA與PCA有很多異同點。LDA是有監督的降維方法，在降維時它會考慮已知的分類關系，通過線性判別式區分出一系列類別間的差異，而PCA是無監督的降維方法，它在降維時不關注數據的分組，目的是找到代表數據集方差最大化方向的一系列正交的主成分軸。可能有點抽象，我們可以根據以下圖像來評估在不同分布的數據中LDA與PCA的表現。

LDA傾向于分類性能最好的投影方向，而PCA選擇樣本點投影具有最大方差的方向。當兩組數據方差大小相近時，LDA的分類性能優于PCA。

在某些方面，如每類數據中涉及的對象數量相對較少或是均值相近時，PCA的性能反而優于LDA。

LDA在R語言中的實現

在使用LDA分析之前，我們得清楚它的幾點假設：

1）樣本量容量：樣本量應該超過自變量的數目。根據經驗，對于少數（4或5）個自變量，樣本量應該超過20。假如樣本容量為n，那自變量數目應小于n-2。雖然這種低樣本量可能有效，但通常不鼓勵這樣做，最好有4~5倍的樣本量。

2）正態分布：測試數據最好符合多元正態分布。你可以用頻率分布的直方圖或者mshapiro.test()函數對測試數據進行檢驗。對于LDA來說，正態分布并不是必須的，如果非正態性并不是由異常值引起的，那么結果仍然是可靠的。

3）方差齊次：LDA對方差-協方差矩陣的異質性非常敏感。在接受一項重要研究的最終結論之前，最好回顧一下組內方差和相關性矩陣。可以用散點圖來檢驗方差齊性，使用數據轉換方式來修正非其次。

我將從機器學習的角度介紹LDA的功能，首先將數據集分為兩部分，一部分作為訓練集構建LDA分類預測模型，一部分作為測試集評估預測模型的精確性。我們使用R中自帶的iris數據集，數據集內包含 3 類共 150 條記錄，每類各 50 個數據，每條記錄都有 4 項特征：花萼(Sepal)長度、花萼寬度、花瓣(Petal)長度、花瓣寬度，可以通過這4個特征預測鳶尾花卉屬于（setosa, versicolour, virginica）中的哪一品種。而LDA可以通過預先提供的品種分類，對特征數據進行降維投影。

library(MASS)
library(ggplot2)
#iris <- scale(iris[,1:4]) #對數據進行標準化
set.seed(1)#設置種子保證(包含隨機函數的)代碼結果可重復
trainset<- sample(rownames(iris),nrow(iris)*0.7) #隨機抽取訓練集
traindata<- subset(iris, rownames(iris) %in% trainset) #區分訓練級與測試集數據
testdata<- subset(iris, !rownames(iris) %in% trainset)

ldamodel<- lda(traindata, Species~.)
ldamodel

①：Coefficents of linear discriminants 是每個分類變量的線性判別系數，可以根據線性函數表達式Z=b₁x₁+b₂x₂+ b₃x₃+ b₄x₄生成得到用于LDA分類決策的線性回歸組合。例如LD1 = 0.828*Sepal.Length + 1.438*Sepal.Width - 2.179*Petal.Length - 2.656*Petal.Width，可在降維后預測訓練集的分類

②：Proportion of trace，類似于PCA中的“方差解釋率”，可用于評估LDA各軸的重要性。

graphset<- cbind(trainset, predict(ldamodel)$x)#通過predict函數獲得數據集通過LDA的投影點坐標并構建繪圖數據集
ggplot(graphset, aes(LD1,LD2)+
geom_point+
theme_bw()+
theme(panel.grid.major = element_blank(),
      panel.grid.minor = element_blank())+
stat_ellipse(level = 0.95)+
xlab("LDA1(99.2%)")+ ylab("LDA2(0.8%)")#Proportion of trace

根據圖片可以看出，LDA投影的第一軸將訓練數據集區分的效果最好，接下來讓我們來檢驗模型對訓練集和測試集分類的精確度。

predictions <- predict(ldamodel, traindata)
mean(predictions$class == traindata$Species)

Predictions<- predict(ldamodel, testdata)
mean(predictions$class == testdata$Species)

在沒有對模型進行優化的情況下，訓練集97%的對象能被分類到正確的類別中，而測試集中所有的對象都匹配到正確的類別中，說明LDA分類模型的精確度是相當可靠的。我們證明了LDA分類的可信度，現在就可以試著用它來對數據集進行降維分類了。

library(tidyverse)
lda<- lda(Species~.,iris) %>% 
    predict()
cbind(iris,lda$x) %>%
    ggplot(.,aes(LD1, LD2,color=Species)) +
    geom_point()+
    theme_bw()+
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())+ 
stat_ellipse(level = 0.95)+
xlab("LDA1(99.12%)")+ ylab("LDA2(0.88%)")

LDA進行預測分類

既然知道了LDA可以根據預先提供的分類信息準確地對數據集進行分類，那我們是否可以用a數據集中的分類特征訓練機器學習模型，再使用模型去預測具有相同分類特征的b數據集呢？

library(mlr)
test<- makeClassifTask(data=iris, target = "Species") #訓練lda機器學習模型
lda <- makeLearner("classif.lda")
ldaModel <- train(lda, test)
LdaModelResult <- getLearnerModel(ldaModel)

#LdaPreds <- predict(LdaModelResult)$x
#head(LdaPreds) #正常的lda降維分
kFold <- makeResampleDesc(method = "RepCV", folds = 10, reps = 50,
                          stratify = TRUE)
ldaCV <- resample(learner = lda, task = test, resampling = kFold,
                  measures = list(mmce, acc))#10倍交叉檢驗，檢驗模型精確度

交叉驗證的結果顯示模型的準確度達到98%

newcase<- tibble(Sepal.Length= runif(50,min=4,max=8), 
                 Sepal.Width= runif(50,min=2,max=4.5), 
                 Petal.Length= runif(50,min=1,max=7), 
                 Petal.Width= runif(50,min=0,max=2.5))#創建新的待測數據集

case<- predict(ldaModel,newdata = newcase)#預測數據集結果
case$data