单因素协方差分析(One-Way ANCOVA)——R软件实现

关键词：R语言; R软件; 协方差分析; 单因素协方差分析; 平行性检验; 平行线检验

一、案例介绍

为研究A、B、C三种饲料对增加小鼠体重的影响，将初始体重相近的45只小鼠随机分成三组，分别喂养A、B、C三种饲料，但在实验设计时未对小鼠的进食量加以限制，现测得三组小鼠的进食量(Food)和所增体重(Weight)，请推断A、B、C三种饲料对小鼠的增重效果是否有差别？部分数据见图1。本文案例可从“附件下载”处下载。

二、问题分析

本案例的分析目的是比较A、B、C三种饲料对增加小鼠体重的影响，即A、B、C三组体重增量是否存在差异。但显然，每个老鼠的进食量会对体重增量产生影响，因此针对这种情况可将进食量作为体重增量的影响因素进行单因素协方差分析(One-Way ANCOVA)。但需要满足9个条件：

条件1：观察变量为连续变量。本研究中观察变量为体重增量，为连续变量，该条件满足。

条件2：自变量为二分类或多分类变量。本研究中自变量为A、B、C三组，该条件满足。

条件3：协变量是连续变量。本研究中协变量为进食量，为连续变量，该条件满足。

条件4：各研究对象之间具有独立的观测值。本研究中各个研究对象均为独立样本，不存在互相干扰的情况，该条件满足。

条件5：观察变量不存在显著异常值，该条件需要通过软件分析后判断。

条件6：各组内因变量的残差服从正态(或近似正态)分布。

条件7：组间因变量的残差方差齐，该条件需要通过软件分析后判断。

条件8：各组内协变量和因变量之间存在线性关系，该条件需要通过软件分析后判断。

条件9：各组内协变量和因变量的回归直线平行，即通过平行性检验，该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 导入数据

mydata <- read.csv("单因素协方差分析.csv")  #导入CSV数据
View(mydata)  #查看数据

在数据栏目中可以查看全部数据情况，数据集中共有3个变量和45个观察数据，3个变量分别代表被三组小鼠的进食量(Food)、所增体重(Weight)及分组(group)。

如果数据集较大也可使用如下命令查看数据框结构：

str(mydata)  #查看数据框结构

(二) 适用条件判断

条件5—条件7需要通过模型残差进行判断，因此先生成模型残差。

1. 生成模型残差及预测值

## 计算数据并提取残差值 ##

library(lmtest) #调用程序包“lmtest”
lmresults<-lm(Weight~Food*group,data=mydata) #将回归分析结果赋值给lmresults
res<-rstandard(lmresults) #提取标准化残差赋值给res
pre<-predict(lmresults) #提取预测值赋值给pre
mydata0<-data.frame(mydata,res,pre) #重新形成数据集
View(mydata0)  #查看数据

1. 条件5判断(异常值判断)

(1) 软件操作

## 绘制箱线图 ##

labels <- c("A","B","C")   #设置分组标签
boxplot(mydata0$res ~ mydata0$group, names = labels, 
        xlab = c("group"), ylab = expression("res"))

(2) 结果解读

图4残差的箱线图未提示任何异常值和极端值，满足条件5。

2. 条件6判断(正态性检验)

(1) 软件操作

## 正态性检验 ##

shapiro.test(mydata0$res)

(2) 结果解读

图5的正态性检验结果显示P值为0.8803 >0.1，可认为残差服从正态分布。此外，本案例也可以绘制Q-Q图，结果也提示残差满足正态性要求(请读者自行操作)。综上，本案例满足条件6。

3. 条件7判断(方差齐性检验)

(1) 软件操作

## 方差齐性检验 ##

library(car) #调用包“car”
mydata0$group <- factor(mydata0$group,labels=c("A","B","C")) #把变量group转换为分类变量
leveneTest(res~group,data=mydata0) #方差齐性检验

(2) 结果解读

由图6的“Levene’s Test for Homogeneity of Variance Test (Levene’s方差齐性检验)”结果可知，F=1.4175，P=0.2537，提示本案例组间因变量的方差齐，满足条件8。

4. 条件8判断(线性关系检验)

(1) 软件操作

## 绘制食物和体重的相关关系图 ##

library(ggplot2)
ggplot(data=mydata0, aes(x=Food, y=Weight,color=group))+geom_point()+stat_smooth(method="lm",se=TRUE)

## 拆分数据 ##

mydata1<-subset(mydata0, mydata0$group=="A")  #按组生成数据子集1
mydata2<-subset(mydata0, mydata0$group=="B")  #按组生成数据子集2
mydata3<-subset(mydata0, mydata0$group=="C")  #按组生成数据子集3

## 相关系数检验  ##

cor.test(mydata1$Weight,mydata1$Food,method = "pearson",alternative = "two.sided") #计算A组相关系数
cor.test(mydata2$Weight,mydata2$Food,method = "pearson",alternative = "two.sided") #计算B组相关系数
cor.test(mydata3$Weight,mydata3$Food,method = "pearson",alternative = "two.sided") #计算C组相关系数

(2) 结果解读

由图7可知，A、B、C三组内“Food”与“Weight”均呈现线性关系；图8—图10“cor.test (相关检验)”分析结果列出了详细的相关系数，均较高(>0.8)，P值均<0.001。综上，本案例各组内协变量和因变量之间均存在线性关系，满足条件8。

5. 条件9判断(平行性检验)

