【R分享實(shí)戰(zhàn)】方差分解分析（Variationpartitionanalysis）及其顯著性檢驗(yàn)

【R分享實(shí)戰(zhàn)】方差分解分析（Variationpartitionanalysis）及其顯著性檢驗(yàn)“

世上所有的成功只是努力的另一個(gè)代名詞?！?/p>

--科白君"R實(shí)戰(zhàn)"專(zhuān)題·第13篇編輯

|

科白維尼

3470字|8分鐘閱讀本期推送內(nèi)容高分期刊中頻頻登場(chǎng)的方差分解分析(VPA)到底是啥？可能很多同學(xué)還沒(méi)接觸過(guò)該分析。目前，群落分析中常見(jiàn)的環(huán)境因子分析包括典范對(duì)應(yīng)分析（canonicalcorrespondenceanalysis）和冗余分析(Redundancyanalysis)，這兩種類(lèi)型分析都是基于降維的思想，將樣本、物種、環(huán)境因子的信息映射到二維平面上，從而判斷三者間的關(guān)系，可用于發(fā)現(xiàn)對(duì)群落結(jié)構(gòu)有影響的環(huán)境變量。而VPA分析可以看做是CCA/RDA分析的一種升級(jí)。用幾組解釋變量（如環(huán)境、氣候、土壤因子等數(shù)據(jù)）來(lái)共同解釋一組響應(yīng)變量的變化（如微生物數(shù)據(jù)），當(dāng)我們需要某個(gè)解釋變量所能夠解釋的方差變異程度信息（即某個(gè)環(huán)境因子對(duì)群落結(jié)構(gòu)變化的貢獻(xiàn)度）時(shí)，就可以采用VPA分析加以補(bǔ)充。01方差分解分析方差分解的數(shù)學(xué)原理，可以參考該鏈接，有比較清晰且易理解的介紹/gustame/constrained-analyses/variation-partitioning。從土壤微生物生態(tài)領(lǐng)域的角度：就是將每個(gè)解釋變量（環(huán)境因子）進(jìn)行獨(dú)立運(yùn)行CCA或RDA，獲得每個(gè)解釋變量（環(huán)境因子）對(duì)于響應(yīng)變量（微生物群落）的方差變異的解釋貢獻(xiàn)度，之后通過(guò)將多組數(shù)據(jù)取交并集的方式獲得每個(gè)解釋變量（環(huán)境因子）的獨(dú)立解釋貢獻(xiàn)度以及環(huán)境因子共同解釋的貢獻(xiàn)度。今天目的是教會(huì)大家利用R語(yǔ)言進(jìn)行VPA分析及其對(duì)應(yīng)解釋度的顯著性檢驗(yàn)。廢話(huà)不多說(shuō)，進(jìn)入代碼講解~02利用R語(yǔ)言進(jìn)行VPA分析這里主要到"vegan"包中的varpart()函數(shù)并使用自帶的數(shù)據(jù)集來(lái)進(jìn)行VPA分析，并使用plot以圖形式輸出結(jié)果。溫馨提示：關(guān)于VPA分析時(shí)所用的環(huán)境因子的分類(lèi)是可以根據(jù)你自己選擇來(lái)劃分的，例如現(xiàn)有氣候因子年均溫、年降雨量；土壤理化因子有pH、SOC、TN、TP、TK、含水率、陽(yáng)離子交換量。響應(yīng)變量為土壤微生物豐度數(shù)據(jù)。因此，我將所有變量劃分為兩個(gè)解釋變量：氣候和土壤理化。加載自帶數(shù)據(jù)集"mite"、"mite.env"、"mite.pcnm"。其中mite為物種豐度數(shù)據(jù)，mite.env和mite.pcnm為環(huán)境因子數(shù)據(jù)。#intall.packages("vegan")#加載vegan包library(vegan)#讀取vegan包自帶數(shù)據(jù)集data(mite)head(mite)data("mite.env")head(mite.env)data("mite.pcnm")head(mite.pcnm)結(jié)果：注意，數(shù)據(jù)格式的模仿：數(shù)據(jù)集的行為樣本編號(hào)，列為物種或者環(huán)境變量。接著，我們利用varpart()函數(shù)進(jìn)行擬合模型。由于mite.pcnm的數(shù)據(jù)集太大了，我們就取了取了前5個(gè)作為第二個(gè)環(huán)境因子。#擬合模型varpart#簡(jiǎn)單了解函數(shù)用法fit<-varpart(mite,mite.env,mite.pcnm[,1:5],transfo="hel")

