R统计绘图-判别分析(LDA/QDA)

R统计绘图-判别分析(LDA/QDA) by EcoEvoPhylo

zhe

点击名片关注我们

判别分析是一种有监督的分类模型。当数据中样本的分类很明确时，判别分析比逻辑斯蒂回归更有效，泛化能力更强。但是不适合自变量与因变量之间具有错综复杂关系的数据集。判别分析包括线性判别分析和二次判别分析。判别分析使用贝叶斯定理确定每个样本/观测属于某个类别的概率，然后选择高概率的类别作为该观测的类别。

一、数据准备

数据是模拟数据，包含3个分类变量，8个定量变量。

# 设置工作路径#knitr::opts_knit$set(root.dir="D:\\EnvStat\\EcoEvoPhylo\\biomarker")# 使用Rmarkdown进行程序运行Sys.setlocale('LC_ALL','C') # Rmarkdown全局设置#options(stringsAsFactors=F)# R中环境变量设置，防止字符型变量转换为因子
library(tidyverse)library(magrittr)# 1.1 数据导入data <- read.csv("data.csv",                 header = TRUE,                 row.names = 1,                 stringsAsFactors = TRUE                 )data %>% head()
# 1.2 数据描述统计data %>%  summary

图1|数据表，data.csv。

二、线性判别分析

线性判别分析假设每种类别中的观测服从多元正态分布，且不同类别之间具有相同的协方差。

# 2.1 基础LDA分析library(MASS)lda.fit0 <- lda(depth ~ .,data = data[-c(2:3)])lda.fit0
# 2.2 提取分析结果## 2.2.1 结果描述lda.fit0$means # 各变量的组内均值lda.fit0$prior # 各组先验概率lda.fit0$counts # 各组样本计数lda.fit0$scaling # 线性判别系数,线性判别输出的系数用于形成LDA判决评分的自变量的线性组合。lda.fit0$svd # 给出线性判别变量的组间和组内标准差之比。它们的平方是典型的F-统计量。
## 2.2.2 提取LDA线性判别系数lda.var <- data.frame(lda.fit0$scaling)write.csv(lda.var,"lda.var.csv",quote = FALSE)lda.var %>% head
## 2.2.3 提取LDA样本判别评分lda.ind <- data.frame(  predict(lda.fit0)$x, # 判别评分  depth = data$depth, # 先验分类  class = predict(lda.fit0)$class # 预测分类                      )write.csv(lda.ind,"lda.ind.csv",quote = FALSE)lda.ind %>% head
## 2.2.4 LDA轴特征值lda.eig <- round(lda.fit0$svd^2/sum(lda.fit0$svd^2)*100,2)lda.eig
## 2.2.5 观测属于某个类别的概率predict(lda.fit0)$posterior
## 2.2.6 预测的观测所属类别predict(lda.fit0)$class
## 2.2.7 判别正确率-混淆矩阵caret::confusionMatrix(data$depth, predict(lda.fit0)$class)

图2|LDA线性判别系数，lda.var.csv。LDA分析只能形成(因变量类别数-1)个成分轴，这里有A，B，C三个类别，因此形成了2个LDA轴。线性判别系数是标准线性组合，用于确定观测的判别评分的特征。

图3|LDA样本判别评分，lda.ind.csv。depth是样本实际类别，class是模型预测分类。

图4|混淆矩阵。包含模型预测准确率等信息。

# 2.3 LDA分析leave-one-out交叉验证lda.fit1 <- lda(depth ~ .,data = data[-c(2:3)],               method = "moment",               #先验概率：默认各分类占总样本数的比例#               prior = c(1,1,1)/3, # 设置顺序必须与因子分类水平一致。               #某变量组内方差<tol^2时，将会停止运行#               tol = 1.0e-4,               #缺失值处理方式#               na.action = na.omit,               #执行leave-one-out交叉验证#               CV = TRUE,               )
lda.fit1$class # 交叉验证结果，每个样本对应的类别。table(lda.fit1$class,data$depth)mean(lda.fit1$class == data$depth) # 模型分类正确率。lda.fit1$posterior # 样本归于各分类的后验概率。

# 2.4 LDA分析结果可视化## 2.4.1 样本判别评分图library(ggplot2)lda.ind %>% ggplot(aes(LD1,LD2,color = depth)) +  geom_point(size=2.5) +  stat_ellipse() +  theme_bw() +  labs(x = paste("LDA1 (",lda.eig[1],"%)",sep = ""),       y = paste("LDA2 (",lda.eig[2],"%)",sep = "")) -> p.ind
ggsave("lda.ind.pdf",p.ind,family = "Times")

图5|LDA样本判别评分图。LDA1能很好的代表depth的组间差异。

## 2.4.2 线性判别系数图lda.var %>% ggplot() +  geom_segment(    aes(      x = 0,xend = LD1,      y = 0,yend = LD2),    linewidth = 1,    color = "blue",    arrow = arrow(      type = "closed",      length = unit(0.2, "cm")),    lineend = "butt",    arrow.fill = "white") +  ggrepel::geom_text_repel(    aes(      LD1,LD2,      label = rownames(lda.var)),    color = "black")+  ggthemes::theme_base() +  geom_hline(yintercept = 0) +  geom_vline(xintercept = 0) +  labs(x = paste("LDA1 (",lda.eig[1],"%)",                 sep = ""),       y = paste("LDA2 (",lda.eig[2],"%)",                 sep = "")) -> p.var
ggsave("lda.var.pdf",p.var,family = "Times")

