


然後計算降維後資料集之間的典範相關係數。

即可得到典範相關係數。如果上述降維為多個正交軸(最大min[p,q]),則得到各軸的典範相關係數。
CCorA(canonical correlation analysis)出現時間較早,早期也稱為CCA。但考慮到後來出現了典範對應分析(canonical correspondence analysis,CCA),為了區分二者,儘量避免使用CCA作為CCorA的簡稱。
CCorA屬於對稱分析(symmetrical analysis),兩組資料集相互之間無響應變數和解釋變數之分,二者地位等同。例如基因共表達、物種間互作、群落演替過程中物種–環境長期相互作用等。因此,CCorA是一種用於描述資料的探索性分析方法,而非確立因果關係的模型,不涉及假設檢驗過程。但如果有必要,仍可以透過置換檢驗確立相關性的顯著性。
CCorA還可用於確定一組典範變數,即每個組內變數的正交線性組合,可以最好地解釋組內和組間的變異性。
CCorA要求兩組資料為符合多元正態分佈的定量資料。且資料集中的變數具有不同的量綱時,需要標準化處理。
本篇簡介R包vegan的CCorA方法。
library(vegan)
#檢視文件,最下方的示例
?CCorA
資料集
vegan包的內建資料集“mite”,記錄了70個土壤觀測地點的35種甲蟎的個體數量,行為觀測地點,列為甲蟎物種。
library(vegan)
#甲蟎資料集
#詳情 https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/mite
data(mite)
head(mite[ ,1:10])

接下來使用該甲蟎資料,展示vegan中的CCorA過程。
按照文件中的示例,應該是期望確立兩組亞群甲蟎物種之間是否存在某種關聯(正相關可表示趨於共存,負相關可表示趨於競爭)。且兩組資料均為物種變數,地位等同。
#對土壤地點設定分組(分組參考自 Legendre 2005, Fig. 4 的 PCA)
group.1 <- c(1, 2, 4:8, 10:15, 17, 19:22, 24, 26:30)
group.2 <- c(3, 9, 16, 18, 23, 25, 31:35)
#按分組將物種多度資料拆分為兩組資料集
#物種多度資料分析,推薦執行 Hellinger 轉化
mite.hel.1 <- decostand(mite[,group.1], 'hel')
mite.hel.2 <- decostand(mite[,group.2], 'hel')
rownames(mite.hel.1) = paste('S', 1:nrow(mite), sep = '')
rownames(mite.hel.2) = paste('S', 1:nrow(mite), sep = '')
多元正態性評估
#Q-Q 圖評估多元正態性
#若資料服從多元正態性,則點將落在直線附近
par(mfrow = c(1, 2))
qqplot(qchisq(ppoints(nrow(mite.hel.1)), df = ncol(mite.hel.1)), mahalanobis(mite.hel.1, colMeans(mite.hel.1), cov(mite.hel.1)))
abline(a = 0, b = 1)
qqplot(qchisq(ppoints(nrow(mite.hel.2)), df = ncol(mite.hel.2)), mahalanobis(mite.hel.2, colMeans(mite.hel.2), cov(mite.hel.2)))
abline(a = 0, b = 1)

“mite.hel.2”符合,但“mite.hel.1”偏離程度比較明顯。
理論上該“mite.hel.1”資料集不能繼續用於CCorA,但由於本文僅為示例演示,所以請允許我們忽略多元正態性的問題,繼續演示。
實際的情況中,不妨試下其它的轉化資料方法可不可行,否則不可使用。
CCorA執行及結果概述
vegan中,CCorA()用於執行CCorA。
#CCorA,詳情 ?CCorA
#Y 和 X 地位等同
#這裡兩組均為物種丰度,量綱一致,且 Hellinger 轉化降低了高丰度物種的權重,故無需對兩組資料集標準化
#若量綱不同(比如不同的環境變數),則需要將對應的資料集標準化
#999 次置換確定相關性的顯著性
out <- CCorA(Y = mite.hel.1, X = mite.hel.2, stand.Y = FALSE, stand.X = FALSE, permutations = 999)
out
Rank顯示了Y和X矩陣中的變數個數。
Pillai's trace是典範相關係數的和。
Significance of Pillai's trace為置換檢驗的結果,其中from F-distribution代表了隨機置換過程中所得的期望Pillai's trace,based on permutations代表了p值,這裡可知相關性是顯著的。
Canonical Correlations代表了各個軸的典範相關係數,它們的加和即等於Pillai's trace。透過前幾個有代表性的軸所示的相關性,兩組資料集整體相關程度很高。
RDA R squares,Y|X和X|Y的RDA的雙變數冗餘係數(R2),adj. RDA R squares是校正後的R2。它們嚴格來說不屬於CCorA的計算部分,但可以評估一組變數的總變差能夠由另一組變數在一個或多個正交軸上的解釋程度。因為即使R2很小時,也可能出現較大的典範相關係數,此時可根據R2評估典範相關係數的合理性。

對於主要資訊的提取。
#根據 summary() 的提示,提取主要資訊
summary(out)
#例如
#Pillai's trace
out$Pillai
#各個軸的典範相關係數
out$Eigenvalues
#校正後的 RDA R2
out$RDA.adj.Rsq
#置換檢驗 p 值
out$p.perm

關於理解變數間相關性
上面給出了兩個資料集整體的相關性。
對於各變數間的相關性,例如我們想知道哪些物種和哪些物種之間存在較大的相關,可透過排序圖觀測。
vegan中提供了一些視覺化方案,例如biplot()。biplot()一次會輸出兩個圖,分別代表矩陣Y和X。由於CCorA是對稱分析,所以兩個排序圖互有關聯,典範軸顯示了兩個矩陣的共同變化趨勢。每個排序圖都是兩組變數線性組合的結果,所以單個變數解讀比較困難,需兩個圖結合起來解讀。
#作圖觀測,如 biplot(),詳情 ?CCorA
#“objects”生成兩個物件圖,第一個表示 Y 矩陣中的物件,第二個表示 X 矩陣中的物件;以展示第 1、2 軸為例
biplot(out, plot.type = 'objects', plot.axes = c(1, 2))

#“variables” 生成兩個變數圖,第一個表示 Y 矩陣中的變數,第二個表示 X 矩陣中的變數;以展示第 1、2 軸為例
biplot(out, plot.type = 'variables', plot.axes = c(1, 2), cex = c(0.7, 0.6))

#“ov”生成四個圖,包含兩個物件和兩個變數;以展示第 1、2 軸為例
biplot(out, plot.type = 'ov', plot.axes = c(1, 2), cex = c(0.7, 0.6))

#“biplots”產生兩個雙序圖,第一個用於 Y 矩陣,第二個用於 X 矩陣;以展示第 1、2 軸為例
biplot(out, plot.type = 'biplots', plot.axes = c(1, 2), cex = c(0.7, 0.6))

好了我們結合上述這兩個雙序圖作個解讀,加深對這種變數間相關性的理解。CCorA中,變數間的相關性是透過在兩個資料集中,變數在物件中的相似分佈特徵,作為呈現的。因此,排序圖中的相關性可以這樣解讀。

Hotelling, H. 1936. Relations between two sets of variates. Biometrika 28: 321-377.



