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

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

一、案例介绍

研究A、B两种镇痛药物联合运用对癌症患者的镇痛效果。药物A取3个剂量：1.0、2.5、5.0 mg；药物B也取3个剂量：5、15、30 μ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.08959，>0.05且较为接近0.1，提示因变量残差近似服从正态分布。本案例满足条件5。

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

(1) 软件操作

## 方差齐性检验 ##

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

(2) 结果解读

图6 Levene’s方差齐性检验结果显示，F=0.3043、P=0.9546；提示数据总体方差相等。本案例满足条件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水平时，镇痛时间分别为：79.00±4.00、94.33±2.52、108.67±4.51min；A药在2.5 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：90.67±4.16、106.67±4.04、116.67±2.52 min；A药在5.0 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：108±4.58、118.33±5.51、132.00±4.00 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("emdata.png",width=8,height=5,dpi=300) #保存图片

2. 结果解读

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

(五) 交互作用判断

由于本案例有两个因素(一个为A药，另一个为B药)，因此需要首先判断两个因素(药物)之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。

1. 软件操作

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

2. 结果解读

图10“ANOVA-Time (方差分析-时间)”显示，Drug A与Drug B之间的交互作用结果为F_{Drug A*Drug B}=0.699，P=0.603，提示Drug A与Drug B之间的交互作用无统计学意义。因此，本案例可直接以主效应分析结果进行判断，如果交互作用有统计学意义，则需要分析单独效应。

(六) A药主效应分析

图10“ANOVA-Time (方差分析-时间)”显示，A药的主效应检验结果为，F_{Drug A}=88.162，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。

(七) B药主效应分析

图10 “ANOVA-Time (方差分析-时间)”显示，B药的主效应检验结果为，F_{Drug B}=95.269，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。

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

上面分析得出了A、B药物在不同药物浓度时，镇痛时间差异均有统计学意义的结论，因此需要进行事后检验，开展两两比较。

1. 软件操作

## 事后检验 ##

AOVa<-aov(Time~a, mydata) #检验药物A不同组之间的差异
TukeyHSD(AOVa, p.adjust.methods="bonferroni") #A药事后两两比较

AOVb<-aov(Time~b,mydata) #检验药物B不同组之间的差异
TukeyHSD(AOVb, p.adjust.methods="bonferroni") #B药事后两两比较

2. 结果解读

图11“TukeyHSD(事后检验比较-Drug A)”结果提供了A药不同浓度间两两比较的“diff (均数差)”、 “lwr (95%置信区间下限)”、“upr (95%置信区间上限)”和“p.adj (调整过后的P值)”。可知，2.5 mg和1.0 mg与5.0 mg相比，均值差逐渐增大。5 mg比2.5 mg时均值增加了14.778 min，比1.0 mg时均值增加了25.444  min，差异有统计学意义(P<0.05)。表明：随着A药药物浓度的增加，镇痛时间呈上升趋势。

图12“TukeyHSD(事后检验比较-Drug A)”结果提供了A药不同浓度间两两比较的“diff (均数差)”、 “lwr (95%置信区间下限)”、“upr (95%置信区间上限)”和“p.adj (调整过后的P值)”。可知，15.0 μg和30.0 μg与5.0 μg相比，均值差逐渐增大。15.0 μg比5.0 μg时均值增加了13.889 min，至30.0 μg时增加了26.556 min，差异均有统计学意义(P<0.05)。表明：随着B药药物浓度的增加，镇痛时间呈上升趋势。

四、结论

本研究采用两因素方差分析探讨A、B两种镇痛药物联合运用在癌症患者的镇痛效果。通过专业知识判断，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过Levene’s检验，提示数据总体方差相等；两药物之间无交互作用(F_{Drug A*Drug B}=0.699，P=0.603)，故进行主效应分析。

A药在1.0 mg水平，B药为5.0、15.0、30.0 μg水平时，镇痛时间分别为：79.00±4.00、94.33±2.52、108.67±4.51min；A药在2.5 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：90.67±4.16、106.67±4.04、116.67±2.52 min；A药在5.0 mg水平，B药为5、15、30 μg水平时，镇痛时间分别为：108±4.58、118.33±5.51、132.00±4.00 min。A药的主效应检验结果为，F_{Drug A}=88.162，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。B药的主效应检验结果为，F_{Drug B}=95.269，P<0.001，认为不同药物浓度时，镇痛时间差异有统计学意义。进一步采用“Tukey's HSD”校正法进行两两比较可知，A药2.5 mg和1.0 mg与5.0 mg相比，均值差逐渐增大，差异均有统计学意义(P<0.05)；表明：随着A药药物浓度的增加，镇痛时间呈上升趋势。B药15.0 μg和30.0 μg与5.0 μg相比，均值差逐渐增大，差异均有统计学意义(P<0.05)。表明：随着B药药物浓度的增加，镇痛时间呈上升趋势。综上可知，不同剂量药物A的镇痛效果不同，不同剂量药物B的镇痛效果不同，A、B两种药物不存在交互作用。

五、分析小技巧

正态性检验：两因素或多因素方差分析时，有两种选择来测试正态性：如果每组有较多观察数，且组别较少时，可使用原始数据检查每个组的正态性。如果有很多组，或每个组的观察数很少，可使用残差检查整体的正态性。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。

交互作用判断：两因素方差分析时，需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。此时，单纯研究某个因素的作用并无意义，应分别探讨另一个因素不同水平时对该因素的作用。当不存在交互作用时，说明两因素的作用彼此独立，逐一分析各因素的主效应即可；计算主效应时，在模型中仍需要保留交互项。