非约束排序3 主成分分析（PCA）-基于Hellinger转化的

2021-05-23 本文已影响0人 fafu生信小蘑菇

第五章非约束排序3 主成分分析（PCA）-基于Hellinger转化的数据以及其他内容

由于前几天电脑被我捣鼓坏了，今天才修好，所以这几天就没有更新，加油。

上次我们讲了主成分分析的概念，如何提取、解读和绘制vegan包输出的PCA结果，如何将新变量和新对象投影到PCA双序图中以及如何组合聚类分析结果和排序结果。不知道你们都消化了吗？那么今天继续学习PCA的部分，今天的主要内容是：

转化后的物种数据PCA分析

PCA的应用

使用PCA.newr（）函数进行PCA分析

PCA中缺失值的估算

1.数据和包准备

setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
# 加载包，函数和数据
library(ade4)
library(vegan)
library(gclus)
library(ape)
library(missMDA)
library(FactoMineR)

#加载本章后面部分所需要的函数
source("cleanplot.pca.R")
source("PCA.newr.R")
source("CA.newr.R")

# 导入Doubs 数据
load("Doubs.RData")
# 删除无物种数据的样方8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]

2. 转化后的物种数据PCA分析

PCA是展示样方之间欧式距离的线性方法，理论上说并不适用物种多度数据的分析（特别是含有许多0值的情形，欧几里得距离对双零问题很敏感）。对于物种组成数据一般不使用PCA这种线性模型去做，而是使用PCoA或NMDS等方法。

如果仍要使用PCA的话，一般需要对数据进行转化，但是并不是说物种多度数据一定要作转化。如果物种多度数据比较均匀，且含有很少的0值（避开了欧几里得距离的“双零”问题），数据不转化也可以的。但是如果高丰度和低丰度物种种群的丰度差别比较大，数据并不“均匀”，可以在运行rda( )时，使用“scale=TRUE”消除高丰度物种与稀有物种在重要性方面的差异，降低欧几里得距离对高数值的敏感性

但是如果0值过多，物种多度数据还是需要进行数据转化，推荐使用Hellinger转化。

2.1 Hellinger转化的鱼类数据PCA分析

首先需要使用decostand（）函数对鱼类数据进行Hellinger转化，然后用rad（）进行PCA分析，当然你也可以试一试别的数据转化方法。

# 物种Hellinger 转化
spe.h <- decostand(spe, "hellinger")
(spe.h.pca <- rda(spe.h))

图2.1-1 Hellinger 转化后数据

图2.1-2 spe.h.pca

# 带断棍模型的碎石图
dev.new(
  title = "PCA特征值-物种数据的碎石图", 
  noRStudioGD = TRUE
)
screeplot(spe.h.pca,
          bstick = TRUE, 
          npcs = length(spe.h.pca$CA$eig)
)

图2.1-3 PCA特征值-物种数据的碎石图

# PCA 双序图
spe.pca.sc1 <- scores(spe.h.pca, display = "species", scaling = 1) #标尺1
spe.pca.sc2 <- scores(spe.h.pca, display = "species", scaling = 2) #标尺2

dev.new(title = "鱼类的主成分分析",
        width = 12,
        height = 6,
        noRStudioGD = TRUE
)
par(mfrow = c(1, 2))
cleanplot.pca(spe.h.pca, scaling = 1, mar.percent = 0.06)
cleanplot.pca(spe.h.pca, scaling = 2, mar.percent = 0.06)

图2.1-4 经过Hellinger转化后鱼类物种数据的PCA分析双序图

这里物种不像环境变量那样能够明显分组，但是同样可以看出物种沿着梯度更替分布。在1型标尺的双序图内，可以观察到有8个物种对于第一、二轴有很大贡献。

那么这些物种与聚类分析中聚类簇的指示种是否重合？我们可以重新运行一下没有转化的物种数据（文件名spe）的PCA，比较一下哪个排序图能更好地展示物种沿着河流路线的梯度分布情况？

