第五章非約束排序3 主成分分析（PCA）-基于Hellinger轉(zhuǎn)化的數(shù)據(jù)以及其他內(nèi)容

由于前幾天電腦被我搗鼓壞了，今天才修好，所以這幾天就沒有更新，加油。

上次我們講了主成分分析的概念，如何提取、解讀和繪制vegan包輸出的PCA結(jié)果，如何將新變量和新對(duì)象投影到PCA雙序圖中以及如何組合聚類分析結(jié)果和排序結(jié)果。不知道你們都消化了嗎？那么今天繼續(xù)學(xué)習(xí)PCA的部分，今天的主要內(nèi)容是：

轉(zhuǎn)化后的物種數(shù)據(jù)PCA分析

PCA的應(yīng)用

使用PCA.newr（）函數(shù)進(jìn)行PCA分析

PCA中缺失值的估算

1.數(shù)據(jù)和包準(zhǔn)備

setwd("E:/Rstudio_working/lai jian shan/DATA/Data")
# 加載包，函數(shù)和數(shù)據(jù)
library(ade4)
library(vegan)
library(gclus)
library(ape)
library(missMDA)
library(FactoMineR)

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

# 導(dǎo)入Doubs 數(shù)據(jù)
load("Doubs.RData")
# 刪除無物種數(shù)據(jù)的樣方8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]

2. 轉(zhuǎn)化后的物種數(shù)據(jù)PCA分析

PCA是展示樣方之間歐式距離的線性方法，理論上說并不適用物種多度數(shù)據(jù)的分析（特別是含有許多0值的情形，歐幾里得距離對(duì)雙零問題很敏感）。對(duì)于物種組成數(shù)據(jù)一般不使用PCA這種線性模型去做，而是使用PCoA或NMDS等方法。

如果仍要使用PCA的話，一般需要對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)化，但是并不是說物種多度數(shù)據(jù)一定要作轉(zhuǎn)化。如果物種多度數(shù)據(jù)比較均勻，且含有很少的0值（避開了歐幾里得距離的“雙零”問題），數(shù)據(jù)不轉(zhuǎn)化也可以的。但是如果高豐度和低豐度物種種群的豐度差別比較大，數(shù)據(jù)并不“均勻”，可以在運(yùn)行rda( )時(shí)，使用“scale=TRUE”消除高豐度物種與稀有物種在重要性方面的差異，降低歐幾里得距離對(duì)高數(shù)值的敏感性

但是如果0值過多，物種多度數(shù)據(jù)還是需要進(jìn)行數(shù)據(jù)轉(zhuǎn)化，推薦使用Hellinger轉(zhuǎn)化。

2.1 Hellinger轉(zhuǎn)化的魚類數(shù)據(jù)PCA分析

首先需要使用decostand（）函數(shù)對(duì)魚類數(shù)據(jù)進(jìn)行Hellinger轉(zhuǎn)化，然后用rad（）進(jìn)行PCA分析，當(dāng)然你也可以試一試別的數(shù)據(jù)轉(zhuǎn)化方法。

# 物種Hellinger 轉(zhuǎn)化
spe.h <- decostand(spe, "hellinger")
(spe.h.pca <- rda(spe.h))

圖2.1-1 Hellinger 轉(zhuǎn)化后數(shù)據(jù)

圖2.1-2 spe.h.pca

# 帶斷棍模型的碎石圖
dev.new(
  title = "PCA特征值-物種數(shù)據(jù)的碎石圖", 
  noRStudioGD = TRUE
)
screeplot(spe.h.pca,
          bstick = TRUE, 
          npcs = length(spe.h.pca$CA$eig)
)

圖2.1-3 PCA特征值-物種數(shù)據(jù)的碎石圖

# PCA 雙序圖
spe.pca.sc1 <- scores(spe.h.pca, display = "species", scaling = 1) #標(biāo)尺1
spe.pca.sc2 <- scores(spe.h.pca, display = "species", scaling = 2) #標(biāo)尺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 經(jīng)過Hellinger轉(zhuǎn)化后魚類物種數(shù)據(jù)的PCA分析雙序圖

這里物種不像環(huán)境變量那樣能夠明顯分組，但是同樣可以看出物種沿著梯度更替分布。在1型標(biāo)尺的雙序圖內(nèi)，可以觀察到有8個(gè)物種對(duì)于第一、二軸有很大貢獻(xiàn)。

那么這些物種與聚類分析中聚類簇的指示種是否重合？我們可以重新運(yùn)行一下沒有轉(zhuǎn)化的物種數(shù)據(jù)（文件名spe）的PCA，比較一下哪個(gè)排序圖能更好地展示物種沿著河流路線的梯度分布情況？

PCA用于物理或化學(xué)屬性變量排序分析，在群落數(shù)據(jù)分析方面的應(yīng)用越來越廣泛。

