看到本筆記系列的名字么？：R在數(shù)量生態(tài)學(xué)中的應(yīng)用--矩陣·度量·聚類·排序·空間哼丈。其實到排序這一部分已經(jīng)算是接近尾聲了关串，因為空間分析哪一部分我打算放棄，目前的生態(tài)數(shù)據(jù)規(guī)模很少有空間數(shù)據(jù)课兄。接下來，搬好小板凳晨继，咱們好好講講排序烟阐。

如果說聚類分析的目的在于尋找數(shù)據(jù)的間斷性，那么排序（ordination）的目的就在于尋找數(shù)據(jù)的連續(xù)性（通過連續(xù)的排序軸展示數(shù)據(jù)的主要趨勢）。

要講排序就要知道我們是在哪里排序蜒茄，這里就引出多維空間的概念唉擂。假設(shè)我有10個樣本的物種豐度表：

行為樣本，列為物種（也可以理解為特征）檀葛。我們要對這10個樣本進(jìn)行排序：

• 假如只有一個物種: 那么根據(jù)這一個物種在每個樣本中的值就可以排序啦玩祟。
• 假如有兩個物種（"CHA"， "TRU"）：我們可以建立二維坐標(biāo)軸屿聋，橫坐標(biāo)是"CHA"卵凑，縱坐標(biāo)為TRU，根據(jù)這兩個物種的值胜臊，我們也可以畫出點(diǎn)來。
• 假如有三個物種：那么就是三維坐標(biāo)來畫點(diǎn)伙判，也是可以畫的象对。

那么大于三個物種的時候呢？那就是n維空間中的點(diǎn)了宴抚，是無法畫出來的勒魔。
所以我們要找到一種方法，將n維空間中的點(diǎn)菇曲，在二維平面內(nèi)展現(xiàn)出來冠绢。
由于這么多的點(diǎn)無法共面，所以要找到一個平面常潮，使空間中的所有點(diǎn)都能投影在這個平面上弟胀，而且投影的點(diǎn)之間的距離，越能反應(yīng)真實距離越好喊式。這個投影過程就是排序運(yùn)算過程孵户。好的排序方法是投影過程信息損失最少。

排序的主要目的之一是生成可視化的排序圖岔留，這就決定了排序過程實際上是講多維空間的數(shù)據(jù)盡可能的數(shù)據(jù)點(diǎn)排列在可視化的低維空間夏哭，也就是使最前面的幾個排序軸盡可能包含數(shù)據(jù)結(jié)構(gòu)變化的主要趨勢。本章講的非約束排序只是描述性方法献联，不存在檢驗評估排序結(jié)果是否顯著性的問題竖配，下一章約束排序則需要對排序結(jié)果進(jìn)行顯著性檢驗。

降維空間內(nèi)的排序

大部分排序方法（NMDS除外）都是基于關(guān)聯(lián)矩陣特征向量的提取里逆。這也產(chǎn)生一個問題：多少個排序軸值得保留和解讀进胯？

PCA

如果一個數(shù)據(jù)矩陣中每個變量都是正態(tài)分布，這樣的矩陣可以稱作多元正態(tài)分布矩陣运悲。PCA排序分析的對象是離差矩陣（dispersion matrix）龄减，即包含方差和協(xié)方差的變量之間的關(guān)聯(lián)矩陣，或不同綱量的變量之間的相關(guān)系數(shù)矩陣班眯∠Ｍ＃可見烁巫，如果變量之間不相關(guān)，PCA的分析也就沒有意義了宠能。PCA致力于分析定量數(shù)據(jù)亚隙，展示對象的歐氏距離，線性關(guān)系违崇。

使用rda()函數(shù)對Doubs環(huán)境數(shù)據(jù)進(jìn)行PCA分析

