某研究對三種環境中的土壤進行了採樣（環境A、B、C，各環境下採集9個土壤樣本），結合二代測序觀測了土壤細菌群落物種組成，並對多種土壤理化指標進行了測量。獲得兩個數據集。
（1）文件「otu_table.txt」為通過16S測序所得的細菌OTU豐度表格，下文統稱為「物種」，儘管OTU並不完全等於物種。
（2）文件「env_table.txt」記錄了多種土壤理化指標信息，即環境數據。
接下來期望通過CCA分析這些數據，用以分析群落物種組成與環境的關係（作為示例，暫且忽略本數據更適用線性模型還是單峰模型）。
vegan包中，通過cca()執行CCA，該函數的輸入格式大致如下（與RDA的rda()函數輸入格式相似）。
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可獲得的約束軸數為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。
通過summary()查看CCA結果概要。關於CCA中I型和II型標尺的基本特徵，可參考前文
#查看統計結果信息，以 I 型標尺為例
otu_cca.scaling1 <- summary(otu_cca, scaling = 1)
otu_cca.scaling1
有關CCA排序圖的表現方式、I型和II型標尺等的基本描述，可參考前文。
對於本示例的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
我們看到本示例的約束軸解釋率還算可以（以前兩軸為例，解釋量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)
vegan中前向選擇過程可通過ordistep()、ordiR2step()等函數完成，它們所依據的原理是不同的。詳情可見前文RDA。
#前向選擇，這裡以 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