lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

R语言中的判别分析方法

已有 2586 次阅读 2020-1-12 10:30 |个人分类:R统计分析|系统分类:科研笔记| _R语言, _机器学习, _分类, _排序分析

R语言中的判别分析方法

 

判别分析(Discriminant analysis)用于基于一个或多个预测变量来预测对象属于某个给定类的概率,它对于连续或分类预测变量均适用。与logistic回归相比,判别分析在处理多分类问题时更稳健。

几种常见的判别分析类型主要如下。

线性判别分析Linear discriminant analysisLDA),使用预测变量的线性组合预测给定的观测类别。LDA假设预测变量服从(多元)正态分布,并且类别间具有相同的方差(对于单变量分析)或相同的协方差矩阵(对于多变量分析)。

二次判别分析Quadratic discriminant analysisQDA),比LDA更灵活,假设预测变量服从(多元)正态分布,但允许各类别的协方差矩阵不等。

混合判别分析(Mixture discriminant analysisMDA),假定每个类别都是子类别的高斯(正态)混合。

弹性判别分析(Flexible Discriminant Analysis FDA),使用预测变量的非线性组合,例如样条曲线(splines)。

正则化判别分析(Regularized discriminant anlysisRDA),在预测变量数量大于训练集对象数量的情况下,正则化可改进协方差矩阵的估计。(注意其与冗余分析RDA)的缩写相同,不要混淆)

接下来,展示这些不同类型的判别分析在R中的实现方法。

 


测试数据集



数据集概要


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

这些鸢尾花来自三种物种,分别为setosan=50)、versicolorn=50)和virginican=50)。

包含四种性状,分别为萼片长度(sepal lengthcm)、萼片宽度(sepal widthcm)、花瓣长度(petal lengthcm)和花瓣宽度(petal widthcm)。

1.png

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

2.png

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

 

探索性分析


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

#变量分布
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')

3.png

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

 

相比直接选取部分趋势明显的变量作为代表,通过降维技术所确定的变量组合特征通常是种更为实用的选择,例如PCA

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

4.png

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

 

接下来接入本篇的方法部分,结合已知的物种分类,使用带监督模式的判别分析实现数据降维,选择代表性的变量组合特征,以及实现分类。

 


线性判别分析(LDA


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)

5.png

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

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

 

变量的标准化也是推荐的一步,对于消除变量间的量纲差异或者较大方差变量的权重时非常有效,可提高降维的精度。

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

 

为了更好地展示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,这里以MASS包中的方法为例。

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

library(MASS)
 
#拟合模型,详情 ?lda
model <- lda(Species~., data = train.data)
model

6.png

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)

7.png

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

有关ggplot2可视化的调整可见“排序图绘制”。

 

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

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

8.png

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

 

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

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

9.png

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

#后验概率也可通过热图展示
heatmap(predictions$posterior, Colv = NA, cexRow = 0.8, cexCol = 1)

10.png

 

实际应用中,还可以通过指定一个后验概率阈值,过滤模糊识别的结果,只保留较为可信的预测分类信息以提升优度等。

 


二次判别分析(QDA



从某种意义上说,QDA由于不假定各类别的协方差相等,因此它比LDA灵活得多。这种优势主要体现在训练集数据量较大、或者观测类别较多时的情况,此时LDA要求的等协方差矩阵的假设经常被拒绝。

当训练集数据较少时,LDA往往比QDA更好。

有关QDALDA的比较

 

MASS包也提供了QDA的方法。

#使用测试集拟合 QDA 模型,详情 ?qda
model <- qda(Species~., data = train.data)
model

11.png

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

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

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

 

随后,构建的分类器将对测试集对象预测分类。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions$class == train.data$Species)

12.png

结果显示,99%以上的训练集对象能够被分类到正确的类别中。

和上述LDA结果(使用相同的训练集数据,精度约~98%)相比,本示例也显示了QDA对训练集对象分类属性的识别有了提升,尽管甚微。

#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

13.png

对于测试集,约97%以的对象能够被分类到正确的类别中,表明QDA分类器的精度也是可靠的。

 


混合判别分析(MDA


无论LDAQDA均要求每个类别的数据服从多元正态(高斯)分布。

对于MDA中存在的多个类别,均被假定为更多子类别的高斯混合,其中每个对象都有属于每个类别的概率。但MDA仍然假设各类别的协方差矩阵相等。

 

MDA可使用mda包实现。

library(mda)
 
#使用测试集拟合模型,详情 ?mda
model <- mda(Species~., data = train.data, subclasses = 3)
model
 
plot(model)

13.png

例如在上式中,定义每个既定类别均由3个子类别组成,通过这些子类别进行拟合分类。Training Misclassification Error代表了对训练集对象预测分类错误率,示例结果显示错误率很低。

散点图展示了降维后的数据特征,MDA第一主轴承载了最能体现类间差异的特征。

 

对于模型预测功能的实现,和上文方法类似。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions == train.data$Species)    #结果等于 1-Training Misclassification Error
 
#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions == test.data$Species)

 

有些情况下,通过识别更小子类的MDA方法,可能会胜过全局识别模式LDAQDA

如下所示,存在3个主要的类别,每个类别又由3个子类别构成,各类别整体不满足正态(高斯)分布,但各子类别满足。由于子类别之间“不连续”,LDAQDA的决策边界无法有效分类,但MDA可以通过正确识别子类别而实现良好的分类。

14.png

 


弹性判别分析(FDA


FDALDA的灵活扩展,它使用样条曲线等预测变量的非线性组合。FDA可用于对每组内变量之间的多元非正态或非线性关系建模,从而实现更准确的分类。

 

FDA也可使用mda包实现。

#使用测试集拟合模型,详情 ?fda
model <- fda(Species~., data = train.data)
model
 
plot(model)

15.png

同样地,Training Misclassification Error代表了对训练集对象预测分类错误率,示例结果显示错误率很低。

散点图展示了降维后的数据特征,FDA第一主轴承载了最能体现类间差异的特征。

 

对于模型预测功能的实现,和上文方法类似。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions == train.data$Species)    #结果等于 1-Training Misclassification Error
 
#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions == test.data$Species)

 


正则化判别分析(RDA


RDA通过正则化组协方差矩阵建立分类规则,从而可以针对数据中的多重共线性提供更稳健的模型。这对于包含高度相关预测变量的大型多元数据集可能非常有用。

RDA可视为LDAQDA之间的一种中间状态。RDAQDA的独立协方差收缩为与LDA类似的等协方差。在预测变量的数量大于训练集对象数量的情况下,可以改善协方差矩阵的估计,潜在地提高了模型的准确性。

备注:正则化判别分析(RDA)与冗余分析RDA)的缩写相同,不要混淆

 

RDA可使用klaR包实现。

library(klaR)
 
#使用测试集拟合模型,详情 ?rda
model <- rda(Species~., data = train.data)
model

16.png

Regularization parametersgammalambda为两种正则化参数。

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

Misclassification rate,代表了对训练集对象预测分类错误率。

#对训练集预测分类
predictions <- predict(model, train.data)
mean(predictions$class == train.data$Species)
 
#对测试集预测分类
predictions <- predict(model, test.data)
mean(predictions$class == test.data$Species)

 


参考资料


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

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

 


更多精彩,可关注个人公众号“生信小白鱼”,感谢大家支持。




http://blog.sciencenet.cn/blog-3406804-1213689.html

上一篇:clusterProfiler的纯无参自定义物种注释的GO、KEGG富集分析及GSEA
下一篇:LEfSe分析的在线+本地运行教程

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-7-14 08:09

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部