fit進(jìn)行VPA分析時(shí)，需要記?。旱谝粋€(gè)數(shù)據(jù)物種數(shù)據(jù)，之后兩個(gè)數(shù)據(jù)分別代表不同的環(huán)境變量，transfo表示對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換，hel為hellinger轉(zhuǎn)換，避免“弓形效應(yīng)”。第一個(gè)結(jié)果：為了讓大家更好的理解VPA分析的結(jié)果，這里舉例解讀如圖所示微生物群落的方差總體解釋量為Y，檢測(cè)了兩個(gè)環(huán)境因子數(shù)據(jù)X1、X2。VPA是將兩個(gè)環(huán)境因子的數(shù)據(jù)各自獨(dú)立運(yùn)行RDA分析（RDA則為標(biāo)準(zhǔn)化后的解釋變量對(duì)響應(yīng)變量逐一進(jìn)行多元回歸分析，獲得擬合值、殘差值，最終整合成為擬合值矩陣，該矩陣進(jìn)行PCA、然后排序所得結(jié)果），獲得每個(gè)環(huán)境因子對(duì)于群落整體變差校正解釋量，再運(yùn)行兩者共同存在時(shí)獲得校正后解釋度數(shù)據(jù)R方，具體得到如下結(jié)果：[獨(dú)立運(yùn)行X1的RDA]：A+C部分解釋貢獻(xiàn)度占比[獨(dú)立運(yùn)行X2的RDA]：B+C部分解釋貢獻(xiàn)度占比[共同運(yùn)行X1和X2的RDA]：A+B+C部分解釋貢獻(xiàn)度占比最后可以得到：[A]=[共同運(yùn)行X1、X2的RDA]-[獨(dú)立運(yùn)行X2的RDA][B]=[共同運(yùn)行X1、X2的RDA]-[獨(dú)立運(yùn)行X1的RDA][C]=[共同運(yùn)行X1、X2的RDA]-[A]-[B][未被這兩種環(huán)境因子解釋到的殘差(D)]=Y-[A]-[B]-[C]第二個(gè)結(jié)果：使用plot函數(shù)可視化得到的VPA分析結(jié)果：plot(fit,bg=c("hotpink","skyblue"))#bg表示兩個(gè)變量的背景顏色結(jié)果：通過(guò)維恩圖的展示就可以清楚的看到每種環(huán)境因子對(duì)于總體變異的解釋程度、共有解釋程度以及殘差。需要注意的事，圖中共有的部分產(chǎn)生的原因在于環(huán)境因子數(shù)據(jù)對(duì)微生物解釋存在著共線(xiàn)性而產(chǎn)生，如果環(huán)境因子完全相互獨(dú)立理論上重疊部分=0；此外，如果解釋貢獻(xiàn)度出現(xiàn)負(fù)數(shù)，則說(shuō)明這組環(huán)境因子數(shù)據(jù)對(duì)群落數(shù)據(jù)方差變化解釋程度比使用隨機(jī)變量的解釋程度還要低，一般當(dāng)做貢獻(xiàn)度為0進(jìn)行解釋?zhuān)ㄗh在選擇環(huán)境因子時(shí)減少共線(xiàn)性程度較高的環(huán)境因子以及貢獻(xiàn)度為負(fù)數(shù)的環(huán)境因子數(shù)據(jù)，以保證結(jié)果準(zhǔn)確。03VPA分析后續(xù)的顯著性檢驗(yàn)實(shí)際上RDA和方差分解的結(jié)果是相同的（這里不做展示），可以發(fā)現(xiàn)環(huán)境X1與X2共同解釋了26%的群落結(jié)構(gòu)變異，其中X1解釋了16%的群落組成變異，氣候解釋了5%的群落組成變異。然而對(duì)于二者是否顯著，是可以做進(jìn)一步檢驗(yàn)。特殊情況，關(guān)于冗余分析（RDA）中環(huán)境因子共同解釋部分出現(xiàn)負(fù)值的說(shuō)明。我推薦賴(lài)江山老師對(duì)于該結(jié)果的深刻解讀。/blog-267448-1187530.html。代碼如下：#首先將兩個(gè)環(huán)境因子合并mite.total<-data.frame(mite.env,mite.pcnm[,1:5])#描述partialRDA的公式#X1formula_X1<-formula(mite~SubsDens+WatrCont+Substrate+Shrub+Topo+Condition(V1)+Condition(V2)+Condition(V3)+Condition(V4)+Condition(V5))#X2formula_X2<-formula(mite~Condition(SubsDens)+Condition(WatrCont)+Condition(Substrate)+Condition(Shrub)+Condition(Topo)+V1+V2+V3+V4+V5)#X1和X2formula_X1X2<-formula(mite~SubsDens:WatrCont:Substrate:Shrub:Topo:V1:V2:V3:V4:V5)#利用partialRDA進(jìn)行999次的置換檢驗(yàn)，最后得出模型的顯著性anova(rda(formula_X1,data=mite.total))anova(rda(formula_X2,data=mite.total))anova(rda(formula_X1X2,data=mite.total))其中：Condition是表示控制一部分變量，只探究非condition部分的變量，即控制mite.pcnm則只探究mite.env。反之，探究另一個(gè)環(huán)境因子。特別強(qiáng)調(diào)：最后一種情況是我理解上的共同解釋部分的顯著性檢驗(yàn)，但是對(duì)錯(cuò)與否目前我也沒(méi)找到正確答案。這里需要更多的統(tǒng)計(jì)學(xué)專(zhuān)業(yè)的同學(xué)給出建議來(lái)完善。結(jié)果如下