平行性检验通过判断group与Food的交互项是否有统计学意义决定。

(1) 软件操作

## 平行性检验  ##

library(multcomp)
fit <- aov(Weight ~ Food * group, data = mydata0)
summary(fit)

(2) 结果解读

“图11中Food*group”的交互项检验结果为，F_group*Food =0.186，P_group*Food=0.831，表示交互项无统计学意义，即满足平行性检验。因此，本案例满足条件9。

(三) 统计描述

由于“group*Food”的交互项无统计学意义，因此交互项不放入模型中。

1. 软件操作

## 计算标准差 ##

library(psych)  #调用包“psych”
describeBy(mydata$Weight, group = mydata$group)  #Weight变量分组描述统计

describeBy(mydata$Food, group = mydata$group)  #Food变量分组描述统计

2. 结果解读

图12中列出了A、B、C三组体重增量分别为36.63±20.37、46.43±21.39、119.03±31.94 g，图13中列出了A、B、C三组的进食量分别为273.22±46.25、274.15±41.66、493.59±73.80 g。

(四) 单因素协方差分析

1. 统计学描述

(1) 软件操作

## 计算边际均值 ##

lmresults1<-lm(Weight~Food+group,data=mydata0) #剔除交互项重新计算回归模型
library(emmeans) #调用程序包"emmeans"
em<-emmeans(lmresults1, ~ Food | group, data=mydata0) #计算边际均值并赋值给
em #查看边际均值

plot(em) #绘制边际均数图

(2) 结果解读

图14“emmeans(估算边际均值)”中提供了三组的“emmean [估算边际均值(偏最小二乘均值))”、“SE (标准误)”、“df(自由度)”、均值的“lower.CL (95%置信区间下限)”及“upper.CL(95%置信区间上限)”，A、B、C三组体重增量的估算边际均值分别为64.6 (95%CI: 55.6~73.7)、74.1 (95%CI: 65.0~83.1)、63.3 (95%CI: 50.1~76.6) g，可见三组体重增量的估算边际均值与三组体重增量的均值存在较大差异。图15的估算边际均值图也绘制了三组体重增量的情况，可见B组最高、A组其次、C组最低。但这种差异是否具有统计学意义，还需要依据统计学推断结果做出判断。

2. 统计学推断

(1) 软件操作

## 协方差分析 ##

fit2<-aov(Weight~Food+group,data=mydata0) #协方差分析
summary(fit2) #查看协方差分析结果

(2)结果解读

图16“aov (方差分析)” 结果显示，Food变量的F_Food=421.111，P_Food<0.001，提示Food的确对Weight具有影响，需要将其列为协变量进行控制。group变量的F_group=1.999，P_group=0.149，提示三组小鼠体重的增重效果差异无统计学意义。

(四) 事后检验(两两比较)

本案例中由于三组间整体比较差异无统计学意义，因此无需再进行事后两两比较。假设本案例三组间整体比较差异有统计学意义，其分析步骤如下。

1. 软件操作

TukeyHSD(fit2) #事后两两比较

2.结果解读

图17 “TukeyHSD (事后检验比较)”中提供了各组两两比较的“diff (均数差)”、 “lwr (95%置信区间下限)”、“upr (95%置信区间上限)”和“p.adj (调整过后的P值)”。可知，三组中任意两组比较差异均无统计学意义(P>0.05)。

四、结论

本研究采用单因素协方差分析，判断在调整进食量后A、B、C三种饲料对小白鼠的增重效果是否有差别。通过对模型残差绘制箱线图提示，数据不存在异常值；通过Shapiro-Wilk检验，提示各组因变量残差满足正态性要求；通过Levene’s检验，提示各组因变量残差满足方差齐性要求；通过绘制散点图和相关性分析，提示各组内协变量和因变量之间存在线性关系；通过平行性检验，发现“group*Food”的交互项无统计学意义(F_group*Food =0.186，P_group*Food=0.831)，提示满足平行性检验要求。

A、B、C三组小白鼠体重增量分别为36.63±20.37、46.43±21.39、119.03±31.94 g。单因素协方差分析显示，在调整进食量后A、B、C三组体重增量的估算边际均值分别为64.6 (95%CI: 55.6~73.7)、74.1 (95%CI: 65.0~83.1)、63.3 (95%CI: 50.1~76.6) g，差异无统计学意义(F_group=1.999，P_group=0.149)。本研究结果提示A、B、C三种饲料对小白鼠体重的增重效果无差异。

五、知识小贴士

• 协方差分析是针对在实验设计阶段难以控制其取值水平，或者无法严格控制的因素，在统计分析阶段对其进行统计控制的一种分析方法，实质为线性回归分析和方差分析的结合。适用于完全随机设计、随机区组(配伍)设计、拉丁方设计、析因设计等类型的方差分析。
• 协方差分析一般要求协变量在组间的观察范围相差不宜太大(分析前最好先对协变量均数间的差别作假设检验)，否则修正后的边际均值的差值可能会落在回归线的延长线上，回归线外推后，是否仍然满足平行线和线性关系尚不可知，其协方差分析的结论可能不一定准确。

六、分析小技巧

• 正态性检验：对协方差分析的正态性条件考察应检验因变量残差的正态性，此时检验多组整体残差的正态性即可，无须检验每组残差的正态性。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。
• 方差齐性检验：对协方差分析的方差齐性条件考察应检验组间因变量残差是否相等。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。