在前面文章中研究的因素可能是多个，但是重复测量因素只有一个，本文介绍存在两个重复测量因素(Two Repeated Measures Factor)但无交互作用的情况及其在R软件中的实现过程。
关键词：R; 重复测量; 重复测量资料; 重复测量方差分析; 两重复测量因素; 球形检验; 交互作用; 主效应; 单独效应
一、案例介绍
评价在某麻醉药物的作用下，不同部位某蛋白酶活性变化情况。在使用该种药物后，检测12只小白鼠用药后10min、20min、30min、40min 4个不同时刻A部位和B部位蛋白酶的活性。部分数据见图1。本案例数据可从“附件下载”处下载。
二、问题分析
本案例资料具有两个重复测量因素，一个是时间因素，有4个水平，即time1~time4，分别代表用药后10min、20min、30min、40min 4个不同时刻。另一个是身体不同部位，有2个水平，即A部位和B部位。因此需要使用两重复测量因素方差分析进行数据分析。需要满足以下5个条件：
条件1：观察变量唯一，且为连续变量。本研究中观察变量只有蛋白酶活性，且为连续变量，该条件满足。
条件2：有两个分析因素。本研究有时间因素time和部位因素两个因素，该条件满足。
条件3：观察变量为重复测量数据，即不满足独立性。本研究中在4个时间点时测量的蛋白酶活性均是针对同一批样本，因此不满足独立性；相同个体不同部位也不满足独立性；该条件满足。
条件4：观察变量不存在显著的异常值，该条件需要通过软件分析后判断。
条件5：各组、各水平(时间点)观察变量为正态(或近似正态)分布，该条件需要通过软件分析后判断。
关于方差齐性的问题，由于本案例不存在分组因素，因此无须检验方差齐同。
三、软件操作及结果解读
(一) 导入数据
语句见下，结果见图2。
mydata = read.csv("./两因素重复测量方差分析一(无交互作用).csv", header = T)
str(mydata)
mydata$ID = factor(mydata$ID)
mydata$Site <- factor(mydata$Site, levels = c(1, 2),
labels = c("A部位", "B部位"))
(二) 适用条件判断
1. 条件4判断(异常值判断)
(1) 软件操作
boxplot(mydata$time1 ~ mydata$Site, boxwex = 0.25, col = c('lightblue'),
xlab = '', ylab = '活性', main = 'time = 10 min')
(2) 结果解读
异常值可通过箱线图和专业知识进行判断。图3为时刻1时A部位和B部位活性的箱线图，提示不存在异常值。依次可判定A、B部位其他时刻数据也不存在异常值(请读者自行操作)。
2. 条件5判断(正态性检验)
(1) 软件操作
mydata1 = subset(mydata, mydata$Site == "A部位")
mydata2 = subset(mydata, mydata$Site == "B部位")
p1 = p2 = {}
for (i in 1:4) {
p1[i] = shapiro.test(mydata1[, i + 2])[["p.value"]]
p2[i] = shapiro.test(mydata2[, i + 2])[["p.value"]]
}
p1; p2
(2) 结果解读
Shapiro-Wilk (夏皮罗-威尔克，S-W)正态性检验的结果见图4。可见大多数数据的P值均>0.05，提示各时间点下两部位的数据整体服从正态分布。关于正态性检验的注意事项详见推文(医学统计学核心概念及重要假设检验的软件实现(2/4)——正态性假设检验的SPSS实现)。
(三) 统计描述
1. 软件操作
##宽数据转换为长数据##
library(reshape2) #调用包“reshape2”
mydata_long <- melt(data = mydata, id.vars = c('ID', 'Site'), #保留不变的变量
measure.vars = c('time1', 'time2', 'time3', 'time4'), #想要转换的变量
variable.name = 'time', #转换后的分类变量名
value.name = 'tvalue') #转换后的数值变量名
## 统计描述 ##
library(plyr)
em <- ddply(mydata_long, c("time", "Site"), summarise,
mean = mean(tvalue, na.rm = TRUE),
sd = sd(tvalue, na.rm = TRUE),
n = sum(!is.na(tvalue)),
se = sd/sqrt(n),
lower.CI = t.test(tvalue)[["conf.int"]][1],
upper.CI = t.test(tvalue)[["conf.int"]][2])
em
## 绘制边际均数图 ##
library(ggplot2)
ggplot(em, aes(x = time, y = mean, color = Site, group = Site)) + #绘制分组图
geom_errorbar(aes(ymin = lower.CI, ymax = upper.CI), width = .05)+ #添加误差线
geom_line() + #添加连线
geom_point() #添加点
ggsave("emplot.png", width = 8, height = 5, dpi = 300) #保存图片
2. 结果解读
图5列出了两部位time1、time2、time3、time4时间点的“估算边际平均值”“标准误差”及均值的“95% 置信区间 (95%CI)”。估算边际均值图(图6)绘制了两部位4个时间点蛋白酶活性的变化情况，可见两组数据均有下降，但下降的幅度较为接近。
(四) 球形检验及方差分析
1. 软件操作
install.packages("ez") #安装包“ez”
library(ez) #调用包“ez”
ezANOVA(data = mydata_long, dv = tvalue, wid = ID,
within = .(time, Site), detailed = T) #重复测量方差分析
2. 结果解读
“莫奇来球形度检验(Mauchly’s Test for Sphericity)”中给出了球形假设检验结果(图7)，可见“时间”“时间*部位”的检验结果均P<0.05，提示不满足球形假设。当违背了球形假设条件时，可采用“多变量检验”结果，或者使用单变量检验校正法的结果，包括格林豪斯-盖斯勒、辛-费德特和下线3种校正方法。单变量检验校正法一般建议采用前两种方法，当epsilon (ε)<0.75时，使用格林豪斯-盖斯勒法，当epsilon (ε)>0.75时，使用辛-费德特法。当违背了球形假设条件时，如果单变量和多变量检验结果不一致，以多变量检验结果为准。
(五) 交互作用判断
由于本案例有两个因素(一个为时间因素，另一个为部位)，因此需要首先判断两个因素之间是否存在交互作用。如果交互作用有统计学意义，则需要分析单独效应。
“ANOVA (方差分析)”结果(图7)显示，time与部位(Site)之间的交互作用结果为F时间*部位=1.283，P=0.296，提示时间与部位之间交互作用无统计学意义。因此，本案例采用主效应分析结果即可。
(六) 时间效应分析
“ANOVA (方差分析)”结果(图7)显示，F时间=369.928，P<0.001，提示不同时间点蛋白酶活性差异有统计学意义。
(七) 部位效应分析
“ANOVA (方差分析)”结果(图7)显示，F部位=0.091，P=0.768，提示不同部位间蛋白酶活性差异无统计学意义。
(八) 事后检验(两两比较)
1. 软件操作
install.packages("bruceR")
library(bruceR) #调用程序包"bruceR"
library(magrittr) #调用程序包"magrittr"
library(dplyr) #调用程序包"dplyr"
MANOVA(data = mydata_long, subID = "ID",
dv = "tvalue", within = c("time")) %>%
EMMEANS("time")
2. 结果解读
由时间效应分析发可知不同时间蛋白酶活性存在差异，则需要进一步两两比较，本案例采用bonferroni法进行事后两两比较。“多重比较”结果(图8)提供了各时间点两两比较的“平均值差值及其标准误”和“显著性(校正P值)”。结果显示，随着时间的延长，time2至time4时刻与time1时刻相比，均值差逐渐增大，且差异有统计学意义。
四、结论
本研究采用两重复测量因素方差分析评价在麻醉药物的作用下，A、B机体部位某蛋白酶活性变化情况。数据不存在异常值，整体数据服从正态分布；时间与部位之间交互作用无统计学意义(F时间*部位=1.283，P>0.05)，故进行主效应分析。
时间效应分析显示，不同时间点蛋白酶活性差异有统计学意义(F时间=369.928，P<0.001)。两两比较发现，随着时间的延长，time2至time4时刻与time1时刻相比，均值差逐渐增大，且具有统计学差异。。部位效应分析显示，不同部位间蛋白酶活性差异无统计学意义(F部位=0.091，P=0.768)。
综上，该麻醉药对A部位B部位的蛋白酶活性影响差异无统计学意义，随着时间的延长，蛋白酶活性开始降低。