R语言重复测量方差分析

一、概念

重复测量指的是对同一观察对象的同一观察指标在不同的时间点上进行多次测量。用于分析观察指标在不同时间上的变化规律。这类资料在医学研究中比较常见。

在实际工作当中，重复测量资料常被误认为是配对设计或随机单位组设计进行分析，不仅损失了重复测量资料所蕴含的信息，还容易的得到错误的结论。由于同一受试对象在不同的时间点的观测值之间往往彼此不独立，存在某种程度的相关，因此不满足常规统计方法所要求的的独立性假定。使得其分析方法有别于一般的统计分析方法。

在临床试验中，对受试者一次用药后，对同一个指标多个时间点测量适用；如果长时间用药，即指标测量和用药交叉进行，虽然也是获得了同一指标不同时间点的测量值，但是主要研究的是不同治疗措施之间差异，不适用重复测量方差分析。

二、使用条件

在对重复测量资料进行方差分析时，要求：

①样本是随机的；

②在处理的同一水平上观测是独立的；即每个时间点各个样本的观测值互不影响；

③每一水平的测定值都来自正态总体；即每个时间点的样本均符合正态分布；

④协方差符合对称性或者球形性；进行重复资料方差分析前需要先进行球形性检验；若满足球形性要求，则直接进行方差分析，；不满足球形性要求，则需要对与时间有关的F统计量分子、分母的自由度进行校正，以减少犯I类错误的的概率，或者直接进行多变量方差分析。

①②点通过数据资料即可判断；③需要进行正态性检验；④需要进行球形性检验。

三、重复测量资料的变异来源

对于重复测量资料的方差分析，主要考虑两个方面因素的影响。一是处理分组，可通过施加干预和随机分组来实现；另一个因素是测量时间，有研究者根据专业知识和要求确定；由此，重复测量资料的变异可分解为处理因素、时间因素、处理因素与时间因素的交互作用、受试对象间的随机误差、重复测量的随机误差五部分。

四、重复测量资料的优缺点

1、优点

①每个个体作为自身对照，克服了个体间的变异，减少了一个差异来源，分析时可更好地集中于处理效应。

②重复测量设计的每一个个体作为自身对照，研究所需的个体相对较少，研究成本更低。

2、缺点

①滞留效应：前面的处理效应可能滞留到下一次的处理；

②潜隐效应：前面的处理效应有可能激活原本以前不活跃的效应；

③学习效应：在研究过程中，研究对象对实验逐步熟悉，其反应能力有可能逐步得到提高。

五、R代码实现及结果解释

1、数据

fenzu	id	subid	time	weight
A药	1	1	0	86.4
A药	2	2	0	107
A药	3	3	0	65.8
A药	4	4	0	88.2
A药	5	5	0	77.6
A药	6	6	0	63.2
A药	7	7	0	69.8
A药	1	1	8	84.2
A药	2	2	8	102.8
A药	3	3	8	64
A药	4	4	8	87.5
A药	5	5	8	75.4
A药	6	6	8	62.4
A药	7	7	8	68
A药	1	1	16	84.2
A药	2	2	16	99.4
A药	3	3	16	63.6
A药	4	4	16	85
A药	5	5	16	76
A药	6	6	16	62.8
A药	7	7	16	65.4
A药	1	1	24	85
A药	2	2	24	90.6
A药	3	3	24	62.4
A药	4	4	24	83.8
A药	5	5	24	75
A药	6	6	24	62.2
A药	7	7	24	65.6
B药	1	8	0	66.4
B药	2	9	0	93
B药	3	10	0	78
B药	4	11	0	73
B药	5	12	0	71.4
B药	6	13	0	91.9
B药	7	14	0	68.8
B药	1	8	8	63.4
B药	2	9	8	90.4
B药	3	10	8	78.2
B药	4	11	8	74
B药	5	12	8	68.6
B药	6	13	8	89.4
B药	7	14	8	65.6
B药	1	8	16	63.8
B药	2	9	16	89.4
B药	3	10	16	74.8
B药	4	11	16	71.8
B药	5	12	16	64.8
B药	6	13	16	94.6
B药	7	14	16	64.6
B药	1	8	24	64
B药	2	9	24	88.6
B药	3	10	24	73.6
B药	4	11	24	70.4
B药	5	12	24	62.8
B药	6	13	24	97.5
B药	7	14	24	63.6

 2、R代码

library(readxl)
data<-read_xlsx("c:/users/bioknow/desktop/manova.xlsx")
#异常值检验  学生化残差图或者箱式图判断
boxplot(data$weight~data$time,boxwex=0.25,col=c("lightblue"),xlab = "time(天)",ylab="weight",sep="")