# 導(dǎo)入本章所需的程序包 
library(ade4)
library(vegan)
library(gclus)  
library(ape)
rm(list = ls())
setwd("D:\\Users\\Administrator\\Desktop\\RStudio\\數(shù)量生態(tài)學(xué)\\DATA")
# 導(dǎo)入CSV文件數(shù)據(jù) 
spe <- read.csv("DoubsSpe.csv", row.names=1)
env <- read.csv("DoubsEnv.csv", row.names=1)
spa <- read.csv("DoubsSpa.csv", row.names=1)
# 刪除沒有數(shù)據(jù)的樣方8
spe <- spe[-8,]
env <- env[-8,]
spa <- spa[-8,]
# 顯示環(huán)境變量數(shù)據(jù)集的內(nèi)容
summary(env)    # 描述性統(tǒng)計
# 全部環(huán)境變量數(shù)據(jù)的PCA分析（基于相關(guān)矩陣：參數(shù)scale=TURE）
# ********************************************************
env.pca <- rda(env, scale=TRUE) # 參數(shù)scale=TRUE 表示對變量進(jìn)行標(biāo)準(zhǔn)化
env.pca

Call: rda(X = env, scale = TRUE)

              Inertia Rank
Total              11     
Unconstrained      11   11
Inertia is correlations 

Eigenvalues for unconstrained axes:
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11 
6.098 2.167 1.038 0.704 0.352 0.319 0.165 0.112 0.023 0.017 0.006 

summary(env.pca) # 默認(rèn)scaling=2
summary(env.pca, scaling=1)

注意函數(shù)summary（）內(nèi)的參數(shù)scaling阿弃，為繪制排序圖所選擇的標(biāo)尺類型，與函數(shù)rda（）內(nèi)數(shù)據(jù)標(biāo)準(zhǔn)化的參數(shù)scale無關(guān)羞延。

Call:
rda(X = env, scale = TRUE) 

Partitioning of correlations:
              Inertia Proportion
Total              11          1
Unconstrained      11          1

Eigenvalues, and their contribution to the correlations 

Importance of components:
                         PC1    PC2     PC3     PC4     PC5     PC6     PC7     PC8      PC9     PC10      PC11
Eigenvalue            6.0980 2.1671 1.03760 0.70351 0.35185 0.31913 0.16455 0.11171 0.023109 0.017361 0.0060618
Proportion Explained  0.5544 0.1970 0.09433 0.06396 0.03199 0.02901 0.01496 0.01016 0.002101 0.001578 0.0005511
Cumulative Proportion 0.5544 0.7514 0.84570 0.90966 0.94164 0.97066 0.98561 0.99577 0.997871 0.999449 1.0000000

Scaling 1 for species and site scores
* Sites are scaled proportional to eigenvalues
* Species are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores:  4.189264 


Species scores

         PC1     PC2      PC3     PC4     PC5     PC6
das  1.45634  1.1597 -0.83818 -0.6394  1.1820 -0.5578
alt -1.40158 -1.3396  0.58576  0.4854  0.7004  0.8234
pen -0.77254 -1.1499 -1.80693 -3.1715  0.1565  1.1780
（......）


Site scores (weighted sums of species scores)

        PC1      PC2       PC3      PC4      PC5      PC6
1  -1.05160 -0.65504 -0.536188 -0.74739  0.04135  0.08372
2  -0.77560 -0.36292  0.104662  0.13751  0.16547 -0.30155
3  -0.70642 -0.21672  0.417881 -0.05503  0.18807 -0.11896
4  -0.65572 -0.13076  0.064532  0.16799 -0.04275 -0.01082
5  -0.31707 -0.29517  0.238423  0.19922  0.11293  0.20052
（······）

PCA 結(jié)果術(shù)語：

• Inertia （慣量）：變量自相關(guān)系數(shù)的總和渣淳，也等于變量的個數(shù)。

• Constrained Unconstrained （約束非約束）：PCA是非約束排序

• Eigenvalue （特征根）：每個排序軸重要性（方差）的指標(biāo)伴箩，可以用特征根數(shù)值入愧，也可以用占總變差的比例表示（每軸特征根除以總慣量）

