單因素協(xié)方差分析（ANCOVA）。當(dāng)方差分析中存在協(xié)變量時(shí)北专，即可稱為協(xié)方差分析竖螃。其中單因素協(xié)方差分析是最常見的，在單因素方差分析中引入了協(xié)變量逗余。
何謂協(xié)變量，舉例來說季惩。對(duì)懷孕母鼠喂食某種藥物录粱，并觀察藥物處理組和對(duì)照組相比，新生小鼠體重是否具有區(qū)別画拾。由于各母鼠的懷孕時(shí)間有所差別啥繁，而這個(gè)懷孕時(shí)間可能會(huì)對(duì)小鼠體重產(chǎn)生影響而不可忽略，就可視作協(xié)變量處理青抛。再舉一例旗闽，重測(cè)序分析中統(tǒng)計(jì)基因的突變頻率，由于各基因長(zhǎng)度是不一樣的，即使堿基是隨機(jī)突變的适室，更長(zhǎng)的基因也可能會(huì)出現(xiàn)更多的突變堿基數(shù)量嫡意。此時(shí)基因的長(zhǎng)度就是一個(gè)協(xié)變量，必須考慮在內(nèi)

加載數(shù)據(jù)

#讀入文件
soil <- read.table('soil.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
group  <- read.table('group.txt', sep = '\t', header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)
soil <- merge(soil, group, by = 'sample')
#以 shannon 指數(shù)為例捣辆，同時(shí)將分組列轉(zhuǎn)換為因子變量
shannon <- soil[ ,c('sample', 'treat', 'shannon', 'days')]
shannon$treat <- factor(shannon$treat)

單因素協(xié)方差分析（ANCOVA）
假設(shè)我們使用同一來源的土壤進(jìn)行盆栽試驗(yàn)（土壤類型蔬螟、植物類型一致），分別添加了三種化學(xué)物質(zhì)（a汽畴、b旧巾、c），探究不同類型的化學(xué)物質(zhì)是如何影響植物根際細(xì)菌群落的我們?cè)谥参锘ㄆ谌佑^察忍些，由于各盆栽中植物開花期時(shí)間并非一致鲁猩，所以我們考慮將植物生長(zhǎng)時(shí)間（days）為協(xié)變量.

評(píng)估檢驗(yàn)的假設(shè)條件

ANCOVA要求數(shù)據(jù)服從正態(tài)分布，以及各組方差相等罢坝，同時(shí)還假定回歸斜率相同廓握。

正態(tài)性檢驗(yàn)

#QQ-plot 檢查數(shù)據(jù)是否符合正態(tài)分布
library(car)
qqPlot(lm(shannon~treat, data = shannon), simulate = TRUE, main = 'QQ Plot', labels = FALSE)

由圖可知，我們的數(shù)據(jù)符合正態(tài)分布模型（盡管似乎有個(gè)離群點(diǎn)炸客，這時(shí)候可以考慮刪除這個(gè)樣本后再繼續(xù)疾棵，本示例演示不再將它刪除，直接進(jìn)入下一步分析）痹仙。

方差齊性檢驗(yàn)

#使用 Bartlett 檢驗(yàn)進(jìn)行方差齊性檢驗(yàn)（p 值大于 0.05 說明方差齊整）
bartlett.test(shannon~treat, data = shannon)
Bartlett test of homogeneity of variances

data:  shannon by treat
Bartlett's K-squared = 4.5383, df = 2, p-value = 0.1034

回歸斜率的同質(zhì)性檢驗(yàn)

這里ANCOVA包含植物生長(zhǎng)時(shí)間×化學(xué)物質(zhì)類型的交互項(xiàng)是尔，需對(duì)回歸斜率的同質(zhì)性進(jìn)行檢驗(yàn)，若交互效應(yīng)顯著开仰，則意味著植物生長(zhǎng)時(shí)間和植物根際菌群的Shannon指數(shù)的關(guān)系依賴于所添加的化學(xué)物質(zhì)類型

#檢驗(yàn)回歸斜率的同質(zhì)性拟枚，交互效應(yīng)不顯著則支持斜率相等的假設(shè)（即 p 值大于 0.05 說明回歸斜率相等）
summary(aov(shannon~days*treat, data = shannon))
  Df Sum Sq Mean Sq F value  Pr(>F)   
days         1 0.0843  0.0843   1.174 0.28728   
treat        2 0.9704  0.4852   6.758 0.00378 **
days:treat   2 0.1120  0.0560   0.780 0.46764   
Residuals   30 2.1539  0.0718                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

根據(jù)aov()公式可知，這實(shí)際上是一個(gè)雙因素方差分析（众弓，根據(jù)雙因素方差分析中的交互項(xiàng)結(jié)果來判斷回歸斜率的同質(zhì)性恩溅。
結(jié)果顯示交互作用不顯著，支持了斜率相等的假設(shè)。如果假設(shè)不成立，可以嘗試變換協(xié)變量或因變量水泉，或使用能對(duì)每個(gè)斜率獨(dú)立解釋的模型酥夭，或使用無需回歸斜率相等的非參數(shù)ANCOVA方法（如sm包sm.ancova()）。

