两因素重复测量方差分析(Two-Way Repeated-Measures ANOVA)一(无交互作用)——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：观察变量为重复测量数据，即不满足独立性。本研究中两个组别在三个时间点时测量的体重均是针对同一批样本，因此不满足独立性，该条件满足。

条件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中的第2和第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.3597、0.7700、0.9191，B组三个时间点的P值分别为0.8127、0.8091、0.1820，均>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.2533，F=1.7413、P=0.2035，F=2e-04、P=0.9876；提示每个时间点的两组之间都满足方差齐性。本案例满足条件6。关于方差齐性检验的更多内容请阅读医学统计学核心概念及重要假设检验的软件实现(4/4)——方差齐性检验及SPSS实现。

(三) 统计描述

1. 软件操作

## 宽数据转换为长数据 ##

mydata$ID<-1:20
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.60±0.89 kg和1.51±0.27、3.86±0.58、5.67±0.91 kg。图7“emmeans(估算边际均值)”中提供了数据的“emmean [估算边际均值(偏最小二乘均值))”、“SE (标准误)”、“df(自由度)”、均值的“lower.CL (95%置信区间下限)”及“upper.CL(95%置信区间上限)”。图8的估算边际均值图绘制了两组三个时间点体重的变化情况，可见两组的体重均有增加，并且增加的幅度基本保持一致。

(四) 球性检验

1. 软件操作

## 球形检验 ##

#install.packages("ez") #安装包“ez”
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.886，P=0.356，表示数据满足球形假设。因此本案例可以直接采用非校正方法分析的结果。

(五) 交互作用判断

由于本案例有两个因素(一个为时间因素time，另一个为分组因素group)，因此需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。

图9“ANOVA (方差分析)”显示，time与group之间的交互作用结果为F_time*group=0.106，P=0.900，提示time与group之间的交互作用无统计学意义，本案例可直接以主效应分析结果进行判断，如果交互作用有统计学意义，则需要分析单独效应。

(六) 时间效应分析

图9“ANOVA (方差分析)”中time的检验结果可知，F_time=437.3，P<0.001(图9)，认为体重变化具有时间变化趋势。

(七) 组间效应分析

图9“ANOVA (方差分析)”中group的检验结果可知，F_group=1.52e-04，P=0.990(图11)，表示两种饲料对家兔的增重效果差异无统计学意义。

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

1. 软件操作

## 事后检验 ##

AOV<-aov(tvalue~time,mydata_long) #检验不同时间之间的差异
TukeyHSD(AOV,p.adjust.methods="bonferroni") #事后两两比较

2. 结果解读

图10 “TukeyHSD (事后检验比较-时间)”表格中提供了各时间点两两比较的“diff (均数差)”、 “lwr (95%置信区间下限)”、“upr (95%置信区间上限)”和“p.adj (调整后的P值)”。可知，随着时间的延长，各时刻与time1时刻相比，均值差逐渐增大。time2比time1时均值增加了2.36 kg，至time3时增加了4.11 kg，time2到time3时增加了1.75 kg，差异均有统计学意义(P<0.001)。

四、结论

本研究采用两因素重复测量方差分析比较A、B两种饲料对家兔的增重效果有无差别。通过专业知识判断，数据不存在异常值；通过Shapiro-Wilk检验，提示各组数据服从正态分布；通过Levene’s检验，提示每个时间点的两组之间都满足方差齐性；球形度检验提示，满足球形假设(W=0.886，P=0.356)；组别与time无交互作用(F_time*group=0.106，P=0.900)，故进行主效应分析。

A、B两组在试验开始时、第1个月和第2个月时的体重分别为1.54±0.33、3.91±0.91、5.60±0.89 kg和1.51±0.27、3.86±0.58、5.67±0.91 kg。两因素重复测量方差分析结果显示，在试验开始时、第1个月和第2个月2组家兔的体重均呈上升趋势(F_time=437.3，P<0.001)，但两种饲料对家兔的增重效果无差异(F_group=1.52e-04，P=0.990)。进一步采用“Bonferroni”校正法进行两两比较可知，随着时间的延长，第1个月和第2个月与试验开始时相比，体重增加均有统计学意义(P<0.001)。综上可知，两种饲料对家兔的增重效果无差异。