线性判别分析（LDA）及其在R中实现

线性判别分析（Linear Discriminant Analysis，LDA）是一种用于数据降维以及分类预测的机器学习方法，由RA Fisher（1936）提出，也称为Fisher判别式。

最初的LDA只能处理两种分类问题，Rao（1948）拓展了“多类别线性判别分析”或称“多元判别分析”，将LDA推广至适用于多分类的情形。

LDA概述

LDA目的是将数据集投影到具有良好分类特征的低维空间上，投影后类内方差最小、类间方差最大，用于表征数据结构以及识别分类，因此它在高维数据集的分析中非常流行。

与其它机器学习算法（例如神经网络、随机森林等）相比，LDA的主要优点是可更好避免过拟合以及计算简单。

LDA数据准备

已知分类

需要已知样本的归属，支持两组或多组情况。已知样本分类将为降维过程提供监督，实现判别的目的。

删除异常值

考虑从数据中删除异常值，它们通常会对LDA的精度产生影响。

正态性与方差齐性

LDA最初是为多元正态分布数据开发的，也就是LDA的前提假设数据满足多元正态性。即便数据集整体不呈正态分布，但也应至少满足在每个分组内的子数据集均应呈正态分布。

LDA对方差齐性没有很严格的要求。这意味着如果以图形展示数据，则每个分组内的数据分布都可描述为一个椭圆体（取决于正态性），并允许椭圆体位置、形状（代表了组间均值、方差）等有所不同。

尽管如此，通常推荐在执行LDA前通过标准化数据的方法，将已服从正态性但方差不等的各变量标准化为等方差（如均值0，标准差1）的结构，即让数据满足方差齐性。数据标准化对于消除变量间的量纲差异或者较大方差变量的权重时非常有效。

LDA原理简述

直观理解，当代表每个分组的椭圆体具有如下所示的相似排列和形状时（表征了数据分布模式），可选择使用LDA；当椭圆体具有明显不同的排列或形状时，一种称为二次判别分析（Quadratic Discriminant Analysis，QDA）的方法更为合适（因此通常QDA更灵活）。

本篇暂不讨论QDA，以下是LDA的工作原理概述。

以包含2个变量，9个对象（其中5个属于分组1，4个属于分组2，即上图示例）的数据矩阵为例。首先将对象投影到变量空间上，一个变量代表一个维度，这里共2个变量所以就是二维变量空间。

考虑图a中的代表两组数据中对象分布的椭圆，它们具有相似的结构，执行LDA。在图上定义任意的直线，并将椭圆中心（代表变量集的均值位置）垂直投影到直线上（图b），称为正交投影。类似的过程，将所有变量也投影到直线上（图c）。

由于各组内的数据均呈正态分布，因此对于每个椭圆的投影数据也均为正态分布（图d），均值表示投影中心，并具有一定的方差，方差反映了离散程度。由于投影直线是任意指定的，因此可存在n种投影情况，具有不同的均值和方差。然后就需从中确定最佳投影直线，所谓最佳，即数据的投影能够最大限度地区分组内和组间特征。

这里假设两组数据在各个变量维度的均值分别为μ₁和μ₂，方差为σ²₁和σ²₂，则变量空间中存在点M，可在各个变量维度通过公式(μ₁-M)²/σ²₁+(μ₂-M)/σ²₂得到一个最小值，即获得最小均值差。对于经过M点的直线，称为分割线（separating line），最佳投影直线与分割线正交，且保证各对象在该直线上的投影具有更小的组内方差，以及更大的组间方差，使组内方差与组间方差之比最大化（图e，两组数据等方差时，M点在两椭圆中心点的中间位置；或f，两组数据不等方差时，M点靠近具有较小方差的椭圆）。

该投影直线也称为超平面（hyperplane），最终实现降维的目的，并可根据对象点在超平面中投影位置进行分类。

如上简单概括了由二维降维至一维的过程。LDA最终可获得的维度k≤n-1，n为变量数量。对于多维情况，各维（LDA轴）计算的原理同上，且各轴之间正交。

