變差分解概述
以兩組解釋變數為例闡述，下圖Venn圖概括了響應變數的總變差組成。
當同時存在兩組解釋變數A和B時，圖中響應變數Y的總變差可以分解為以下幾個部分：
[a+b]變數A能夠解釋的總變差（變數A的邊際效應，marginal effect，即當約束模型中僅存在變數A這一個解釋變數時，約束模型所能夠被解釋的總變差）；
[b+c]變數B能夠解釋的總變差；
[a]僅能被變數A所單獨解釋的變差（變數A的區域性效應，partial effect，即將變數B作為協變數處理時，變數A所能解釋的變差部分）；
[c]僅能被變數B所單獨解釋的變差；
[b]變數A和B共同解釋的變差（由於變數A和B具有相關性，無法歸其所屬）；
[d]變數A和B均無法解釋的變差，即殘差。
常規RDA的變差分解
在RDA中響應變數的總變差=總方差，且早期變差分解也僅能在RDA中使用，因此RDA的變差分解也常稱為方差分解（variance partitioning）。
以上述兩組解釋變數的變差分解為例，若量化其中 [a]、[b]、[c]、[d]四部分的變差，必須執行3次RDA分析。變差分解的概念步驟如下。
（1）如果有必要，執行變數選擇（例如對所有解釋變數執行前向選擇），只保留顯著的變數。
（2）執行RDA分析，獲得各解釋變數所能解釋的變差。對應於上圖，即：
單獨以A作為解釋變數進行Y的RDA分析，可以獲得[a+b]部分的值；
單獨以B作為解釋變數進行Y的RDA分析，可以獲得[b+c]部分的值；
將A和B一起作為解釋變數進行Y的RDA分析，可以獲得[a+b+c]部分的值。
（3）R2校正也是不可或缺的。
如果每組解釋變數（變數集）包含相同數量的變數（例如，兩個土壤和兩個氣候變數），則組間解釋的變差是可以直接比較的；若變數數量不一致，則需要在比較前校正R2。同時由於隨機相關的存在，R2一般會伴隨解釋變數數量的增加或樣方數的減少而增大，解釋變數的累加會造成被解釋變差表面上膨脹。考慮到這些原因，我們通常會計算上述3個RDA分析的校正後R2。
（4）透過減法計算獲得各部分校正後的變差：
[a]adj = [a+b+c]adj – [b+c]adj；
[c]adj = [a+b+c]adj – [a+b]adj；
[b]adj = [a+b]adj – [a]adj = [b+c]adj – [c]adj；
[d]adj = 1 – [a+b+c]adj。
上述3個RDA分析都可以應用置換檢驗。
同上文所述，[a]為變數A的區域性效應（partial effect），即將變數B作為協變數處理時，變數A所能解釋的變差部分；[c]的理解同理。因此對於[a]和[c]部分的計算和檢驗也可以透過偏RDA分析進行。由於目前沒有公式可以直接計算偏RDA的校正R2，但上面描述的減法過程可以避免這個問題，Peres-Neto等（2006）也是用這種方法獲得各部分變差的無偏估計。
同上文所述，[b]由變數A、B所共享，歸因於變數A和B存在某種相關性，無法明確[b]的所屬，因此[b]部分不能直接估計和檢驗。同時可知，若變數A和B互相獨立（不相關、正交），那麼二者將不存在共同解釋的部分，即此時[b]將為0。
在某些情況下，獲得的校正R2可能是負值，表明解釋變數對響應變數的解釋程度甚至不及隨機變數的解釋程度，該解釋變數的存在意義不大。在對分析結果作生態解釋時，負的校正R2可以忽略（當作0看待）。
在R中，RDA的變差分解可使用vegan包中varpart()函式來實現。
##常規 RDA 的變差分解
#資料的網盤連結：https://pan.baidu.com/s/18xTqAZPa2Al2sDQpS1ysFA
#讀取物種資料，細菌門水平丰度表
phylum <- read.delim('phylum_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
#推薦在對純物種多度資料執行 PCA、RDA 等線性排序方法前，執行 Hellinger 轉化
phylum_hel <- decostand(phylum, method = 'hellinger')
#讀取環境變數
env <- read.delim('env_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
library(vegan)
#執行變差分解，詳情 ?vegdist
#以兩組（pH 為一組，5 種碳氮磷鉀指標為另一組）環境變數為例
rda_vp <- varpart(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')])
rda_vp
plot(rda_vp, digits = 2, Xnames = c('pH', 'CNPK'), bg = c('blue', 'red'))
這裡選擇了6個環境變數（前文RDA分析中，前向選擇後保留了這6種環境變數，此處拿來展示），根據環境變數的屬性，將它們劃分為兩組：一組為土壤化學屬性，pH；另一組則包括了土壤碳氮磷鉀含量，DOC、SOM、AP、AK、NH4。對這兩組環境變數執行變差分解後，結果如下。
其中，pH即對應於上圖中的“a+b”部分，解釋變差約13.2%，其中單獨解釋部分為“a”，解釋變差約1%；同理，DOC、SOM、AP、AK、NH4這5個的合集即對應於上圖中的“b+c”部分，解釋變差約57.1%，其中單獨解釋部分為“c”，解釋變差約44.9%；兩組環境變數共同解釋部分即為“b”部分，解釋變差約12.2%；兩組環境變數均不能解釋的部分為“d”，約41.9%。
#再以四組為例，碳、氮、磷、鉀指標間比較
rda_vp <- varpart(phylum_hel, env[c('DOC', 'SOM')], env['NH4'], env['AP'], env['AK'])
rda_vp
plot(rda_vp, digits = 2, Xnames = c('C', 'N', 'P', 'K'), bg = c('blue', 'red', 'green', 'yellow'))
在前文RDA中，我們發現在經過變數的前向選擇後，“TC”這個環境變數被剔除了。通常由於有些解釋變數之間可能存在較強的線性相關，即共線性問題，可能會造成迴歸係數不穩定，所以前向選擇程式一般會將這些變數剔除。那麼，對於“TC”這個環境變數來講，是否是因為這個原因呢？
#結合前文 RDA 中的示例分析
#檢視前向選擇中被剔除的環境變數“TC”，與這 6 個被保留的環境變數之間解釋變差的“共享程度”
rda_vp <- varpart(phylum_hel, env['TC'], env[c('pH', 'DOC', 'SOM', 'AP', 'AK', 'NH4')])
plot(rda_vp, digits = 2, Xnames = c('TC', 'forward_env'), bg = c('blue', 'red'))
事實的確如此，TC所解釋的變差，均與這6個被保留的（前向選擇後所保留的）環境變數之間產生交集，即存在非常強的共線性問題。而對於TC本身來講，它所能單獨解釋的變差部分為0，所以前向選擇程式自動將它剔除了。
還可以透過置換檢驗，檢驗各解釋變差部分的顯著性。例如，對pH所解釋變差部分的置換檢驗如下所示。
##置換檢驗確定解釋變差的顯著性
#需結合 RDA 執行式，vegan 中為 rda()
#解釋變差的置換檢驗，以 pH 所能解釋的全部變差為例；999 次置換，詳情 ?anova.cca
anova.cca(rda(phylum_hel, env['pH']), permutations = 999)
#若考慮 pH 單獨解釋的變差部分，需將其它變數作為協變數；999 次置換
anova.cca(rda(phylum_hel, env['pH'], env[c('DOC', 'SOM', 'AP', 'AK', 'NH4')]), permutations = 999)
注意觀察，存在不顯著的解釋變差部分嗎？對於這部分，應當如何處理呢？
db-RDA（CAP）的變差分解
##db-RDA 的變差分解
#資料的網盤連結：https://pan.baidu.com/s/1itLEWnV_sLpKohHleMwggA
#讀入物種資料，以細菌 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)
#計算樣方距離，以 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))
library(vegan)
#以兩組環境變數（DOC、AP+AK）為例，執行變差分解，詳情 ?varpart
#輸入 varpart() 的為已經計算好的 Bray-curtis 距離測度
#引數 add=T 意為校正 PCoA 的負特徵根，這個最好和先前計算 db-RDA 的執行式保持一致（即二者要麼都校正，要麼都不校正）
db_rda_vp <- varpart(dis_bray, env['DOC'], env[c('AP', 'AK')], add = TRUE)
db_rda_vp
plot(db_rda_vp, digits = 2, Xnames = c('DOC', 'AP+AK'), bg = c('blue', 'red'))
對於該示例，DOC（“a+b”部分）解釋了總變差的14.5%，其中單獨（“a”部分）解釋的比例為4.2%；AP+AK單獨（“c”部分）解釋了9.4%，和DOC共同（“b”部分）解釋了10.3%；DOC+AP+AK未能（“d”部分）解釋的變差佔比76.0%。
##置換檢驗確定解釋變差的顯著性
#需結合 db-RDA 執行式，vegan 中提供了 capscale()、dbrda() 等函式可直接執行 db-RDA
#下式均以 capscale() 為例
#解釋變差的置換檢驗，以 DOC 所能解釋的全部變差為例；999 次置換，詳情 ?anova.cca
anova.cca(capscale(dis_bray~DOC, env, add = TRUE), permutations = 999)
#若考慮 DOC 單獨解釋的變差部分，需將其它變數作為協變數；999 次置換
anova.cca(capscale(dis_bray~DOC+Condition(AP+AK), env, add = TRUE), permutations = 999)
CCA的變差分解
由於在CCA中，響應變數的總變差透過慣量（而非方差）表徵，因此“方差分解”這個名稱最好不要在CCA中使用，還是使用“變差分解”為好。
早期由於計算方法的原因，CCA的變差分解難以實現。例如在RDA中常用的Ezekiel校正法在CCA中不適用；且RDA變數選擇中的常用方法，如Borcard終止準則（Blanchet等，2008）也不能很好地用於CCA。這些阻礙了CCA分析中正確估計解釋變數能夠解釋的變差比例（即無法準確計算校正後R2），也過度誇大了I類錯誤，導致了無法對CCA模型執行變差分解。
後來隨著技術的更新，一系列的方法被引入。同樣以R包vegan為例，RsquareAdj()已能夠對CCA的R2實現校正（早期是不能的，僅能用於提取CCA的R2但不能校正），ordiR2step()也已將Borcard終止準則引入至CCA的變數前向選擇中，類似地，varpart()函式也可用於執行CCA的變差分解。
##CCA 的變差分解
#資料的網盤連結：https://pan.baidu.com/s/1UWlvIUWrnl4gyoEhpCGVtw
#讀入物種資料，以細菌 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)
library(vegan)
#以兩組環境變數（DOC、AP+AK）為例，執行變差分解，詳情 ?varpart
#引數 chisquare = TRUE，執行 CCA 的變差分解；預設情況下 chisquare = FALSE，即執行 RDA 的變差分解
#注：如果 varpart() 不支援 CCA，請更新 R 版本（如 R3.6 的 vegan）
cca_vp <- varpart(otu, env['DOC'], env[c('AP', 'AK')], chisquare = TRUE)
cca_vp
plot(cca_vp, digits = 2, Xnames = c('DOC', 'AP+AK'), bg = c('blue', 'red'))
對於該示例，DOC（“a+b”部分）解釋了總變差的12.4%，其中單獨（“a”部分）解釋的比例為3.8%；AP+AK單獨（“c”部分）解釋了8.9%，和DOC共同（“b”部分）解釋了8.6%；DOC+AP+AK未能（“d”部分）解釋的變差佔比78.7%。
##置換檢驗確定解釋變差的顯著性
#需結合 CCA 執行式，vegan 中為 cca()
#解釋變差的置換檢驗，以 DOC 所能解釋的全部變差為例；999 次置換，詳情 ?anova.cca
anova.cca(cca(otu~DOC, env), permutations = 999)
#若考慮 DOC 單獨解釋的變差部分，需將其它變數作為協變數；999 次置換
anova.cca(cca(otu~DOC+Condition(AP+AK), env), permutations = 999)
注意事項
最後，還需注意幾點。
變差分解的一個主要目的是計算共同解釋的變差，但在解讀時必須非常謹慎，因為共同部分很難確定與哪組變數有關。
變差分解共同解釋部分（對應於前圖“b”部分）與方差分析中因子的互動作用（interaction）不同。互動作用是雙因素方差分析中，一個因子不同水平與其它因子不同水平之間的協同作用。變差分解是由於解釋變數共線性引起解釋變差的重疊，不是協同作用的結果。換句話說，如果兩個（組）變數相互獨立（不相關，正交），則共同解釋部分“[b]”將為0。
參考資料
DanielBorcard, FranoisGillet, PierreLegendre, et al. 數量生態學:R語言的應用（賴江山譯）. 高等教育出版社, 2014.
Blanchet F G , Legendre P , Borcard D . Forward selection of explanatory variables. Ecology, 2008, 89(9):2623-2632.
Borcard D , Drapeau L P . Partialling out the Spatial Component of Ecological Variation. Ecology, 1992, 73(3):1045-1055.
David Zelený博士：https://www.davidzeleny.net/anadat-r/doku.php/en:varpart
Peres-Neto P R , Legendre P , Dray, Stéphane, et al. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology, 2006, 87(10):2614-2625.
