R语言统计分析——单因素协方差分析

参考资料：R语言实战【第2版】

        单因素协方差分析（ANCONA）扩展了单因素方差分析（ANOVA），包含一个或多个定量的协变量。下面使用multcomp包中的litter数据集进行操作：

# 加载数据集
data(litter,package="multcomp")
# 查看数据集
head(litter)

attach(litter)
# 查看剂量处理数据
table(dose)
# 查看不同剂量处理下的均值
aggregate(weight~dose,FUN=mean)
# 单因素协方差分析
fit<-aov(weight~gesttime+dose)
# 查看分析结果
summary(fit)

        ANCOVA的F检验表明：①小鼠怀孕时间与幼崽出生体重相关；②在控制怀孕时间的情况下，药物剂量与出生体重相关：即在控制怀孕时间下，药物的不同剂量下幼崽出生体重均值显著不同。

        由于使用了协变量，我们需要获取调整后的组均值，即去除协变量效应后的组均值。可以使用effects包中的effect()函数来计算调整的均值：

library(effects)
# 计算调整后均值（不同剂量下小鼠幼崽出生体重调整后均值）
effect("dose",fit)

        本例中，调整的均值与aggregate()函数得出的未调整的均值类似，但并非所有的情况都是如此。总之，effects包为复杂的研究设计提供了强大的计算调整均值的方法。

        下面我们还需要了解哪些剂量处理间差异显著，即多重比较。我们使用multcomp包来对所有均值进行成对比较。另外，multcomp包还可以用来检验用户自定义的均值假设。假设我们仅对未用药条件与其他三种用药剂量影响结果是否不同感兴趣，我们可以如下操作：

#加载multcomp包
library(multcomp)
contrast<-rbind("no drug vs. drug"=c(3,-1,-1,-1))
summary(glht(fit,linfct = mcp(dose=contrast)))

         对照c(9,-1,-1,-1)设定第一组与其他三组的均值进行比较。上面结果显示：p值小于0.05，即未用药组与其他用药条件相比，小鼠出生体重显著的高。

1、评估检验的假设条件

        ANCOVA与ANOVA相同，都需要正态性和同方差性假设，可利用car包中的qqPlot()函数和bartlett.test()函数进行操作，具体参照：R语言统计分析——单因素方差分析。另外，ANCOVA还假定了回归斜率相同，本例中就是假设四个处理组通过怀孕时间来预测出生体重的回归斜率都是相同的。当ANCOVA模型包含怀孕时间×剂量的交互项时，可对回归斜率的同质性进行检验：交互效应若显著，则意味着时间和幼崽出生体重间的关系依赖于药物剂量的水平。

## 检验回归斜率的同质性
# 加载multcomp包
library(multcomp)
# 单因素协方差分析（含交互作用项）
fit2<-aov(weight~gesttime*dose,data=litter)
# 查看结果
summary(fit2)

        含交互项的协方差分析结果显示：交互效应不显著，支持斜率相等的假设。倘若假设不成立，我们可以尝试变换协变量或因变量，或使用能对每个斜率独立解释的模型，或使用不需要假设回归斜率同质性的非参数ANCOVA方法，如sm包中的sm.ancova()函数。

2、结果可视化

        HH包中的ancova()函数可以绘制因变量、协变量和因子之间的关系图。如下：

# 加载HH包
library(HH)
# 作图
ancova(weight~gesttime+dose,data=litter)

        从上图中，我们可以看出，用怀孕时间预测出生体重的回归直线是相互平行的，只是截距项不同。随着怀孕时间增加，幼崽出生体重也会增加。另外，还可以看到0剂量组截距项最大，5剂量组截距项最小。若用ancova(weight~gesttime*dose)，生成的图形将允许斜率和截距项依据组别而发生变化（如下图），这对可视化哪些违背回归斜率同质性的实例展示非常有用。