再补充一点关于数据的正态性问题

由上述过程也不难看出，对于违反正态性的数据，理论上无法根据公式准确计算M点的位置，因此最终分类结果也不准确。

尽管如此，违反正态性假设的数据仍可尝试执行LDA，有些情况下仍可以得到较好的效果。例如Duda等（2001）指出，LDA在人脸和物体识别的任务中经常取得良好的表现，尽管经常与正态性假设违背。

LDA与PCA

LDA和主成分分析（PCA）均是线性降维技术，通用的LDA方法与PCA比较相似，但是除了找到使数据方差最大的成分轴外，还对最大化多个类之间距离的轴感兴趣。

因此二者的主要区别在于，PCA属于“无监督”算法，它在降维时不关注数据的分组（类）关系，目标是找到代表数据集最大化方差方向的一系列正交的主成分（特征向量）；LDA则是“监督”算法，考虑已知的类别关系，通过线性判别式定义一系列表示最大化类别之间距离的正交轴。

尽管听起来很直观，对于已知类别标签的数据而言，LDA似乎优于PCA，但并非总是如此。Martinez等（2011）在比较PCA和LDA排序图中展示的分组关系后表明，如果每类数据中涉及的对象数量相对较少，PCA的性能反而优于LDA。

R包MASS的LDA

接下来，展示R语言执行LDA的过程，以MASS包的方法为例。

数据集

iris数据集，记录了150朵鸢尾花的花朵性状测量值。

这些鸢尾花来自三种物种，分别为setosa（n=50）、versicolor（n=50）和virginica（n=50）。

包含四种性状，分别为萼片长度（sepal length，cm）、萼片宽度（sepal width，cm）、花瓣长度（petal length，cm）和花瓣宽度（petal width，cm）。

#数据集，详情 ?iris
data(iris)
head(iris)

接下来期望从中找到合适的“变量组合”，作为区分不同鸢尾花的代表特征。

* 探索性分析

首先不妨查看一下各变量的数值分布，哪些花朵性状在不同物种之间具有明显的区别。

#变量分布
library(ggplot2)
 
ggplot(reshape2::melt(iris, id = 'Species'), aes(x = value)) +
geom_histogram(aes(fill = Species), color = 'gray') +
facet_wrap(~variable, ncol = 2, scales = 'free')

直方图清晰地表明，花瓣的长度和宽度似乎可以作为潜在特征，用以区分三种鸢尾花物种。

相比直接选取部分趋势明显的变量作为代表，通过降维技术所确定的变量组合特征通常是种更为实用的选择。最首先想到的降维方法，可能就是PCA。

#通过 PCA 降维数据，观测数据特征
pca <- princomp(iris[1:4])
plot(pca$scores[ ,1], pca$scores[ ,2], col = rep(c('red', 'green', 'blue'), summary(iris$Species)))

PCA显示，第一主成分可作为区分不同鸢尾花物种花朵性状的潜在特征。

接下来接入本篇的方法部分，使用LDA实现数据降维，选择代表性的变量组合特征，以及实现分类。

正态性评估及变量标准化

如上文概述中提到，LDA要求输入数据满足（多元）正态性，可通过QQ图评估。

#QQ-plot 检验多元正态性
qqplot(qchisq(ppoints(nrow(iris[1:4])), df = ncol(iris[1:4])), mahalanobis(iris[1:4], colMeans(iris[1:4]), cov(iris[1:4])))
abline(a = 0, b = 1)

QQ图显示了示例数据集满足多元正态性。

如果正态性假设被拒绝，可尝试转化数据的方式（如log转化，但要保证这种转化方式是合理的），获得正态分布的数据。

尽管也可直接使用非正态性的数据直接执行LDA，如概述中提到的。

变量的标准化也是推荐的一步。

使标准化后的数据集保证方差齐性，对于消除变量间的量纲差异或者较大方差变量的权重时非常有效，可提高降维的精度。

