方差分析二重复测量资料方差分析沉毅

1、方差分析（二）：重复测量资料的方差分析沈 毅2005.31 重复测量（repeated measure）是指对同一观察对象的同一观察指标在不同时间点上进行的多次测量，用于分析该观察指标在不同时间上的变化特点。这类测量资料在临床和流行病学研究中比较常见。例如，为研究某种药物对高血压病人的治疗效果，需要定时多次测定受试者的血压，以分析其血压的变动情况。一些传统的统计方法，如t检验、方差分析、线性回归模型等，都要求各次观察是相互独立的。而重复测量资料由于是对同一受试者的某项观察指标进行的多次测量，在同一受试者的多次测量之间可能存在某种相关性，用通常的统计方法就不能充分揭示出其内在的特点，有时甚至会得

2、出错误的结论。因此，有关重复测量资料的分析方法是近代统计学研究的热点之一。2重复测量数据结构与独立数据结构的区别及其优缺点 为了解重复测量数据与独立观察数据之间的区别，下面用一个完全随机设计的独立数 据结构与具有一个受试者内因素(时间)的重复测量数据结构进行比较。 完全随机设计的独立数据结构：从正常人、可疑硅沉着病者及一期硅沉着病病人中各随机抽取5人，测量他们的血清黏蛋白含量(mgL)，结果列于下表中第一部分。 重复测量的数据结构：对5名粉尘作业工人的血清黏蛋白含量(mg几)连续3年的测量 结果列于下表中第二部分。3 从上表的例子可以看出，独立数据结构的各个观察值是彼此独立的，它适宜于用通常的

3、方差分析方法做统计分析。重复测量数据结构是对每一受试者的同一观察指标(血清黏蛋白含量)进行的多次测量。由于这种多次测量之间可能存在相关性，就需要用特殊的统计方法进行分析。 重复测量设计的主要优点是可以减少样本含量，其次是能够控制个体变异，即个体差异。例如在单因素实验中，可以用随机区组(或称配伍组)设计方法来缩小随机误差。而重复测量设计是以同一个受试者作为一个区组，故可以把它看成为是随机区组设计的一种极端形式。但在随机区组设计下的每一测量都是在不同受试者身上进行的，它们对某种处理因素的反应是独立的，符合独立性的假定。而在重复测量设计下的测量是在同一受试者身上进行的，它们对同一处理因素在不同时间上

4、的反应可能是不独立的，后一次的测量结果可能受到前一次测量结果的影响。因此，对同一个体在不同时间上的测量值之间就可能存在相关关系。这给分析工作带来了一定的复杂性。 在实际工作中，重复测量资料比独立观察资料往往更为多见。如在临床研究中，需要观察病人在不同时间的某些生理、生化或病理指标的变化趋势，不同时间或疗程的治疗效果。在流行病学研究中观察队列人群在不同时间上的发病情况。在卫生学研究中，纵向观察儿童的生长发育规律等。重复测量资料在自然科学和社会科学的很多领域内都可见到。4 第一节 重复测量资料方差分析对协方差阵的要求 在对重复测量资料进行方差分析时，除要求样本是随机的、在处理的同一个水平上的观察是

5、独立的以及每一水平上的测定值都来自正态总体外，特别强调协方差阵（covariance matrix）的球形性（sphericity）或为符合对称性（compound symmetry）。Box（1954）指出，若球形性质得不到满足，则方差分析的F值是有偏的，这会造成过多的拒绝本来是真的无效假设(增加型错误)。 5一、协方差阵的球形性检验 方差是指在某一时点上测定值变异性的大小，而协方差是指在两个不同时点上测定值相互变异性的大小。如果在某个时点上的取值不影响其他时点上的取值，则协方差为0，反之，则不为0。由方差协方差构成的矩阵称协方差阵。 6 设 k、l为两个测定时点， 代表协方差阵中的元素。当

6、kl时为方差，kl时为协方差。共有a个测定时点，将这a个方差和（a-1）/2个协方差排成协方差阵V为：7 协方差阵的球形性质是指该矩阵主对角线元素（方差）相等、非主对角线元素（协方差）为零。用Mauchly氏法检验协方差阵的球形性质。Mauchly氏检验的 P值若大于研究者所选择的显著性水准时，说明协方差阵的球形性质得到满足。否则，必须对与时间有关的F统计量的分子、分母自由度进行调整，以便减少犯类错误的概率。调整系数为（epsilon）。 8二、自由度调整方法 1. 调整系数的计算 有两种调整系数。 （1）Greenhouse-Geisser调整系数 为： 式（10-2）中的 是矩阵（10-1

