单因素协方差分析（ANCOVA）。当方差分析中存在协变量时，即可称为协方差分析。其中单因素协方差分析是最常见的，在单因素方差分析中引入了协变量。
何谓协变量，举例来说。对怀孕母鼠喂食某种药物，并观察药物处理组和对照组相比，新生小鼠体重是否具有区别。由于各母鼠的怀孕时间有所差别，而这个怀孕时间可能会对小鼠体重产生影响而不可忽略，就可视作协变量处理。再举一例，重测序分析中统计基因的突变频率，由于各基因长度是不一样的，即使碱基是随机突变的，更长的基因也可能会出现更多的突变碱基数量。此时基因的长度就是一个协变量，必须考虑在内
加载数据
#读入文件
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 指数为例，同时将分组列转换为因子变量
shannon <- soil[ ,c('sample', 'treat', 'shannon', 'days')]
shannon$treat <- factor(shannon$treat)
单因素协方差分析（ANCOVA）
假设我们使用同一来源的土壤进行盆栽试验（土壤类型、植物类型一致），分别添加了三种化学物质（a、b、c），探究不同类型的化学物质是如何影响植物根际细菌群落的我们在植物花期取样观察，由于各盆栽中植物开花期时间并非一致，所以我们考虑将植物生长时间（days）为协变量.
评估检验的假设条件
ANCOVA要求数据服从正态分布，以及各组方差相等，同时还假定回归斜率相同。
正态性检验
#QQ-plot 检查数据是否符合正态分布
library(car)
qqPlot(lm(shannon~treat, data = shannon), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
由图可知，我们的数据符合正态分布模型（尽管似乎有个离群点，这时候可以考虑删除这个样本后再继续，本示例演示不再将它删除，直接进入下一步分析）。
方差齐性检验
#使用 Bartlett 检验进行方差齐性检验（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
回归斜率的同质性检验
这里ANCOVA包含植物生长时间×化学物质类型的交互项，需对回归斜率的同质性进行检验，若交互效应显著，则意味着植物生长时间和植物根际菌群的Shannon指数的关系依赖于所添加的化学物质类型
#检验回归斜率的同质性，交互效应不显著则支持斜率相等的假设（即 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
根据aov()公式可知，这实际上是一个双因素方差分析（，根据双因素方差分析中的交互项结果来判断回归斜率的同质性。
结果显示交互作用不显著，支持了斜率相等的假设。如果假设不成立，可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用无需回归斜率相等的非参数ANCOVA方法（如sm包sm.ancova()）。
单因素协方差分析（ANCOVA）
单因素ANCOVA
上述检验均已通过，接下来进行单因素ANCOVA。
如果不满足上述前提假设，一是可以考虑转化数据（当然，我们需要确保转换后的数据能够被合理解释，否则将无意义），二是可以考虑使用非参数的检验方法。上述提及了一个对应单因素协方差分析的非参数替代方法，sm包sm.ancova()。
对于带协变量的项，以单因素协方差为例，aov()函
数书写为aov(y~x+A)的样式，其中x为协变量，A为因子变量，注意需要将协变量写在因子前面，如上所示，协变量植物生长时间（days）在前，因子化学物质类型（treat）在后，顺序不可颠倒。
#满足假设，单因素协方差分析，详情使用?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结果表明：（1）植物生长时间相隔几天的差异并未导致其根际菌群产生较大的变异（p值不显著）；（2）控制植物生长时间，不同类型的化学物质处理下的植物根际细菌群落的Shannon指数存在显著不同（p值显著）。
#查看各组均值及标准差
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
#由于使用了协变量，若想获取去除协变量效应后的组均值（调整的组均值）
library(effects)
effect('treat', fit)
treat effect
treat
a b c
9.294176 9.446374 9.045950
因变量、协变量和因子之间的关系图
#HH 包 ancova() 可绘制因变量、协变量和因子之间的关系图
#详情使用 ?ancova 查看帮助
library(HH)
ancova(shannon~days+treat, data = shannon)
如上示例，通过“ancova(shannon~days+treat, data = shannon)”，我们又执行了一次ANCOVA，结果屏幕输出，和上文结果一致。同时通过关系图可知，3条回归线相互平行，只是截距项不同，b组截距项最大，c组截距项最小；回归线拟合效果并不理想，也对应了我们先前的结论，在“数天”这么一个短期时间内，植物根际菌群的Shannon指数没有发生剧烈的改变。
ancova(shannon~days*treat, data = shannon)
我们还通过“ancova(shannon~days*treat, data = shannon)”，执行了一次双因素ANOVA。不过在这里我们并不是想做个双因素ANOVA分析，而是在更改了函数公式后，意在可视化允许斜率和截距项依据组别而发生改变，从而体现那些违背回归斜率同质性的实例。（上文已知，回归斜率的同质性检验是通过双因素ANOVA中的交互作用判断的，本示例中的回归斜率的同质性检验已通过，大家可以尝试自行寻找一例回归斜率不相等的数据，然后使用ancova()作图查看其与回归斜率相等的数据的差异）
暂时不太懂。。。
来源：小白鱼的生统笔记
网友评论