PCA用于物理或化学属性变量排序分析，在群落数据分析方面的应用越来越广泛。

数据预转化并不会改变PCA线性排序的属性。物种数据预转化保证物种数据本身特异性，即并没有过分强调双零的重要性。

物种数据1型标尺PCA双序图揭示的是群落的潜在梯度结构，样方也是沿着这些潜在梯度（排序轴）依次排列，双序图平衡贡献圆可以识别物种变量对于排序轴的贡献程度。

2型标尺PCA双序图可以揭示物种之间的相关性，但如果是用转化后的物种数据进行排序，排序图上的物种之间相关性并不完全等同于原始数据的Pearson相关系数。

2.2 环境因子被动加入物种数据PCA分析

在排序过程直接加入环境变量的约束分析（即典范排序）经常被使用，但是如果想在物种数据简单排序后，再被动加入环境变量呢？这样我们可以总览物种数据的结构的同时，也可以推测这种结构与环境因子的关系。

通过vegan包里envfit（）函数实现，它也适用于CA、PCoA和NMDS。

“envfit 函数在排序图中确定环境变量向量或因子的平均值。[...]排序图上点到环境变量向量上的投影距离代表与该环境变量相关性，因子的点是因子水平的平均值”。

envfit（）函数的输出结果包含代表因子水平的点或定量环境变量的箭头的坐标，这些坐标可以被动投影到排序图。此外，用户可以通过p1ot（）函数只显示通过envfit（）函数计算的置换检验获得的相应的显著性水平的环境变量。

dev.new(
  title = "带有环境因子的PCA双序图-标尺2",
  noRStudioGD = TRUE
)
biplot(spe.h.pca, main = "鱼类数据PCA - 2型标尺")
(spe.h.pca.env <-
    envfit(spe.h.pca, env, scaling = 2)) 
# 默认2型标尺

图2.2-1

图2.2-2

# 显示显著环境变量，用户设定变量箭头颜色
plot(spe.h.pca.env, p.max = 0.05, col = 3)

#这行命令是在最后版的双序图中增加了显著环境变量
#注意：必须给envfit提供与绘图相同的标尺

图2.2-3

在约束排序中，envfit（）函数也可以用于解释变量与排序轴回归 $r^2$ 的显著性的置换检验。但envfit（）并不是检验解释变量对响应变量是否有显著效应最好的方法。

3. PCA的应用领域

PCA在生态学应用主要是用于基于定量环境变量的样方排序，以及经过适当转化后的群落组成数据的样方排序。PCA原始假设要求数据符合多元正态分布。但在生态学领域，只要偏离正态分布不要太离谱，PCA对于数据是否正态分布并不敏感。

PCA的主要计算步骤是基于离散矩阵（线性协方差或相关）的特征分解过程。虽然协方差或相关系数必须用定量变量才能算出，但二元（1-0）数据也可以用于PCA分析。

PCA使用条件：

PCA分析必须从相同量纲的变量表格开始。因为需要将变量总方差分配给特征根，所以变量必须有相同的物理单位，方差和才有意义（方差的单位是变量单位的平方），或者变量是无量纲的数据，例如标准化或对数转化后的数据。

矩阵数据不能倒置，因为对象（样方）之间的协方差或相关没有意义。

协方差或相关系数必须用定量变量才能算出。但由于半定量变量之间的Pearson相关系数等同于Spearman秩相关系数，所以基于半定量数据的PCA排序结果中变量之间关系是反映Spearman秩相关。

二元的数据也可以进行PCA分析。排序图对象之间的距离是1减去简单匹配系数 $s_1$ 后的平方根（即 $\sqrt{1-S_1}$ ）乘以变量数量的平方根。

物种有-无数据进行PCA分析前，可以用Hellinger转化或弦转化数据进行预处理。因为有-无数据的Hellinger距离或弦距离均等于 $\sqrt2\sqrt{1-Ochiai\ similarity}$ ，因此，有-无数据的Hellinger转化或弦转化后PCA分析的排序图（1型标尺）内样方之间的距离是Ochiai距离。 $\sqrt2\sqrt{1-Ochiai\ similarity}$ 是一种距离测度，这种测度也适用于有一无物种数据的分析。

