两因素方差分析(Two-Way ANOVA)二(有交互作用)——R软件实现

关键词：R语言; R软件; 两因素方差分析; 交互作用; 主效应; 单独效应; 简单效应

一、案例介绍

观察A、B两种镇痛药物联合运用在产妇分娩时的镇痛效果。A药(Drug_A)取3个剂量：1.0、2.5、5.0；B药(Drug_B)也取3个剂量：5.0、15.0、30.0μg，共9个处理组。将27名产妇随机分成9组，每组3名产妇，记录每名产妇分娩时的镇痛时间Time (min)。试分析A、B两药联合运用的镇痛效果。部分数据见图1。本文案例可从“附件下载”处下载。

二、问题分析

本案例的分析目的是分析A、B两种药物联合运用的镇痛效果。临床上，药物之间联合运用往往会相互影响，这种影响可能为正向的增强效应，也可能为反向的拮抗作用。针对这种情况，可以使用多因素方差分析，由于本案例为药物A和药物B两个因素，因此可以使用两因素方差分析。但需要满足6个条件：

条件1：观察变量唯一，且为连续变量。本研究中观察变量只有镇痛时间，且为连续变量，该条件满足。

条件2：有两个因素，且都为分类变量。本研究中有A、药物B物两个因素，都为分类变量，该条件满足。

条件3：观测值之间相互独立。本研究中各研究对象的观测值都是独立的，不存在互相干扰的情况，该条件满足。

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

条件5：各组、各水平观测值为正态(或近似正态)分布，该条件需要通过软件分析后判断。

条件6：相互比较的各处理水平(组别)的总体方差相等，即方差齐同，该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 导入数据

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

在数据栏目中可以查看全部数据情况，数据集中共有3个变量和27个观察数据，3个变量分别代表A药物(Drug_A)、B药物(Drug_B)和镇痛时间(Time)。 用如下命令查看数据框结构：

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

(二) 适用条件判断

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

1. 生成模型残差

##将用药量重新定义为分组变量##

mydata$a <-factor(mydata$Drug_A,labels = c("1mg","2.5mg","5mg"))
mydata$b <-factor(mydata$Drug_B,labels = c("5μg","15μg","30μg"))

##数据整理及残差值提取 ##

AOV<-aov(Time~a*b,mydata) #计算方差分析模型
res<-rstandard(AOV) #提取标准化残差赋值给res
pre<-predict(AOV) #提取预测值赋值给pre
mydata<-data.frame(mydata,res,pre) #重新形成数据集
View(mydata)  #查看数据

2. 条件4判断(异常值判断)

(1) 软件操作

## 绘制箱线图 ##

library(ggplot2)  #调用包“ggplot2”
boxp<- ggplot(mydata,aes(x=b, y=res, color=a))+ #定义颜色与分组
  stat_boxplot(geom='errorbar')+ #添加误差线
  geom_boxplot() #绘制箱线图
ggsave("boxp.png", width=8, height=5, dpi=300) #保存图片

(2) 结果解读

残差的箱线图将被保存在R数据分析路径下，图4残差的箱线图未提示任何异常值和极端值，满足条件4。

3. 条件5判断(正态性检验)

(1) 软件操作

##正态性检验 ##

shapiro.test(mydata$res)

(2) 结果解读

图5正态性检验结果显示P= 0.1582，> 0.1，提示因变量残差近似服从正态分布。本案例满足条件5。

4. 条件6判断(方差齐性检验)

(1) 软件操作

##方差齐性检验 ##

library(car) #调用包“car”
leveneTest(res~a*b,mydata,center="mean") #方差齐性检验

(2) 结果解读

图6 Levene’s方差齐性检验显示，F=0.6957、P=0.6912；提示数据总体方差相等。本案例满足条件6。

(三) 统计描述

1. 软件操作

##数据分组描述  ##

library(psych)  #调用包“psych”
describeBy(mydata$Time, group = list(mydata$Drug_A,mydata$Drug_B))

##统计描述 ##

lm<-lm(Time~a*b,data=mydata) #计算回归模型
library(emmeans) #调用程序包"emmeans"
em<-emmeans(lm, ~ a*b, data=mydata) #计算边际均值并赋值给em
em #查看边际均值

2. 结果解读

图7-1至图7-3列出了各组的均值和标准差，可知A药在1.0 mg水平，B药为5.0、15.0、30.0 μg水平时，镇痛时间分别为：83±20.21、100±18.03、85±10.00 min；A药在2.5 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：90±21.80、115±21.80、135±15.00 min；A药在5.0 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：110±21.80、95±27.84、177±15.28 min。

 (四) 两因素方差分析