數(shù)據(jù)預(yù)轉(zhuǎn)化并不會(huì)改變PCA線性排序的屬性。物種數(shù)據(jù)預(yù)轉(zhuǎn)化保證物種數(shù)據(jù)本身特異性，即并沒有過分強(qiáng)調(diào)雙零的重要性。

物種數(shù)據(jù)1型標(biāo)尺PCA雙序圖揭示的是群落的潛在梯度結(jié)構(gòu)，樣方也是沿著這些潛在梯度（排序軸）依次排列，雙序圖平衡貢獻(xiàn)圓可以識(shí)別物種變量對(duì)于排序軸的貢獻(xiàn)程度。

2型標(biāo)尺PCA雙序圖可以揭示物種之間的相關(guān)性，但如果是用轉(zhuǎn)化后的物種數(shù)據(jù)進(jìn)行排序，排序圖上的物種之間相關(guān)性并不完全等同于原始數(shù)據(jù)的Pearson相關(guān)系數(shù)。

2.2 環(huán)境因子被動(dòng)加入物種數(shù)據(jù)PCA分析

在排序過程直接加入環(huán)境變量的約束分析（即典范排序）經(jīng)常被使用，但是如果想在物種數(shù)據(jù)簡(jiǎn)單排序后，再被動(dòng)加入環(huán)境變量呢？這樣我們可以總覽物種數(shù)據(jù)的結(jié)構(gòu)的同時(shí)，也可以推測(cè)這種結(jié)構(gòu)與環(huán)境因子的關(guān)系。

通過vegan包里envfit（）函數(shù)實(shí)現(xiàn)，它也適用于CA、PCoA和NMDS。

“envfit 函數(shù)在排序圖中確定環(huán)境變量向量或因子的平均值。[...]排序圖上點(diǎn)到環(huán)境變量向量上的投影距離代表與該環(huán)境變量相關(guān)性，因子的點(diǎn)是因子水平的平均值”。

envfit（）函數(shù)的輸出結(jié)果包含代表因子水平的點(diǎn)或定量環(huán)境變量的箭頭的坐標(biāo)，這些坐標(biāo)可以被動(dòng)投影到排序圖。此外，用戶可以通過p1ot（）函數(shù)只顯示通過envfit（）函數(shù)計(jì)算的置換檢驗(yàn)獲得的相應(yīng)的顯著性水平的環(huán)境變量。

dev.new(
  title = "帶有環(huán)境因子的PCA雙序圖-標(biāo)尺2",
  noRStudioGD = TRUE
)
biplot(spe.h.pca, main = "魚類數(shù)據(jù)PCA - 2型標(biāo)尺")
(spe.h.pca.env <-
    envfit(spe.h.pca, env, scaling = 2)) 
# 默認(rèn)2型標(biāo)尺

圖2.2-1

圖2.2-2

# 顯示顯著環(huán)境變量，用戶設(shè)定變量箭頭顏色
plot(spe.h.pca.env, p.max = 0.05, col = 3)

#這行命令是在最后版的雙序圖中增加了顯著環(huán)境變量
#注意：必須給envfit提供與繪圖相同的標(biāo)尺

圖2.2-3

在約束排序中，envfit（）函數(shù)也可以用于解釋變量與排序軸回歸 $r^2$ 的顯著性的置換檢驗(yàn)。但envfit（）并不是檢驗(yàn)解釋變量對(duì)響應(yīng)變量是否有顯著效應(yīng)最好的方法。

3. PCA的應(yīng)用領(lǐng)域

PCA在生態(tài)學(xué)應(yīng)用主要是用于基于定量環(huán)境變量的樣方排序，以及經(jīng)過適當(dāng)轉(zhuǎn)化后的群落組成數(shù)據(jù)的樣方排序。PCA原始假設(shè)要求數(shù)據(jù)符合多元正態(tài)分布。但在生態(tài)學(xué)領(lǐng)域，只要偏離正態(tài)分布不要太離譜，PCA對(duì)于數(shù)據(jù)是否正態(tài)分布并不敏感。

PCA的主要計(jì)算步驟是基于離散矩陣（線性協(xié)方差或相關(guān)）的特征分解過程。雖然協(xié)方差或相關(guān)系數(shù)必須用定量變量才能算出，但二元（1-0）數(shù)據(jù)也可以用于PCA分析。

PCA使用條件：

PCA分析必須從相同量綱的變量表格開始。因?yàn)樾枰獙⒆兞靠偡讲罘峙浣o特征根，所以變量必須有相同的物理單位，方差和才有意義（方差的單位是變量單位的平方），或者變量是無量綱的數(shù)據(jù)，例如標(biāo)準(zhǔn)化或?qū)?shù)轉(zhuǎn)化后的數(shù)據(jù)。