7、）中第k行第l列元素， 是所有元素的总平均值， 是主对角线元素的平均值， 是第k行的平均值。 的取值在1.0与1/（a-1）之间。 9 （2）HuynhFeldt调整系数 研究表明，当真值在0.7以上时，用 进行自由度凋整后的统计学结论偏于保守，故Huynh和Feldt提出用平均调整值 值进行调整。 值的计算公式为 式（10-3）中的g是对受试对象的某种特征（如性别或年龄）进行分组的组数，n是每组的观察例数。 10 2调整规则 只对具有重复测定性质的时间效应的F值的自由度，和处理时间交互作用的F值的自由度进行调整。由于F值有两个自由度v1和v2，调整的分子自由度 ，分母自由度 具体计算时可用

8、或 代替。用调整所得的 及 的值查 F界值表，得 。由于1，所以调整后的F临界值要大于调整前的F临界值。 11第二节 单因素重复测量资料的方差分析 一、单因素重复测量资料的例子12 例10-1 一项关于不同药物治疗心律失常效果的对比研究。对9例经常出现心室早搏的病人于用药前测定其心率后进行随机化给药。一部分病人按A药安慰剂（C药）B药的顺序给药，另一部分病人按B药安慰剂（C药）A药的顺序给药。安慰剂（C药）持续一周作为药物后效的清除期。比较用药前与各种药物及A药与B药之间的心律差别。表10-1列出9名受试病人在用药前、安慰剂（C药）期及用药（A与B）期的心率。 上表的下部及右侧所列的数据均从表

9、体部分的原始测量值派生而来，是为方差分析而准备的。本例药物水平数a4，每组观察例数n9，观察值总个数Nan36。 13二、方差分析的步骤 1提出检验假设 检验假设为： H0：1234； H1：ih，至少有一个不等式成立。14 2计算离均差平方和、自由度及均方 有总离均差平方和、处理因素离均差平方和、受试对象间离均差平方和及受试对象内离均差平方和等。计算公式为： （1）总离均差平方和SS总及总自由度v总的计算 本例，SS总201647-（2665）2/（49）4362.97，v总36-135。 （2）处理因素的离均差平方和SS处理及自由度v处理的计算15 （3）受试对象间离均差平方和SS对象间及

10、自由度v对象间的计算 （4）受试对象内离均差平方和SS对象内及自由度v对象内的计算16 （5）误差的离均差平方和SS误差与自由度v误差的计算 SS误差SS总-SS处理-SS对象间，v误差（n-1）（a-1） （10-8） 本例，SS误差4362.97-1185.42-2023.721153.83，v误差（9-1）（4-1）24 根据以上4种离均差平方和与自由度计算所得的均方见表 10-2。 3. 计算F值 由于主要是处理因素的统计学检验，故只计算处理因素的F值。 F处理MS处理/MS误差，F处理服从v1v处理与v2v误差的F分布 本例，F处理395.14/48.088.22；v13，v224。

11、所有计算结果均列于表10-2中。 17 SAS程序给出本例的协方差阵Mauchly球形性检验的结果为P0.1628，故不必进行自由度的调整。查F界值表得：F0.01（3，24）3.01，F0.01（3，24）4.72。本例处理因素的F值为8.22，大于F0.01（3，24），故拒绝无效假设，说明处理因素间的差别具有统计学意义。 由表10-2看出，受试对象内离均差平方和等于处理因素的离均差平方和与误差的离均差平方和两项之和。 SAS程序给出 。用 调整的处理因素的分子自由度为 ；分母自由度为 。查F界值表得调整自由度后的F临界值F0.05（2，19）3.52，比未调整的F临界值大。未调整的概率

12、P0.0006， 调整概率P0.0020。 18三、平均值之间的多重比较 以上用单因素重复测量方差分析方法对心率资料进行分析之后所得到的统计学结论是：拒绝无效假设，即在治疗药物的四个水平中，至少有一个水平的总体平均值不同于其他水平的总体平均值。为了确定这个特殊总体，必须进行平均值之间的多重比较。但此处不能采用上一讲中介绍的多重比较方法，因为那些方法都是建立在独立样本基础之上的。这里可采用配对样本的差值t检验，因为配对样本就是重复测量试验中一种最简单的对比研究设计。其检验步骤如下： 19 1. 计算每一个病人在不同给药情况下的差值：di（j-h）Yij-Yih，i为病人号，j，h为药物水平号，各