• Scaling(標(biāo)尺）：排序結(jié)果投影到排序空間的可視化方式。

  • Scaling1（1型標(biāo)尺）=（distance biplot嗤谚，距離雙序圖）：特征向量被標(biāo)準(zhǔn)化為單位長度棺蛛，關(guān)注的是對象之間的關(guān)系。雙序圖中對象之間的距離近似于歐氏距離巩步，代表變量的箭頭之間的夾角沒有意義旁赊。
  • Scaling2 （2型標(biāo)尺）=（correlation，biplot 相關(guān)雙序圖）：特征向量被標(biāo)準(zhǔn)化為特征根的平方根椅野，關(guān)注的是變量之間的關(guān)系终畅。雙序圖中對象之間的距離不再近似于歐氏距離，代表變量的箭頭之間的夾角反映變量之間的相關(guān)性鳄橘。

• Species scores：變量箭頭在排序圖中的坐標(biāo)

• Site scores ： 對象在排序圖中的坐標(biāo)声离。

提取、解讀和繪制vegan程序輸出的PCA結(jié)果

PCA不是統(tǒng)計檢驗瘫怜，而是一種探索分析方法术徊，其目的是在低維空間盡可能多滴展示數(shù)據(jù)的主要趨勢特征。選擇多少個排序軸并沒有統(tǒng)一的標(biāo)準(zhǔn)鲸湃。不過也有兩種判別方法：

• Kaiser-Guttman準(zhǔn)則
• 斷棍模型（broken stick model）

# 查看和繪制PCA輸出的部分結(jié)果 
# ****************************
?cca.object  # 解釋vegan包輸出的排序結(jié)果對象結(jié)構(gòu)和如何提取部分結(jié)果
# 特征根
(ev <- env.pca$CA$eig)
        PC1         PC2         PC3         PC4         PC5         PC6         PC7         PC8         PC9        PC10        PC11 
6.097995220 2.167125814 1.037602930 0.703508127 0.351852003 0.319125849 0.164550682 0.111707052 0.023109109 0.017361390 0.006061823 
# 應(yīng)用Kaiser-Guttman準(zhǔn)則選取排序軸 
 PC1      PC2      PC3 
6.097995 2.167126 1.037603 

# 斷棍模型
n <- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) {
    bsm$p[i] = bsm$p[i-1] + (1/(n + 1 - i))
}
bsm$p <- 100*bsm$p/n
bsm
    j          p
1   1  0.8264463
2   2  1.7355372
3   3  2.7456382
4   4  3.8820018
5   5  5.1807031
6   6  6.6958547
7   7  8.5140365
8   8 10.7867637
9   9 13.8170668
10 10 18.3625213
11 11 27.4534304
# 繪制每軸的特征根和方差百分比 
par(mfrow=c(2,1))
barplot(ev, main="特征根", col="bisque", las=2)
abline(h=mean(ev), col="red")   # 特征根平均值
legend("topright", "平均特征根", lwd=1, col=2, bty="n")
barplot(t(cbind(100*ev/sum(ev),bsm$p[n:1])), beside=TRUE, 
    main="% 變差", col=c("bisque",2), las=2)
legend("topright", c("% 特征根", "斷棍模型"), 
    pch=15, col=c("bisque",2), bty="n")
    #兩種方法是否保留相同的軸數(shù)赠涮？

樣方和變量的雙序圖

 兩種PCA雙序圖：1型標(biāo)尺和2型標(biāo)尺 
#********************************
# 使用biplot（）函數(shù)繪制排序圖 
par(mfrow=c(1,2))
biplot(env.pca, scaling=1, main="PCA-1型標(biāo)尺")
biplot(env.pca, main="PCA-2型標(biāo)尺")  # 默認(rèn) scaling = 2

# 使用cleanplot.pca（）函數(shù)繪圖
source("cleanplot.pca.R")
cleanplot.pca(env.pca, point=TRUE) 
# point=TRUE表示樣方用點(diǎn)表示，變量用箭頭表示

cleanplot.pca(env.pca)              
# 默認(rèn)樣方僅用序號標(biāo)識（同vegan包的標(biāo)準(zhǔn)）

cleanplot.pca(env.pca, ahead=0)    
# ahead=0表示變量用無箭頭的線表示
#左圖中圓圈代表什么意義暗挑？看下面解釋內(nèi)容

• 1型標(biāo)尺雙序圖中的平衡貢獻(xiàn)圓：
  其半徑代表變量的向量長度對排序的平均貢獻(xiàn)率笋除。因此在任何二維排序圖中，如果某個變量的箭頭長度長于圓的半徑炸裆，代表他對這個排序空間的貢獻(xiàn)大于所有變量的平均長度垃它。

組合聚類結(jié)果和排序結(jié)果

• 在排序圖內(nèi)用不同顏色去區(qū)分不同樣方組
• 將聚類樹添加到排序圖上

