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

关键词：R语言; R软件; 单因素方差分析; F检验; Welch检验; 韦尔奇检验; 事后检验; 两两比较

一、案例介绍

某医生用A、B、C三种方案治疗血红蛋白低下的贫血患者，治疗两个月后，记录每名受试者血红蛋白的上升克数。问3种治疗方案对患者贫血的疗效是否有差别？部分数据见图1。本文案例可从“附件下载”处下载。

二、问题分析

本案例的分析目的是比较3种治疗方案对患者贫血的疗效是否有差别，即判断3种治疗方案患者的血红蛋白上升克数是否存在差异。针对这种情况可以使用单因素方差分析。但需要满足6个条件：

条件1：观察变量为连续变量。本研究中的血红蛋白上升克数为连续变量，该条件满足。

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

条件3：观测值可分为多组(≥3)。本研究中分为A、B、C三组，该条件满足。

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

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

条件6：多组观测值的整体方差相等，该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 导入数据

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

在数据栏目中可以查看全部数据情况，数据集中共有3个变量和78个观察数据，3个变量分别代表被调查者的编号(ID)、治疗分组(group)及血红蛋白的上升克数(Hb)。
如果数据集较大也可使用如下命令查看数据框结构：

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

(二) 适用条件判断

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

(1) 软件操作

## 拆分数据##

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

##分组描述数据子集##

summary(mydata1$Hb)  #描述子集1基本情况
summary(mydata2$Hb)  #描述子集2基本情况
summary(mydata3$Hb)  #描述子集3基本情况

##分组查看缺失值情况##

is.na(mydata1$Hb)  #查看子集1有无缺失值
is.na(mydata2$Hb)  #查看子集2有无缺失值
is.na(mydata3$Hb)  #查看子集3有无缺失值

## 分组绘制箱线图##

labels <- c("A","B","C")   #设置分组标签
boxplot(mydata$Hb ~ mydata$group, names = labels, 
        xlab = c("group"), ylab = expression("Hb"))  #分组绘制箱线图

(2) 结果解读

图3“summary (描述性分析)”命令运行结果，列出了观察变量的“Min.(最小值)、“1st QU.(第一四分位数)”、“Median(中位数)”、“Mean(平均值)”、“3st QU.(第三四分位数)和“Max.(最大值)”，依据专业可判断血红蛋白的上升克数均可能存在0.1 g和3.7 g的情况；此外，图5中的箱线图也未提示任何异常值。综上，本案例未发现需要删除的异常值，满足条件4。

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

(1) 软件操作

##分组绘制Q-Q图##

par(mfrow = c(1, 3))  #设置1行3个图
qqnorm(mydata1$Hb, ylab="Hb", main="A")  #绘制子集1的qq图
qqline(mydata1$Hb)  #增加趋势线
qqnorm(mydata2$Hb, ylab="Hb", main="B")  #绘制子集2的qq图
qqline(mydata2$Hb)  #增加趋势线
qqnorm(mydata3$Hb, ylab="Hb", main="C")  #绘制子集3的qq图
qqline(mydata3$Hb)  #增加趋势线

##正态性检验 ##

tapply(mydata$Hb, mydata$group, shapiro.test) #分组进行Shapiro-Wilk正态性检验

(2) 结果解读

图6 Q-Q图上三组散点基本围绕对角线分布，提示三组数据呈正态分布；图7的“Shapiro-Wilk normality test (S-W正态性检验)”结果分别显示三组的P=0.296、0.486、0.435，均>0.1，也提示三组数据服从正态分布。综上，本案例满足条件5。

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

(1) 软件操作

##计算标准差##

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

## levene法方差齐性检验##

library(dplyr)  #调用软件包“dplyr”
mydata4 <- mydata %>% mutate(group1=factor(group,labels=c("A","B","C")))  #把变量group转换为分类变量

library(car)  #调用软件包“car”
leveneTest(Hb~group1,data=mydata4,center=mean)  #levene法方差齐性检验

##分组描述性分析——均值标准差图##

install.packages(“ggpubr”)  #安装“ggpubr”包
library(ggpubr)  #调用包“ggpubr”
ggerrorplot(mydata, x = "group", y = "Hb", desc_stat = "mean_sd") #画平均值标准差图

(2) 结果解读

图8“describeBy (描述性分析)”运行结果提供了各组观察变量的统计量，包括：“vars(能够返回变量的数量)”、“n (样本量)”、“Mean (均数)”、“sd (标准差)”、“median (中位数)”、“terimmed(截尾均值)”、“绝对中位差(mad)”、“min(最小值)”、“max(最大值)”、range(全距)、“skew(偏度)”、“kurtosis(丰度)”和“se (标准误)”。结果可知，A、B、C三组的标准差分别为0.87、0.74、0.54 g，三组的标准差数值存在差异，但还需要依据统计学检验的结果进行判断。图9“LeveneTest(Levene检验)”结果显示，F=2.6102，P=0.08019<0.1，提示三组数据的方差不齐，不满足条件6。图10为三组人群血红蛋白指标含量的分布图。

