某研究對三種環境中的土壤進行了取樣（環境A、B、C，各環境下采集9個土壤樣本），結合二代測序觀測了土壤細菌群落物種組成，並對多種土壤理化指標進行了測量。獲得兩個資料集。
（1）檔案“otu_table.txt”為透過16S測序所得的細菌OTU丰度表格，下文統稱為“物種”，儘管OTU並不完全等於物種。
（2）檔案“env_table.txt”記錄了多種土壤理化指標資訊，即環境資料。
接下來期望透過db-RDA分析這些資料，用以分析群落物種組成與環境梯度的關係。
vegan包中db-RDA的初步計算
首先展示db-RDA的初步計算過程，對於內容詳解見下文。
手動分步計算
按照上述db-RDA原理，不妨先從頭走一遍，以便加深理解。
先計算PCoA，再將結果作為RDA的輸入，最後被動新增物種得分。方法大致如下。
library(vegan)
#讀入物種資料，以細菌 OTU 水平丰度表為例
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
#讀取環境資料
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
##根據原理一步步計算 db-RDA
#計算樣方距離，以 Bray-curtis 距離為例，詳情 ?vegdist
dis_bray <- vegdist(otu, method = 'bray')
#或者直接使用現有的距離矩陣，這裡同樣為 Bray-curtis 距離
dis_bray <- as.dist(read.delim('bray_distance.txt', row.names = 1, sep = '\t', check.names = FALSE))
#PCoA 排序，這裡透過 add = TRUE校正負特徵值，詳情 ?cmdscale
pcoa <- cmdscale(dis_bray, k = nrow(otu) – 1, eig = TRUE, add = TRUE)
#提取 PCoA 樣方得分（座標）
pcoa_site <- pcoa$point
#db-RDA，環境變數與 PCoA 軸的多元迴歸
#透過 vegan 包的 RDA 函式 rda() 執行，詳情 ?rda
db_rda <- rda(pcoa_site, env, scale = FALSE)
#被動擬合物種得分
v.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) – 1)
db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) – 1)
db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
#先作圖展示下，詳情 ?plot.cca
#樣方展示為點，物種暫且展示為“+”，環境變數為向量
par(mfrow = c(1, 2))
plot(db_rda, type = 'n', display = c('wa', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，雙序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
plot(db_rda, type = 'n', display = c('wa', 'sp', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，三序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'sp', pch = 3, col = 'gray', cex = 1)
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
capscale()函式
即上述分步過程相比，capscale()是vegan中打包好的db-RDA執行函式，可以直接一步出結果。平常說的“CAP（Canonical analysis of principal coordinates）”就是指它，就是db-RDA。
#或者，capscale() 提供了直接執行的方法，詳情 ?capscale
#輸入原始資料，指定距離型別，例如樣方距離以 Bray-curtis 距離為例
#計算結果中包含物種得分
db_rda <- capscale(otu~., env, distance = 'bray', add = TRUE)
#或者直接輸入已經計算/或讀入的距離測度
db_rda <- capscale(dis_bray~., env, add = TRUE)
#但是這種情況下，物種得分丟失，需要被動新增
v.eig <- t(otu) %*% db_rda$CCA$u/sqrt(nrow(otu) – 1)
db_rda$CCA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
v.eig <- t(otu) %*% db_rda$CA$u/sqrt(nrow(otu) – 1)
db_rda$CA$v <- decostand(v.eig, 'normalize', MARGIN = 2)
#先作圖展示下
#樣方展示為點，物種暫且展示為“+”，環境變數為向量
par(mfrow = c(1, 2))
plot(db_rda, type = 'n', display = c('wa', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，雙序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
plot(db_rda, type = 'n', display = c('wa', 'sp', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，三序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'sp', pch = 3, col = 'gray', cex = 1)
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
dbrda()函式
dbrda()是vegan中打包好的db-RDA另一執行函式，直接一步出結果。dbrda()和capscale()在計算方法上存在小差異，但返回值基本是一致的。
#以及 dbrda()，詳情 ?dbrda
#輸入原始資料，指定距離型別，例如樣方距離以 Bray-curtis 距離為例
db_rda <- dbrda(otu~., env, distance = 'bray', add = TRUE)
#或者直接輸入已經計算/或讀入的距離測度
db_rda <- dbrda(dis_bray~., env, add = TRUE)
#兩種情況下，dbrda() 均預設不計算物種得分，可透過 sppscores() 新增，詳情 ?sppscores
sppscores(db_rda) <- otu
#先作圖展示下
#樣方展示為點，物種暫且展示為“+”，環境變數為向量
par(mfrow = c(1, 2))
plot(db_rda, type = 'n', display = c('wa', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，雙序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
plot(db_rda, type = 'n', display = c('wa', 'sp', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，三序圖')
points(db_rda, choices = 1:2, scaling = 1, display = 'sp', pch = 3, col = 'gray', cex = 1)
points(db_rda, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(db_rda, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
關於在db-RDA中使用歐氏距離或卡方距離
上述均以群落分析中常用的Bray-curtis距離作演示。
基於歐式距離的db-RDA與常規RDA的比較
由概念還可知，若在db-RDA中使用歐式距離，將得到和常規RDA一樣的結果。
#使用原始物種多度資料執行 RDA，推薦先作個 Hellinger 轉化
hel <- decostand(otu, method = 'hellinger')
#分別執行基於 Hellinger 轉化後物種資料的 tb-RDA 以及與使用歐氏距離的 db-RDA
tb_rda_test <- rda(hel~., env, scale = FALSE)
db_rda_test <- capscale(hel~., env, distance = 'euclidean')
par(mfrow = c(1, 2))
plot(tb_rda_test, scaling = 1, main = 'tb-RDA三序圖')
plot(db_rda_test, scaling = 1, main = 'db-RDA三序圖')
可看到二者結果是一致的。
注：不要在意兩張圖“方向相反”，排序圖的方向是沒有任何意義的，以排序空間中樣方、物種和環境變數的相對位置為準。
基於卡方距離的db-RDA與CCA的比較
順便再比較一下使用卡方（χ2）距離的db-RDA和CCA的區別。
#但是基於卡方距離的 db-RDA，和 CCA 並不完全相同
cca_test <- cca(otu~., env, scale = FALSE)
chi <- decostand(otu, method = 'chi.square') #先做卡方標準化，再計算歐氏距離，即可得卡方距離
db_rda_test <- capscale(chi~., env, distance = 'euclidean')
par(mfrow = c(2, 2))
plot(cca_test, scaling = 1, main = 'db-RDA三序圖')
plot(db_rda_test, scaling = 1, main = 'CCA三序圖')
plot(cca_test, display = c('wa', 'cn'), scaling = 1, main = 'db-RDA雙序圖')
plot(db_rda_test, display = c('wa', 'cn'), scaling = 1, main = 'CCA雙序圖')
二者的樣方得分和環境變數得分是一致的，區別在於物種得分。
從計算過程來看，CCA是CA與多元迴歸的結合，因此它和基於χ2距離的db-RDA（基於χ2距離的PCoA與多元迴歸的結合）結果的相似性是可以預料到的。對於物種得分的不同，CCA在計算過程中用加權平均法對物種排序，而db-RDA中物種排序是在排序結束後被動加入的，方法的不同帶來了差異。
db-RDA結果的初步解讀
其中可能capscale()（如上所述，俗稱CAP）最為常用，下文就以它的結果為例作簡介。
db-RDA結果總概
#db-RDA 結果解讀，以 capscale() 函式結果為例，簡介
db_rda <- capscale(otu~., env, distance = 'bray', add = TRUE)
#檢視統計結果資訊，以 I 型標尺為例
db_rda.scaling1 <- summary(db_rda, scaling = 1)
db_rda.scaling1
db-RDA排序圖
關於更美觀的視覺化方案，可參考上文示例，或者文末的ggplot2方法。
#作圖檢視排序結果，詳情 ?plot.cca
#三序圖，包含 I 型標尺和 II 型標尺，樣方座標展示為使用物種加權計算的樣方得分
par(mfrow = c(1, 2))
plot(db_rda, scaling = 1, main = 'I 型標尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling1 <- scores(db_rda, 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(db_rda, scaling = 2, main = 'II 型標尺', display = c('wa', 'sp', 'cn'))
rda_sp.scaling2 <- scores(db_rda, choices = 1:2, scaling = 2, display = 'sp')
#arrows(0, 0, rda_sp.scaling2[ ,1], rda_sp.scaling2[ ,2], length = 0, lty = 1, col = 'red')
#隱藏物種，以 I 型標尺為例展示雙序圖
#比較分別使用物種加權計算的樣方座標以及擬合的樣方座標的差異
par(mfrow = c(1, 2))
plot(db_rda, scaling = 1, main = 'I 型標尺，加權', display = c('wa', 'cn'))
plot(db_rda, scaling = 1, main = 'I 型標尺，擬合', display = c('lc', 'cn'))
主要結果內容提取
若有需要，考慮在結果中將需要的資訊提取出來。
#scores() 提取排序得分（座標），以 I 型標尺為例，前四軸為例
#使用物種加權和計算的樣方得分
db_rda_site.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'wa')
#物種變數（響應變數）得分
db_rda_sp.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'sp')
#環境變數（解釋變數）得分
db_rda_env.scaling1 <- scores(db_rda, choices = 1:4, scaling = 1, display = 'bp')
#或者在 summary() 後提取，以 I 型標尺為例，前四軸為例
db_rda.scaling1 <- summary(db_rda, scaling = 1)
#使用物種加權和計算的樣方得分
db_rda_site.scaling1 <- db_rda.scaling1$site[ ,1:4]
#物種
db_rda_sp.scaling1 <- db_rda.scaling1$species[ ,1:4]
#環境
db_rda_env.scaling1 <- db_rda.scaling1$biplot[ ,1:4]
#若需要輸出在本地
#樣方
write.table(data.frame(db_rda_site.scaling1), 'db_rda_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種
write.table(data.frame(db_rda_sp.scaling1), 'db_rda_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境
write.table(data.frame(db_rda_env.scaling1), 'db_rda_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#不建議直接在原始資料集中提取，因為這樣提取的座標數值未經標尺縮放處理，不利於反映生物學問題
#db_rda$CCA$u[ ,1:4]
#db_rda$CCA$v[ ,1:4]
#db_rda$CCA$biplot[ ,1:4]
#RDA 基於多元迴歸，若對典範特徵係數（即每個解釋變數與每個約束軸之間的迴歸係數）感興趣
#coef() 提取典範係數
rda_coef <- coef(db_rda)
關於RDA的總變差（即響應變數矩陣的總變差能夠被解釋變數解釋的部分，如果用比例表示，即等同於多元迴歸的R2，也就是約束軸的總解釋量）以及各約束軸解釋量等，不建議直接在原始結果中獲取。初始所得的RDA模型需要經過校正和檢驗後才可使用，否則可能會帶來較大的偏差，原因可參考前文。
R2經過校正後通常會降低，這是必然的。
#RsquareAdj() 提取 R2，詳情 ?RsquareAdj()
r2 <- RsquareAdj(db_rda)
db_rda_noadj <- r2$r.squared #原始 R2
db_rda_adj <- r2$adj.r.squared #校正後的 R2
對於示例資料的db-RDA排序結果，透過上文，我們已經得知了響應變數矩陣的總變差為4.111，未校正前的約束軸解釋的變差總和為2.476（佔比0.6024，即未校正前的R2）；那麼根據校正後的R2（0.39186），可還原校正R2後的約束軸解釋的變差總和約為“4.111×0.39186=1.6109”。
對於各約束軸校正後的解釋量。例如，已知RDA1的解釋變差佔約束軸總解釋變差的0.4427，由於校正後的R2為0.39186，那麼可推斷校正R2後的RDA1解釋率為“0.4427×0.39186=0.1735”；同理，校正R2後RDA2解釋率“0.2532×0.39186=0.0992”；其餘約束軸的解釋量以此類推。
#概念如上所述，如下為程式碼實現
#關於約束軸承載的特徵值或解釋率，應當在 R2 校正後重新計算
db_rda_exp_adj <- db_rda_adj * db_rda$CCA$eig/sum(db_rda$CCA$eig)
db_rda_eig_adj <- db_rda_exp_adj * db_rda$tot.chi
約束軸的置換檢驗及p值校正
我們看到本示例的約束軸解釋率還算可以（以前兩軸為例，解釋量0.1735+0.0992=0.2727，不算很低），那麼，這些約束軸是否有效呢？
對於置換檢驗，首先檢驗全模型，通過後再檢驗各約束軸。
#置換檢驗
#所有約束軸的置換檢驗，即全域性檢驗，基於 999 次置換，詳情 ?anova.cca
db_rda_test <- anova.cca(db_rda, permutations = 999)
#各約束軸逐一檢驗，基於 999 次置換
db_rda_test_axis <- anova.cca(db_rda, by = 'axis', permutations = 999)
#p 值校正（Bonferroni 為例）
db_rda_test_axis$`Pr(>F)` <- p.adjust(db_rda_test_axis$`Pr(>F)`, method = 'bonferroni')
本示例中僅前兩軸顯著，也就是說，第三及之後的約束軸理論上不能再用作分析了。
注：如果檢驗未透過，在考慮更換其它方法之前，不妨先試試剔除一些環境變數，因為也有可能是某些環境變數與群落特徵之間缺乏較好的線性關係所致。
非約束殘差軸的重要性確定
此外，當db-RDA約束軸承載的變差比例較低時，可能需要審視非約束殘差軸。
#斷棍模型和 Kaiser-Guttman 準則幫助確定重要的殘差軸
pcoa_eig <- db_rda$CA$eig
n <- length(pcoa_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
barplot(t(cbind(100 * pcoa_eig/sum(pcoa_eig), bsm$p[n:1])), beside = TRUE, main = '% 變差', col = c('orange', 'bisque'), las = 2)
abline(h = mean(100 * pcoa_eig/sum(pcoa_eig)), col = 'red')
legend('topright', c('% 特徵根', '斷棍模型', '平均特徵根'), pch = c(15, 15, NA), col = c('orange', 'bisque', 'red'), lwd = c(NA, NA, 1), bty = 'n')
結果指示可能還有1-2個殘差非約束軸是比較重要的。若有必要，可繪製前1-2個非約束軸的排序圖，檢視為何還有相當一部分變差未能解釋。
db-RDA的變數選擇（前向選擇為例）
RDA中，有時解釋變數太多，需要想辦法減少解釋變數的數量。減少解釋變數的數量可能有兩個並不衝突的原因：第一是尋求簡約的模型，利於我們對模型的解讀；第二是當解釋變數過多時會導致模型混亂，例如有些解釋變數之間可能存在較強的線性相關，即共線性問題，可能會造成迴歸係數不穩定。
方差膨脹因子（Variance Inflation Factor，VIF）顯示環境變數間共線性很嚴重。有幾個變數的VIF值特別大（超過10甚至20），所以有必要剔除一些變數。
#計算方差膨脹因子，詳情 ?vif.cca
vif.cca(db_rda)
#前向選擇，以 ordiR2step() 的方法為例，基於 999 次置換檢驗，詳情 ?ordiR2step
db_rda_forward_pr <- ordiR2step(capscale(otu~1, env, distance = 'bray', add = TRUE), scope = formula(db_rda), R2scope = TRUE, direction = 'forward', permutations = 999)
#以 db_rda 和 db_rda_forward_pr 為例，簡要繪製雙序圖比較變數選擇前後結果
par(mfrow = c(1, 2))
plot(db_rda, scaling = 1, main = '原始模型，I 型標尺', display = c('wa', 'cn'))
plot(db_rda_forward_pr, scaling = 1, main = '前向選擇後，I 型標尺', display = c('wa', 'cn'))
#細節部分檢視
summary(db_rda_forward_pr, scaling = 1)
#比較選擇前後校正後 R2 的差異，詳情 ?RsquareAdj
#可以看到變數選擇後，儘管去除了很多環境變數，但總 R2 並未損失很多
RsquareAdj(db_rda)$adj.r.squared
RsquareAdj(db_rda_forward_pr)$adj.r.squared
#所有約束軸的全域性檢驗，999 次置換，詳情 ?anova.cca
db_rda_forward_pr_test <- anova.cca(db_rda_forward_pr, permutations = 999)
#各約束軸逐一檢驗，999 次置換
db_rda_forward_pr_test_axis <- anova.cca(db_rda_forward_pr, by = 'axis', permutations = 999)
#p 值校正（Bonferroni 為例）
db_rda_forward_pr_test_axis$`Pr(>F)` <- p.adjust(db_rda_forward_pr_test_axis$`Pr(>F)`, method = 'bonferroni')
結果顯示前兩軸顯著，我們考慮將簡約db-RDA模型的前兩軸座標提取出，以便用於後續分析或視覺化。
#提取或輸出變數選擇後的排序座標，以 I 型標尺為例，前兩軸為例
#提取方式可參考上文，這裡透過 scores() 提取
#使用物種加權和計算的樣方得分
db_rda_forward_pr_site.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'wa')
write.table(data.frame(db_rda_forward_pr_site.scaling1), 'db_rda_forward_pr_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種變數（響應變數）得分
db_rda_forward_pr_sp.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'sp')
write.table(data.frame(db_rda_forward_pr_sp.scaling1), 'db_rda_forward_pr_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境變數（解釋變數）得分
db_rda_forward_pr_env.scaling1 <- scores(db_rda_forward_pr, choices = 1:2, scaling = 1, display = 'bp')
write.table(data.frame(db_rda_forward_pr_env.scaling1), 'db_rda_forward_pr_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
此外，如前所述，儘管變數選擇策略具有明顯的優勢，但一定切記不應盲目依賴自動選擇程式在迴歸模型中選擇相關的環境變數，即僅透過統計學手段得到的最優變數集合可能並非具有很好的生物學意義。
附帶一個ggplot2視覺化示例
上文已經展示了透過內建的plot()函式，可繪製排序圖。
同樣地，也可在匯出座標後，透過其它作圖包（如ggplot2）繪製。
#以前向選擇後的簡約模型 db_rda_forward_pr 為例作圖展示前兩軸
#plot(db_rda_forward_pr) 方法可參考上文
#下面是 ggplot2
#提取樣方和環境因子排序座標，前兩軸，I 型標尺
db_rda_forward_pr.scaling1 <- summary(db_rda_forward_pr, scaling = 1)
db_rda_forward_pr.site <- data.frame(db_rda_forward_pr.scaling1$sites)[1:2]
db_rda_forward_pr.env <- data.frame(db_rda_forward_pr.scaling1$biplot)[1:2]
#手動新增分組
db_rda_forward_pr.env$name <- rownames(db_rda_forward_pr.env)
db_rda_forward_pr.site$name <- rownames(db_rda_forward_pr.site)
db_rda_forward_pr.site$group <- c(rep('A', 9), rep('B', 9), rep('C', 9))
#計算校正 R2 後的約束軸解釋率
exp_adj <- RsquareAdj(db_rda_forward_pr)$adj.r.squared * db_rda_forward_pr$CCA$eig/sum(db_rda_forward_pr$CCA$eig)
rda1_exp <- paste('RDA1:', round(exp_adj[1]*100, 2), '%')
rda2_exp <- paste('RDA1:', round(exp_adj[2]*100, 2), '%')
#ggplot2 作圖
library(ggplot2)
p <- ggplot(db_rda_forward_pr.site, aes(CAP1, CAP2)) +
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'), legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) +
labs(x = rda1_exp, y = rda2_exp) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = db_rda_forward_pr.env, aes(x = 0, y = 0, xend = CAP1, yend = CAP2), arrow = arrow(length = unit(0.2, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = db_rda_forward_pr.env, aes(CAP1 * 1.1, CAP2 * 1.1, label = name), color = 'blue', size = 3)
p
