某研究對三種環境中的土壤進行了取樣（環境A、B、C，各環境下采集9個土壤樣本），結合二代測序觀測了土壤細菌群落物種組成，並對多種土壤理化指標進行了測量。獲得兩個資料集。
（1）檔案“phylum_table.txt”為透過16S測序所得的細菌類群丰度表格，此處統計至細菌“門水平”（即界門綱目科屬種的“門”這一水平）。（為了方便演示，我沒使用OTU水平的）
（2）檔案“env_table.txt”記錄了多種土壤理化指標資訊，即環境資料。
接下來期望透過RDA分析這些資料，用以分析群落物種組成與環境梯度的關係。
RDA的初步執行
vegan包是生態學統計分析中最常用的R包之一，已經打包好了關於RDA計算、檢驗等一些列命令，方便我們快速有效地執行RDA。
vegan包中，透過rda()執行RDA，該函式的輸入格式大致如下。
library(vegan)
#詳情 ?rda
#呼叫格式 1
rda(Y, X, W)
#或者格式 2
rda(Y~var1+var2+var3+factorA+var2*var3+Condition(var4))
在第1種呼叫格式中，Y為響應變數矩陣，X為解釋變數矩陣，W為偏RDA分析需要輸入的協變數矩陣。當不考慮協變數時，可直接輸入響應變數矩陣Y和解釋變數矩陣X，即“rda(Y, X)”。這種寫法雖然簡單，但要求解釋變數矩陣或協變數矩陣中不能含有因子變數（定性變數）。（當僅存在響應變數矩陣Y時，則執行PCA排序）
因此，通常建議使用第2種寫法。其中Y為響應變數矩陣，var1、var2、var3為定量解釋變數，factorA為因子變數，var2*var3為變數2和變數3的互動作用項，var4為協變數。與上述第1種寫法相比，不僅能夠識別因子變數，還能夠對變數間的互動作用加以考慮。下文的具體示例中，均以第2種寫法作展示。
RDA包含RDA、基於轉化的冗餘分析（tb-RDA）、偏冗餘分析（偏RDA）、基於距離的冗餘分析（db-RDA）、非線性RDA等多種模式。實際的資料分析中，需根據情況選擇合適的分析模式。關於這些RDA的基本概念描述，可參考前文。
以下分別簡介在R中執行這幾種RDA模式的命令方法。
常規RDA
即提供響應變數矩陣Y，以及解釋變數矩陣X，將原始資料直接透過rda()函式分析。
##讀取資料
#讀取物種資料，細菌門水平丰度表
phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#讀取環境變數
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#使用全部的環境資料
rda_result <- rda(phylum~., env, scale = FALSE)
使用物種資料和環境變數資料，其中原始物種組成資料（資料框phylum）作為響應變數，原始環境因子資料（資料框env）作為解釋變數。式中scale引數用於對響應變數矩陣（與解釋變數矩陣無關）執行標準化處理。scale = FALSE，不執行標準化，使用擬合值的協方差矩陣運算RDA；若為TRUE，則執行標準化，基於相關矩陣運算RDA。響應變數的量綱不一致時（不過，物種多度資料不存在這種情況），或者存在極端值影響時（也可以透過額外的資料轉化來解決），推薦scale = TRUE。
“~.”表示將資料框env中所有的列變數（環境因子）作為解釋變數帶入RDA排序，是一種變數全選的省略模式（在不考慮協變數以及變數間的互動作用項的情況下，可以直接使用），能夠自動識別定量解釋變數以及定性解釋變數型別。對於特殊情況，可以寫作如下樣式。
#若只關注區域性環境資料，除了在原始表格中修改變數個數外，還可直接在 rda() 中指定
#例如只考慮 pH、TC、TN、AP、AK 這 5 種環境變數
rda_result <- rda(phylum~pH+TC+TN+AP+AK, env, scale = FALSE)
#預設情況下，rda(phylum~., env)，不包含環境變數間的互動作用
#若想關注某兩個環境變數間的互動作用，例如新增 TC 和 TN 的互動作用項
rda_result <- rda(phylum~pH+TC+DOC+SOM+TN+NO3+NH4+AP+AK+TC*TN, env, scale = FALSE)
#當存在協變數時（即偏 RDA，可參見下文）
#例如控制土壤 pH 影響後（pH 作為協變數），觀測其它環境因素的影響
rda_result <- rda(phylum~TC+DOC+SOM+TN+NO3+NH4+AP+AK+Condition(pH), data = env, scale = FALSE)
此時我們將RDA結果儲存在物件rda_result中，內容說明參見下文“RDA結果的初步解讀”。
tb-RDA（本篇重點展示的是它）
對應於tb-PCA，通常情況下，對於物種多度資料，不建議直接用於RDA排序。特別是當物種丰度表中出現很多“0”值的情況下，應當首先進行Hellinger轉化後，再輸入至RDA中。這種首先對原始資料作一定的轉化後在執行RDA的方法稱為基於轉化的冗餘分析（Transformation-based redundancy analysis，tb-RDA）。除了第一步的資料轉化外，其與過程均和常規的RDA相同。
因此我們首先對示例物種丰度資料進行Hellinger轉化，並使用轉化後的資料執行RDA。
##tb-RDA
#物種資料 Hellinger 預轉化（處理包含很多 0 值的群落物種資料時，推薦使用）
phylum_hel <- decostand(phylum, method = 'hellinger')
#使用全部的環境資料
rda_tb <- rda(phylum_hel~., env, scale = FALSE)
#參考下文中對 RDA 結果解讀的說明，可分別使用 plot(rda_result) 和 plot(rda_tb) 檢視比較 Hellinger 轉化前後二者的差異
rda()命令中，使用Hellinger轉化後的物種資料（資料框phylum_hel）作為響應變數，原始環境因子資料（資料框env）作為解釋變數。式中scale引數用於對響應變數矩陣（與解釋變數矩陣無關）執行標準化處理。scale = FALSE，不執行標準化，使用擬合值的協方差矩陣運算RDA；若為TRUE，則執行標準化，基於相關矩陣運算RDA。Hellinger轉化後的物種多度資料，通常無需標準化，首先物種資料的量綱均相同，其次Hellinger轉化消除了極端值的影響，降低了高丰度物種的權重。
此時我們將tb-RDA結果儲存在物件rda_tb中，內容說明參見下文“RDA結果的初步解讀”。
偏RDA
有時，我們期望在控制變數矩陣W對響應變數Y的影響的前提下，將關注點集中在另一組變數矩陣X對響應變數Y的影響上，應用到RDA中即為偏冗餘分析（Partial canonical ordination ，偏RDA）。
例如對於示例資料，我們期望控制土壤pH影響後（pH作為協變數），觀測其它環境因素的影響。
##偏 RDA
#控制土壤 pH 影響後（pH 作為協變數），觀測其它環境因素的影響；物種資料 Hellinger 預轉化
rda_part <- rda(phylum_hel~TC+DOC+SOM+TN+NO3+NH4+AP+AK+Condition(pH), data = env, scale = FALSE)
rda()命令中，使用Hellinger轉化後的物種資料（資料框phylum_hel）作為響應變數，原始環境因子資料（資料框env）作為解釋變數。由於加入了協變數，此時“~”後不再直接使用“.”表示全部的環境因子，而是透過指定的方式，將所考慮的解釋變數的名稱以“+”連線，並將協變數（示例為pH）新增至Condition()裡。
此時我們將偏RDA結果儲存在物件rda_part中，內容說明參見下文“RDA結果的初步解讀”。
db-RDA
基於距離的冗餘分析（Distance-based redundancy analysis，db-RDA）。
考慮到它比較特殊，且又非常的重要，本篇暫且跳過這種模式，我把它單獨整理了一篇，放在下篇展示。
非線性關係的RDA
當響應變數與解釋變數明顯為非線性關係（例如單峰關係，當然，若響應變數與解釋變數存在明顯的單峰響應模式，可使用CCA來解決），常規的線性RDA模型不適合時，可考慮使用這種方法來解決。
由於我並沒有充分理由認為示例資料中響應變數與解釋變數存在明顯的非線性關係（鑑定過程繁瑣，懶得觀測……），故這裡不再使用示例資料演示這種非線性關係的RDA過程了。大家有興趣可參考賴江山老師譯作“數量生態學:R語言的應用”170-174頁內容。
RDA結果的初步解讀和資訊提取
上文簡介了使用R的vegan包中的rda()命令，初步獲得幾種不同模式RDA排序結果的方法。在初步得到排序結果後，我們首先需要對結果進行解讀。
下文僅對rda()執行結果中的各部分內容所反映的資訊（意義）及其提取方法作簡要描述，並初步展示排序圖的繪製方法，即主要以軟體所得執行結果為例展開敘述以幫助快速上手軟體使用。對於RDA排序結果解讀的更詳細說明，即主要的基本概念部分，如變差（方差）解釋程度、模型初步評估、排序圖的詳解、I型或II型標尺的縮放原理等，本文不作介紹，可參考前文。
RDA結果總概
以上述tb-RDA排序結果“rda_tb”為例，對vegan所得RDA結果的主要結構內容進行簡要的說明。
##RDA 結果解讀，以下以 tb-RDA 結果為例
#檢視統計結果資訊，以 I 型標尺為例
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
rda_tb.scaling1
summary()展示結果時，根據所重要要關注的內容，可選根據I型標尺或II型標尺縮放排序座標後展示出。如果對排序樣方之間的距離更感興趣，或者大多數解釋變數為因子型別，則考慮I型標尺；如果對變數之間的相關關係更感興趣，則考慮II型標尺。關於I型標尺和II型標尺，可參考前文。這裡我們展示了按I型標尺縮放後排序結果（引數scaling = 1；同理，scaling = 2展示II型標尺；引數中還有scaling = 0和scaling = 3的選項，這個我就不清楚了），部分內容如下所示。
僅透過初步結果，我們發現該RDA模型承載的總方差比例達到了70%以上（0.7338），特別是RDA1和RDA2的解釋量總和為“0.53445+0.118470=0.65292”（即承載了響應變數矩陣65.29%的方差），非常可觀的一個數字。做到這裡是不是很開心地去整理資料去了？別急，仔細看下紅字部分的標註，這些結果並不是很可信，有待於校正和檢驗通過後才能被接受。詳見下文“RDA的R2校正和約束軸的置換檢驗”。
對於偏RDA結果，解讀方式也大致同上。但是有這麼幾個地方存在區別。
RDA排序圖
先不急介紹RDA模型的驗證方法，還得先把基礎部分介紹完……
排序結果最終以排序圖作為呈現。對於vegan所得rda排序結果，可直接使用plot()繪製排序圖。若想提前瞭解更多，可使用?plot.cca()檢視詳細說明（注意不是?plot()）。
此處同樣以上述tb-RDA排序結果“rda_tb”為例，繪製簡單樣式的RDA三序圖作為展示。圖中對於排序物件（樣方），展示為“標籤點”；響應變數（物種變數）和解釋變數（環境變數）均以向量表示。
#作圖檢視排序結果，三序圖，包含 I 型標尺和 II 型標尺
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = 'I 型標尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling1 <- scores(rda_tb, choices = 1:2, scaling = 1, display = 'sp')
arrows(0, 0, rda_sp.scaling1[ ,1], rda_sp.scaling1[ ,2], length = 0, lty = 1, col = 'red')
plot(rda_tb, scaling = 2, main = 'II 型標尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling2 <- scores(rda_tb, choices = 1:2, scaling = 2, display = 'sp')
arrows(0, 0, rda_sp.scaling2[ ,1], rda_sp.scaling2[ ,2], length = 0, lty = 1, col = 'red')
plot()命令中，display引數用於定義在排序圖中展示哪些資訊。其中，“sp”代表物種，“wa”代表使用物種加權和計算的樣方座標，“lc”代表以擬合值（解釋變數的線性組合）計算樣方的座標，“cn”代表約束成分（即解釋變數）。可知對於樣方座標有兩種可選展示方式，上述示例使用了“wa”。選擇哪種排序座標在排序圖繪製樣方的點，至今還有爭議（參見“數量生態學:R語言的應用”148頁內容），但無論如何，排序圖的解讀都需要依據統計顯著性檢驗結果，與多元迴歸一樣，不顯著的迴歸結果不能被解讀，必須丟棄（RDA的統計檢驗部分參見下文）。
plot()預設作圖結果中，物種同樣以“標籤點”來表示，因此我們需要額外繪製箭頭以展示物種向量。scores()用於提取物種、樣方或環境變數的排序座標，同樣以display引數定義。這裡使用scores()提取了物種變數的排序座標，並使用arrows()在排序圖中為物種向量新增箭頭，以(0,0)出發，指向對應的座標位置處。
由於物種資訊較多，影響觀測。如果僅考慮透過排序圖展示環境變數對樣方的約束資訊，那麼通常可以在圖中忽略物種資訊，即展示為雙序圖樣式。
下圖展示了僅包含樣方和環境變數資訊的雙序圖，同時分別使用display = c('wa', 'cn')和display = c('lc', 'cn')引數，在排序圖中分別展示了使用物種加權和計算的樣方座標以及以擬合值（解釋變數的線性組合）計算樣方的座標。（可觀察比較二者的區別）
#隱藏物種資訊，以 I 型標尺為例展示雙序圖，並檢視分別使用物種加權計算的樣方座標（wa）以及擬合的樣方座標（lc）的差異
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = 'I 型標尺，加權', display = c('wa', 'cn'))
plot(rda_tb, scaling = 1, main = 'I 型標尺，擬合', display = c('lc', 'cn'))
以上暫時僅對排序圖的視覺化做個簡單介紹。對於更詳細的視覺化細節部分，可參見本篇末尾。
備註：對於vegan的排序結果使用plot()繪製排序圖時，所得排序圖中解釋變數（環境變數）向量的箭頭頂點座標可能並非我們透過summary()所觀測到的實際座標位置，而更像是在原座標的基礎上同比例放大（或縮小）了一定的數值後展示的。這是因為在某些情況下，若直接展示原始座標，則向量箭頭的長度相較於排序圖將會顯得比例失調（過長或過短），嚴重影響我們對圖的解讀，此時通常考慮放大（或縮小）一定數值後展示出（可參見Legendre和Legendre，1998，404頁）。vegan預設的做圖方式中即考慮了這一點。
RDA結果中的主要資訊提取
若有需要，考慮在結果中將需要的資訊提取出來。以下同樣以上述tb-RDA排序結果“rda_tb”為例展示RDA重要資訊的提取方法。
上述在初步繪製排序圖觀測時，提到了可使用scores()提取物種、樣方或環境變數的排序座標，並且已經使用scores()提取了物種變數的排序座標。同樣地，修改display引數後即可提取樣方和環境變數的排序座標。
#scores() 提取排序得分（座標），以 I 型標尺為例，前四軸為例
#使用物種加權和計算的樣方得分
rda_site.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'wa')
#物種變數（響應變數）得分
rda_sp.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'sp')
#環境變數（解釋變數）得分
rda_env.scaling1 <- scores(rda_tb, choices = 1:4, scaling = 1, display = 'bp')
#或者在 summary() 後提取，以 I 型標尺為例，前四軸為例
rda_tb.scaling1 <- summary(rda_tb, scaling = 1)
names(rda_tb.scaling1)
#使用物種加權和計算的樣方得分
rda_site.scaling1 <- rda_tb.scaling1$site[ ,1:4]
#物種
rda_sp.scaling1 <- rda_tb.scaling1$species[ ,1:4]
#環境
rda_env.scaling1 <- rda_tb.scaling1$biplot[ ,1:4]
#若需要輸出在本地
#樣方
write.table(data.frame(rda_site.scaling1), 'rda_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種
write.table(data.frame(rda_sp.scaling1), 'rda_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境
write.table(data.frame(rda_env.scaling1), 'rda_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#不建議直接在原始資料集中提取，因為這樣提取的座標數值未經標尺縮放處理，不利於反映生物學問題
#rda_tb$CCA$u[ ,1:4]
#rda_tb$CCA$v[ ,1:4]
#rda_tb$CCA$biplot[ ,1:4]
RDA基於多元迴歸，若對典範特徵係數（即每個解釋變數與每個約束軸之間的迴歸係數）感興趣，可使用coef()提取（使用?coef()以瞭解更多）。
#coef() 提取典範係數
rda_coef <- coef(rda_tb)
關於約束模型解釋的總變差（即響應變數矩陣的總方差能夠被解釋變數解釋的部分，如果用比例表示，即等同於多元迴歸的R2，也就是約束軸的總解釋量）以及各約束軸解釋量等，不建議直接在原始結果中獲取。原因很簡單，正如前文所提到，初始所得的RDA模型需要經過校正和檢驗後才可使用，否則可能會帶來較大的偏差。
因此，下文的內容很重要，對於嚴謹的RDA模型的構建，請不要忽略。
上文中提到，與多元迴歸的未校正R2一樣，RDA的R2也有偏差，需要進行校正（即原始R2不可直接使用）。RDA結果中所展示的各約束軸的解釋量也是基於未校正的R2得來的，也需要根據校正後的R2重新計算（即原始所得約束軸的解釋量是不可信的，同樣不能直接使用）。
對於校正前後R2的提取，可使用RsquareAdj()完成。
#RsquareAdj() 提取 R2，詳情 ?RsquareAdj()
r2 <- RsquareAdj(rda_tb)
rda_noadj <- r2$r.squared #原始 R2
rda_adj <- r2$adj.r.squared #校正後的 R2
R2經過校正後通常會降低，這是必然的。
對於示例資料的tb-RDA排序結果，透過上文，我們已經得知了響應變數矩陣的總方差為0.02686，未校正前的約束軸解釋的方差總和為0.01971（佔比0.7338，即未校正前的R2）；那麼根據校正後的R2（0.5928），可還原校正R2後的約束軸解釋的方差總和約為“0.02686×0.5928=0.01592”。
同樣地，我們也可獲得各約束軸校正後的解釋量。例如，已知RDA第一個約束軸（按解釋量或特徵值的大小由高往低排序，各約束軸記為RDA1、RDA2、RDA3……）的解釋方差佔約束模型總解釋方差的0.72835（特徵值佔比0.72835，同樣地可知，該軸的解釋量佔約束模型總解釋量的比例也是0.72835），由於校正R2後的約束模型的總解釋量降至0.5928，那麼可推斷校正R2後的第一個約束軸的解釋量（第一軸的解釋方差佔響應變數矩陣總方差的比例）即為“0.5928×0.72835=0.43177”。同理，校正R2後RDA第二個約束軸的解釋量“0.5928×0.16145=0.09571”，其餘約束軸的解釋量以此類推。
對於我們的示例資料的tb-RDA結果，儘管在R2校正後約束模型的總體解釋量有所下降（由0.7338下降至0.5928），但看起來似乎仍很可觀。特別是的前兩個約束軸的解釋量仍然非常好：RDA1和RDA2的解釋量總和為“0.43177+0.09571=0.52748”，即承載了響應變數矩陣52.75%的方差，達到了一半以上。
約束軸的置換檢驗及p值校正
首先需對全模型執行置換檢驗檢視顯著性，僅當全模型透過檢驗時，才能依次對每一個約束軸執行檢驗判斷是否能夠保留。否則，RDA模型則從全域性上來講就不合理，本次的RDA結果不能再被使用，需要考慮更改分析方案（更換排序模型，或者再嘗試對資料進行有效的轉化後使用等）。
vegan包中執行置換檢驗的命令是anova()（很尷尬它就叫這個名字，R語言中方差分析的命令也是anova()……二者名稱一樣但用法完全不同，注意區分二者的用途）。關於此命令的詳情，可使用?anova.cca()檢視幫助（不要使用?anova()，不然彈出的將會是方差分析的命令幫助）。
以下以前文tb-RDA結果為例，展示RDA全模型及各約束軸的置換檢驗過程。
##置換檢驗
#所有約束軸的置換檢驗，基於 999 次置換檢驗
rda_tb_test <- anova(rda_tb, permutations = 999)
#或者使用
rda_tb_test <- anova.cca(rda_tb, step = 1000)
#各約束軸逐一檢驗，基於 999 次置換檢驗
rda_tb_test_axis <- anova(rda_tb, by = 'axis', permutations = 999)
#或者使用
rda_tb_test_axis <- anova.cca(rda_tb, by = 'axis', step = 1000)
若使用anova()，使用permutations引數指定置換次數；若使用anova.cca()，使用step引數設定產生參照分佈的樣方總數，step數值減1即為實際的置換次數。對於RDA全模型的置換檢驗，直接將前文所得tb-RDA結果“rda_tb”輸入至命令內即可；對於各約束軸，則還需在命令中新增by = 'axis'。此外，anova()還可以透過by = 'margin'檢驗協變數等，這裡不再多說，有需要請?anova.cca()參閱該命令的幫助文件。
結果就很明瞭了。首先全模型是顯著的，因此我們繼續檢驗每個約束軸，發現僅第一軸（RDA1）和第二軸（RDA2）是顯著的（預設顯著性閾值p<0.05），因此第三軸及以後的結果需要被剔除。
當然，必要時可以考慮做下p值校正。
對於上述示例結果，我們加入Bonferroni校正，使用p.adjust()完成。
#p 值校正（Bonferroni 為例）
rda_tb_test_axis$`Pr(>F)` <- p.adjust(rda_tb_test_axis$`Pr(>F)`, method = 'bonferroni')
結果如下所示。對於RDA1和RDA2兩個約束軸，結果仍然是可信的。
對於示例資料來講，根據檢驗結果，我們僅能在全模型中提取前兩個約束軸的資訊出來以進行更進一步的分析。不過，這樣的結果似乎已經是挺不錯的了。如上文所述，在校正了約束模型的R2後，RDA1和RDA2的解釋量總和為“0.43177+0.09571=0.52748”（即承載了響應變數矩陣52.75%的方差，達到了一半以上；再結合置換檢驗結果，即便不考慮其它的約束軸，約束模型所能解釋的方差佔比也是蠻可觀的），表明給定的環境變數對群落物種組成分佈具有相當高的解釋度。特別是第一軸，提示我們在解讀排序圖時，可能需要主要觀察樣方或物種沿第一約束軸的分佈規律，以及相關環境變數與約束軸、樣方及物種的相關性程度等。
事實上，即便後面的約束軸通過了檢驗，我們一般也不將它們列為考慮的範疇，畢竟它們的解釋量很低，強制將它們加入生態學模型解釋物種組成成因也可能會帶來混亂的結果。我們更期望的還是具有較高解釋量的約束軸能夠被保留，以及樣方或物種的排序分佈具有明顯趨勢的模型不要被捨棄等。
注意，非約束模型中排序軸取捨的常用方法不適用於約束模型。例如斷棍模型、凱撒格特曼規則（Kaiser-Guttman rule）等，我們可能在PCA分析中常用它們作為排序軸的選擇標準，但它們在RDA約束軸的取捨中意義不大。對於約束軸的取捨，還需透過置換檢驗獲得。然而，對於RDA的非約束殘差軸，斷棍模型等仍然適用。
非約束殘差軸的重要性確定
當我們發現RDA模型的承載方差比例較低，存在相當一部分比例的殘差未參與解釋時，可能需要審視非約束殘差軸。以下接著以上文tb-RDA結果為例，提取其中非約束殘差軸的特徵值，並分別依據斷棍模型和Kaiser-Guttman準則確定殘差軸。
##斷棍模型和 Kaiser-Guttman 準則確定殘差軸
#提取殘差特徵值
pca_eig <- rda_tb$CA$eig
#Kaiser-Guttman 準則
pca_eig[pca_eig > mean(pca_eig)]
#斷棍模型
n <- length(pca_eig)
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
# 繪製每軸的特徵根和方差百分比
par(mfrow = c(2, 1))
barplot(pca_eig, main = '特徵根', col = 'bisque', las = 2)
abline(h = mean(pca_eig), col = 'red')
legend('topright', '平均特徵根', lwd = 1, col = 2, bty = 'n')
barplot(t(cbind(100 * pca_eig/sum(pca_eig), bsm$p[n:1])), beside = TRUE, main = '% 變差', col = c('bisque', 2), las = 2)
legend('topright', c('% 特徵根', '斷棍模型'), pch = 15, col = c('bisque', 2), bty = 'n')
結果如下。綜合考慮，大約3個殘差非約束軸是比較重要的。若有必要，我們可能需要繪製前2-3個非約束軸的排序圖，檢視為何還有相當一部分方差未能解釋。
變數選擇（前向選擇）
約束排序中，有時解釋變數太多，需要想辦法減少解釋變數的數量。一般來講，減少解釋變數的數量可能有兩個並不衝突的原因：第一是尋求簡約的模型，利於我們對模型的解讀；第二是當解釋變數過多時會導致模型混亂，例如有些解釋變數之間可能存在較強的線性相關，即共線性問題，可能會造成迴歸係數不穩定。
方差膨脹因子
每個變數的共線性程度可以使用變數的方差膨脹因子（Variance Inflation Factor，VIF）度量，VIF是衡量一個解釋變數的迴歸係數的方差由共線性引起的膨脹比例。如果VIF值超過20，表示共線性很嚴重。實際上，VIF值超過10則可能就會有共線性的問題，需要處理。
vegan包中計算約束模型中變數VIF值的命令為vif.cca()。以上文tb-RDA結果為例計算VIF。
#計算方差膨脹因子
vif.cca(rda_tb)
可以看出有幾個變數的VIF值特別大（超過10甚至20），所以有必要剔除一些變數。
多元迴歸變數篩選通常有三種模式：前向選擇（forward selection）、後向選擇（backward selection）以及逐步選擇（stepwise selection，也稱雙向選擇，forward-backward selection）。其中前向選擇在RDA分析中最為常用，以下將以前向選擇為例簡介RDA模型中的變數選擇方法。
前向選擇
在R中，解釋變數的前向選擇通常使用vegan包中ordistep()函式、ordiR2step()函式，或者packfor包中的forward.sel()函式來完成。不同的命令執行前向選擇所依據的原理是不同的，以下分別簡介用法。
ordistep()
vegan包中ordistep()函式執行前向選擇過程大致如下（除了前向選擇，ordistep()同樣可以執行後向選擇和逐步選擇，本篇不再細說，詳細請使用?ordistep()檢視說明）。
（1）依次分別執行每個解釋變數與響應變數的RDA排序。
（2）選擇“最好”的解釋變數。對於ordistep()，根據置換檢驗中F統計量的顯著性水平（p值）選擇最好的變數，即最小的p值。如果遇到p值相等的情況，擁有最小AIC資訊準則（Akaike information criterion，AIC）的變數應該入選。最好的變數也就是最顯著的變數。
（3）尋找模型中第二個、第三個、第四個……解釋變數。上一步選取了只含有一個最好變數的模型，現在再加入某一個剩下的解釋變數。對於ordistep()，選擇偏RDA置換檢驗最小p值的模型，即表示新加入的變數對模型的貢獻最顯著。如果遇到p值相等的情況，具有最小AIC值的模型應該入選。
（4）這個過程一直持續，直到無顯著的解釋變數為止。即如果加入新變數的偏RDA置換檢驗顯著性p值大於或等於α（例如，以α=0.05為準），選擇過程即被終止。
在上文，透過計算計算方差膨脹因子，我們已知tb-RDA示例結果中的某些解釋變數最好加以剔除。以下使用ordistep()，展示在示例tb-RDA模型中執行前向選擇的過程。ordistep()可以直接將rda()的輸出作為輸入（此外，ordistep()也同樣可用於CCA模型，vegan包cca()的輸出）。由於原始的tb-RDA模型就已經是包含所有解釋變數的RDA全模型了，我們可直接在它基礎上繼續執行前向選擇。
#vegan 包 ordistep() 前向選擇，基於 999 次置換檢驗，詳情 ?ordistep()
rda_tb_forward_p <- ordistep(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), direction = 'forward', permutations = 999)
#簡要繪製雙序圖比較變數選擇前後結果
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = '原始模型，I 型標尺', display = c('wa', 'cn'))
plot(rda_tb_forward_p, scaling = 1, main = '前向選擇後，I 型標尺', display = c('wa', 'cn'))
#細節部分檢視
summary(rda_tb_forward_p, scaling = 1)
#比較選擇前後校正後 R2 的差異
RsquareAdj(rda_tb)$adj.r.squared
RsquareAdj(rda_tb_forward_p)$adj.r.squared
ordistep()中，rda(phylum_hel~1, env)意為從env（環境解釋變數資料集矩陣）中的第一個變數開始執行選擇；direction = 'forward'，前向選擇；permutations = 999，999次置換。
全模型中共計9個環境變數，經過ordistep()前向選擇後保留了其中的6個。儘管解釋變數的數量少了1/3，但是我們發現模型的總解釋量幾乎沒有改變，即變數選擇後在提供了簡約的模型的同時，並沒有犧牲解釋能力。
然後再根據排序圖，或者summary()統計，重新解讀模型，觀察比較變數選擇前後的差異；並重新校正模型R2以及檢驗約束軸，必要時增添p值校正，最終選擇合適的約束軸用於解釋生態學現象等。此處不再更多的展示，請大家自行測試。
但是ordistep()的選擇標準過於寬鬆，有時會選擇顯著但不包含任何變數的模型（誇大I類錯誤），或選擇包括過多變數的模型（誇大被解釋方差的量）。為了解決這個問題，可將Borcard終止準則（Blanchet等，2008）引入作為前向選擇的第二個終止原則，加強版的前向選擇程式，也就是下文將要介紹的ordiR2step()和forward.sel()。
在使用ordistep()時，若是發現結果似乎異常，例如下圖中的這種結果，變數選擇後的校正R2高於了全模型的校正R2？則推薦使用ordiR2step()或forward.sel()執行變數選擇過程（其實也推薦在一開始就直接使用這種加強版的前向選擇程式ordiR2step()或forward.sel()，忽略ordistep()）。
ordiR2step()
ordiR2step()的前向選擇原理與ordistep()類似，但引入全模型R2adj作為第二個終止原則（即增添了Borcard終止準則）：如果當前所選模型的R2adj達到或超過全模型的R2adj，或備選變數的偏RDA置換檢驗p值不顯著，則變數選擇將停止。
與上文一致，以下仍然以示例tb-RDA模型為例做演示。和ordistep()用法一致，除了前向選擇，ordiR2step()也可以執行後向選擇和逐步選擇，並且可以直接將rda()的輸出作為輸入（此外，ordiR2step()也同樣可用於CCA模型，vegan包cca()的輸出）。本文不再細說，詳細請使用?ordiR2step()檢視說明。
#vegan 包 ordiR2step() 前向選擇，基於 999 次置換檢驗，詳情 ?ordiR2step()
rda_tb_forward_r <- ordiR2step(rda(phylum_hel~1, env, scale = FALSE), scope = formula(rda_tb), R2scope = rda_adj, direction = 'forward', permutations = 999)
#簡要繪製雙序圖比較變數選擇前後結果
par(mfrow = c(1, 2))
plot(rda_tb, scaling = 1, main = '原始模型，I 型標尺', display = c('wa', 'cn'))
plot(rda_tb_forward_r, scaling = 1, main = '前向選擇後，I 型標尺', display = c('wa', 'cn'))
#細節部分檢視
summary(rda_tb_forward_r, scaling = 1)
#比較選擇前後校正後 R2 的差異
RsquareAdj(rda_tb)$adj.r.squared
RsquareAdj(rda_tb_forward_r)$adj.r.squared
ordiR2step()中，rda(phylum_hel~1, env)意為從env（環境解釋變數資料集矩陣）中的第一個變數開始執行選擇；rda_adj即為RDA全模型的校正R2，在上文中已經透過RsquareAdj()獲得，引數scope = formula(rda_tb)即指定全模型的校正R2作為判斷的依據；direction = 'forward'，前向選擇；permutations = 999，999次置換。
全模型中共計9個環境變數，經過ordistep()前向選擇後保留了其中的6個。儘管解釋變數的數量少了1/3，但是我們發現模型的總解釋量幾乎沒有改變，即變數選擇後在提供了簡約的模型的同時，並沒有犧牲解釋能力。
然後再根據排序圖，或者summary統計，重新解讀模型，觀察比較變數選擇前後的差異；並重新校正模型R2以及檢驗約束軸，必要時增添p值校正，最終選擇合適的約束軸用於解釋生態學現象等。此處不再更多的展示，請大家自行測試。
備註：在本示例中，ordiR2step()和ordistep()的所得結果是相同的，巧合了而已（本身資料集不大，解釋變數的數量也不多）。但在實際的資料分析中，二者的結果有時差別還是挺大的。
提取變數選擇後的排序座標，以此為例。
#提取或輸出變數選擇後的排序座標
#以 ordiR2step() 結果為例，以 I 型標尺為例，前四軸為例
#提取方式可參考上文，這裡透過 scores() 提取
#使用物種加權和計算的樣方得分
rda_tb_forward_r_site.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'wa')
write.table(data.frame(rda_tb_forward_r_site.scaling1), 'rda_tb_forward_r_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種變數（響應變數）得分
rda_tb_forward_r_sp.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'sp')
write.table(data.frame(rda_tb_forward_r_sp.scaling1), 'rda_tb_forward_r_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境變數（解釋變數）得分
rda_tb_forward_r_env.scaling1 <- scores(rda_tb_forward_r, choices = 1:4, scaling = 1, display = 'bp')
write.table(data.frame(rda_tb_forward_r_env.scaling1), 'rda_tb_forward_r_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
forward.sel()
由於forward.sel()和ordiR2step()的選擇原理相似，都是引入了Borcard終止準則後的加強版的選擇程式，所以一般情況下二者前向選擇的結果是一致的（儘管具有不同的引數設定邏輯）。
但是與ordiR2step()相比，forward.sel()目前僅適用於RDA（改良後的ordiR2step()既可用於RDA，vegan包rda()的輸出；也可用於CCA，vegan包cca()的輸出），且僅可用前向選擇（ordiR2step()同時支援前向、後向、逐步選擇），所以看起來似乎並沒有ordiR2step()好用……
而且forward.sel()所在的packfor包僅支援特定的R版本，所以本篇就不展示它了，大家有興趣可自行研究。
變數選擇中的注意事項
儘管變數選擇策略具有明顯的優勢，但一定切記不應盲目依賴自動選擇程式在迴歸模型中選擇相關的環境變數，因為可能得到生態學上不相關的模型，或者被忽視的其它變數組合（當前保留變數集中的部分解釋變數+某些被剔除的解釋變數）在某種程度上也可能比當前模型更具代表性（Legendre和Legendre，1998）。即僅透過統計學手段得到的最優變數集合可能並非代表了最具生態學意義的模型。
排序圖的視覺化示例
到這裡，對RDA分析的介紹基本就結束了。其實到了這一步，我們期間透過R2校正、置換檢驗、前向選擇步驟等，再結合自己的專業知識進行合理的解讀和評估後，也差不多得到了一個合理的RDA模型了。當然，我的意思可不是指這裡的示例資料的RDA分析已經得到了一個理想模型了，而是指在實際情況的分析中，我們透過多方考慮所得的最終模型，這個模型除了在統計學上合理外，在解釋生物學意義中也是較為理想的。這個最終的模型和最初直接由rda()所得模型相比，可能差異巨大。不過到了這裡，想必大家也能理解了為什麼最初所得的RDA模型最好不要直接使用的原因了。
那麼最後，再對RDA排序圖的視覺化方案做個簡介。
plot()作圖
在前文，已經簡單地使用plot()函式（實為plot.cca()）繪製了RDA排序圖用於展示。配合plot()函式使用，除了前述scores()（提取座標）、arrows()（繪製箭頭）等命令外，還有point()（繪製散點）、text()（控制標籤字型）等。通常，這幾個命令相互配合使用，一個美觀的RDA排序圖即可呈現。
我們使用plot()，對上文中經過解釋變數前向選擇後所得的簡約tb-RDA模型“rda_tb_forward_r”（僅使用它作為示例展示視覺化方法，並沒有說它是最終的模型，儘管它是透過多步計算所得，但本文中並沒有結合生態學意義去評估或解釋它），繪製雙序圖展示其中的前兩個約束軸。不過在作圖前，首先需要確認以下資訊。
（1）根據置換檢驗，我們已知前兩個約束軸是顯著的，可以使用。
（2）根據“RsquareAdj(rda_tb_forward_r)$adj.r.squared”可知簡約模型的校正後R2為“0.5814848”，然後再根據“summary(rda_tb_forward_r)”的結果可知約束軸RDA1和RDA2的承載方差分別佔約束模型總方差的“0.73786”和“0.168480”，那麼由此推斷約束軸RDA1和RDA2的解釋量分別約為“0.5814848×0.73786 = 42.91%”和“0.5814848×0.168480 = 9.80%”。解釋量可觀。
（3）根據所關注的重點，確定使用哪種標尺。例如這次想重點觀測樣方間的關係及其環境梯度，那麼就以I型標尺為例展示。
（4）可首先使用“plot(rda_tb_forward_r, scaling = 1, display = c('wa', 'cn'))”檢視簡約作圖結果，排序空間內的樣方分佈是否區分明顯，環境變數對排序空間的貢獻是否達到了預期等。
確認都沒問題後，就可以作圖了，一個簡單的示例如下。
#plot() 作圖示例，詳情 ?plot.cca
#以前向選擇後的簡約模型 rda_tb_forward_r 為例，展示前兩軸，I 型標尺，雙序圖，預設使用物種加權和計算的樣方座標
#注意 RDA 軸的解釋率，需要手動結合校正後的 R2 計算下
png('rda_test1.png', width = 600, height = 600, res = 100, units = 'px')
plot(rda_tb_forward_r, display = c('wa', 'cn'), choices = 1:2, scaling = 1, type = 'n', main = 'RDA雙序圖，I型標尺', xlab = 'RDA1 (42.91%)', ylab = 'RDA2 (9.80%)')
text(rda_tb_forward_r, display = 'cn', choices = 1:2, scaling = 1, col = 'blue', cex = 0.8)
points(rda_tb_forward_r, display = 'wa', choices = 1:2, scaling = 1, pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
dev.off()
其中，plot()控制全域性；text()定義了環境變數向量的屬性；points()定義了樣方座標散點的屬性，由於示例資料中27個樣方共分為3組，我們使用3種顏色表示出。
ggplot2作圖
對於習慣ggplot2做圖方式的同學們來說，我們只需要將樣方排序座標、物種排序座標以及環境因子的排序座標分別提取出構建資料集後，即可使用ggplot2繪製排序圖了。
#以前向選擇後的簡約模型 rda_tb_forward_r 為例，展示前兩軸，I 型標尺，雙序圖，預設使用物種加權和計算的樣方座標
#提取樣方和環境因子排序座標，前兩軸，I 型標尺
rda_tb_forward_r.scaling1 <- summary(rda_tb_forward_r, scaling = 1)
rda_tb_forward_r.site <- data.frame(rda_tb_forward_r.scaling1$sites)[1:2]
rda_tb_forward_r.env <- data.frame(rda_tb_forward_r.scaling1$biplot)[1:2]
#讀取樣本分組資料（附件“group.txt”）
group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#合併樣本分組資訊，構建 ggplot2 作圖資料集
rda_tb_forward_r.site$sample <- rownames(rda_tb_forward_r.site)
rda_tb_forward_r.site <- merge(rda_tb_forward_r.site, group, by = 'sample')
rda_tb_forward_r.env$sample <- rownames(rda_tb_forward_r.env)
#ggplot2 作圖
library(ggplot2)
p <- ggplot(rda_tb_forward_r.site, aes(RDA1, RDA2)) +
geom_point(aes(color = group)) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5), legend.key = element_rect(fill = 'transparent')) +
labs(x = 'RDA1 (42.91%)', y = 'RDA2 (9.80%)', title = 'RDA雙序圖，I型標尺', color = '') +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = rda_tb_forward_r.env, aes(x = 0, y = 0, xend = RDA1,yend = RDA2), arrow = arrow(length = unit(0.1, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = rda_tb_forward_r.env, aes(RDA1 * 1.1, RDA2 * 1.1, label = sample), color = 'blue', size = 3)
p
#ggsave('rda_test2.pdf', p, width = 5, height = 4)
ggsave('rda_test2.png', p, width = 5, height = 4)
備註：若覺得箭頭長度太長或太短，影響觀測，可以同比例縮放所有箭頭的長度後再繪製出來，這種做法是被允許的（上文中在介紹排序圖時也提到了這種做法的可行性；而且能看到上述plot()函式的默認出圖，也是這樣做的，函式預設使用一個內部常數對結果重新標定，這時變數向量的箭頭長度不完全等於原始值，而是等比例原始值）。