#正态性检验 shapiro.test()
shapiro.test(data$weight[data$time=="0S"])
shapiro.test(data$weight[data$time=="8S"])
shapiro.test(data$weight[data$time=="16S"])
shapiro.test(data$weight[data$time=="24S"])

#球形性检验前数据变形状
library(reshape2)
data$col<-paste(data$fenzu,data$time,sep = ".")
data.dcast<-dcast(data[,c(2,6,5)],id~col)
data.dcast<-data.dcast[,-1]
data.dcast<-as.matrix(data.dcast)
#进行多元线性回归
mlmfit<-lm(data.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)

#方差分析
data$time<-as.factor(data$time)
data$subid<-as.factor(data$subid)
data$fenzu<-as.factor(data$fenzu)
fit<-aov(weight~time*fenzu+Error(subid/(time)),data = data)
summary(fit)

3、结果解释

①正态性检验结果

> shapiro.test(data$weight[data$time=="0"])

	Shapiro-Wilk normality test

data:  data$weight[data$time == "0"]
W = 0.91782, p-value = 0.2045

> shapiro.test(data$weight[data$time=="8"])

	Shapiro-Wilk normality test

data:  data$weight[data$time == "8"]
W = 0.9192, p-value = 0.2141

> shapiro.test(data$weight[data$time=="16"])

	Shapiro-Wilk normality test

data:  data$weight[data$time == "16"]
W = 0.87601, p-value = 0.051

> shapiro.test(data$weight[data$time=="24"])

	Shapiro-Wilk normality test

data:  data$weight[data$time == "24"]
W = 0.8778, p-value = 0.05407

正态性检验结果可知，四个时间点的效应量都是满足正态性分布要求的P>0.05。

②球形性检验

> data.dcast
     A药.0 A药.16 A药.24 A药.8 B药.0 B药.16 B药.24 B药.8
[1,]  86.4   84.2   85.0  84.2  66.4   63.8   64.0  63.4
[2,] 107.0   99.4   90.6 102.8  93.0   89.4   88.6  90.4
[3,]  65.8   63.6   62.4  64.0  78.0   74.8   73.6  78.2
[4,]  88.2   85.0   83.8  87.5  73.0   71.8   70.4  74.0
[5,]  77.6   76.0   75.0  75.4  71.4   64.8   62.8  68.6
[6,]  63.2   62.8   62.2  62.4  91.9   94.6   97.5  89.4
[7,]  69.8   65.4   65.6  68.0  68.8   64.6   63.6  65.6


> #球形性检验
> mauchly.test(mlmfit,M=~group+time,data=idata)

	Mauchly's test of sphericity

data:  SSD matrix from lm(formula = data.dcast ~ 1)
W = 1.1907e-07, p-value = 8.909e-08

球形性检验需要使用款格式数据，这里展现出来了，原始数据是长格式数据。供大家参考，格式别整错了。

遗憾的本例数据球形性检验P<0.05,协方差矩阵不满足球形性。理论上我们得采用多变量方差分析或者广义估计方程来分析了。这里为了后续结果的说明，姑且认为满足球形性要求吧。

③fang

> summary(fit)

Error: subid
          Df Sum Sq Mean Sq F value Pr(>F)
fenzu      1     39    39.3   0.059  0.812
Residuals 12   7960   663.3               

Error: subid:time
           Df Sum Sq Mean Sq F value   Pr(>F)    
time        3 118.68   39.56   7.246 0.000634 ***
time:fenzu  3   5.86    1.95   0.358 0.783914    
Residuals  36 196.55    5.46                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fenzu和time的交互作用没有统计学意义F=0.358，p=0.784。只需解读两个内因素fenzu和time的主效应即可。

fenzu对weight的主效应无统计学意义，F=0.059,P=0.812。即不同fenzu对weight没影响。

time对weight的主效应有统计学意义，F=7.246，P=0.0006。即不同time对weight有影响。如要评价不同时间点之前的差异，需进行两两比较。

④fenzu和time对weight的交互影响图

 从图可以看出，两条线近乎平行，说明fenzu和time两因素间不存在交互作用。

目前，根据自己使用软件的情况，相对R语言，SPSS、SAS两种软件重复测量方差分析的输出结果比较全面。在实际应用当中可采用这两款软件分析获取更为全面的结果。

 

参考文献：

1、李晓松.医学统计学（第二版）[M].北京：高等教育出版社。

2、汪海波,罗莉,吴为,孟玲,杨世宏,汪海玲.SAS统计分析与应用从入门到精通（第二版）[M].北京:人民邮电出版社。

3、两因素重复测量方差分析，史上最详细SPSS教程！

4、R语言：两因素重复测量方差分析 - 知乎