数量生态学笔记||线性判别式分析（LDA）
2018-10-14 本文已影响296人
周运来就是我
线性判别式分析（linear discriminant analysis , LDA）跟PCA非常相似、唯一不同的是LDA的结果是将数据投影到不同分类、PCA的结果是将数据投影到最高相似分组，而且过程无一例外的都基于特征值与特性向量实现降维处理。PCA变换基于在原数据与调整之后估算降维的数据之间最小均方错误，PCA趋向提取数据最大相同特征、而忽视数据之间微小不同特征、所以如果在OCR识别上使用PCA的方法就很难分辨Q与O个英文字母、而LDA基于最大类间方差与最小类内方差，其目的是减少分类内部之间差异，扩大不同分类之间差异。所以LDA在一些应用场景中有比PCA更好的表现。
LDA与RDA和CCA的不同之处是其响应变量是样方的分组情况。LDA的目的是确定一套独立的定量解释变量能够能够多大程度上解释当前样方的分组结果。LDA需要从标准化的解释变量计算判别函数。判别函数值用于量化的解释变量对于对象分组的相对贡献。另外，通过原始的解释变量计算的判别函数称为识别函数，目的是判断新的对象应该归于哪一组。
运行LDA是必须保证解释变量组内协方差矩阵齐性，但生态数据通常无法满足这个条件。
核心：向最大化类间差异、最小化类内差异的方向线性投影
线性鉴别分析的基本思想是通过线性投影来最小化同类样本间的差异，最大化不同类样本间的差异。具体做法是寻找一个向低维空间的投影矩阵W，样本的特征向量x经过投影之后得到的新向量：
y = Wx
同一类样投影后的结果向量差异尽可能小，不同类的样本差异尽可能大。直观来看，就是经过这个投影之后同一类的样本进来聚集在一起，不同类的样本尽可能离得远。下图是这种投影的示意图：
上图中特征向量是二维的，我们向一维空间即直线投影，投影后这些点位于直线上。在上面的图中有两类样本，通过向右上方的直线投影，两类样本被有效的分开了。绿色的样本投影之后位于直线的下半部分，红色的样本投影之后位于直线的上半部分。
训练时的优化目标是类间差异与类内差异的比值：
最后归结于求解矩阵的特征值与特征向量：
LDA是有监督的机器学习算法，在计算过程中利用了样本标签值。这是一种判别模型，也是线性模型。LDA也不能直接用于分类和回归问题，要对降维后的向量进行分类还需要借助其他算法，如kNN。
library(ade4)
library(vegan)
#library(packfor)
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\数量生态学\\DATA")
# 此程序包可以从 https://r-forge.r-project.org/R/?group_id=195 下载
# 如果是MacOS X系统，packfor程序包内forward.sel函数的运行需要加载
# gfortran程序包。用户必须从"cran.r-project.org"网站内选择"MacOS X"，
# 然后选择"tools"安装gfortran程序包。
library(MASS)
library(ellipse)
library(FactoMineR)
# 附加函数
source("evplot.R")
source("hcoplot.R")
# 导入CSV数据文件
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)
# 删除没有数据的样方8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]
# 提取环境变量das（离源头距离）以备用
das <- env[, 1]
# 从环境变量矩阵剔除das变量
env <- env[, -1]
spe.hel <- decostand(spe, "hellinger")
> # 线性判别式分析（LDA）
> # *******************
> # 基于Hellinger转化鱼类数据的样方Ward聚类分析（分4组）
> gr <- cutree(hclust(vegdist(spe.hel, "euc"), "ward.D"), 4) #如果之前没有保存gr
> # 线性判别式分析（LDA）
> # *******************
> # 基于Hellinger转化鱼类数据的样方Ward聚类分析（分4组）
> gr <- cutree(hclust(vegdist(spe.hel, "euc"), "ward.D"), 4) #如果之前没有保存gr
> env.pars2 <- as.matrix(env[, c(1, 9, 10)])
> # 使用函数betadisper（）（vegan包）（Marti Anderson检验）验证解释变量组内# 协方差矩阵的齐性
> env.pars2.d1 <- dist(env.pars2)
> (env.MHV <- betadisper(env.pars2.d1, gr))
Homogeneity of multivariate dispersions
Call: betadisper(d = env.pars2.d1, group = gr)
No. of Positive Eigenvalues: 3
No. of Negative Eigenvalues: 0
Average distance to median:
1 2 3 4
197.458 215.135 32.645 5.845
Eigenvalues for PCoA axes:
PCoA1 PCoA2 PCoA3 <NA> <NA> <NA> <NA> <NA>
2036227.3390 442.1009 31.4911 NA NA NA NA NA
> anova(env.MHV)
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 3 226555 75518 4.9758 0.007641 **
Residuals 25 379427 15177
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> permutest(env.MHV) #置换检验
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 226555 75518 4.9758 999 0.008 **
Residuals 25 379427 15177
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> #组内协方差矩阵不齐，需要对alt和dbo两个变量进行对数化
> env.pars3 <- cbind(log(env$alt), env$oxy, log(env$dbo))
> colnames(env.pars3) <- c("alt.ln", "oxy", "dbo.ln")
> row.names(env.pars3) <- row.names(env)
> env.pars3.d1 <- dist(env.pars3)
> (env.MHV2 <- betadisper(env.pars3.d1, gr))
Homogeneity of multivariate dispersions
Call: betadisper(d = env.pars3.d1, group = gr)
No. of Positive Eigenvalues: 3
No. of Negative Eigenvalues: 0
Average distance to median:
1 2 3 4
0.9153 1.1817 1.0133 0.7599
Eigenvalues for PCoA axes:
PCoA1 PCoA2 PCoA3 <NA> <NA> <NA> <NA> <NA>
146.6682 6.6185 2.4819 NA NA NA NA NA
> permutest(env.MHV2)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 3 0.5337 0.17790 0.2795 999 0.836
Residuals 25 15.9136 0.63654
> #组内协方差矩阵显示齐性，可以继续分析
> # LDA计算（判别函数）
> env.pars3.df <- as.data.frame(env.pars3)
> gr
1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
1 1 1 2 2 2 1 2 1 1 1 1 1 2 2 2 2 2 3 3 3 4 4 4 3 3 3 3 3
> head(env.pars3.df)
alt.ln oxy dbo.ln
1 6.839476 12.2 0.9932518
2 6.837333 10.3 0.6418539
3 6.817831 10.5 1.2527630
4 6.749931 11.0 0.2623643
5 6.744059 8.0 1.8245493
6 6.740519 10.2 1.6677068
> (spe.lda <- lda(gr ~ alt.ln + oxy + dbo.ln, data=env.pars3.df))
Call:
lda(gr ~ alt.ln + oxy + dbo.ln, data = env.pars3.df)
Prior probabilities of groups:
1 2 3 4
0.3103448 0.3103448 0.2758621 0.1034483
Group means:
alt.ln oxy dbo.ln
1 6.464881 11.388889 1.048586
2 6.245151 9.944444 1.209785
3 5.385688 8.387500 1.557074
4 5.477515 5.200000 2.707430
Coefficients of linear discriminants:
LD1 LD2 LD3
alt.ln -2.5129208 -1.7224562 -1.0771922
oxy -0.7730690 -0.2117976 0.8197343
dbo.ln -0.3348384 -2.7147468 2.3012850
Proportion of trace:
LD1 LD2 LD3
0.9032 0.0876 0.0092
> # 此结果对象内包含大量解读LDA的必要信息
> summary(spe.lda)
Length Class Mode
prior 4 -none- numeric
counts 4 -none- numeric
means 12 -none- numeric
scaling 9 -none- numeric
lev 4 -none- character
svd 3 -none- numeric
N 1 -none- numeric
call 3 -none- call
terms 3 terms call
xlevels 0 -none- list
> # 显示三个变量分组平均值
> spe.lda$means
alt.ln oxy dbo.ln
1 6.464881 11.388889 1.048586
2 6.245151 9.944444 1.209785
3 5.385688 8.387500 1.557074
4 5.477515 5.200000 2.707430
> #计算范数标准化特征向量（matrix C，eq.11.33），即标准化判别函数系数
> (Cs <- spe.lda$scaling)
LD1 LD2 LD3
alt.ln -2.5129208 -1.7224562 -1.0771922
oxy -0.7730690 -0.2117976 0.8197343
dbo.ln -0.3348384 -2.7147468 2.3012850
> # 计算典范特征根
> spe.lda$svd^2
[1] 53.6898562 5.2088449 0.5478896
> # 计算在典范变量空间内样方的坐标
> (Fp <- predict(spe.lda)$x)
LD1 LD2 LD3
1 -4.08638573 -0.89640510 0.368028660
2 -2.49450633 0.46365901 -1.995824032
3 -2.80466835 -1.20357224 -0.404993616
4 -2.68895361 1.49616437 -2.201175466
5 -0.87807014 -2.09926524 -1.059020040
6 -2.51740981 -2.13333526 0.387269196
7 -2.90386963 0.07319673 -0.891988271
9 0.10415477 -1.24335518 -1.988893958
10 -1.49957974 -0.97965055 0.082158618
11 -1.88806718 0.38774388 0.504579631
12 -2.43308232 -0.02501061 1.334323268
13 -2.36655395 0.63877375 1.047520033
14 -2.35214071 -0.52520225 2.062059098
15 -1.57722535 1.28900168 0.253631503
16 -0.32438764 1.07783858 -0.206474231
17 -0.23770978 -0.21870301 1.018179031
18 -0.03051363 1.18888940 0.008410463
19 -0.14515674 0.79740513 0.706294065
20 0.34427132 1.44578197 0.169066314
21 1.44181530 0.83677120 0.075460197
22 1.38965445 0.44108254 0.553586732
23 3.22326372 -2.24627627 1.120313152
24 4.22156847 -1.19694492 -0.421313261
25 5.07604329 -1.72116640 -0.573615673
26 3.85542329 -0.32572787 -0.218162145
27 3.29378337 0.46604929 -0.152484763
28 2.84858565 1.28339090 -0.129929770
29 2.33552915 1.38947026 0.517474958
30 3.09418785 1.53939626 0.035520310
> # 替代方式: Fp <- scale(env.pars3.df, center=TRUE, scale=FALSE) %*% C
> # 对象的分类
> (spe.class <- predict(spe.lda)$class)
[1] 1 2 1 2 2 1 1 2 2 1 1 1 1 2 2 2 2 2 2 3 3 4 4 4 3 3 3 3 3
Levels: 1 2 3 4
> # 对象属于分类组的后验概率
> (spe.post <- predict(spe.lda)$posterior)
1 2 3 4
1 9.858670e-01 1.413296e-02 8.589109e-10 1.420163e-15
2 4.940452e-01 5.059487e-01 6.085899e-06 6.418503e-12
3 8.588135e-01 1.411862e-01 2.982838e-07 1.233446e-11
4 4.792624e-01 5.207302e-01 7.439251e-06 3.764956e-13
5 1.881617e-01 8.115272e-01 3.063205e-04 4.798814e-06
6 8.836761e-01 1.163235e-01 3.978965e-07 2.998536e-10
7 7.973581e-01 2.026411e-01 7.545732e-07 9.620716e-13
9 2.171690e-02 9.658988e-01 1.224973e-02 1.345408e-04
10 4.808279e-01 5.190573e-01 1.147554e-04 2.862482e-08
11 6.122120e-01 3.876919e-01 9.607204e-05 3.223522e-10
12 8.718114e-01 1.281831e-01 5.490190e-06 1.869581e-11
13 8.146431e-01 1.853415e-01 1.540601e-05 1.071606e-11
14 9.112721e-01 8.872342e-02 4.439650e-06 6.168551e-11
15 3.980760e-01 6.011260e-01 7.980082e-04 4.615970e-10
16 6.299642e-02 8.902208e-01 4.678228e-02 4.895917e-07
17 1.417366e-01 8.258918e-01 3.236084e-02 1.077488e-05
18 3.993984e-02 8.237538e-01 1.363046e-01 1.813186e-06
19 7.959866e-02 8.283919e-01 9.200680e-02 2.625695e-06
20 1.533174e-02 5.625259e-01 4.221367e-01 5.688134e-06
21 3.189551e-04 6.254116e-02 9.366679e-01 4.720092e-04
22 6.638353e-04 8.237900e-02 9.158176e-01 1.139555e-03
23 5.212073e-08 4.643126e-05 2.492989e-02 9.750236e-01
24 2.218889e-10 3.122087e-06 5.333234e-02 9.466645e-01
25 3.359000e-13 1.801844e-08 2.924798e-03 9.970752e-01
26 4.997460e-09 4.293353e-05 5.209165e-01 4.790406e-01
27 5.899119e-08 2.395345e-04 9.602048e-01 3.955563e-02
28 2.153179e-07 5.145082e-04 9.973440e-01 2.141263e-03
29 2.518020e-06 1.857210e-03 9.975367e-01 6.035869e-04
30 5.013311e-08 1.724659e-04 9.981144e-01 1.713044e-03
> # 预设分类和预测分类的比较表格
> (spe.table <- table(gr, spe.class))
spe.class
gr 1 2 3 4
1 7 2 0 0
2 1 8 0 0
3 0 1 7 0
4 0 0 0 3
> # 正确分类的比例
> diag(prop.table(spe.table, 1))
1 2 3 4
0.7777778 0.8888889 0.8750000 1.0000000
> # 绘制典范变量空间内样方位置图，样方颜色代表不同的组
> plot(Fp[, 1], Fp[, 2], type="n")
> text(Fp[, 1], Fp[, 2], row.names(env), col=c(as.numeric(spe.class)+1))
> abline(v=0, lty="dotted")
> abline(h=0, lty="dotted")
> # 绘制95%的椭圆图
> for(i in 1:length(levels(as.factor(gr)))){
+ cov <- cov(Fp[gr==i, ])
+ centre <- apply(Fp[gr==i, ], 2, mean)
+ lines(ellipse(cov, centre= centre, level=0.95))
+ }
> # 新样方的分类（识别）
> # 生成一个新的样方（ln(alt)=6.8, oxygen=90 和 ln(dbo)=3.2）
> newo <- c(6.8, 90, 3.2)
> newo <- as.data.frame(t(newo)) # 必须在同一行
> colnames(newo) <- colnames(env.pars3)
> (predict.new <- predict(spe.lda, newdata=newo))
$`class`
[1] 1
Levels: 1 2 3 4
$posterior
1 2 3 4
1 1 7.578586e-65 8.588537e-153 6.08234e-184
$x
LD1 LD2 LD3
1 -64.87086 -23.29703 69.26423
> # 基于刀切法（jackknife）分类的LDA（即留一法交叉验证)
> spe.lda.jac <- lda(gr ~ alt.ln + oxy + dbo.ln, data=env.pars3.df, CV=TRUE)
> summary(spe.lda.jac)
Length Class Mode
class 29 factor numeric
posterior 116 -none- numeric
terms 3 terms call
call 4 -none- call
xlevels 0 -none- list
> # 正确分类的数量和比例
> spe.jac.class <- spe.lda.jac$class
> spe.jac.table <- table(gr, spe.jac.class)
> diag(prop.table(spe.jac.table, 1))
1 2 3 4
0.7777778 0.6666667 0.7500000 1.0000000
参考
线性判别分析（Linear Discriminant Analysis，LDA）
lda.pdf
Linear Discriminant Analysis in R: An Introduction
lda
2014_python_lda
用一句话总结常用的机器学习算法
PCA与LDA比较