1. 软件操作

##绘制边际均数图 ##

emdata<-data.frame(em)
emplot<-ggplot(emdata,aes(x = a, y=emmean, color=b,group=b))+ #绘制分组图
  geom_errorbar(aes(ymin=lower.CL,ymax=upper.CL),width=.05)+ #添加误差线
  geom_line()+  #添加连线
  geom_point()  #添加点
ggsave("emplot.png",width=8,height=5,dpi=300) #保存图片

2. 结果解读

图9的估算边际均值图绘制了各组镇痛时间的变化情况，可见两药物的变化曲线趋势并不平行，存在明显的交互作用。图8“emmeans(估算边际均值)”中提供了数据的“emmean (估算边际均值(偏最小二乘均值))”、“SE (标准误)”、“df(自由度)”、均值的“lower.CL (95%置信区间下限)”及“upper.CL(95%置信区间上限)”。

(五) 交互作用判断

由于本案例有两个因素(一个为A药，另一个为B药)，因此需要首先判断两个因素(药物)之间是否存在交互作用。尽管图9的变化趋势已经提示两药物之间存在明显的交互作用，但仍需要统计学推断结果的支持。如果交互作用有统计学意义，则需要分析单独效应。

1. 软件操作

summary(AOV) #查看方差检验结果

2. 结果解读

图10“ANOVA-Time (方差分析-时间)”显示，Drug_A与Drug_B之间的交互作用结果为F_{Drug_A*Drug_B}=5.073，P=0.00647，提示Drug_A与Drug_B之间的交互作用有统计学意义，本案例需要分析单独效应。

(六) 单独效应分析

1. 软件操作

##单独效应分析 ##

install.packages("bruceR") #安装程序包"bruceR"
library(bruceR) #调用程序包"bruceR"
library(magrittr) #调用程序包"magrittr"
library(dplyr) #调用程序包"dplyr"
MANOVA(data=mydata,dv="Time",between=c("Drug_A","Drug_B")) %>%
  EMMEANS("Drug_A",by="Drug_B") %>%
  EMMEANS("Drug_B",by="Drug_A")

2. 结果解读

(1) A药的单独效应

由图11“Simple Effects of ‘Drug_A’ (A药的单独效应)”可知，在B药为5.0 μg水平时，A药各水平之间的镇镇痛时间差异无统计学意义(F=1.489，P=0.252)；在B药为15.0 μg水平时，A药各水平之间的镇痛时间差异无统计学意义(F=0.838，P=0.449)；但在B药为30.0 μg水平时，A药各水平之间的镇痛时间差异有统计学意义(F=16.289，P<0.001)。图12的结果和图7结果一致，为各组的“Mean (估算边际均值(偏最小二乘均值))”、“S.E. (标准误)”及均值的“95% Confidence Interval (95%CI)”。

(1) B药的单独效应

由图13“Simple Effects of ‘Drug_B’ (B药的单独效应)”可知，在A药为1.0 mg水平时，B药各水平之间的镇痛时间差异无统计学意义(F=0.652，P=0.533)；在A药为2.5 mg水平时，B药各水平之间的镇痛时间差异有统计学意义(F=3.931，P=0.038)；但在A药为5.0 mg水平时，B药各水平之间的镇痛时间差异有统计学意义(F=14.613，P<0.001)。图14的结果和图7结果一致，为各组的“Mean (估算边际均值(偏最小二乘均值))”、“SE (标准误)”及均值的“95% Confidence Interval (95%CI)”。

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

1. A药的单独效应事后检验

上面分析得出了“在B药为30.0 μg水平时，A药各水平之间的镇痛时间差异有统计学意义”的结论。由于A药在B药为5.0 μg和15.0 μg水平下的整体检验差异无统计学意义，此时便不再进行事后两两比较。现对A药在B药为30.0 μg水平时进行事后检验，结果如图15所示。

图15“Pairwise Comparisons of ‘Drug_A’ (事后检验比较‘Drug_A’)”表格中提供了在B药在不同水平时A药不同浓度间通过“Bonferroni (邦弗罗尼)”法两两比较的“均数差”、“SE (标准误)”、“df (自由度)”、“t (统计量t值)”、“P* (校正P值)”。可知，A药浓度为2.5 mg和5.0 mg与1.0 mg相比，均值差逐渐增大。2.5 mg比1.0 mg时均值增加了50.00 min，至5.0 mg时增加了91.67 min，差异均有统计学意义(P<0.001)；2.5 mg到5.0 mg时增加了41.67 min，差异无统计学意义(P=0.055)。表明，在B药为30.0 μg水平时，随着A药浓度的增加，镇痛时间呈上升趋势。

