系列19我们已经介绍随机区组设计和析因设计的方差分析，今天接着介绍重复测量资料的方差分析及其SAS实现。
图9-4 定量数据假设检验方法选择
2. 方差分析
2.1 完全随机设计的方差分析
2.2 随机区组设计的方差分析
2.3 析因设计的方差分析
2.4 重复测量设计的方差分析
重复测量（Repeated Measure）是指对同一观察对象的同一观察指标在不同时间或不同环境下进行的多次测量，用于分析观察指标的变化趋势及有关的影响因素。重复测量资料在不同科学领域内都很常见。在医学中，为了解病人的病情变化对相关疫情指标进行连续的监测所得到资料；在毒物试验中，将某种化合物给予受孕雌鼠以观察每胎中胎鼠的畸胎数或死亡数等资料。
前面介绍的t检验、完全随机设计方差分析、随机区间设计的方差分析要求各次观察是相互独立的，而多重复测量资料同一受试对象的某观察指标进行的多次测量，那么对于同一个受试对象来说，多次测量间可能不是独立的，存在着某种相关性。此时，应使用多重复测量的统计方法以充分揭示出数据内在的特点，否则有可能会得出错误的结论。重复测量设计方差分析的前提条件包括：处理因素的各处理水平的受试对象个体间是相互独立的随机样本，总体均数服从正态分布；相互比较的各处理水平的总体方差齐；各时间点组成的协方差阵具有球形性特征。
重复测量资料最常见的情况是同一个受试对象前后有两次测量结果，如饭前饭后血糖的测量值。前后测量设计资料包括未设立平行对照和未设立平行对照的。这类资料的特点、及其与配对设计资料的区别详见表9-5。
表9-5 前后测量设计与配对设计的区别
重复测量次数大于2时，即多重复测量设计。多重复测量设计与前面讲的随机区组设计很相似，但两者间是有区别的，多重复测量设计的“处理”是在区组受试对象间随机分配，区组内的各时间点是固定的，不能随机分配；随机区组设计的每个区组内受试对象相互独立，处理在区组内随机分配，每个受试对象授受的处理是不同的。由于多重复测量设计的区组内受试对象可能存在相关性，我们不能直接使用随机区组的方差分析，如果重复测量数据满足“球对称”假设，可用随机区组方差分析方法比较处理间差异；如果不满足Mauchly “球对称”假设，仍应用随机区组设计方差分析会增大I类错误，这时需要对区组内效应的F界值进行校正，常用的校正方差有Greenhoouse-Geisser、Huynh-Feldt和Lower-bound三种方法。当样本含量较小时，Mauchly“球对称”检验效能较低，易于发生II类错误，这时应对处理组内效应F值的自由度进行校正，以得到可信的统计结论。
重复测量数据的方差分析比前面介绍的方差分析要稍复杂一点，我们将按步骤给大家介绍，完成SAS实现。
（1）单组重复测量资料的方差分析
*===导入数据，数据来源孙振球、徐勇勇主编医学统计学表12-3;PROC IMPORT DATAFILE="E:\Jindingtongji\SAS\DATA\ANOVA_R1.CSV" OUT=REPEATED DBMS=CSV REPLACE;*====第一步：数据的概貌分析;*===数据转置;PROC TRANSPOSE DATA=REPEATED OUT=REPEATED_T(RENAME=(COL1=BG)); VAR MIN0 MIN45 MIN90 MIN135; BY CODE;RUN;*===每一个受试对象血糖浓度变化趋势;PROC SGPLOT DATA=REPEATED_T; VLINE _NAME_/RESPONSE=BG GROUP=CODE; XAXIS DISCRETEORDER=DATA LABEL="TIME"; YAXIS VALUES=(4 TO 7 BY 0.2) LABEL="BLOOD GLUCOSE(mmol/L)";RUN;图9-27 各受试对象不同时间血糖浓度变化
概貌分析有利于对数据的轮廓有一个大致的描述。从图9-27中可以看出8个受试对象随着放置时间的延长血糖浓度均有一个下降趋势。其中5号受试对象的变化趋势是下降——平缓——下降的过程；8号受试对象从MIN0到MIN90的下降过程缓慢，从MIN90到MIN135的下降幅度较快。
*===受试对象平均血糖变化趋势;PROC MEANS DATA=REPEATED MEAN; VAR MIN0 MIN45 MIN90 MIN135; OUTPUT OUT=RE_STAT MEAN(MIN0 MIN45 MIN90 MIN135)=M_MIN0 M_MIN45 M_MIN90 M_MIN135;RUN;PROC TRANSPOSE DATA=RE_STAT OUT=RE_STAT_T(RENAME=(COL1=BG)); VAR M_MIN0 M_MIN45 M_MIN90 M_MIN135;RUN;PROC SGPLOT DATA=RE_STAT_T; VLINE _NAME_/RESPONSE=BG; XAXIS DISCRETEORDER=DATA LABEL="TIME"; YAXIS VALUES=(4 TO 7 BY 0.2) LABEL="BLOOD GLUCOSE(mmol/L)";RUN;*====第二步：重复测量数据的方差分析;PROC ANOVA DATA=REPEATED; MODEL MIN0 MIN45 MIN90 MIN135=/NOUNI; REPEATED TIME 4/PRINTE;RUN;图9-29 球性检验结果
图9-29的“球对称”假设检验结果，显示：本例不符合“球对称”假设P=0.0073提示在后续分析中要应用校正系数进行校正。
图9-30 多变量方差分析结果
图9-30的多变量方差分析结果显示：本例的四种多变量分析方法的P=0.0009，均小于0.05，提示不同时间的的血糖浓度存在差异。
图9-31 重复测量方差分析结果
由于“球对称”假设不满足，重复测量方差分析结果要采用G-G（Greenhouse-Geisser Epsilon）的结果，P<0.0001，即可认为不同时间点的血糖浓度有差异，与之前的多变量结果一致。
*====第三步：两两比较;PROC TTEST DATA=REPEATED ALPHA=0.0167; PAIRED MIN0*(MIN45 MIN90 MIN135);RUN;图9-32 两两比较结果
经方差分析显示，受试对象在各时间点的血糖浓度是不相同的，至少有一个时间点的总体平均值不同于其他时间点的总体平均值。为了找出差异，需要进行两两比较。平均值之间两两比较的方法很多，有研究表明：当校正系数ε>0.70时，即资料的协方差矩阵基本符合球性条件时以Tukey法为优；当ε<0.70时，即资料的协方差矩阵远离球性条件较远时以Bonferroni法为优。
本例的ε=0.5283，因此采用Bonferroni法。Bonferroni法是一种适用于多重比较的配对t检验法，根据Bonferroni不等式对t分布的临界值进行校正，即名义水准α’=α/m。假设本例事先设计时重点关注受试对象不同放置时间的血糖浓度与最初血糖浓度比较，有MIN0:MIN45，MIN0:MIN90，MIN0:MIN135共3次比较，名义水准α’=0.05/3=0.0167。两两比较结果显示：MIN0时间点与MIN45时间点的差异无统计学意义（P=0.0176），MIN0时间点与MIN90时间点的差异有统计学意义（P=0.0001），MIN0时间点与MIN135时间点的差异有统计学意义（P<0.0001），详见图9-32。
*====第四步：趋势分析;PROC ANOVA DATA=REPEATED; MODEL MIN0 MIN45 MIN90 MIN135=/NOUNI; REPEATED TIME 4(0 45 90 135) POLYNOMIAL/SUMMARY PRINTE;RUN;图9-33 趋势分析结果
当对同一受试对象的某项观察指标在多个时间点上进行重复测量时，有时想要建立一个数学模型，将反应变量表示为时间的函数，以便用更为精炼的方式表达反应变量随时间变化的趋势，也可以比较不同处理组在趋势上的差异。统计学上有很多描述反应变量与时间关系的模型，如指数模型、对数模型等。
对数模型的选择，有的出于理论考虑，有的出于经验近似。通常多项式函数配合方法是一种经验近似，且不要求数据的协方差矩阵满足球性条件。图9-33中TIME_1至TIME_3表示配合1到4阶正交多项式模型，图中结果显示一阶和二阶的P<0.0001，即血糖浓度随时间变化成一阶（线性）或二阶关系，由于一阶的F=117.12大于二阶的F=109.84，在此我们选择线性关系模型。
（2）多组重复测量资料的方差分析
多组重复测量资料是指将受试对象按处理因素的不同水平分成多个组，对这些组内的每个受试对象在不同时间点上的反应变量进行测量。分组因素称为受试者间因素，时间因素是重复测量因素，称为受试者内因素。本次以两组重复测量资料为例进行讲解。
*===导入数据，数据来源余松林重复测量资料分析方法与SAS程序的表3-2; PROC IMPORT DATAFILE="E:\Jindingtongji\SAS\DATA\ANOVA4_3.CSV" OUT=REPEATED2 DBMS=CSV REPLACE;RUN;*===第一步：概貌分析;PROC MEANS DATA=REPEATED2 MEAN; CLASS GROUP; VAR DAY0 DAY10 DAY20; OUTPUT OUT=RE_STAT2 MEAN(DAY0 DAY10 DAY20)=M_DAY0 M_DAY10 M_DAY20;RUN;DATA RE_STAT2; SET RE_STAT2; IF GROUP=" " THEN DELETE;RUN;PROC TRANSPOSE DATA=RE_STAT2 OUT=RE_STAT2_T(RENAME=(COL1=WEIGHT)); VAR M_DAY0 M_DAY10 M_DAY20; BY GROUP;RUN;PROC SGPLOT DATA=RE_STAT2_T; VLINE _NAME_/RESPONSE=WEIGHT GROUP=GROUP; XAXIS DISCRETEORDER=DATA LABEL="TIME"; YAXIS VALUES=(250 TO 400 BY 40) LABEL="WEIGHT(kg)";RUN;图9-34 三组不同时间点体重变化
图9-34的结果显示：用药后到用药10天内三组体重变化相当，用药10天后到20天三组体重变化较大，其中对照组体重变化较大。
*===第二步：方差分析;PROC ANOVA DATA=REPEATED2; CLASS GROUP; MODEL DAY0 DAY10 DAY20=GROUP/NOUNI; REPEATED TIME 3/PRINTE;RUN;图9-35 球性检验结果
图9-35的“球对称”假设检验结果，显示：本例不符合“球对称”假设P<0.0001提示在后续分析中要应用校正系数进行校正。
图9-36 不同时间点的多变量方差分析结果
图9-36的多变量方差分析结果显示：本例的四种多变量分析方法的P值均小于0.0001，提示不同时间的的体重存在差异。
图9-37 时间与处理组交互作用的多变量方差分析结果
图9-37的多变量方差分析结果显示：本例的四种多变量分析方法的P值均小于0.05，提示时间与处理组间存在在交互作用。
图9-38 重复测量方差分析结果（一）
图9-38的重复测量方差分析结果显示：处理组间的F=1.24，P=0.2997，提示处理组间的差异无统计学意义。
图9-39 重复测量结果方差分析（二）
图3-39的重复测量方差分析结果显示：本例“球对称”假设不满足，时间点的结果采用G-G结果P<0.0001，提示不同时间点的体重平均值不同；处理组与时间的交互作用P=0.0042，提示不同处理组在不同时间点的体重平均值不同。
*===第三步：趋势分析;PROC ANOVA DATA=REPEATED2; CLASS GROUP; MODEL DAY0 DAY10 DAY20=GROUP/NOUNI; REPEATED TIME 3(0 10 20) POLYNOMIAL/SUMMARY PRINTE;RUN;如果不同处理组与时间存在交互作用，则它们在时间上的变化趋势是不同的。因此本例配合正交多项式做趋势分析，进一步验证不同处理组增加的体重是否存在时间趋势上的差异。正交多项式对比的方法是将不同处理组的数据配合不同的正交多项式曲线，检验这些曲线的参数是否来自同一总体。从低阶到高阶，每次检验一个参数，是一种单自由度比较，用于分析这些曲线在哪些阶次上发生变化。本例有3个时间点，因此最高只能配合到2阶多项式。
图9-40 趋势分析结果
趋势对比分析首先是检验最高阶次的参数在各组之间是否具有统计意义。如果在某一阶上的组间差异具有统计学意义，由可以认为曲线间具有不同的趋势。否则，应继续对次高阶的参数做评价。如果在任何阶次上都不存在有统计学意义有差异，说明这两条曲线的变化趋势是一致的。图9-40的趋势分析结果显示：体重平均值在1阶（P<0.0001）和2阶（P=0.0016）上的P值均小于0.05，提示三个处理组的体重在各时间点上的平均水平是不同的。处理组效应在1阶（P=0.0052）和2阶（P=0.0205）的P值也均小于0.05，提示三个处理组的体重的时间变化趋势不一致，即三个处理的时间变化曲线彼此不平行，这一点与我们之前的方差分析结果一致。
整理不易，欢迎点亮再看哦！
参考文献：
[1] 冯国双. 白话统计[M]. 北京：电子工业出版社，2018.
[2] 夏庄坤, 徐唯, 潘红莲, 等. 深入解析SAS——数据处理、分析优化与商业应用[M]. 北京: 机械工业出版社,2014.
[3] Marfio F. Triola. ElementaryStatistics[M]. New York: Christine Stavrou, 2010.
[4] 余松林, 向惠云. 重复测量资料分析方法与SAS程序[M]. 北京: 科学出版社,2004.
[5] 高惠璇. SAS系统Base SAS软件使用手册[M]. 北京：中国统计出版社，1997.
----------------------------------------------
----------------------------------------------
精鼎特邀
联系客服