單因素協(xié)方差分析（ANCOVA）

單因素ANCOVA

上述檢驗(yàn)均已通過，接下來進(jìn)行單因素ANCOVA。
如果不滿足上述前提假設(shè)，一是可以考慮轉(zhuǎn)化數(shù)據(jù)（當(dāng)然锌订，我們需要確保轉(zhuǎn)換后的數(shù)據(jù)能夠被合理解釋，否則將無意義）画株，二是可以考慮使用非參數(shù)的檢驗(yàn)方法辆飘。上述提及了一個(gè)對(duì)應(yīng)單因素協(xié)方差分析的非參數(shù)替代方法啦辐，sm包sm.ancova()。
對(duì)于帶協(xié)變量的項(xiàng)蜈项，以單因素協(xié)方差為例芹关，aov()函
數(shù)書寫為aov(y~x+A)的樣式，其中x為協(xié)變量战得，A為因子變量充边，注意需要將協(xié)變量寫在因子前面，如上所示常侦，協(xié)變量植物生長(zhǎng)時(shí)間（days）在前浇冰，因子化學(xué)物質(zhì)類型（treat）在后，順序不可顛倒聋亡。

#滿足假設(shè)肘习，單因素協(xié)方差分析，詳情使用?aov查看幫助
fit <- aov(shannon~days+treat, data = shannon)
summary(fit)
  Df Sum Sq Mean Sq F value  Pr(>F)   
days         1 0.0843  0.0843   1.190 0.28346   
treat        2 0.9704  0.4852   6.852 0.00333 **
Residuals   32 2.2659  0.0708                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANCOVA結(jié)果表明：（1）植物生長(zhǎng)時(shí)間相隔幾天的差異并未導(dǎo)致其根際菌群產(chǎn)生較大的變異（p值不顯著）坡倔；（2）控制植物生長(zhǎng)時(shí)間漂佩，不同類型的化學(xué)物質(zhì)處理下的植物根際細(xì)菌群落的Shannon指數(shù)存在顯著不同（p值顯著）。

#查看各組均值及標(biāo)準(zhǔn)差
aggregate(shannon$shannon, by = list(shannon$treat), FUN = mean)
 Group.1       x
1       a 9.30400
2       b 9.44425
3       c 9.03825
aggregate(shannon$shannon, by = list(shannon$treat), FUN = sd)
 Group.1         x
1       a 0.1632561
2       b 0.3136513
3       c 0.2899420
#由于使用了協(xié)變量罪塔，若想獲取去除協(xié)變量效應(yīng)后的組均值（調(diào)整的組均值）
library(effects)
effect('treat', fit)
treat effect
treat
       a        b        c 
9.294176 9.446374 9.045950 

因變量投蝉、協(xié)變量和因子之間的關(guān)系圖

#HH 包 ancova() 可繪制因變量、協(xié)變量和因子之間的關(guān)系圖
#詳情使用 ?ancova 查看幫助
library(HH)
ancova(shannon~days+treat, data = shannon)

如上示例征堪，通過“ancova(shannon~days+treat, data = shannon)”瘩缆，我們又執(zhí)行了一次ANCOVA，結(jié)果屏幕輸出佃蚜，和上文結(jié)果一致庸娱。同時(shí)通過關(guān)系圖可知，3條回歸線相互平行谐算，只是截距項(xiàng)不同熟尉，b組截距項(xiàng)最大，c組截距項(xiàng)最兄拗斤儿；回歸線擬合效果并不理想，也對(duì)應(yīng)了我們先前的結(jié)論恐锦，在“數(shù)天”這么一個(gè)短期時(shí)間內(nèi)雇毫，植物根際菌群的Shannon指數(shù)沒有發(fā)生劇烈的改變。

ancova(shannon~days*treat, data = shannon)

我們還通過“ancova(shannon~days*treat, data = shannon)”踩蔚，執(zhí)行了一次雙因素ANOVA。不過在這里我們并不是想做個(gè)雙因素ANOVA分析枚粘，而是在更改了函數(shù)公式后馅闽，意在可視化允許斜率和截距項(xiàng)依據(jù)組別而發(fā)生改變，從而體現(xiàn)那些違背回歸斜率同質(zhì)性的實(shí)例。（上文已知福也，回歸斜率的同質(zhì)性檢驗(yàn)是通過雙因素ANOVA中的交互作用判斷的局骤，本示例中的回歸斜率的同質(zhì)性檢驗(yàn)已通過，大家可以嘗試自行尋找一例回歸斜率不相等的數(shù)據(jù)暴凑，然后使用ancova()作圖查看其與回歸斜率相等的數(shù)據(jù)的差異）
暫時(shí)不太懂峦甩。。现喳。
來源：小白魚的生統(tǒng)筆記