应该避免用变量箭头顶点的距离来评估变量之间的相关性，而是用箭头之间的角度来判断相关性。

4.使用PCA.newr（）函数进行PCA分析

为了方便能够快速地用自己的数据进行PCA分析，本书作者编写了PCA.newr（）和biplot.PCA.newr（）两个函数（在自编函数PCA.newr.R里面）。下面用Doubs环境数据为例演示这两个函数的使用流程。

#PCA，这里默认是1型标尺双序图
# PCA; 这里默认是1型标尺双序图
env.PCA.PL <- PCA.newr(env, stand = TRUE)
dev.new(
  title = "PCA on environmental variables - scaling 1",
  noRStudioGD = TRUE
)
biplot.PCA.newr(env.PCA.PL)

# PCA;生成2型标尺双序图
dev.new(
  title = "PCA on environmental variables - scaling 2",
  noRStudioGD = TRUE
)
biplot.PCA.newr(env.PCA.PL, scaling = 2)

图4-1

图4-2

图4-3

这里主成分轴正负方向是随机的，可能与vegan包输出的排序图成镜像关系。但没有关系，因为对象或变量之间的相对位置没有变化。

5. PCA中缺失值的估算

如果有缺失值该怎么办?

是删除缺失值所在的行和列，把珍贵的同行或列数据丢弃不用？或者，为了避免这种情况，用一些方法找到有意义的替代值进行运算（统计语时中称为“插值"）。例如：利用相应变量的平均值或通过其他相关变量的国归分析获得的估计插值。当这两种解决方案出现时，我们不会考虑估计缺失值不确定性及其在进一步分析中的用途。结果，补充后的变量的标准误被低估，导致接下来的运算一些假设统计检验失效。

为了解决这些缺点，Josse和Huson（2012）提出了“一个正则化（regu-larized）的迭代PCA算法，为主坐标分析和主成分分析提供点估计”。这个算法与vegan包里的算法不一样，PCA是一轴一轴逐步迭代出来。当碰到缺失值的时候，会以该变量的平均值来替代。但不断迭代过程，缺失值会根据新的计算结果不断获得新的估计值，直到收敛为止，也就是让缺失值获得稳定的补充值。

实施复杂程序的目的在于：

①避免过度拟合估算模型（即当数据里面有比较多的缺失值的时候，有太多参数需要调整的问题）

②克服低估内插值方差的问题

注意：用推算得到的PCA的解决方案在不同PCA轴之间并不是嵌套关系，即基于s轴解决方案不等于与（s+1）或（s+2）解决方案。所以，做出所需的轴数适当的先验决定非常重要。该决定可能是经验性的（例如，使用具有特征值之和大于断棍模型预测的轴），或者在PCA重建过程，它可以通过数据本身的交又验证程序告知数据哪里删除，以及计算重建误差。我们的目的是保留最小化重建误差的轴数。

missMDA包，汇总了PCA分析所有可能通过选代方式插值缺失值的方法。这个差值的函数为imputePcA（）。接下来将以Doubs环境数据做两个实验，第一次只有3个缺失值（即所有319个数值的1%的缺失值），第二次是32个缺失值（10%）。

我们的第一个实验包括删除选择代表各种情况的三个值。

第一个值将靠近带相对对称分布的变量的平均值附近。pH平均值=8.04。让我们删除样方2中pH值（8.0）。

第二个将接近高度不对称分布的变量的均值，pho平均值=0.57。让我们删除样方19中的pho值（第18行，0.60）。

第三个将从不对称分布变量中删除，但删除后远离均值。bod平均值为5.01。让我们删除样方23的bod值（第22行，16.4，是第二高的值）。

看一下imputePcA（）如何执行迭代？

