單因素協方差分析(ANCOVA)在R中實作
前文簡介了單因素方差分析(單因素ANOVA)在R語言中的實作方法，本篇繼續簡介單因素協方差分析(ANCOVA)。當方差分析中存在協變量時，即可稱為協方差分析。其中單因素協方差分析是最常見的，在單因素方差分析中引入了協變量。
何謂協變量，舉例來說。對懷孕母鼠喂食某種藥物，并觀察藥物處理組和對照組相比，新生小鼠體重是否具有差別。由于各母鼠的懷孕時間有所差别，而這個懷孕時間可能會對小鼠體重産生影響而不可忽略，就可視作協變量處理。再舉一例，重測序分析中統計基因的突變頻率，由于各基因長度是不一樣的，即使堿基是随機突變的，更長的基因也可能會出現更多的突變堿基數量。此時基因的長度就是一個協變量，必須考慮在内。好了，下面接正文。
本文使用的作圖資料的網盤連結(提取碼z4w4)：
https://pan.baidu.com/s/1J-9GsmoHuQ_CEpxeWyEQsA
資料預處理
示例資料說明
我們首先将示例資料讀到R中，并從中挑選部分資料作為示範。
#讀入檔案
soil group soil
#以 shannon 指數為例，同時将分組列轉換為因子變量
shannon shannon$treat str(shannon)
head(shannon)
假設存在這麼一個研究：
我們使用來源于同一環境中的土壤進行盆栽試驗(土壤類型一緻)，并種植了同種植株(植物類型一緻)。我們将盆栽(1植株/1盆栽)分為了3組，分别在土中添加了三種化學物質(a、b、c)；然後等待植物到達花期後，收集每個植株的根際土，并通過16S測序，獲得了植物根際細菌群落的Alpha多樣性指數，意在探究不同類型的化學物質是如何影響植物根際細菌群落的。
由于各盆栽中植株到達開花期時間并非一緻，無法保證所有植株均在同一天取樣觀察。盡管期間并未相隔很長時間(天)，理論上單因素ANOVA就可以滿足需求，但我們仍然想要将各盆栽中植株到達花期的時間作為協變量去考慮(如果覺得植物生長時間相隔幾天的差異可能會導緻其根際菌群産生較大的變異時，盡管實際上可能有些不妥，但作為示例，請忽略這點)，嘗試使用單因素協方差分析(ANCOVA)。
對應于上述挑選出的測試資料“shannon”：sample，試驗樣本名稱，每一個樣本即對應一個盆栽，各盆栽中土壤類型、植物類型完全相同，而且均是1植株/1盆栽；treat，在土壤中添加的三種化學物質(a、b、c)，這列作為分組列，需要轉換為因子變量類型，各組處理間互相獨立；shannon，Alpha多樣性指數中的Shannon指數，數值變量；days，各盆栽中植株的開花時間(即生長天數)，這裡為數值變量。
這裡，植物根際菌群的Shannon指數為因變量，植物生長天數(days)為協變量，考慮使用ANCOVA，探究不同類型的化學物質處理下的植物根際細菌群落的Shannon指數是否存在顯著不同。
評估檢驗的假設條件
ANCOVA要求資料服從正态分布，以及各組方差相等，同時還假定回歸斜率相同。
正态性檢驗
同單因素ANOVA的方法，可使用Q-Q圖來檢驗正态性假設。
#QQ-plot 檢查資料是否符合正态分布
library(car)
qqPlot(lm(shannon~treat, data = shannon), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
由圖可知，我們的資料符合正态分布模型(盡管似乎有個離群點，這時候可以考慮删除這個樣本後再繼續，本示例示範不再将它删除，直接進入下一步分析)。
方差齊性檢驗 同單因素ANOVA的方法，可使用Bartlett檢驗進行方差齊性檢驗。
#使用 Bartlett 檢驗進行方差齊性檢驗(p 值大于 0.05 說明方差齊整)
bartlett.test(shannon~treat, data = shannon)
結果顯示，我們的資料各組方差相等。
回歸斜率的同質性檢驗
這裡ANCOVA包含植物生長時間×化學物質類型的互動項，需對回歸斜率的同質性進行檢驗，若互動效應顯著，則意味着植物生長時間和植物根際菌群的Shannon指數的關系依賴于所添加的化學物質類型。
#檢驗回歸斜率的同質性，互動效應不顯著則支援斜率相等的假設(即 p 值大于 0.05 說明回歸斜率相等)
summary(aov(shannon~days*treat, data = shannon))
根據aov()公式可知，這實際上是一個雙因素方差分析(關于雙因素ANOVA，下篇再簡介)，根據雙因素方差分析中的互動項結果來判斷回歸斜率的同質性。 結果顯示互動作用不顯著，支援了斜率相等的假設。 如果假設不成立，可以嘗試變換協變量或因變量，或使用能對每個斜率獨立解釋的模型，或使用無需回歸斜率相等的非參數ANCOVA方法(如sm包sm.ancova())。
單因素協方差分析(ANCOVA)
單因素ANCOVA
上述檢驗均已認證，接下來進行單因素ANCOVA。
如果不滿足上述前提假設，一是可以考慮轉化資料(當然，我們需要確定轉換後的資料能夠被合了解釋，否則将無意義)，二是可以考慮使用非參數的檢驗方法。上述提及了一個對應單因素協方差分析的非參數替代方法，sm包sm.ancova()。
對于帶協變量的項，以單因素協方差為例，aov()函數書寫為aov(y~x+A)的樣式，其中x為協變量，A為因子變量，注意需要将協變量寫在因子前面，如上所示，協變量植物生長時間(days)在前，因子化學物質類型(treat)在後，順序不可颠倒。
#滿足假設，單因素協方差分析，詳情使用?aov檢視幫助
fit summary(fit)
ANCOVA結果表明：(1)植物生長時間相隔幾天的差異并未導緻其根際菌群産生較大的變異(p值不顯著)；(2)控制植物生長時間，不同類型的化學物質處理下的植物根際細菌群落的Shannon指數存在顯著不同(p值顯著)。
對于各組均值的獲得方式。
#檢視各組均值及标準差
aggregate(shannon$shannon, by = list(shannon$treat), FUN = mean)
aggregate(shannon$shannon, by = list(shannon$treat), FUN = sd)
#由于使用了協變量，若想擷取去除協變量效應後的組均值(調整的組均值)
library(effects)
effect('treat', fit)
因變量、協變量和因子之間的關系
我們可以使用HH包ancova()函數，繪制因變量、協變量和因子之間的關系圖檢視。 ancova()函數既可根據輸入的公式執行對應的方差分析，并在螢幕輸出方差分析結果； 同時又可生成一張關系圖，由兩部分組成，左側面闆取決于因子的水準groups，右側面闆是所有groups的疊加。
#HH 包 ancova() 可繪制因變量、協變量和因子之間的關系圖
#詳情使用 ?ancova 檢視幫助
library(HH)
ancova(shannon~days+treat, data = shannon)
ancova(shannon~days*treat, data = shannon)
如上示例，通過“ancova(shannon~days+treat, data = shannon)”，我們又執行了一次ANCOVA，結果螢幕輸出，和上文結果一緻。同時通過關系圖可知，3條回歸線互相平行，隻是截距項不同，b組截距項最大，c組截距項最小；回歸線拟合效果并不理想，也對應了我們先前的結論，在“數天”這麼一個短期時間内，植物根際菌群的Shannon指數沒有發生劇烈的改變。
我們還通過“ancova(shannon~days*treat, data = shannon)”，執行了一次雙因素ANOVA。不過在這裡我們并不是想做個雙因素ANOVA分析，而是在更改了函數公式後，意在可視化允許斜率和截距項依據組别而發生改變，進而展現那些違背回歸斜率同質性的執行個體。(上文已知，回歸斜率的同質性檢驗是通過雙因素ANOVA中的互動作用判斷的，本示例中的回歸斜率的同質性檢驗已認證，大家可以嘗試自行尋找一例回歸斜率不相等的資料，然後使用ancova()作圖檢視其與回歸斜率相等的資料的差異)
友情連結
R語言執行單因素方差分析及多重比較
R語言執行兩組間差異分析Wilcox秩和檢驗
R語言執行兩組間差異分析T檢驗
葉綠體基因注釋工具PGA
葉綠體/線粒體線上注釋網站GeSeq
線粒體線上注釋網站MITOS
R語言繪制蝴蝶(柱狀)圖
R語言繪制雙向柱狀圖
R語言繪制分組柱狀圖
R語言繪制堆疊柱狀圖
R語言繪制餅圖(扇形圖)
R語言繪制花瓣圖