图6|LDA线性判别系数图。

## 2.4.3 样本-变量双序图### ggord包可快速绘制LDA双序图，也可直接使用ggplot2绘制。 # devtools::install_github('fawda123/ggord')library(ggord)ggord(  #自行提供分组信息、样本得分与变量得分数据进行绘图#  #此方法不能自行计算轴特征值，后续也无法修改轴标题#  grp_in = data$depth,  obs = lda.ind,  #为了图绘制美观，变量系数全部乘以一个常数。#  vecs = lda.var*3,    #直接提供排序对象进行绘图，此方法会直接计算轴特征值，#  # ord_in = lda.fit0,    xlims = c(-5.5,5.5),  ylims = c(-3,3),  size = 2.5, # 点大小。  veclsz = 0.8, # 变量系数箭头粗细。  repel = FALSE) +  geom_hline(yintercept = 0) +  geom_vline(xintercept = 0) +  theme(  plot.background = element_blank()  )  -> p.biplot  ggsave("lda.biplot.pdf",p.biplot,family = "Times")

图7|LDA样本-变量双序图。此数据中V1的LDA系数较大，其他变量的系数都较小，等比例乘以常数之后，多数变量箭头还是很短。绘图时可根据自己数据，挑选变量，选择性绘图。

## 2.4.4 lefse分析中的LDA得分### 获取第一特征向量w <- lda.fit0$scaling[,1]w
### 标准化第一特征向量w_unit <- w/sqrt(sum(w^2))w_unit
### 计算观测新坐标xy_matrix <- as.matrix(data[-c(1:3)])LD <- xy_matrix %*% w_unit
### 计算LDA得分###计算组间距离，作为效应系数###如果进行了自举检验，则计算重复的均值；###类别超过2组时，两两计算LDA得分，然后选择最大值。pairs <- combn(unique(as.character(data$depth)), 2)pairs res_lda_pair <- list()for (p in seq(ncol(pairs))) {effect_size <- abs(  mean(LD[data[,"depth"] == pairs[,p][1]]) - mean(LD[data[,"depth"] == pairs[,p][2]]));coeff <- abs(w_unit * effect_size); # 标准化后特征向量*效应系数(组间距离) rres <- lda.fit0$means %>% as.data.frame;rres1 <- apply(rres[rownames(rres) %in% pairs[,p],], 2, function(x) abs(x[1] - x[2]));res_lda_pair[[p]] <- (rres1 + coeff[names(rres1)]) * 0.5;                }res_lda_pair res <- sapply(colnames(data[-c(1:3)]), function(k){  lapply(seq(ncol(pairs)), function(p){    res_lda_pair[[p]][k]   }) %>% unlist %>% max})res <- sapply(res, function(x) {log10(1 + abs(x)) * ifelse(x > 0, 1, -1)})res
### 检索到的另一种计算多分类lefse中LDA得分的简单方法data.frame(t(lda.fit0$means)) %>%  mutate(name = rownames(.)) %>%   rowwise() %>%   mutate(    max = max(A, B, C),    min = min(A, B, C),    ) %>%  mutate(LDA = log10(1+(max-min)/2)) %>% # 趋势是一致的。  ungroup() -> lefse.ldalefse.lda$class <- NA # 各自变量在哪个类别中富集。for(i in 1:nrow(lefse.lda))  {lefse.lda[i, "class"] <- names(which.max(lefse.lda[i, c("A", "B","C")]))}lefse.lda                       lefse.lda$res_lda <- reswrite.csv(lefse.lda,"lefse.lda.csv",quote = FALSE)lefse.lda

图8|lefse分析中的LDA得分，lefse.lda.csv。两种计算方法的结果有差异。计算过程参考：https://rdrr.io/cran/microeco/src/R/trans_diff.R和https://zhuanlan.zhihu.com/p/265681612。

三、lefse分析结果可视化

二次判别分析(QDA)假设每个类别的观测服从多元正态分布，但允许每个类别具有各自的协方差，意味着判别分数的计算允许引入二次项。

# 3.1 基础QDA分析library(MASS)qda.fit0 <- qda(depth ~ .,data = data[-c(2:3)])qda.fit0
# 3.2 结果提取## 3.2.1 结果描述-因为是二次函数，结果中没有线性判别系数。qda.fit0$means # 各变量的组内均值qda.fit0$prior # 各组先验概率
## 3.2.2 分类正确性caret::confusionMatrix(data$depth,predict(qda.fit0)$class) # 预测正确率只有些微提升。

图9|QDA混淆矩阵。QDA分类器使用预测数据的二次函数，而不是线性函数，因此输出结果没有线性判别系数。

# 3.3 QDA分析leave-one-out交叉验证qda.fit1 <- qda(depth ~ .,data = data[-c(2:3)],               method = "moment",               #先验概率：默认各分类占总样本数的比例#               prior = c(1,1,1)/3, # 设置顺序必须与因子分类水平一致。               #某变量组内方差<tol^2时，将会停止运行#               tol = 1.0e-4,               #缺失值处理方式#               na.action = na.omit,               #执行leave-one-out交叉验证#               CV = TRUE,               )
qda.fit1$class # 交叉验证结果，每个样本对应的类别。table(qda.fit1$class,data$depth)mean(qda.fit1$class == data$depth) # 模型分类正确率。qda.fit1$posterior # 样本归于各分类的后验概率。