# PCA缺失值估计实验一:3个缺失值

# 用 NA 替换3个备选的值
env.miss3 <- env
env.miss3[2, 5] <- NA    # pH
env.miss3[18, 7] <- NA   # pho
env.miss3[22, 11] <- NA  # dbo

#三个变量的新均值（无缺失值）
mean(env.miss3[, 5], na.rm = TRUE)
mean(env.miss3[, 7], na.rm = TRUE)
mean(env.miss3[, 11], na.rm = TRUE)

图5-1

图5-2

图5-3

# 估算
env.imp <- imputePCA(env.miss3)

# 估算值

env.imp$completeObs[2, 5]    # 原始值: 8.0
env.imp$completeObs[18, 7]   # 原始值: 0.60
env.imp$completeObs[22, 11]  # 原始值 16.4

# 带估算值的PCA
env.imp3 <- env.imp$completeObs
env.imp3.pca <- rda(env.imp3, scale = TRUE)

# 原始PCA和带估算值的PCA的Procrustes比较
env.pca <- rda(env,scale=TRUE)
pca.proc <- procrustes(env.pca, env.imp3.pca, scaling = 1)

图5-4

图5-5

图5-6

dev.new(title = "Imputation in PCA",
        width = 14,
        height = 7,
        noRStudioGD = TRUE
)

plot(pca.proc, 
     main = "原始和带估计值PCA的procrustes比较\n3缺失值")
points(pca.proc, display = "target", col = "red")
text(
  pca.proc,
  display = "target",
  col = "red",
  pos = 4,
  cex = 0.6
)

图5-7

图5-7 PCA中缺失数据的估算。Procrustes比较原始环境变量PCA和对三个缺失值进行估算的PCA。1型标尺，这里只绘制样方图。红色是原始PCA的样方；蓝色是估算值的PCA中的样方。3个缺失值，即1%

上面这三个估算值，唯一比较好的估算是pH，它恰好落在原始值上。请记住，pH在数据中是对称分布的。另外两个估算远离原始值，这两个变量都有强不对称分布。显然，这比缺失值接近或不接近变量的平均值更重要（pho属于这种情况，但bod不是这种情况）。

那么如何影响排序结果（前2轴）？

为了观察这个结果，在带有估算值矩阵上计算新的PCA，并将其与原始PCA进行比较。通过Procrustes旋转新环境变量PCA轴，使其与原始PCA一致。Procrustes旋转是最小化两个排序结构同一标号的对象图上的距离平方和，目的在于获得最多的重叠区域。Procrustes旋转可以通过vegan包的函数procrustes（）来执行。

显然，实际值和重建值之间的差异对排序的影响很小。唯一可见的区别是在样方23，其中实际值与重建值之间的差异最大。

#缺失值估计实验二:32个缺失值
#从319个值中随机抽取32个用NA替代
rnd <- matrix(sample(c(rep(1, 32), rep(0, 287))), 29, 11)
env.miss32 <- env
env.miss32[rnd == 1] <- NA
# 每个样方里面多少NA?
summary(t(env.miss32))
# 显示NA数量的替代方法：
# sapply(as.data.frame(t(env.miss32)), function(x) sum(is.na(x)))

# 内插估计
env.imp2 <- imputePCA(env.miss32)

# 带估计值的PCA
env.imp32 <- env.imp2$completeObs
env.imp32.pca <- rda(env.imp32, scale = TRUE)

# 原始PCA和带估计值的PCA的procrustes比较
pca.proc32 <- procrustes(env.pca, env.imp32.pca, scaling = 1)

图5-8

图5-9

dev.new(title = "Imputation in PCA",
        width = 14,
        height = 7,
        noRStudioGD = TRUE
)

plot(pca.proc32, 
     main = "原始和带估计值PCA的procrustes比较\n32缺失值")
points(pca.proc32, display = "target", col = "red")
text(
  pca.proc32,
  display = "target",
  col = "red",
  pos = 4,
  cex = 0.6
)