【原】【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