矩陣數(shù)據(jù)不能倒置，因?yàn)閷?duì)象（樣方）之間的協(xié)方差或相關(guān)沒有意義。

協(xié)方差或相關(guān)系數(shù)必須用定量變量才能算出。但由于半定量變量之間的Pearson相關(guān)系數(shù)等同于Spearman秩相關(guān)系數(shù)，所以基于半定量數(shù)據(jù)的PCA排序結(jié)果中變量之間關(guān)系是反映Spearman秩相關(guān)。

二元的數(shù)據(jù)也可以進(jìn)行PCA分析。排序圖對(duì)象之間的距離是1減去簡(jiǎn)單匹配系數(shù) $s_1$ 后的平方根（即 $\sqrt{1-S_1}$ ）乘以變量數(shù)量的平方根。

物種有-無數(shù)據(jù)進(jìn)行PCA分析前，可以用Hellinger轉(zhuǎn)化或弦轉(zhuǎn)化數(shù)據(jù)進(jìn)行預(yù)處理。因?yàn)橛?無數(shù)據(jù)的Hellinger距離或弦距離均等于 $\sqrt2\sqrt{1-Ochiai\ similarity}$ ，因此，有-無數(shù)據(jù)的Hellinger轉(zhuǎn)化或弦轉(zhuǎn)化后PCA分析的排序圖（1型標(biāo)尺）內(nèi)樣方之間的距離是Ochiai距離。 $\sqrt2\sqrt{1-Ochiai\ similarity}$ 是一種距離測(cè)度，這種測(cè)度也適用于有一無物種數(shù)據(jù)的分析。

應(yīng)該避免用變量箭頭頂點(diǎn)的距離來評(píng)估變量之間的相關(guān)性，而是用箭頭之間的角度來判斷相關(guān)性。

4.使用PCA.newr（）函數(shù)進(jìn)行PCA分析

為了方便能夠快速地用自己的數(shù)據(jù)進(jìn)行PCA分析，本書作者編寫了PCA.newr（）和biplot.PCA.newr（）兩個(gè)函數(shù)（在自編函數(shù)PCA.newr.R里面）。下面用Doubs環(huán)境數(shù)據(jù)為例演示這兩個(gè)函數(shù)的使用流程。

#PCA，這里默認(rèn)是1型標(biāo)尺雙序圖
# PCA; 這里默認(rèn)是1型標(biāo)尺雙序圖
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型標(biāo)尺雙序圖
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

這里主成分軸正負(fù)方向是隨機(jī)的，可能與vegan包輸出的排序圖成鏡像關(guān)系。但沒有關(guān)系，因?yàn)閷?duì)象或變量之間的相對(duì)位置沒有變化。

5. PCA中缺失值的估算

如果有缺失值該怎么辦?

是刪除缺失值所在的行和列，把珍貴的同行或列數(shù)據(jù)丟棄不用？或者，為了避免這種情況，用一些方法找到有意義的替代值進(jìn)行運(yùn)算（統(tǒng)計(jì)語時(shí)中稱為“插值"）。例如：利用相應(yīng)變量的平均值或通過其他相關(guān)變量的國(guó)歸分析獲得的估計(jì)插值。當(dāng)這兩種解決方案出現(xiàn)時(shí)，我們不會(huì)考慮估計(jì)缺失值不確定性及其在進(jìn)一步分析中的用途。結(jié)果，補(bǔ)充后的變量的標(biāo)準(zhǔn)誤被低估，導(dǎo)致接下來的運(yùn)算一些假設(shè)統(tǒng)計(jì)檢驗(yàn)失效。

為了解決這些缺點(diǎn)，Josse和Huson（2012）提出了“一個(gè)正則化（regu-larized）的迭代PCA算法，為主坐標(biāo)分析和主成分分析提供點(diǎn)估計(jì)”。這個(gè)算法與vegan包里的算法不一樣，PCA是一軸一軸逐步迭代出來。當(dāng)碰到缺失值的時(shí)候，會(huì)以該變量的平均值來替代。但不斷迭代過程，缺失值會(huì)根據(jù)新的計(jì)算結(jié)果不斷獲得新的估計(jì)值，直到收斂為止，也就是讓缺失值獲得穩(wěn)定的補(bǔ)充值。

實(shí)施復(fù)雜程序的目的在于：

①避免過度擬合估算模型（即當(dāng)數(shù)據(jù)里面有比較多的缺失值的時(shí)候，有太多參數(shù)需要調(diào)整的問題）

②克服低估內(nèi)插值方差的問題

