随机区组设计(Randomized Block Design)资料的方差分析—R软件实现

关键词：R; 随机区组设计; 方差分析; 仅研究主效应的实验设计

随机区组设计是配对设计的扩展。其具体做法是先按影响实验结果的非处理因素(如性别、年龄、体重、病情程度等)将实验对象配成区组(block)，再分别将各区组内的实验对象随机分配到各处理组。若将区组视为一种处理因素的不同水平，则随机区组设计等同于两因素设计。

一、案例介绍

欲研究A (drug=1)、B (drug=2)、C (drug=3)3种药物(drug)对贫血患儿的治疗效果。将24例贫血患儿按照贫血程度分成8个区组(block)，每个区组的3名患儿分别随机接受A、B、C3种药物治疗。记录治疗30天后每名患儿血红蛋白(Hb)的增量(g/L)。部分数据见图1。本案例数据可从“附件下载”处下载。

二、问题分析

本案例是随机区组设计的方差分析，可以按照两因素方差分析进行数据分析。但需要满足6个条件：

条件1：观察变量唯一，且为连续变量。本研究中观察变量只有血红蛋白的增量，且为连续变量，该条件满足。

条件2：有两个因素，且都为分类变量。本研究中有处理效应(药物)及区组(贫血程度)两个因素，都为分类变量，该条件满足。

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

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

条件5：各组、各水平观测值为正态(或近似正态)分布或因变量残差整体服从正态分布，该条件需要通过软件分析后判断。

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

三、软件操作及结果解读

(一) 导入数据

语句如下，结果见图2。

mydata = read.csv("随机区组设计.csv", header = T)
str(mydata)

(二) 适用条件判断

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

1. 生成模型残差

##重新定义分组变量##

mydata$Block <- factor(mydata$block)  #将变量“block”设置为因子变量
mydata$Drug <- factor(mydata$drug, labels = c("A", "B", "C"))  #将变量“drug”设置为因子变量并为各水平赋值标签

见图3。

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

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

见图4。

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

(1) 软件操作

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

(2) 结果解读

方差齐性检验结果见图5。可知F=3.3026，P=0.05661，提示各处理水平的总体方差相等。满足条件4。

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

(1) 软件操作

##正态性检验##

tapply(mydata$res, mydata$drug, shapiro.test) 

(2) 结果解读

正态性检验结果见图6。可知A，B，C三种药物的p值分别为0.1655、0.7382、0.6628，均>0.1，提示因变量残差整体服从正态分布。满足条件5。

4. 条件6判断(异常值检测)

(1) 软件操作

##绘制箱线图##

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

(2) 结果解读

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

(三) 统计描述

1. 软件操作

##数据分组描述(图8)##

library(psych) #调用包“psych”
describeBy(mydata$Hb, group = mydata$Drug)

2. 结果解读

图8列出了处理组各组的均数和标准差，可知经A、B、C药治疗后贫血患儿血红蛋白的平均增量分别为9.16±1.20、10.28±1.13、11.58±0.94 (g/L)。

(四) 两因素方差分析

1. 软件操作

##查看方差分析结果(图9)##

summary(AOV)

##拟合线性回归模型##

blrfit = lm(Hb ~ Block + Drug, data = mydata)
summary(blrfit)

##重新定义参照组##

mydata$Drug <- relevel(mydata$Drug, ref = "C") #将药物C设置为参照组
blrfit1 = lm(Hb ~ Block + Drug, data = mydata)
summary(blrfit1)

2. 结果解读

随机区组设计的方差分析结果见图9，可知区组因素间与处理因素间的差异均有统计学意义(P<0.001)。

将以上两因素一同纳入线性回归模型分析，可知以药物A为参照(图10)，药物B的疗效平均比药物A高1.112 g/L (β=1.112)和药物C的疗效平均比药物A高2.412 g/L (β=2.412)，差异均有统计学意义(P<0.001)。若以药物C为参照(图11)，药物C的疗效平均比药物B高1.300 g/L(β=-1.300)，药物A和药物C疗效比较结果同前，差异均有统计学意义(P<0.001) 。

四、结论

本案例基于R软件实现随机区组设计的方差分析。案例欲研究A、B、C3种药物对贫血患儿的治疗效果，采用随机区组设计。经A、B、C药治疗后贫血患儿血红蛋白的平均增量分别为9.16±1.20、10.28±1.13、11.58±0.94 g/L。区组因素间与处理因素间的差异均有统计学意义(P<0.001)。药物B的疗效平均比药物A高1.112 g/L，药物C的疗效平均比药物A高2.412 g/L，药物C的疗效平均比药物B高1.300 g/L(β=-1.300)，差异均有统计学意义(P<0.001)。