13、di（j-h）列于表10-3的右边各列中。由于设计时只考虑用药前与各种药物及A药与B药之间的差纵情况，表中只列出了di（1-2）、di（1-3）、di（1-4）及di（2-4）四种组合，而不是所有可能的六种组合。 2. 进行差值t检验（参阅成对测量值之差的比较），可分别得到t值为：t（1：2）4.41，t（1：3）0.03，t（1：4）3.19，t（2：4）-0.96 3计算校正临界值t 由于是对同一份资料进行多重比较，为克服累积类错误对结果判断所造成的影响，根据Bonferroni不等式的原理对临界t值进行调整。 20 首先确定比较的次数 c。因该研究已事先确定只作 4次比较，故c4。若在方

14、差分析之后再作多重比较，则只能取所有可能的比较次数。例如本例在方差分析之后再进行比较时，则比较的次数应为 c4（41）/26。 其次是选择累积 类错误的概率0.10。采用双侧检验，每次检验所用的I类错误概率水准为0.10/40.0125，自由度vn-19-18，查附表2：t界值表得t0.0125/2（8）3.21。与前面计算出的t值相比较，可见用药前心率与服用A药后心率之差具有统计学意义。用药后心率平均降低12.44次/分，而用药前心率与服安慰剂后心率之间以及A药与B药之间心率之差无统计学意义。用药前心率与用B药后心率之差接近显著性水平。 21第三节 两因素重复测定资料的方差分析 两因素重复测

15、定资料中的两因素是指一个组间因素（处理因素）和一个组内因素（时间因素）。组间因素是指分组或分类变量，它把所有受试对象按分类变量的水平分为几个组。组内因素是指重复测定的时间变量。例10-1只有组内因素，没有组间因素。 22一、例子 例10-2 一项药物代谢动力学研究，目的是对比某种药物的不同剂型在体内的代谢速度。剂型分胶囊型和片剂型。将16名受试对象随机分为两组，每组8名。一组给予胶囊，另一组给予片剂，分别在服药后1、2、4、6及8小时测定血中的药物浓度。测定结果见表10-4。 S（9.73）2+（54.61）2+（30.30）2+（70.20）2238012.38 本例的组间因素是药物剂型，组

16、内因素是测定时间。各下标的意义是：i（i1，2，g）为组间因素的分组号，j（j1，2，p）为测定时间点的序号，k（k1，2，ni）为组间因素第i水平的受试对象号，受试对象总数为n1n2ng。当各ni相等时，则用 n代替ni。测量值总个数Ngnp本例g2；各组受试对象数n8，p5，受试对象总数为2816例，测量值总个数N80。23 二、方差分析模型 一个组间因素，一个组内因素的方差分析模型为： Yijk+i+i+（）ij+（i）k+ijk 模型中各参数的意义是：为总体平均值；i为处理组i的效应；i为第j个测定时间点的效应；（）ij为第i组在第j个测定时点上的效应，属交互作用，为固定效应；（i）k

17、为第i组第k个观察对象的效应，属随机效应；ijk为误差项。给定限制条件为： 模型中的参数估计值与平均值之间的关系见 表 10-5。24三、方差分析的步骤 1.离均差平方和、自由度及均方的计算 令 为观察值总和， 为观察值平方总和， 为第i组观察值之和， 第j时点观察值之和， 为第i组第k个受试对象的观察值之和， 为在（ij）水平上的观察值之和。观察值总个数为Ngpn。表10-4中同时列出了Tij，Tik，Ti，Tj，T和S的值。将各离均差平方和、自由度及均方的计算公式列于表 10-6中。 用表10-6中所列的计算公式计算出各离均差平方和、自由度及均方见表10-7。 2计算F值本例的各F值见表1

