RDA的基本方法描述
冗餘分析（RDA）和基於轉化的冗餘分析（tb-RDA）
Rao（1964）首次提出冗餘分析（Redundancy analysis，RDA），從概念上講，RDA是響應變數矩陣與解釋變數矩陣之間多元多重線性迴歸的擬合值矩陣的PCA分析，也是多響應變數（multi-response）迴歸分析的拓展。在群落分析中常使用RDA，將物種多度的變化分解為與環境變數相關的變差（variation；或稱方差，variance，因為RDA中變差=方差；由約束/典範軸承載），用以探索群落物種組成受環境變數約束的關係。
包含很多零值的物種多度資料在執行多元迴歸或其它基於歐式距離的分析方法之前必須被轉化，Legendre和Gallagher（2001）提出的基於轉化的RDA（Transformation-based redundancy analysis，tb-RDA）用於解決這個問題。tb-RDA在分析前首先對原始資料做一定的轉化（例如Hellinger預轉化包含很多零值的群落物種資料），並使用轉化後的資料執行RDA。即除了第一步增添了資料轉化外，其餘過程均和常規的RDA相同，只是在原始資料本身做了改動，RDA演算法本質未變。
RDA演算法可以簡要總結如下（詳細過程可參閱Legendre和Legendre（1998）“Numerical Ecology”，579-584頁的內容）。其中矩陣Y是中心化的響應變數矩陣，X矩陣是中心化（或標準化）的解釋變數矩陣。RDA中通常使用標準化後的解釋變數，因為在很多情況下解釋變數具有不同的量綱，解釋變數標準化的意義在於使典範係數的絕對值（即模型的迴歸係數）能夠度量解釋變數對約束軸的貢獻，解釋變數的標準化不會改變回歸的擬合值和約束排序的結果。在群落分析中，響應變數矩陣一般即為物種多度資料，解釋變數矩陣即為環境變數資料。
（1）先將矩陣Y中的每個響應變數分別與矩陣X中的所有解釋變數進行多元迴歸，透過迴歸模型獲得每個響應變數的擬合值（fitted values，即在迴歸線上對應的值）以及殘差（residuals，響應變數的觀測值和擬合值之間的差值），最終得到包含所有響應變數擬合值及殘差的擬合值矩陣Ŷ以及殘差矩陣Yres）。
（2）對擬合值矩陣Ŷ執行PCA，得到典範特徵向量（eigenvectors）矩陣U。使用矩陣U計算兩套樣方排序得分（座標）：一套使用中心化的原始資料矩陣Y獲得在原始變數Y空間內的樣方排序座標（即計算YU，所獲得的座標稱為“樣方得分”，即物種得分的加權和）；另一套使用擬合值矩陣Ŷ獲得在解釋變數X空間內的樣方排序座標（即計算ŶU，所獲得的座標稱為“樣方約束”，即約束變數的線性組合）。
（3）一般來講，RDA過程執行到上步就算完成了。但一般情況下我們會同時對殘差矩陣Yres執行PCA，獲得殘差非約束排序。非約束軸即代表瞭解釋變數未能對響應變數作出解釋的部分，嚴格地來說不屬於RDA的範疇，但能夠幫助我們獲取更多資訊。
Zelený博士使用僅包含一個解釋變數（環境變數）的資料形象化地展示了RDA過程（原文：https://www.davidzeleny.net/anadat-r/doku.php/en:rda_cca）。
（1）執行物種spe1與環境變數env1的線性迴歸（由於此處示例中僅存在一個環境解釋變數，故此迴歸為一元線性迴歸；當存在多解釋變數時，即為多元線性迴歸），將回歸模型擬合的物種丰度值儲存在擬合值矩陣，物種丰度的殘差儲存在殘差矩陣。見下圖1中所示的過程。
（2）如此對物種組成矩陣中的所有物種重複相同的操作，最終獲得包含所有物種丰度擬合值及殘差的兩個矩陣。見下圖2中所示的兩個矩陣。（1）（2）過程即形象化地展示了RDA中的迴歸細節部分。
（3）迴歸過程執行完畢後，使用PCA，在擬合值矩陣中提取約束的排序軸，並在殘差值矩陣中提取非約束的軸。見下圖2中所示的過程，在該示例中，由於僅有一個解釋變數（環境變數env1），因此僅得到一個約束的排序軸（排序圖中的垂直軸是第一個非約束軸）。
RDA排序結果產生的約束軸的數量為min[p, m, n – 1]；如果同時獲得非約束排序結果（即PCA），則非約束軸數量為min[p, n – 1]。其中，p為響應變數數量；m為定量解釋變數數量以及定性解釋變數（因子變數）的因子水平的自由度（即該變數因子水平數減1）；n為排序物件數量。
偏冗餘分析（偏RDA）
偏冗餘分析（Partial canonical ordination，偏RDA）相當於多元偏線性迴歸分析（Davies和Tso，1982），在實際應用中同樣廣泛。
與偏線性迴歸相似，解釋變數被分為兩組：X和W，其中X表示模型中將被考慮的解釋變數，W表示對響應變數Y的影響被控制的協變數。若變數矩陣W對響應變數Y的影響已知，在這種情況下，我們常期望在控制W對Y的影響的前提下，將關注點集中在另一組變數矩陣X對響應變數Y的影響上。例如，可以以氣候變數X作為解釋變數，土壤因子變數W作為協變數，對群落物種資料進行RDA分析。這樣分析的目的是在控制土壤因子影響後，展示單獨能夠被氣候變數線性模型解釋的物種格局。
基於距離的冗餘分析（db-RDA）
儘管tb-RDA的應用拓展了RDA的適用範圍，但無論常規的RDA或tb-RDA，樣方或物種的降維過程實質上均以歐氏距離為基礎。有時我們可能期望關注非歐式距離樣方或物種關係的RDA。
Legendre和Anderson（1999）提出的基於距離的冗餘分析（Distance-based redundancy analysis，db-RDA）用於解決這個問題，並且證明RDA能夠以方差分析方式分析由用於選擇的任何距離矩陣。db-RDA將主座標分析（PCoA）計算的樣方得分矩陣應用在RDA中，其好處是可以基於任意一種距離測度（例如Bray-curtis距離等，而不再僅僅侷限於歐氏距離）進行RDA排序，因此db-RDA在生態學統計分析中被廣泛使用。當然，若是我們直接計算樣方群落間的歐式距離矩陣並將其輸入至db-RDA中計算時，也將會得到和常規的RDA一致的結果。
db-RDA首先基於物種多度資料（根據實際需求，選擇使用原始物種資料或經過某種轉化後的物種資料）計算相異矩陣（如Bray-curtis距離矩陣），作為PCoA的輸入，之後將所有PCoA排序軸上的樣方得分矩陣用於執行RDA，而不再使用原始的物種資料以及解釋變數（環境因子資料）直接作為RDA的輸入。由於在PCoA中可能會產生負特徵值，必要時需要引入一些有效的校正方法。
儘管物種資訊在相異矩陣的計算過程中丟失，但主座標矩陣依然可以視為表徵資料總變差的距離矩陣，因此db-RDA結果反映瞭解釋變數對從整個響應資料中得出的樣方相似性（相異矩陣）的影響。若期望補充物種得分，方法與PCoA中過程一致，物種得分可以透過與它們所在樣方得分的多度加權平均與PCoA軸建立關聯而投影到最終的排序圖中，用以表明響應變數（即物種資訊）對PCoA排序的貢獻程度。
非線性關係的冗餘分析
RDA是透過多元線性迴歸分析後獲得的擬合值矩陣的PCA分析，因此，很多用於多元迴歸的技術都可以在RDA中使用。上面所有的RDA模型都僅使用一階的解釋變數。然而，物種分佈對環境梯度通常是單峰響應（如下圖所示），即物種通常有最適合的生態區域。因此，對於這種情況一階線性模型可能不太適用。
圖注：單一物種沿環境梯度的單峰響應曲線。D為該物種的最適環境梯度，E為較適環境，A、B、C則為位於該物種生態位外（對於該物種來講此時環境很極端，難以生存或增長）。
繪製所有響應變數與單個解釋變數之間的散點圖檢驗非線性關係非常繁瑣，識別和擬合單峰響應最簡單模式是在一階函式基礎上加入二階解釋變數（即二次項）並執行變數的前向選擇，以保留某些一階或二階變數。然而，含有二階解釋變數的RDA結果很難解讀，所以只有充分理由認為是非線性關係時才能使用非線性的RDA。三階的解釋變數也可以用於擬合單峰響應關係，但需要高偏態分佈的響應變數。原始的多項式通常高度相關，因此建議使用正交的多項式。
需注意，所有的非線性RDA僅使用非轉化的響應變數。
RDA結果解讀
從上述解釋計算RDA的步驟中即可看出，RDA的排序軸實際上是解釋變數的線性組合（即線性模型擬合值的排序）。換句話說，RDA的目的是尋找能最大程度解釋響應變數矩陣變差的一些列的解釋變數的線性組合，因此RDA是被解釋變數約束的排序。約束排序與非約束排序的區別很明顯：約束排序過程中解釋變數矩陣控制排序軸的權重（特徵根）、正交性和方向。在RDA中，排序軸解釋或模擬（從統計意義上講）依賴矩陣（響應變數）的變差，並可以檢驗響應變數矩陣Y與解釋變數矩陣X的線性相關顯著性；非約束排序PCA分析則不存在這種情況。儘管在非約束模型中，可以透過在排序後被動地加入解釋變數以達到解釋排序軸的目的（詳見前文），但此舉與約束排序相比具有本質區別。
在群落分析中，對於非約束排序模型（如PCA），我們感興趣的資訊主要是排序圖中樣方和物種變數得分的相對位置、部分排序軸的相對重要性（根據特徵值判斷）以及排序軸的生態解釋等；而對於約束排序模型（如RDA），我們通常更關注環境變數對物種組成的影響（即環境變數所能解釋的變差，以及解釋程度的顯著性）、哪些環境變數對於群落結構的解釋更為重要（變數選擇）以及獲知各變數或變數集解釋的變差（變差分解）等。這些相關的延伸（也很重要）內容不在本文中介紹，若有需要可點選對應連結閱讀。
變差（方差）解釋程度
通常情況下我們在執行RDA時（如使用R語言vegan包的rda()函式執行RDA），能夠同時獲得約束軸（即解釋變數能夠解釋的部分，以約束軸呈現）和非約束軸（即解釋變數未能解釋的部分，多元迴歸的殘差部分，該部分以非約束軸呈現）兩部分資訊，原始響應變數矩陣的總變差為約束軸解釋變差和非約束軸解釋變差的加和。
RDA的約束軸
約束模型解釋變差反映了響應變數變化量的多少與解釋變數有關，如果用比例表示，其值相當於多元迴歸的R2，這個解釋比例值也稱作雙多元冗餘統計（Bimultivariate redundancy statistic）。然而，類似於多元迴歸的未校正R2，RDA的R2也是有偏差的，需要進行校正。
同時，並非每一個約束軸都是合理有效的，還需依據置換的原理檢驗各約束軸的顯著性，對約束軸進行取捨（詳見前文）。因此，與非約束模型PCA等的非約束軸等不同，RDA約束軸的評判方法比較嚴格，若約束軸未透過檢驗，則不應被選擇。（PCA只是探索性分析方法，非約束軸的選擇並無嚴格的標準；RDA已經涉及了統計檢驗的過程，顯著性透過p值衡量）
RDA的非約束軸
如上所述，也就是約束軸未能解釋的，多元迴歸的殘差部分，額外以非約束PCA軸作為呈現。對RDA進行解讀時，最好同時結合約束軸和非約束軸中的資訊，儘管非約束部分嚴格來講不屬於RDA範疇，但很多情形中仍具參考價值。
群落分析中，常透過RDA描述環境變數（解釋變數）解釋樣方物種組成（響應變數）的差異。如果約束軸解釋的變差大於非約束軸解釋的變差，表明響應資料的大部分變化量均可透過解釋變數作出解釋，群落物種組成分佈真實地由給定環境因子所影響（對於RDA結果，即二者呈現出較好的線性梯度）；如果約束軸解釋變差低於非約束軸解釋變差，或者約束軸解釋變差僅佔總變差的較小比例，此時應謹慎對待，因為模型並未顯示出給定環境因子能夠對群落物種的組成作出有效的解釋。
約束模型解釋量偏低的原因可能是還有重要的解釋變數尚未考慮，或是解釋變數之間存在互動作用，或者歸因於實際群落中物種和環境的複雜關係通常很難僅透過簡單的模型有效描述出等（例如常規的RDA基於一階線性模型，但物種和環境的關係多數情況下並非一階線性關係，這種情況下，物種分佈可能並非不受這些環境因子的約束，僅僅歸因於簡單的一階線性模型無法有效描述其關係）。
排序軸特徵值和解釋率
RDA中每一個約束軸的特徵值（eigenvalue）與特徵值總和（約束軸和非約束軸特徵值總和）的比例即為該軸的解釋率。所有約束軸解釋率總和即R2。因此，對於合理的RDA模型來講，選定軸（通常選取特徵值最高的前2-3軸用來觀測）的解釋量不能太低。
少數情況下，殘差之間的排序或相關性（非約束軸）可能比具有良好特徵的約束軸更具生物學意義。如上所述，對於RDA的殘差，即額外以PCA軸的形式呈現。如果有必要，透過觀測非約束空間中的樣方和物種的相對位置可以幫助解讀這些殘差的特徵。
RDA排序圖
維度選擇
對於排序物件、解釋變數以及響應變數的相互關係，最終透過排序圖直觀呈現。一般而言，我們僅選擇前2-3個特徵值較高（且顯著）的約束軸用於觀測（並嘗試對其做出解釋），並表示為二維/三維散點圖的樣式，少數情況下也會根據實際情況選擇特定的排序軸（例如第二軸的趨勢不明顯，第三軸反而明顯，因此跳過了第二軸，使用二維點圖對第一、三軸視覺化並做出解釋；有時也會選擇使用兩個二維點圖，分別展示並解釋第一、二和三、四軸等）。有一點需要切記，就是不要試圖解釋太多的軸，太多的生態維度反而意義不大，正如McCune和Grace（2002）所說：“Very few ecologists have dared to venture into the uncertain waters of four or more dimensions”。
排序圖表示
以R語言vegan包分析群落資料為例，排序物件（樣方）、響應變數（物種）以及解釋變數（環境變數）在各個約束軸中的排序結果分別報告為樣方得分（site scores）、物種得分（species scores）以及解釋變數得分（explanatory variable scores），投影到排序圖中即表示為座標軸上對應位置處的座標。根據是否展示物種向量，排序圖可分為雙序圖（僅展示樣方和環境變數二者關係）和三序圖（展示樣方、物種及環境變數三者關係）。
RDA排序圖中，樣方直接在對應座標處繪製為點。物種變數則呈現為向量，由原點(0,0)起始，指向物種得分的對應座標處，向量的方向表示了該物種丰度增加的方向。解釋變數得分（explanatory variable scores）同樣以向量的形式表示在RDA排序圖中，環境向量的長度表示樣方物種的分佈與該環境因子相關性的大小；向量與約束軸夾角的大小表示環境因子與約束軸相關性的大小，夾角小說明關係密切，若正交則不相關。
圖注：二維散點圖展示某RDA結果的前兩個約束軸。其中A圖為雙序圖，展示樣方和環境變數；B圖為三序圖，展示樣方、物種和環境變數。
PCA中有I、II型標尺的選擇，RDA中也是如此。
無論RDA雙序圖或三序圖，均需要同時展示物件和變數，但沒有同時視覺化物件和變數的最最佳化方法。通常根據關注側重點的不同，考慮使用不同的標度（或稱尺度、標尺）展示排序圖，以重點反映關鍵的資訊。得分集作為RDA輸出的典型特徵，將根據使用的標度而變化。一般有兩種標尺模式，不同模式的排序圖有不同的解讀方式。如果對排序樣方之間的距離更感興趣，或者大多數解釋變數為因子型別，則考慮I型標尺；如果對變數之間的相關關係更感興趣，則考慮II型標尺。
在兩種標尺的RDA排序圖內，對於雙序圖，由於未展示物種資訊故無需考慮樣方和物種關係的解讀方式；對於三序圖，樣方和物種箭頭的解讀同PCA，將樣方物件點垂直投影到物種變數箭頭上位置表示該物種在該樣方內數值在所有樣方內的排序位置。但是，解釋變數的箭頭或因子變數的形心的解讀另有規則。
以下是對RDA I型標尺和II型標尺的簡要描述（詳細過程可參閱Legendre和Legendre（1998）“Numerical Ecology”，585-587頁的內容）。
RDA的I型標尺
I型標尺關注的是物件之間的關係。矩陣U的特徵向量（代表了約束軸中的響應變數得分，參見開頭“RDA演算法簡述”，下述X、Y、Ŷ等同述），被標準化為單位長度。X空間內的樣方得分由公式Z=ŶU獲得，這些向量具有相同的方差λk。Y空間內的樣方得分由公式F=YU獲得，這些向量的方差通常顯著高於λk，因為Y同時涵括了擬合值和殘差值因此其總方差高於Ŷ（僅為擬合值）的方差。透過雙序圖表示矩陣Z和U，或F和U時，特徵向量與樣方得分矩陣的結果較好地還原了原始矩陣：ZU’=Ŷ或FU’=Y。對於定量解釋變數x與樣方得分擬合值的相關性，乘以(λk/Y的總方差)1/2後在各排序軸中呈現，其中λk為約束軸k的特徵值。在這種標尺中，各約束軸中的樣方得分具有不同的方差。
I型標尺的RDA排序圖主要展示如下資訊：
（1）樣方點垂直投影到響應變數或定量解釋變數的箭頭或延長線上，投影點近似於該樣方內該響應變數或解釋變數的數值沿著變數的位置；
（2）響應變數與解釋變數箭頭之間的夾角反映了它們之間的相關性（但響應變數之間的夾角無此含義）；
（3）定性解釋變數的形心與響應變數（物種）箭頭之間的解讀如同樣方點與響應變數之間的解讀（因為定性解釋變數的形心也是一組樣方的形心）；
（4）定性解釋變數的形心之間或形心與樣方點之間的距離近似他們之間的歐式距離。
RDA的II型標尺
II型標尺中，每個特徵向量被標準化為特徵根的平方根，關注的是變數之間的關係。
透過UɅ1/2轉化，將矩陣U的特徵向量（代表了約束軸中的響應變數得分，參見開頭“RDA演算法簡述”，下述X、Y、Ŷ等同述）標準化為(λk)1/2。X空間內的樣方得分透過Z=ŶU獲得後，再經ZɅ–1/2轉化縮放為單位方差。同樣地，Y空間內的樣方得分透過F=YU獲得後，再經FɅ–1/2轉化縮放，並且所得向量的方差通常顯著高於由X所得向量的方差，原因同上述“I型標尺”中的描述。透過雙序圖表示矩陣Z和U，或F和U時，特徵向量與樣方得分矩陣的結果較好地還原了原始矩陣：ZU’=Ŷ或FU’=Y。定量解釋變數x與樣方得分擬合值的相關性直接呈現在各排序軸中。
II型標尺的RDA排序圖主要展示如下資訊：
（1）將樣方點垂直投影到響應變數或定量解釋變數的箭頭或延長線上，投影點位置近似於該響應變數或解釋變數在該樣方內的數值；
（2）響應變數與解釋變數箭頭之間的夾角反映它們之間的相關性，響應變數之間和解釋變數之間也同樣解讀；
（3）定性解釋變數的形心與響應變數箭頭之間的解讀如同樣方點與響應變數之間的解讀，可以透過投影點判斷；
（4）定性解釋變數的形心之間或形心與樣方之間的距離不再近似歐式距離。
看完描述是不是有點懵？接下來透過兩個圖來簡要說明下對RDA排序圖的解讀吧。
（1）a圖，存在排序樣方（樣本）I和II，解釋變數（環境變數）1，探究I或II與1的關係時，將I或II垂直投影在1的向量（箭頭）上，根據交叉點的位置判斷變數1在I或II中的值。交叉點越靠近該變數向量的正方向，則表明所對應的樣方中，該變數的數值越大。例如，假設變數1為土壤碳含量，樣方I投影在1的正方向，樣方II投影在1的負方向上（圖中紅色虛線反向延長線部分），兩個交叉點相比較，I與1的交叉點更位於1延伸方向，因此可知I中的土壤碳含量要比II中的土壤碳含量要高。
若1為響應變數（物種變數），觀察方法同樣適用。例如變數1為物種species1，同樣據此可判斷物種species1在I中的丰度高於在II中的丰度。
注：無論I型標尺或II型標尺，均可據此判斷變數在樣方中的相對數值大小。
（2）b圖，根據向量（箭頭）夾角判斷變數間的相關性。∠a接近90°，即接近正交，表明變數1和2之間的相關性很小，二者相互之間幾乎不存在影響。∠b小於90°，夾角為銳角，表明變數2和3之間存在正相關；銳角角度越小，則正相關性越大。∠c大於90°，夾角為鈍角，表明變數3和4之間存在負相關；鈍角角度越大，則負相關性越大。
注：對於I型標尺，僅可據此觀測解釋變數與響應變數間的相關性；對於II型標尺，既可以據此觀測解釋變數與響應變數間的相關性，也可以觀測解釋變數之間、或響應變數之間的相關性。
（3）對於因子型別的解釋變數5（定性變數，非數值型變數），在圖中以點表示而非以向量表示，探究變數5與其它變數間的相關性時需要根據投影判斷。例如，變數5垂直投影在變數4的正方向，表明與變數4存在正相關；投影在變數2的負方向，表明與變數2存在負相關；相關性的大小，可以透過垂線交叉點與原點(0,0)的距離來表示。
注：對於I型標尺，僅能據此觀測定性解釋變數與響應變數間的相關性；對於II型標尺，既可以據此觀測定性解釋變數與響應變數間的相關性，也可以觀測其與定量解釋變數之間的相關性。
（4）若為I型標尺，還可根據圖中樣方點之間的距離判斷樣方群落之間的相似性。兩個樣方距離越近，則群落相似性越大；反之越低。
（5）此外，還可透過比較解釋變數（環境變數）向量在約束軸上投影的相對長度，判斷環境變數對群落特徵的貢獻度。例如在圖b中，將變數2和變數1均投影至RDA2軸，此時變數2的投影長度相對更長，表明變數2比變數1對RDA2軸形成的貢獻更大。
解釋變數向量與約束軸夾角的大小同樣具有意義，表示解釋變數與約束軸相關性的大小，夾角小說明關係密切，若正交則不相關。例如在圖b中，變數2的向量與RDA2軸的夾角比與RDA1軸的夾角更小，表明變數2與RDA2的關聯程度比與RDA1的關聯程度要高，即相較之下變數2更貢獻於RDA2軸。
注：無論I型標尺或II型標尺，均可據此判斷。
座標位置的特殊調整
此外，在繪製排序圖時，通常會考慮將排序物件或變數的得分（即座標）乘以一個常量，即對排序圖中的點或向量座標同比例放大或縮小一定的數值後展示，以產生易於解讀的生態排序結果（可參見Legendre和Legendre，1998，404頁）。這種做法可能與嚴格的數學定義相違背，但是對於解答生物學問題通常是沒有影響的。
也可參考前述博文的末尾所述。
其它內容
本篇關於RDA的基礎部分，大致就先介紹這些吧，對於初步理解絕對夠用了。
RDA作為多元統計分析中的常見方法，其使用範圍廣泛，因此目前可用於執行RDA的軟體也有很多可供選擇。例如，生態學統計分析的常用R包vegan中，提供了關於RDA模型計算、I型或II型標尺選擇、R2校正、置換檢驗、變數選擇、變差分解等一系列的方法，因此下一篇將簡介R包vegan的RDA分析。
參考資料
DanielBorcard, FranoisGillet, PierreLegendre, et al. 數量生態學:R語言的應用（賴江山譯）. 高等教育出版社, 2014.
David Zelený博士：https://www.davidzeleny.net/anadat-r/doku.php/en:rda_cca_examples
Davies P T, Tso K S. Procedures for Reduced-Rank Regression. Journal of the Royal Statistical Society. Series C (Applied Statistics), 1982, 31(3):244-255.
GUSTA ME Blog：https://mb3is.megx.net/gustame/constrained-analyses/rda
Legendre P, Anderson M J . Distance-Based Redundancy Analysis: Testing Multispecies Responses in Multifactorial Ecological Experiments. Ecological Monographs, 1999, 69(1):1-24.
Legendre P, Gallagher E D . Ecologically Meaningful Transformations for Ordination of Species Data. Oecologia, 2001, 129(2):271-280.
Legendre P, Legendre L. Numerical Ecology. Second English edition. Developments in Environmental Modelling, 1998, 20, Elsevier
Mccune B P, Grace J B. Analysis of Ecological Communities. Journal of Experimental Marine Biology and Ecology, 2002, 289(2).
Rao C R. The use and interpretation of principal component analysis in applied research. Sankhya, A, 1964, 26(4):329-358.