(三) 统计描述及推断

1. 软件操作

##单因素方差分析 ##

library(rstatix)
welch_anova_test(data = mydata,Hb~group)  #校正的单因素方差分析

2. 结果解读

由于本案例满足条件1-5，但不满足方差齐性，所以采用校正单因素方差分析(Welch’s)法的结果，F=5.48，P=0.007；提示各组均数不全相等(至少有两组均数不相同)。

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

上面分析得出了“各组均数不全相等”的结论，但是具体是哪些组别之间存在差异尚不清楚，因此需要进行事后检验，开展两两比较。

1. 软件操作

## 两两比较 ##

games_howell_test (data=mydata4,Hb~group1,conf.level = 0.95, detailed = FALSE)

#注：当不满足方差齐性时采用“Games-Howell”法进行事后检验

2. 结果解读

图12“Games-Howell Post-Hoc Test-Hb (Games-Howell法事后检验)”表格中提供了各组两两比较的“estimate (均数差:group2-group1)”、“conf.low (95%置信区间下限)”、“conf.high (95%置信区间上限)”和“p.adj (调整过后的P值)”。可知，A组和B组之间的差异(B-A)为-0.438 (95%CI：-0.978~0.101) g，但差异无统计学意义(P=0.132)；A组和C组之间的差异(C-A)为-0.665 (95%CI：-1.15~-0.179) g，差异有统计学意义(P=0.005)；B组和C组之间的差异(C-B)为-0.227 (95%CI：-0.661~0.207) g，差异无统计学意义(P=0.421)。因此，可知A组血红蛋白的上升克数比C组高0.665 g，差异有统计学意义(P=0.005)。

四、结论

本研究采用单因素方差分析判断3种方案治疗患者贫血的血红蛋白上升克数是否存在差异。通过专业知识判断，三组数据不存在需要删除的异常值；通过绘制Q-Q图和Shapiro-Wilk检验，提示三组数据满足正态分布；通过Levene’s检验，提示三组数据间方差不齐，采用校正单因素方差分析(Welch’s)法分析数据。

分析结果显示，A、B、C三组血红蛋白的上升克数分别为1.69±0.87、1.25±0.74、1.02±0.54 g，三组均数不全相等(F=5.48，P=0.007)。进一步采用Games-Howell法进行两两比较，显示A组血红蛋白的上升克数比C组高0.665 g，差异有统计学意义(P=0.005)；A组血红蛋白的上升克数比B组高0.438 g，但差异无统计学意义(P=0.132)；B组血红蛋白的上升克数比C组高0.227 g，但差异无统计学意义(P=0.421)。因此，A、B、C三种治疗方案对治疗贫血患者疗效不一样，A方案疗效最好。

五、分析小技巧

(一) 正态性检验

• 严格来讲，单因素方差分析时，需要分别对每一组数据的正态性进行检验。但方差分析对数据的分布具有一定的耐受力，如果数据不是严重偏态或者只有部分组别数据不满足正态性要求，鉴于参数检验的统计学效能优于非参数检验的角度，还是可以使用单因素方差分析方法，而不使用非参数检验。
• 单因素方差分析时，也可考察整体数据的正态性情况。如上所述，单因素方差分析对数据分布具有一定的耐受力，如果数据满足整体正态分布，也是可以使用该种分析方法。但我们建议尽可能分组别检验数据的正态性。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。
  

(二) 方差齐性检验

• 对于不是特别严重的方差不齐，单因素方差分析提供了校正检验方法(Welch’s)，是考虑了方差差异之后的更为稳健的分析结果。但当组间方差差异较大时，校正结果也不一定可信，建议使用非参数检验(Kruskal-Wallis H检验)。如果数据正态性和方差齐性都不满足，最好使用非参数检验(Kruskal-Wallis H检验)。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。

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

• 多重比较一般分为事前检验(Prior tests)和事后检验(Post hoc tests)。事前检验是指在数据收集之前便决定了要通过多重比较来考察多个组与某个特定组之间的差别，多根据专业意义设定比较的策略。如果是事前检验，不论整体分析的结果如何，均可进行比较，并且一般不需要对检验水准进行太多修正。事后检验只有在方差分析得到有统计学意义的F值后才有必要进行，是一种探索性分析。对于事先未计划的多重比较(即事后检验)，各组间的差别只是一种提示，要确认这种差别最好重新设计实验。
• 事前检验(Prior tests)最常用的两两比较方法有：LSD法和Dunnett-t检验。事后检验(Post hoc tests)最常用的方法有：SNK法、Duncan法、Turkey’b和Scheffe法。
• LSD法是最灵敏的方法，因此也容易犯假阳性错误。Dunnett-t法多用于多个实验组与一个对照组比较，此时应指定对照组，多用于确证性研究，少用于探索性研究。Bonferroni法是对LSD法的严格校正，结果更加保守，但当组数较多时，较难发现组间差异，因此如果各组例数相差不大且组数不多时可采用。
• 当不满足方差齐性时常采用“Games-Howell”法进行事后检验。