2. B药的单独效应事后检验

上面分析得出了“在A药为2.5 mg和5.0 mg水平时，B药各水平之间的镇痛时间差异均有统计学意义”的结论。由于B药在A药为1.0 mg水平下的整体检验差异无统计学意义，此时便不再进行事后两两比较。现对B药在A药为2.5 mg和5.0 mg水平时进行事后检验，结果如图16所示。

图16 “Pairwise Comparisons of ‘Drug_B’ (事后检验比较‘Drug_B’)”表格中提供了在A药在不同水平时B药不同浓度间通过“Bonferroni (邦弗罗尼)”法两两比较的“均数差”、“SE (标准误)”、“df (自由度)”、“t (统计量t值)”、“P* (校正P值)”。 可知，当A药为2.5 mg时，B药浓度30.0 μg与5.0 μg相比，镇痛时间差异有统计学意义(P=0.036)，30.0 μg与15.0 μg相比以及15.0 μg与5.0 μg相比差异均无统计学意义(P>0.05)。当A药为5.0 mg时，B药浓度30.0 μg与15.0 μg相比，镇痛时间差异有统计学意义(P<0.001)，30.0 μg与5.0 μg相比，镇痛时间差异有统计学意义(P=0.002)，但15.0 μg与5.0 μg相比，镇痛时间差异无统计学意义(P=1.000)。

四、结论

本研究采用两因素方差分析探讨A、B两种镇痛药物联合运用在产妇分娩时的镇痛效果。通过对模型残差绘制箱线图提示，数据不存在异常值；通过Shapiro-Wilk检验，提示残差服从正态分布；通过Levene’s检验，提示数据总体方差相等；两药物之间存在交互作用(F_{Drug_A*Drug_B}=5.073，P=0.00647)，故进行简单效应分析。

A药在1.0 mg水平，B药为5.0、15.0、30.0 μg水平时，镇痛时间分别为：83±20.21、100±18.03、85±10.00 min；A药在2.5 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：90±21.79、115±21.79、135±15.00 min；A药在5.0 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：110±21.79、95±27.84、177±15.28 min。

A药的简单效应分析显示，在B药为30.0 μg水平时，A药各水平之间的镇痛时间差异有统计学意义(F=16.289，P<0.001)。但在B药为5.0 μg和15.0 μg水平时，A药各水平之间的镇痛时间差异均无统计学意义(P>0.05)。对A药在B药为30.0 μg水平时进行“Bonferroni (邦弗罗尼)”法两两比较显示，1.0 mg至2.5 mg时均值增加了50.00 min，至5.0 mg时增加了91.67 min，差异均有统计学意义(P<0.001)；表明，在B药为30.0 μg水平时，随着A药浓度的增加，镇痛时间呈上升趋势。

B药的简单效应分析显示，在A药为1.0 mg水平时，B药各水平之间的镇痛时间差异无统计学意义(F=0.652，P=0.533)；但在A药为2.5 mg和5.0 mg水平时，B药各水平之间的镇痛时间差异均有统计学意义(P<0.05)。“Bonferroni (邦弗罗尼)”法两两比较显示，当A药为2.5 mg时，B药浓度30.0 μg与5.0 μg相比，镇痛时间差异有统计学意义(P=0.036)，30.0 μg与15.0 μg相比以及15.0 μg与5.0 μg相比差异均无统计学意义(P>0.05)。当A药为5.0 mg时，B药浓度30.0 μg与15.0 μg相比以及30.0 μg与5.0 μg相比，镇痛时间差异均有统计学意义(P<0.01)，但15.0 μg与5.0 μg相比，镇痛时间差异无统计学意义(P=1.000)。

综上可知，对产妇分娩时镇痛，A、B两种药物联合运用时会相互影响，随着药物浓度的增加，联合作用趋于复杂和不稳定。

五、分析小技巧

• 正态性检验：两因素或多因素方差分析时，有两种选择来测试正态性：如果每组有较多观察数，且组别较少时，可使用原始数据检查每个组的正态性。如果有很多组，或每个组的观察数很少，可使用残差检查整体的正态性。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。
• 方差齐性检验：两因素方差分析时，使用残差检验方差齐性。可绘制残差图(“rvfplot”命令)或使用BP(Breusch-Pagan)法(“estat hettest”命令)检验方差齐性。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。
• 交互作用判断：两因素方差分析时，需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。当不存在交互作用时，说明两因素的作用彼此独立，逐一分析各因素的主效应即可；计算主效应时，在模型中仍需要保留交互项。