某研究對三種環境中的土壤進行了取樣（環境A、B、C，各環境下采集9個土壤樣本），結合二代測序觀測了土壤細菌群落物種組成，並對多種土壤理化指標進行了測量。獲得兩個資料集。
（1）檔案“otu_table.txt”為透過16S測序所得的細菌OTU丰度表格，下文統稱為“物種”，儘管OTU並不完全等於物種。
（2）檔案“env_table.txt”記錄了多種土壤理化指標資訊，即環境資料。
接下來期望透過CCA分析這些資料，用以分析群落物種組成與環境的關係（作為示例，暫且忽略本資料更適用線性模型還是單峰模型）。
vegan包執行CCA
library(vegan)
#詳情 ?cca
#呼叫格式 1
cca(Y, X, W)
#或者格式 2
cca(Y~var1+var2+var3+factorA+var2*var3+Condition(var4))
在第1種呼叫格式中，Y為響應變數矩陣，X為解釋變數矩陣，W為偏CCA的協變數矩陣。當不考慮協變數時，可直接輸入響應變數矩陣Y和解釋變數矩陣X，即“cca(Y, X)”。這種寫法雖然簡單，但要求解釋變數矩陣或協變數矩陣中不能含有因子變數（定性變數）。（當僅存在響應變數矩陣Y時，則執行CA排序）
因此，通常建議使用第2種寫法。其中Y為響應變數矩陣，var1、var2、var3為定量解釋變數，factorA為因子變數，var2*var3為變數2和變數3的互動作用項，var4為協變數。與上述第1種寫法相比，不僅能夠識別因子變數，還能夠對變數間的互動作用加以考慮。
本示例中，環境變數全部為定量解釋變數，且暫不考慮解釋變數間的互動作用以及協變數的情況。因此兩種寫法分別為：
#讀入物種資料，以細菌 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)
#CCA，詳情 ?cca
#呼叫格式 1；這裡無協變數矩陣，所以直接輸入響應變數矩陣和解釋變數矩陣
otu_cca <- cca(otu, env)
#或者格式 2；Y~. 是 Y~var1+var2+…+varn 的簡寫，不涉及互動作用及協變數
otu_cca <- cca(otu~., env)
CCA結果的初步解讀
CCA可獲得的約束軸數為min[p–1, m, n–1]。其中，p為響應變數（物種）數量；m為定量解釋變數數量以及定性解釋變數（因子變數）的因子水平的自由度（即該變數因子水平數減1）；n為排序物件（樣方）數量。
CCA的總變差≠總方差，而是透過一個叫均方列聯絡數（mean squared contingency coefficient）或稱總慣量（total inertia）的指標表徵。CCA中的R2即代表了總慣量（而非總方差）被環境變數所解釋的程度，約束軸承載了被成功解釋的慣量部分。由於CCA與CA共享一套基礎演算法，CCA是在CA的基礎上新增約束過程發展而來，其約束演算法源自RDA中使用的多元迴歸。因此其很多特徵與CA（體現在樣方與物種的關係）或RDA（體現在環境變數的解釋規則）相似，可分別參考前文CA或RDA。
CCA結果總概
透過summary()檢視CCA結果概要。關於CCA中I型和II型標尺的基本特徵，可參考前文
#檢視統計結果資訊，以 I 型標尺為例
otu_cca.scaling1 <- summary(otu_cca, scaling = 1)
otu_cca.scaling1
CCA排序圖
對於本示例的CCA結果，如下為一些簡略的排序圖展示，便於觀測資料。關於更美觀的視覺化方案，可參考本文末。
#作圖檢視排序結果，詳情 ?plot.cca
#三序圖，包含 I 型標尺和 II 型標尺，樣方座標展示為使用物種加權計算的樣方得分
par(mfrow = c(1, 2))
plot(otu_cca, scaling = 1, main = 'I 型標尺', display = c('wa', 'sp', 'cn'))
plot(otu_cca, scaling = 2, main = 'II 型標尺', display = c('wa', 'sp', 'cn'))
在這裡，I型標尺中，樣方得分是物種得分加權平均，故排序圖中物種通常分佈在樣方範圍之外；II型標尺中，物種得分是樣方得分的平均值，並按出現該物種的所有樣方中該物種丰度加權，故排序圖中樣方通常分佈在物種範圍之外。
#比較分別使用物種加權計算的樣方座標以及擬合的樣方座標的差異
#隱藏物種，以 I 型標尺為例展示雙序圖
par(mfrow = c(1, 2))
plot(otu_cca, scaling = 1, main = 'I 型標尺，加權', display = c('wa', 'cn'))
plot(otu_cca, scaling = 1, main = 'I 型標尺，擬合', display = c('lc', 'cn'))
主要資訊提取
若有需要，考慮在結果中將需要的資訊提取出來。
#scores() 提取排序得分（座標），以 I 型標尺為例，前四軸為例
#使用物種加權和計算的樣方得分
otu_cca_site.scaling1 <- scores(otu_cca, choices = 1:4, scaling = 1, display = 'wa')
#物種變數（響應變數）得分
otu_cca_sp.scaling1 <- scores(otu_cca, choices = 1:4, scaling = 1, display = 'sp')
#環境變數（解釋變數）得分
otu_cca_env.scaling1 <- scores(otu_cca, choices = 1:4, scaling = 1, display = 'bp')
#或者在 summary() 後提取，以 I 型標尺為例，前四軸為例
otu_cca.scaling1 <- summary(otu_cca, scaling = 1)
#使用物種加權和計算的樣方得分
otu_cca_site.scaling1 <- otu_cca.scaling1$site[ ,1:4]
#物種
otu_cca_sp.scaling1 <- otu_cca.scaling1$species[ ,1:4]
#環境
otu_cca_env.scaling1 <- otu_cca.scaling1$biplot[ ,1:4]
#若需要輸出在本地
#樣方
write.table(data.frame(otu_cca_site.scaling1), 'otu_cca_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種
write.table(data.frame(otu_cca_sp.scaling1), 'otu_cca_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境
write.table(data.frame(otu_cca_env.scaling1), 'otu_cca_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#不建議直接在原始資料集中提取，因為這樣提取的座標數值未經標尺縮放處理，不利於反映生物學問題
#otu_cca$CCA$u[ ,1:4]
#otu_cca$CCA$v[ ,1:4]
#otu_cca$CCA$biplot[ ,1:4]
#若對典範特徵係數（即每個解釋變數與每個約束軸之間的迴歸係數）感興趣
#coef() 提取 CCA 典範係數
cca_coef <- coef(otu_cca)
和多元迴歸的R2一樣，CCA約束軸承載的總變差（響應變數矩陣的總變差能夠被解釋變數解釋的部分，如果用比例表示，即等同於多元迴歸的R2，也就是約束軸的總解釋量）以及各約束軸解釋量等，有待校正。初始所得的CCA模型需要經過校正和檢驗後才可使用，否則可能會帶來較大的偏差，原因可參考前文。
對於校正前後R2的提取，可使用RsquareAdj()完成。R2經過校正後通常會降低，這是必然的。
#RsquareAdj() 提取 R2，詳情 ?RsquareAdj()
r2 <- RsquareAdj(otu_cca)
otu_cca_noadj <- r2$r.squared #原始 R2
otu_cca_adj <- r2$adj.r.squared #校正後的 R2
對於示例資料的CCA排序結果，透過上文，我們已經得知了響應變數矩陣的總變差為1.3043，未校正前的約束軸解釋的變差總和為0.7402（佔比0.5675，即未校正前的R2）；那麼根據校正後的R2（0.33739），可還原校正R2後的約束軸解釋的變差總和約為“1.3043×0.33739=0.4401”。
對於各約束軸校正後的解釋量。例如，已知CCA1的解釋變差佔約束軸總解釋變差的0.3910，由於校正後的R2為0.33739，那麼可推斷校正R2後的CCA1解釋率為“0.3910×0.33739=0.1319”；同理，校正R2後CCA2解釋率“0.2517×0.33739=0.0849”；其餘約束軸的解釋量以此類推。
#概念如上所述，如下為程式碼實現
#關於約束軸承載的特徵值或解釋率，應當在 R2 校正後重新計算
otu_cca_exp_adj <- otu_cca_adj * otu_cca$CCA$eig/sum(otu_cca$CCA$eig)
otu_cca_eig_adj <- otu_cca_exp_adj * otu_cca$tot.chi
約束軸的置換檢驗及p值校正
我們看到本示例的約束軸解釋率還算可以（以前兩軸為例，解釋量0.1319+0.0849=0.2168，不算很低），那麼，這些約束軸是否有效呢？
對於置換檢驗，首先檢驗全模型，通過後再檢驗各約束軸。
#置換檢驗
#所有約束軸的置換檢驗，即全域性檢驗，基於 999 次置換，詳情 ?anova.cca
otu_cca_test <- anova.cca(otu_cca, permutations = 999)
#各約束軸逐一檢驗，基於 999 次置換
otu_cca_test_axis <- anova.cca(otu_cca, by = 'axis', permutations = 999)
#p 值校正（Bonferroni 為例）
otu_cca_test_axis$`Pr(>F)` <- p.adjust(otu_cca_test_axis$`Pr(>F)`, method = 'bonferroni')
本示例中僅前兩軸顯著，也就是說，第三及之後的約束軸理論上不能再用作分析了。
注：如果檢驗未透過，在考慮更換其它方法之前，不妨先試試剔除一些環境變數，因為也有可能是某些環境變數與群落特徵之間缺乏較好的線性關係所致。
CCA的變數選擇（前向選擇為例）
有時解釋變數太多，需要想辦法減少解釋變數的數量。減少解釋變數的數量可能有兩個並不衝突的原因：第一是尋求簡約的模型，利於我們對模型的解讀；第二是當解釋變數過多時會導致模型混亂，例如有些解釋變數之間可能存在較強的線性相關，即共線性問題，可能會造成迴歸係數不穩定。
方差膨脹因子（Variance Inflation Factor，VIF）顯示環境變數間共線性很嚴重。有幾個變數的VIF值特別大（超過10甚至20），所以有必要剔除一些變數。
#計算方差膨脹因子，詳情 ?vif.cca
vif.cca(otu_cca)
#前向選擇，這裡以 ordiR2step() 的方法為例，基於 999 次置換檢驗，詳情 ?ordiR2step
otu_cca_forward_pr <- ordiR2step(cca(otu~1, env), scope = formula(otu_cca), R2scope = TRUE, direction = 'forward', permutations = 999)
#以 otu_cca 和 otu_cca_forward_pr 為例，簡要繪製雙序圖比較變數選擇前後結果
par(mfrow = c(1, 2))
plot(otu_cca, scaling = 1, main = '原始模型，I 型標尺', display = c('wa', 'cn'))
plot(otu_cca_forward_pr, scaling = 1, main = '前向選擇後，I 型標尺', display = c('wa', 'cn'))
#細節部分檢視
summary(otu_cca_forward_pr, scaling = 1)
#比較選擇前後校正後 R2 的差異，詳情 ?RsquareAdj
#可以看到變數選擇後，儘管去除了很多環境變數，但總 R2 並未損失很多
RsquareAdj(otu_cca)$adj.r.squared
RsquareAdj(otu_cca_forward_pr)$adj.r.squared
#所有約束軸的全域性檢驗，999 次置換，詳情 ?anova.cca
otu_cca_forward_pr_test <- anova.cca(otu_cca_forward_pr, permutations = 999)
#各約束軸逐一檢驗，999 次置換
otu_cca_forward_pr_test_axis <- anova.cca(otu_cca_forward_pr, by = 'axis', permutations = 999)
#p 值校正（Bonferroni 為例）
otu_cca_forward_pr_test_axis$`Pr(>F)` <- p.adjust(otu_cca_forward_pr_test_axis$`Pr(>F)`, method = 'bonferroni')
結果顯示前兩軸顯著，我們考慮將簡約CCA模型的前兩軸座標提取出，以便用於後續分析或視覺化。
##提取或輸出變數選擇後的排序座標
#提取方式和上文一致，這裡透過 scores() 提取，以 I 型標尺為例，前兩軸為例
#使用物種加權和計算的樣方得分
otu_cca_forward_pr_site.scaling1 <- scores(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'wa')
write.table(data.frame(otu_cca_forward_pr_site.scaling1), 'otu_cca_forward_pr_site.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#物種變數（響應變數）得分
otu_cca_forward_pr_sp.scaling1 <- scores(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'sp')
write.table(data.frame(otu_cca_forward_pr_sp.scaling1), 'otu_cca_forward_pr_sp.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
#環境變數（解釋變數）得分
otu_cca_forward_pr_env.scaling1 <- scores(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'bp')
write.table(data.frame(otu_cca_forward_pr_env.scaling1), 'otu_cca_forward_pr_env.scaling1.txt', sep = '\t', col.names = NA, quote = FALSE)
一些視覺化示例
上文已經展示了透過內建的plot()函式，繪製簡單的排序圖用於觀測資料。現在新增些引數在裡面，讓它美觀一些。
#以上述前向選擇後的簡約模型 otu_cca_forward_pr 為例作圖展示前兩軸
#計算校正 R2 後的約束軸解釋率
exp_adj <- RsquareAdj(otu_cca_forward_pr)$adj.r.squared * otu_cca_forward_pr$CCA$eig/sum(otu_cca_forward_pr$CCA$eig)
cca1_exp <- paste('CCA1:', round(exp_adj[1]*100, 2), '%')
cca2_exp <- paste('CCA2:', round(exp_adj[2]*100, 2), '%')
#plot() 作圖，詳情 ?plot.cca
#樣方展示為點，物種暫且展示為“+”，環境變數為向量
par(mfrow = c(1, 2))
plot(otu_cca_forward_pr, type = 'n', display = c('wa', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，雙序圖', xlab = cca1_exp, ylab = cca2_exp)
points(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
plot(otu_cca_forward_pr, type = 'n', display = c('wa', 'sp', 'cn'), choices = 1:2, scaling = 1, main = 'I型標尺，三序圖', xlab = cca1_exp, ylab = cca2_exp)
points(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'sp', pch = 3, col = 'gray', cex = 1)
points(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'wa', pch = 19, col = c(rep('red', 9), rep('orange', 9), rep('green3', 9)), cex = 1)
text(otu_cca_forward_pr, choices = 1:2, scaling = 1, display = 'cn', col = 'blue', cex = 0.8)
或者可在匯出座標後，透過其它作圖包（如ggplot2）繪製。
#提取樣方和環境因子排序座標，前兩軸，I 型標尺
otu_cca_forward_pr.scaling1 <- summary(otu_cca_forward_pr, scaling = 1)
otu_cca_forward_pr.site <- data.frame(otu_cca_forward_pr.scaling1$sites)[1:2]
otu_cca_forward_pr.env <- data.frame(otu_cca_forward_pr.scaling1$biplot)[1:2]
#手動新增分組
otu_cca_forward_pr.env$name <- rownames(otu_cca_forward_pr.env)
otu_cca_forward_pr.site$name <- rownames(otu_cca_forward_pr.site)
otu_cca_forward_pr.site$group <- c(rep('A', 9), rep('B', 9), rep('C', 9))
#ggplot2 作圖
library(ggplot2)
library(ggrepel) #用於 geom_label_repel() 新增標籤
p <- ggplot(otu_cca_forward_pr.site, aes(CCA1, CCA2)) +
geom_point(aes(color = group)) +
stat_ellipse(aes(color = group), level = 0.95, show.legend = FALSE, linetype = 2) +
scale_color_manual(values = c('red', 'orange', 'green3')) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
labs(x = cca1_exp, y = cca2_exp, title = 'I型標尺，雙序圖') +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = otu_cca_forward_pr.env, aes(x = 0, y = 0, xend = CCA1, yend = CCA2), arrow = arrow(length = unit(0.2, 'cm')), size = 0.3, color = 'blue') +
geom_text(data = otu_cca_forward_pr.env, aes(CCA1 * 1.2, CCA2 * 1.2, label = name), color = 'blue', size = 3) +
geom_label_repel(aes(label = name, color = group), size = 3, box.padding = unit(0.5, 'lines'), show.legend = FALSE)
p