18、0-7。25 3确定P值并作出统计推断 在查阅F临界值表前，要对组内变量（时间）及组间 组内（剂型时间）交互作用的自由度作调整。调整时要应用 的估计值。G-G调整系数 ； H-F调整系数 。如用 计算的F组内自由度调整值为： F组间*组内的分子自由度的调整值为： 对组间（剂型）效应仍用原来的自由度 查F临界值表得：F0.05（1，14）4.60。对组内（时间）效应及组间组内交互效应的F临界值要用调整自由度得F0.05（2，29）3.33。 统计学推断为：不同剂型药物的血液浓度间的差别无统计学意义（ P0.0645），不同时间的血药浓度间的差别具有统计学意义（ P0.0001），剂型与时间之间无

19、明显的交互作用（ P0.3413）。26第四节 趋势分析 趋势分析（trend analysis）用于分析某项观察指标在整个时间区间内的变化趋势，或者对比不同处理组在趋势上的差异。例如分析不同类型糖尿病病人在一天中的血糖变化规律，分析不同性别儿童的生长发育趋势等。本节所介绍的趋势分析方法是把应变量作为时间的函数来配合多项式曲线，不要求资料的协方差阵满足球形性条件。 27 一、多项式曲线模型 一般的多项式回归模型可表示为： 本节所介绍的多项式曲线是采用正交多项式系数来配合模型。其模型结构为： 模型中的 是在（r）阶下对第j时点平均值 的预报值，Cj(1)，Cj(r)是从第（1）阶起直到（r）阶的

20、第j个正交多项式系数。0，1，r为模型参数，从实际资料中求得它们的估计值，b0，b1，br。设测定的时间点有p个，从理论上讲，幂的阶次可以达到rp。在这种情况下，理论模型与观察资料完全一致。但在实际工作中，一般取幂的阶次 r（p-1）。当 r1时为线形模型。利用多项式曲线作趋势分析是将两个组（如治疗组与对照组或不同的治疗组）之间的比较转换为两个组参数之间的比较。如果任何一对参数之间的差别具有统计学意义，就可认为这两条曲线的趋势不一致。 28二、例 子 例10-3 一项药物代谢动力学研究，目的是比较某种药物的两种剂型（胶囊型和片剂型）在体内的代谢过程。每组8名受试对象，给药后等时距地每隔1小时测

21、定一次血药浓度（g/ml，用曲线下面积表示），共测定5次。总测定值的个数N40。资料见表10-8。 为了解用正交多项式作趋势分析的过程，先对胶囊型组进行单组资料的趋势分析，然后再把片剂型资料包括进来，作两组趋势的对比分析。 29三、单组资料的趋势分析 用表10-8中胶囊型资料作单组资料趋势分析。计算所需的数据经整理后列于表10-9。由于在模型配合过程中需要反复使用误差均方进行假设检验，故首先进行单向方差分析，得到误差均方的估计值为：MS误差899.98，v误差35。 趋势分析的过程可概括为六个步骤，各步骤的计算结果见表 10-11。 1.从正交多项式系数表10-10中选择一组系数C（r）（c1

22、（r） c2（r）cp（r）。系数C（r）的下标（r）是多项式曲线的阶数。本例 p5。从表中水平数为 5的四种类型的多项式中选择一种，得到一组多项式系 数C（r）。配合过程通常先由低阶向高阶逐步过渡，故首先选择一组线性多项式系数为C（1）（-2 -1 0 1 2）。第二步配合二次曲线时选择一组系数为C（2）（2 -1 -2 -1 2），余类推。 30 2计算r次幂下的（r）值： 3计算离均差平方和SS（r） 式（10-11）的分子中的n为观察例数，分母为系数的平方和（见表 10-10 的最后一列）。 4计算F值，确定概率P并作出统计推断 F值的计算公式为 统计量F服从v11及v2=v误差的F分

23、布 5计算该次幂下的回归系数r的估计值br 31 6确定回归方程并求出预报值回归方程为； 本例资料一直配合到四次幂的多项式曲线时仍然有统计学意义，说明该回归曲线的弯曲度较大。最后得到的回归方程如下： 用式（10-15）对本例资料的预报结果为： 同理得到 利用测量的平均值和预报的平均值制作的趋势图（ 图 10-1）中，可见两条曲线几乎完全重叠。 经配合多项式曲线后所得的总的回归平方和 SS回归为各次回归平方和的合计，即： SS回归113028.67+14366.44+10841.03+4005.44142241.53 即四次多项式曲线可以解释总变异（ SS总173734.51）的约 82。32四、两组重复测定资料的趋势对比 对例10-3的实