# 組合聚類分析結(jié)果和排序結(jié)果
# ***************************
# 使用環(huán)境變量數(shù)據(jù)對樣方進(jìn)行基于變量標(biāo)準(zhǔn)化后歐氏距離的Ward聚類分析
env.w <- hclust(dist(scale(env)), "ward")
# 裁剪聚類樹鸿摇，只保留4個聚類簇
gr <- cutree(env.w, k=4)
grl <- levels(factor(gr))
# 提取樣方坐標(biāo)渐逃，1型標(biāo)尺
sit.sc1 <- scores(env.pca, display="wa", scaling=1)
# 按照聚類分析的結(jié)果對樣方進(jìn)行標(biāo)識和標(biāo)色（1型標(biāo)尺）
p <- plot(env.pca, display="wa", scaling=1, type="n", 
    main="PCA（基于相關(guān)矩陣）+聚類簇")
abline(v=0, lty="dotted")
abline(h=0, lty="dotted")
for (i in 1:length(grl)) {
    points(sit.sc1[gr==i,], pch=(14+i), cex=2, col=i+1)
    }
text(sit.sc1, row.names(env), cex=.7, pos=3)
# 在排序圖內(nèi)添加聚類樹
ordicluster(p, env.w, col="dark grey")
legend(locator(1), paste("組",c(1:length(grl))), pch=14+c(1:length(grl)), 
  col=1+c(1:length(grl)), pt.cex=2)
  

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

PCA是一種展示樣方之間的歐氏距離的線性方法，理論上并不適用于物種多度分析。然而我們可以通過對數(shù)據(jù)進(jìn)行轉(zhuǎn)化來做PCA椎例。

# 魚類多度數(shù)據(jù)的PCA分析
# **********************
# 物種數(shù)據(jù)Hellinger轉(zhuǎn)化
spe.h <- decostand(spe, "hellinger")
spe.h.pca <- rda(spe.h)
spe.h.pca
Call: rda(X = spe.h)

              Inertia Rank
Total          0.5025     
Unconstrained  0.5025   27
Inertia is variance 

Eigenvalues for unconstrained axes:
    PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
0.25796 0.06424 0.04632 0.03850 0.02197 0.01675 0.01472 0.01156 
(Showed only 8 of all 27 unconstrained eigenvalues)
#繪制每軸的特征根和方差百分比
ev <- spe.h.pca$CA$eig
evplot(ev)

# PCA雙序圖
cleanplot.pca(spe.h.pca, ahead=0)
#這里物種不像環(huán)境變量那樣能夠明顯分組锌半。但也能看出物種沿著梯度更替
#分布昌粤。在1型標(biāo)尺的雙序圖內(nèi)窍蓝，可以觀察到有8個物種對于第一、二軸有很大
#貢獻(xiàn)务热∫涫龋可以比較一下，這些物種與4.10.3節(jié)聚類分析中聚類簇的指示種是否
#重合崎岂？

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

# 使用PCA（）和biplot.PCA（）兩個函數(shù)進(jìn)行環(huán)境數(shù)據(jù)的PCA分析
# **********************************************************
source("PCA.R")  #導(dǎo)入PCA.R腳本捆毫，此腳本必須在當(dāng)前工作目錄下或給路徑
# PCA，這里默認(rèn)是1型標(biāo)尺雙序圖
env.PCA.PL1 <- PCA(env, stand=TRUE)
biplot.PCA(env.PCA.PL1)
abline(h=0, lty=3)
abline(v=0, lty=3)
# PCA冲甘，生成2型標(biāo)尺雙序圖
env.PCA.PL2 <- PCA(env, stand=TRUE)
biplot.PCA(env.PCA.PL2 ,scaling=2)
abline(h=0, lty=3)
abline(v=0, lty=3)
#這里主成分軸正負(fù)方向是隨機(jī)的冻璃，可能與vegan包輸出的排序圖成鏡像關(guān)
#系。但沒有關(guān)系损合，因為對象或變量之間的相對位置沒有變化。

參考：
相應(yīng)分析的R包c(diǎn)a和mca娘纷，cca嫁审，RDA的R實現(xiàn)整理
使用Vegan包進(jìn)行生態(tài)學(xué)數(shù)據(jù)排序分析的學(xué)習(xí)（一）
原來你是這樣的排序分析