#标准化数据，如标准化为均值 0，标准差 1 的结构
iris[1:4] <- scale(iris[1:4], center = TRUE, scale = TRUE)

LDA 降维及分类

为了更好地展示LDA的分类器功能，将示例数据集分为两部分，一部分作为训练集用于LDA降维及分类器构建，另一部分作为测试集进一步评估LDA预测分类的功效。

#将数据集随机分为训练集（80%）和测试集（20%）
set.seed(123)
training <- sample(rownames(iris), nrow(iris)*0.8)
 
train.data <- subset(iris, rownames(iris) %in% training)
test.data <- subset(iris, ! rownames(iris) %in% training)

经过几个准备步骤后，执行LDA。

LDA算法首先查找使类别之间距离最大化的方向（即LDA轴），作为变量响应分类的最佳线性组合，并通过这种组合进一步预测分类。

library(MASS)

#拟合模型
model <- lda(Species~., data = train.data)
model

lda()确定各组数据的平均值并计算对象属于不同组的概率，将对象划分到概率最高的组中。

Prior probabilities of groups，各组的先验概率，即已知分组中所含对象数量占总数量的比例。例如在本示例中，随机抽取的训练集的setosa组中共含有40个对象（40个鸢尾花观测个体），占训练集所有对象（总计120个鸢尾花观测个体）的33.3%。

Group means，组均值，展示了每个分组中变量的平均值。

Coefficients of linear discriminants，线性判别系数，用于形成LDA决策规则的预测变量的线性组合。例如，LD1 = 0.828*Sepal.Length + 1.438*Sepal.Width - 2.179*Petal.Length - 2.656*Petal.Width。可在降维后根据线性判别系数所得表达式，预测训练集的分类。

Proportion of trace，可以将它理解为各LDA轴的“方差解释率”，以评估各LDA轴的重要程度，示例结果显示LDA第一主轴是非常重要的。

#作图观测对象在前两 LDA 轴中的投影
plot(model, col = rep(c('red', 'green', 'blue'), summary(train.data$Species)), dimen = 2)
 
#或者 ggplot2 作图
ggplot(cbind(train.data, predict(model)$x), aes(LD1, LD2, color = Species)) +
geom_point() +
stat_ellipse(level = 0.95, show.legend = FALSE)

结果显而易见了，LDA第一主轴承载了最能体现类间差异的特征。

对比训练集中对象的既定分组属性和由LDA判别的分组属性的一致性，结果可表征LDA模型的拟合精度。

#模型拟合精度评估，训练集
predictions <- predict(model, train.data)
 
#查看训练集对象的后验概率，即根据概率划分到高概率的类中
head(predictions$posterior)
#查看对训练集对象预测的分类
head(predictions$class)
#比较预测的分类和已知先验分类属性的差异，结果反映了准确度信息
mean(predictions$class == train.data$Species)

结果显示，98%以上的对象能够被分类到正确的类别中。

现在更改为测试集，进一步评估LDA分类器精度。

#模型拟合精度评估，测试集
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

结果显示，约97%以的对象能够被分类到正确的类别中，LDA分类器的精度是相对可靠的。

参考资料

https://online.stat.psu.edu/stat555/node/101/

http://sebastianraschka.com/Articles/2014_python_lda.html

https://www.writeawriting.com/research/types-discriminant-analysis/

http://www.sthda.com/english/articles/36-classification-methods-essentials/146-discriminant-analysis-essentials-in-r/

Duda R O, Hart P E, Stork D G. Pattern Classification. Wiley, 2001.

Fisher R A. The Use of Multiple Measurements in Taxonomic Problems. Annals of Human Genetics, 1936, 7(7):179-188.

Martinez A M, Kak A C. PCA versus LDA. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2001, 23(2):228-233.

Rao C R. The utilization of multiple measurements in problems of biological classification. Journal of the Royal Statistical Society, 1948, 10(2):159-203.

友情链接