在前面文章中我们介绍了两因素重复测量方差分析(Two-way Repeated-Measures ANOVA) ——不存在交互作用时在R软件中的操作步骤，本篇文章将使用实例演示在R软件中实现两因素重复测量方差分析——存在交互作用时的操作步骤。
关键词：R语言; R软件; 重复测量; 重复测量资料; 重复测量方差分析; 两因素重复测量方差分析; 球形检验; 交互作用; 主效应; 单独效应
一、案例介绍
研究A、B两种饲料对家兔的增重效果，选择20只家兔，随机分成两组，第一组用饲料A饲养，第二组用饲料B饲养，并于试验开始后第一个月(time1)、第二个月(time2)、第三个月(time3)分别测量2组家兔体重(kg)，试比较A、B两种饲料对家兔的增重效果有无差别？部分数据见图1。本文案例可从“附件下载”处下载。
二、问题分析
本案例的分析目的是比较A、B两种饲料对家兔的增重效果有无差别。由于3个时间点的数据属于重复测量数据，且有两个组别，可以使用两因素(时间因素time和分组因素group)重复测量方差分析。但需要满足以下6个条件：
条件1：观察变量唯一，且为连续变量。本研究中观察变量只有体重，且为连续变量，该条件满足。
条件2：有两个分析因素。本研究有时间因素time和分组因素group两个因素，该条件满足。
条件3：观察变量为重复测量数据，即不满足独立性。本研究中两个组别在3个时间点时测量的体重均是针对同一批样本，因此不满足独立性，该条件满足。
条件4：观察变量不存在显著的异常值，该条件需要通过软件分析后判断。
条件5：各组、各水平(时间点)观察变量为正态(或近似正态)分布，该条件需要通过软件分析后判断。
条件6：相互比较的各处理水平(组别)的总体方差齐，该条件需要通过软件分析后判断。
三、软件操作及结果解读
(一) 导入数据
mydata <- read.csv("两因素重复测量方差分析二.csv") #导入CSV数据
View(mydata) #查看数据
在数据栏目中可以查看全部数据情况，数据集中共有4个变量和20个观察数据，4个变量分别代表饲料种类(group)、试验开始后第一个月(time1)、第二个月(time2)、第三个月(time3)的家兔体重(kg)。
如果数据集较大也可使用如下命令查看数据框结构：
str(mydata) #查看数据框结构
(二) 适用条件判断
1. 条件4判断(异常值判断)
(1) 软件操作
## 分组描述数据 ##
library(psych) #调用包“psych” describeBy(mydata$time1, group = mydata$group) #time1基本信息分组描述 describeBy(mydata$time2, group = mydata$group) #time2基本信息分组描述 describeBy(mydata$time3, group = mydata$group) #time3基本信息分组描述
## 绘制箱线图 ##
mydata$group <-factor(mydata$group,labels = c("A","B")) #转换为分组变量
## 分组绘制箱线图 ##
par(mfrow = c(1, 3)) #设置画1行3个图片
boxplot(time1~group, mydata, xlab = expression("time1")) #绘制T1箱线图
boxplot(time2~group, mydata, xlab = expression("time2")) #绘制T2箱线图
boxplot(time3~group, mydata, xlab = expression("time3")) #绘制T3箱线图
(2) 结果解读
图3-1至图3-3“describeBy (描述性分析)”结果中，列出了各组观察变量的最小值和最大值，依据专业尚不能认为存在异常值；此外，图4中的第3个箱线图中存在两个离群值，但根据专业尚不能判定为异常值。综上，本案例未发现需要处理的异常值，满足条件4。
2. 条件5判断(正态性检验)
(1) 软件操作
## 正态性检验 ##
tapply(mydata$time1,mydata$group,shapiro.test) #检验time1的正态性 tapply(mydata$time2,mydata$group,shapiro.test) #检验time2的正态性 tapply(mydata$time3,mydata$group,shapiro.test) #检验time3的正态性
(2) 结果解读
正态性检验结果显示A组三个时间点的P值分别为0.360、0.770、0.919，B组三个时间点的P值分别为0.813、0.809、0.128，均>0.1，提示各组数据均服从正态分布。此外，本案例也可以绘制Q-Q图，结果也提示各组数据均服从正态分布(请读者自行操作)。综上，本案例满足条件5。关于正态性检验的注意事项详见文章医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现。
3. 条件6判断(方差齐性检验)
(1) 软件操作
library(car) #调用软件包“car” leveneTest(time1~group,data=mydata,center=mean) leveneTest(time2~group,data=mydata,center=mean) leveneTest(time3~group,data=mydata,center=mean)
(2) 结果解读
图6 Levene’s方差齐性检验显示，time1、time2、time3时，两组之间的方差齐性检验结果分别为F=1.3929、P=0.253，F=1.9389、P=0.181，F=0.952、P=0.342；提示每个时间点的两组之间都满足方差齐性。本案例满足条件6。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。
(三) 统计描述
1. 软件操作
## 宽数据转换为长数据 ##
mydata$ID<-1:20 #给数据增加ID项
library(reshape2) # 调用包“reshape2”
mydata_long<-melt(data=mydata,id.vars=c('ID','group'),#保留不变的变量
measure.vars=c('time1','time2','time3'), #想要转换的变量
variable.name='time', #转换后的分类变量名
value.name='tvalue') #转换后的数值变量名
## 统计描述 ##
lm<-lm(tvalue~group*time,data=mydata_long) #计算回归模型 library(emmeans) #调用程序包"emmeans" em<-emmeans(lm, ~ group*time, data=mydata_long) #计算边际均值并赋值给em em #查看边际均值
## 绘制边际均数图 ##
emdata<-data.frame(em)
library(ggplot2)
ggplot(emdata,aes(x = time, y=emmean, color=group,group=group))+ #绘制分组图
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. 结果解读
图3-1至图3-3列出了各组的均值和标准差，可知A、B两组time1、time2、time3时间点的均值分别为1.54±0.33、3.91±0.91、5.6±0.89 kg和1.51±0.27、2.73±0.56、3.81±0.61 kg。图7“emmeans(估算边际均值)”中提供了数据的“emmean [估算边际均值(偏最小二乘均值)]”、“SE (标准误)”、“df(自由度)”、均值的“lower.CL (95%置信区间下限)”及“upper.CL(95%置信区间上限)”。图8的估算边际均值图绘制了两组三个时间点体重的变化情况，可见两组的体重均有增加，并且增加的幅度不太一致。
(四) 球形检验
1. 软件操作
## 球形检验 ##
library(ez) #调用包“ez” ezANOVA(data=mydata_long,dv=tvalue,wid=ID,within=.(time),between=group,detailed=T) #重复测量方差分析#
2. 结果解读
由图9中的“Mauchly’s Test for Sphericity” (莫奇来球形检验)可知，W=0.952，P=0.661，表示数据满足球形假设。因此本案例可以直接采用非校正方法分析的结果。
(五) 交互作用判断
由于本案例有两个因素(一个为时间因素time，另一个为分组因素group)，因此需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。
图9中的“ANOVA-time (方差分析-时间)”显示，time与group之间的交互作用结果为Ftime*group=31.06，P<0.001，提示time与group之间的交互作用有统计学意义，本案例需要分析单独效应。
(六) 单独效应分析
1. 软件操作
## 单独效应分析 ##
install.packages("bruceR")
library(bruceR) #调用程序包"bruceR"
library(magrittr) #调用程序包"magrittr"
library(dplyr) #调用程序包"dplyr"
MANOVA(data=mydata,dvs="time1:time3", dvs.pattern = "time(.)",
between="group",within = c("time")) %>%
EMMEANS("group",by="time") %>%
EMMEANS("time",by="group")
2. 结果解读
(1) time的单独效应
由图10“Joint Tests of ‘time’ (time的单独效应)”可知，A组内不同时间点的体重增量差异有统计学意义(F=272.57，P<0.001)；B组内不同时间点的体重增量差异也有统计学意义(F=88.16，P<0.001)。图11的结果和图7结果一致，为各组的“Mean (估算边际均值(偏最小二乘均值))”、“S.E. (标准误)”及均值的“95% Confidence Interval (95%CI)”。
(2) group的单独效应
由图12“Joint Tests of ‘group’ (group的单独效应)”可知，time1时两组间的体重增量差异无统计学意义(F=0.049，P=0.827)；time2时两组间的体重增量差异有统计学意义(F=12.290，P=0.003)；time3时两组间的体重增量差异有统计学意义(F=27.493，P<0.001)。图13的结果和图7结果一致，为各组的“Mean (估算边际均值(偏最小二乘均值))”、“SE (标准误)”及均值的“95% Confidence Interval (95%CI)”。
(七) 事后检验(两两比较)
1. time的单独效应事后检验
上面分析得出了“A、B两组内不同时间点的体重增量差异有均有统计学意义”的结论。现通过“Bonferroni (邦弗罗尼)”法分别对A、B两组内三个时间点数据进行两两比较，结果如图14所示。
图14“Pairwise Comparisons of ‘time’ (time成对比较)”结果可知，A组内，相比time1，time2和time3时的体重增量逐渐增加，且差异有统计学意义(P<0.001)；B组内，相比time1，time2和time3时的体重增量逐渐增加，且差异也有统计学意义(P<0.001)。
2. group的单独效应事后检验
上面分析得出了“time2和time3时，两组之间的体重增量差异均有统计学意义”的结论。由于group只有两组，因此无须再进行事后检验。但由图15“Pairwise Comparisons of ‘group’ (group成对比较)”结果可知，time1时，A组的体重增量与B组之间差异无统计学意义(t=0.221，P=0.827)，time2时，A组的体重增量比B组高1.18 kg，差异有统计学意义(t=3.507，P=0.003)，time3时，A组体重与B组之间的差异进一步加大，A组的体重增量比B组高1.79 kg，差异有统计学意义(t=5.243，P<0.001)。
四、结论
本研究采用两因素重复测量方差分析比较A、B两种饲料对家兔的增重效果有无差别。通过专业知识判断，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过Levene’s检验，提示每个时间点两组之间都方差齐；球形度检验提示，满足球形假设条件(W=0.952，P=0.661)；组别与time之间存在交互作用(Ftime*group=31.06，P<0.001)，故进行单独效应分析。
A、B两组在试验开始后第1个月、第2个月和第3个月时的体重分别为1.54±0.33、3.91±0.91、5.6±0.89 kg和1.51±0.27、2.73±0.56、3.81±0.61 kg。时间的单独效应分析显示，A、B两组组内不同时间点的体重增量差异均有统计学意义(P<0.001)。通过“Bonferroni (邦弗罗尼)”法分别对A、B两组内三个时间点数据进行两两比较显示，A、B两组内不同时间点之间的体重增量差异均有统计学意义(P<0.001)。组别的单独效应和事后检验分析显示，time1时，A组的体重增量与B组之间差异无统计学意义(t=0.221，P=0.827)，time2时，A组的体重增量比B组高1.18 kg，差异有统计学意义(t=3.507，P=0.003)，time3时，A组体重与B组之间的差异进一步加大，A组的体重增量比B组高1.79 kg，差异有统计学意义(t=5.243，P<0.001)。综上，A饲料对家兔的增重效果更好。