在前面文章中介绍了随机区组设计方差分析(Randomized Block Design ANOVA)——R软件实现，本文进一步介绍仅研究主效应的实验设计中的两阶段交叉设计方差分析(Two-stage cross-over design ANOVA)在R软件中的实现步骤。
关键词：R; 交叉设计; 两阶段交叉设计; 方差分析; 仅研究主效应的实验设计
在医学研究中，若将A、B两种处理先后施加于同一批实验对象，随机地使半数实验对象先接受A后接受B，另半数对象先接受B后接受A，这种设计称为两阶段交叉设计(crossover design)。应用方差分析的变异分解思想，此种设计的总变异可以分解为处理效应、阶段效应、受试者间的个体差异及随机误差共4项。需要注意的是，该设计有一个限制条件，即前一个阶段的处理效应不能持续到下一个阶段。因此在实际研究设计时，往往在两阶段间设置一个洗脱阶段。
一、案例介绍
假设有一项两阶段交叉试验，分别使用A、B两种药物治疗同一批共10例患者，检测某指标的血浆浓度。对10例患者随机编号，第一阶段样本编号(testid)为奇数号的5例患者用A药物治疗，样本编号(testid)为偶数号的5例患者用B药物治疗。第二阶段对调治疗药物，即样本编号(testid)为奇数号的患者用B药物，样本编号(testid)为偶数号的药物用A药物。试对该两阶段交叉实验结果进行方差分析，评估药物效果。数据见图1。本案例数据可从“附件下载”处下载。
二、问题分析
在两阶段交叉数据分析，可以按照三因素方差分析思路进行数据处理。但需要满足6个条件：
条件1：观察变量唯一，且为连续变量。本研究中观察变量只有某指标的血浆浓度，且为连续变量，该条件满足。
条件2：有3个因素，且都为分类变量。本研究中有处理效应(A、B两种药物)、阶段效应及受试者间的个体差异，都为分类变量，该条件满足。
条件3：观测值相互独立。本研究中各研究对象的观测值都是独立的，不存在互相干扰的情况，该条件满足。
条件4：相互比较的各处理水平(组别)的总体方差相等，即方差齐同，可采用方差齐性检验。该条件需要通过软件分析后判断。
条件5：各组、各水平观测值为正态(或近似正态)分布.。该条件需要通过软件分析后判断。
条件6：观察变量不存在显著的异常值。该条件需要通过软件分析后判断。
三、软件操作及结果解读
(一) 导入数据(图2)
mydata = read.csv("交叉设计方差分析.csv", header = T)
str(mydata)
##设置变量类型，添加标签##
mydata$testid <- factor(mydata$testid)
mydata$drug <- factor(mydata$drug, levels = c(1, 2), labels = c("drug_A", "drug_B"))
mydata$phase <- factor(mydata$phase, levels = c(1, 2), labels = c("第一阶段", "第二阶段"))
(二) 适用条件判断
条件4~6需要通过模型残差进行判断，因此先生成模型残差。
1. 生成模型残差(图3)
##数据整理及残差值提取## AOV <- aov(contain ~ testid + drug + phase, mydata) #计算方差分析模型 res <- rstandard(AOV) #提取标准化残差赋值给res pre <- predict(AOV) #提取预测值赋值给pre mydata <- data.frame(mydata, res, pre) #重新形成数据集 View(mydata) #查看数据
2. 条件4判断(方差齐性检验)
(1) 软件操作
##方差齐性检验## library(car) #调用软件包“car” leveneTest(contain ~ drug, data = mydata, center = mean) #levene法方差齐性检验 t.test(res ~ drug, data = mydata) #T检验残差差异
(2) 结果解读
图4通过levene法对原始数据进行方差齐性检验，结果显示P=0.8622>0.1，提示相互比较的各处理水平的总体方差相等。若直接比较两处理组的残差是否有差异，两独立样本t检验结果显示P=1.000 (图5)，也提示相互比较的两处理水平的总体方差相等。本案例满足条件4。
3. 条件5判断(正态性检验)
(1) 软件操作
##正态性检验##
shapiro.test(mydata$res)
(2) 结果解读
图6正态性检验结果显示P=0.9928>0.05，提示因变量残差近似服从正态分布。本案例满足条件5。
4. 条件6判断(异常值检测)
(1) 软件操作
##绘制箱线图##
library(ggplot2) #调用包“ggplot2”
boxp <- ggplot(mydata, aes(x = drug, y = res))+ #定义颜色与分组
stat_boxplot(geom = 'errorbar')+ #添加误差线
geom_boxplot() #绘制箱线图
boxp
ggsave("boxp.png", width=8, height=5, dpi=300) #保存图片
(2) 结果解读
残差的箱线图将被保存在R数据分析路径下，图7残差的箱线图未提示任何异常值和极端值，满足条件6。
(三) 统计描述
1. 软件操作
##数据分组描述## library(psych) #调用包“psych” describeBy(contain ~ drug + phase, data = mydata) describeBy(contain ~ drug, data = mydata) describeBy(contain ~ phase, data = mydata)
2. 结果解读
由图8可知，第一阶段A、B药物组的血浆浓度分别为69.02±17.46、76.40±19.15 mmol/L，第二阶段A、B药物组的血浆浓度分别为77.20±19.02、70.20±16.77 mmol/L。
由图9可知，A、B药物组总体的血浆浓度分别为73.11±17.75、73.30±17.28 mmol/L。
由图10可知，第一阶段、第二阶段总体的血浆浓度分别为72.71±17.71、73.70±17.30 mmol/L。
(四) 三因素方差分析
1. 软件操作
##查看方差分析结果##
summary(AOV)
2. 结果解读
图11为两阶段交叉设计的方差分析结果。drug的F=0.251，P=0.630；表明尚不能认为两种药物的血浆浓度有差别。phase的F=6.813，P=0.031，表明两阶段血浆浓度差异有统计学意义。testid的F=851.367，P<0.001，表明个体间血浆浓度差异有统计学意义。注意，方差分析并未给出顺序效应(sequential)的结果，因此可以考虑使用混合效应模型进行两阶段交叉设计分析。
(五) 两阶段交叉设计分析(混合效应模型)
重新导入数据“交叉设计方差分析_顺序效应.csv”，其中变量sequential的水平1表示先用drug_A、再用drug_B，水平2表示先用drug_B、再用drug_A。其数据导入、变量类型设置和标签添加步骤不再赘述。直接演示使用混合效应模型进行两阶段交叉设计分析。
1. 软件操作
library (nlme) #加载nlme包 m0.nlme <- lme (contain ~ drug + phase, random = ~ 1 | testid, data = mydata) #分析不含sequential的结果 summary (m0.nlme) #查看m0.nlme模型结果 m1.nlme <- lme (contain ~ drug + phase + sequential, random = ~ 1 | testid, data = mydata) #分析含sequential的结果 summary (m1.nlme) #查看m1.nlme模型结果
2. 结果解读
图12为不添加sequential的混合效应分析结果。drug的t=0.501，P=0.630；表明尚不能认为两种药物的血浆浓度有差别。phase的t=2.610，P=0.031，表明两阶段血浆浓度差异有统计学意义。
图13为添加sequential的固定效应的混合效应分析结果。结果与图12一致。可见顺序效应并无统计学意义(t=0.628，P=0.548)。
四、结论
本研究采用交叉设计分析A、B两种药物对某指标的血浆浓度的影响。通过残差判断，2组数据不存在需要特殊处理的异常值；通过对残差进行正态性检验，提示数据服从正态分布；通过对残差进行统计检验，提示2组药物组数据间方差齐。
分析结果显示，A、B药物组总体的血浆浓度分别为73.11±17.75、73.30±17.28 mmol/L，第一阶段、第二阶段总体的血浆浓度分别为72.71±17.71、73.70±17.30 mmol/L。两种药物的血浆浓度差异无统计学意义(F=0.251，P=0.630)；两阶段血浆浓度差异有统计学意义(F=6.813，P=0.031)。