当时间序列数据的类型是经过不同处理后得到多组观测值,每一个观测值为同一观察对象的同一观察指标在不同的时间点上进行多次测量所得到的的数据值,然后需要分析观察指标在不同时间上的变化规律,例如不同处理方式的结果是否有差异,或者同一处理方式下各个时间点是否有差异时,可进行重复测量方差分析。下面介绍一个实例。
如下表所示,欲分析A、B两种药的减肥效果,每组选3个人,每个人服药后第0、1、2、3个月分别记录下体重值。问题有2个:
(1)A、B两种减肥药服用后有无减肥效果。
(2)A、B两种减肥药效果有无差异。

1. 原始数据result.txt。
原始数据需要转化为长数据格式,内容如下:

2.进行重复测量方差分析。
a.导入数据并作图进行初步观察。
library(reshape2)
mydata <- read.table('result.txt', header=T, sep="\t", stringsAsFactors=F)
with(mydata, interaction.plot(time,group,result,type='b',col=c('red','blue'),pch=1:2))
处理和时间对体重值的交互影响图如下:

可以看到,服用A、B两药后体重是随时间逐步下降的。
b.分析数据是否满足随机、独立、来自正态分布总体及协方差是否符合球形性,只有满足这四个条件才能进行重复测量方差分析。
(1)分析实验设计可以判断样本是随机的,在处理的同一水平上观测是独立的(即每个时间点各个样本的观测值互不影响)。
(2)进行组内正态性检验。采用shapiro-wilk对观测时间各水平进行正态检验,P值大于0.05,则数据符合正态分布,如果小于0.05,则不符合正态分布。
shapiro.test(mydata$result[mydata$time=='0'])
shapiro.test(mydata$result[mydata$time=='1'])
shapiro.test(mydata$result[mydata$time=='2'])
shapiro.test(mydata$result[mydata$time=='3'])
结果P值均大于0.05,数据是来源于正态分布总体。如下图所示:

(3)进行球形性检验。P值大于0.05则说明协方差矩阵满足球形性。
mydata$col <- paste(mydata$group, mydata$time, sep='.')
mydata.dcast <- dcast(mydata[,c(3,6,5)], group_no~col)
mydata.dcast <- mydata.dcast[,-1]
mydata.dcast<-as.matrix(mydata.dcast)
mlmfit=lm(mydata.dcast~1)
group=factor(rep(c("A","B"), c(4, 4)))
time=ordered(rep(1:4,2))
idata=data.frame(group,time)
mauchly.test(mlmfit,M=~group+time,data=idata)
结果P值大于0.05,说明协方差矩阵满足球形性。

c.进行重复测量方差分析。
mydata$time <- as.factor(mydata$time)
mydata$no <- as.factor(mydata$no)
mydata$group <- as.factor(mydata$group)
fit<-aov(result~time*group+Error(no/(time)),data = mydata)
summary(fit)
结果如下:

可以看到
(1)group对result的P=0.0796,说明group效应无统计学意义,即A、B两种减肥药服用后的减肥效果差异无统计学意义。
(2)time对result的P=1.53e-10,说明time效应有统计学意义,即A、B两种减肥药服用后随时间变化体重变化差异有统计学意义,亦即A、B两种减肥药有减肥效果。如要评价不同时间点之前的差异,需进行两两比较。
(3)time和group的交互作用有统计学意义,p=0.00042。