注意：用推算得到的PCA的解決方案在不同PCA軸之間并不是嵌套關(guān)系，即基于s軸解決方案不等于與（s+1）或（s+2）解決方案。所以，做出所需的軸數(shù)適當(dāng)?shù)南闰?yàn)決定非常重要。該決定可能是經(jīng)驗(yàn)性的（例如，使用具有特征值之和大于斷棍模型預(yù)測(cè)的軸），或者在PCA重建過程，它可以通過數(shù)據(jù)本身的交又驗(yàn)證程序告知數(shù)據(jù)哪里刪除，以及計(jì)算重建誤差。我們的目的是保留最小化重建誤差的軸數(shù)。

missMDA包，匯總了PCA分析所有可能通過選代方式插值缺失值的方法。這個(gè)差值的函數(shù)為imputePcA（）。接下來將以Doubs環(huán)境數(shù)據(jù)做兩個(gè)實(shí)驗(yàn)，第一次只有3個(gè)缺失值（即所有319個(gè)數(shù)值的1%的缺失值），第二次是32個(gè)缺失值（10%）。

我們的第一個(gè)實(shí)驗(yàn)包括刪除選擇代表各種情況的三個(gè)值。

第一個(gè)值將靠近帶相對(duì)對(duì)稱分布的變量的平均值附近。pH平均值=8.04。讓我們刪除樣方2中pH值（8.0）。

第二個(gè)將接近高度不對(duì)稱分布的變量的均值，pho平均值=0.57。讓我們刪除樣方19中的pho值（第18行，0.60）。

第三個(gè)將從不對(duì)稱分布變量中刪除，但刪除后遠(yuǎn)離均值。bod平均值為5.01。讓我們刪除樣方23的bod值（第22行，16.4，是第二高的值）。

看一下imputePcA（）如何執(zhí)行迭代？

# PCA缺失值估計(jì)實(shí)驗(yàn)一:3個(gè)缺失值

# 用 NA 替換3個(gè)備選的值
env.miss3 <- env
env.miss3[2, 5] <- NA    # pH
env.miss3[18, 7] <- NA   # pho
env.miss3[22, 11] <- NA  # dbo

#三個(gè)變量的新均值（無缺失值）
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 = "原始和帶估計(jì)值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中缺失數(shù)據(jù)的估算。Procrustes比較原始環(huán)境變量PCA和對(duì)三個(gè)缺失值進(jìn)行估算的PCA。1型標(biāo)尺，這里只繪制樣方圖。紅色是原始PCA的樣方；藍(lán)色是估算值的PCA中的樣方。3個(gè)缺失值，即1%

上面這三個(gè)估算值，唯一比較好的估算是pH，它恰好落在原始值上。請(qǐng)記住，pH在數(shù)據(jù)中是對(duì)稱分布的。另外兩個(gè)估算遠(yuǎn)離原始值，這兩個(gè)變量都有強(qiáng)不對(duì)稱分布。顯然，這比缺失值接近或不接近變量的平均值更重要（pho屬于這種情況，但bod不是這種情況）。

那么如何影響排序結(jié)果（前2軸）？

為了觀察這個(gè)結(jié)果，在帶有估算值矩陣上計(jì)算新的PCA，并將其與原始PCA進(jìn)行比較。通過Procrustes旋轉(zhuǎn)新環(huán)境變量PCA軸，使其與原始PCA一致。Procrustes旋轉(zhuǎn)是最小化兩個(gè)排序結(jié)構(gòu)同一標(biāo)號(hào)的對(duì)象圖上的距離平方和，目的在于獲得最多的重疊區(qū)域。Procrustes旋轉(zhuǎn)可以通過vegan包的函數(shù)procrustes（）來執(zhí)行。

顯然，實(shí)際值和重建值之間的差異對(duì)排序的影響很小。唯一可見的區(qū)別是在樣方23，其中實(shí)際值與重建值之間的差異最大。

#缺失值估計(jì)實(shí)驗(yàn)二:32個(gè)缺失值
#從319個(gè)值中隨機(jī)抽取32個(gè)用NA替代
rnd <- matrix(sample(c(rep(1, 32), rep(0, 287))), 29, 11)
env.miss32 <- env
env.miss32[rnd == 1] <- NA
# 每個(gè)樣方里面多少NA?
summary(t(env.miss32))
# 顯示NA數(shù)量的替代方法：
# sapply(as.data.frame(t(env.miss32)), function(x) sum(is.na(x)))

# 內(nèi)插估計(jì)
env.imp2 <- imputePCA(env.miss32)

# 帶估計(jì)值的PCA
env.imp32 <- env.imp2$completeObs
env.imp32.pca <- rda(env.imp32, scale = TRUE)

# 原始PCA和帶估計(jì)值的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 = "原始和帶估計(jì)值PCA的procrustes比較\n32缺失值")
points(pca.proc32, display = "target", col = "red")
text(
  pca.proc32,
  display = "target",
  col = "red",
  pos = 4,
  cex = 0.6
)