4.11 重复测量设计一元定量资料方差分析

4.11.1 问题与数据

【例4-12】某课题组研究益髓生血颗粒治疗_SEA/ααCS基因型患者血红蛋白H（Hemoglobin H，HbH）病的临床疗效，选取13名_SEA/ααCS基因型HbH病患者，给以益髓生血颗粒治疗，疗程为3个月。分别记录患者治疗前、服药1个月、服药2个月及服药3个月时血红蛋白的含量，结果如表4-14所示。请进行合适的统计分析。

【例4-13】某研究者欲观察隔日禁食结合耐饥食饵对小鼠生理指标的影响，将20只成年ICR小鼠随机地等分为对照组（自由摄食标准饲料）、食饵组（自由摄食耐饥食饵）、隔日禁食组（隔日禁食标准饲料）、辟谷食饵组（隔日禁食耐饥食饵），实验8周，自由饮水，观察4组小鼠体重的变化，如表4-15所示。请进行适当的统计分析。

【例4-14】某研究者用贲门癌患者的标本制成液体，在3种不同处理条件下观察鸡胚背根神经节与鸡胚交感神经节中长出突起的神经节比例。现有贲门癌患者10例，将每人的标本均分成3份，分别给予3种不同的处理（因素A），即A1（加入100ng/ml神经生长因子）、A2（加入200ng/ml神经生长因子）和A3（单用贲门癌培养液）；并对每种处理后的标本中的2种类型的神经节（因素B），即B1（背根神经节）与B2（交感神经节），观测长出突起的神经节的比例（Y）。实验结果如表4-16所示。请进行适当的统计分析。

4.11.2 对数据结构的分析

在例4-12中，对每一个患者在4个时间点上分别测量其血红蛋白含量，说明“时间”因素是一个重复测量因素，并且所有患者均接受同一种治疗方法—服用益髓生血颗粒，因而这是具有一个重复测量的单因素设计定量资料。

在例4-13中，对每一只小鼠在4个时间点上分别测量其体重，说明“时间”因素是一个重复测量的因素。此外，还有2个实验因素“摄食方式”（自由摄食或隔日禁食）和“饲料种类”（标准饲料或耐饥食饵），因而这是具有一个重复测量的三因素设计定量资料。

在例4-14中，涉及2个实验因素，即“处理”和“神经节类型”，前者有3个水平，后者有2个水平。由于是对同一个患者的标本采用不同的处理，同时在不同类型神经节中重复获得指标“长出突起的神经节的比例”的6个观测值，因而2个实验因素都是重复测量因素。该资料所对应的实验设计类型为具有2个重复测量的两因素设计。

4.11.3 分析目的与统计分析方法的选择

分析目的都是比较不同实验条件下定量观测指标的平均值之间的差别是否具有统计学意义。

例4-12，应选用具有一个重复测量的单因素设计一元定量资料方差分析；例4-13，应选用具有一个重复测量的三因素设计一元定量资料方差分析；例4-14，应选用具有2个重复测量的两因素设计一元定量资料方差分析。

4.11.4 SAS程序

对例4-12进行具有一个重复测量的单因素设计一元定量资料方差分析，SAS程序名为SASTJFX4_12.SAS。

SAS程序中第1步为建立数据集，person代表“患者编号”，month代表“观测时间”，y代表观测指标“血红蛋白含量”。第2、3、4、5、6步分别调用MIXED过程，采用VC、CS、UN、AR(1)、SP(POW)5种协方差结构模型对资料进行方差分析。其中，SP(POW)模型在使用时要求将“repeated/type=sp(pow) (c-list)”语句中的“c-list”替换成重复测量因素（可以是多个）的名称，但这个或这些重复测量因素应为定量变量且其水平数不能太少。在本例中，重复测量因素为观测时间，其包含4个水平（0、1、2、3），可近似视为定量变量，故可选用SP(POW)模型。若某重复测量因素为定性变量（如部位等）或虽可视为定量变量但水平数较少（如小于等于3）时，是不适宜选用此模型的。此外，定量的重复测量因素在赋值时需保证赋值数值的间隔与实际观测时间间隔成比例，较好且最简单的方法是以其真实水平代入。如本例，可以0、1、2、3代表治疗前、服药1个月、服药2个月、服药3个月，而不宜随意赋值，如赋为1、3、5、9则不妥。

对例4-13进行具有一个重复测量的三因素设计一元定量资料方差分析，SAS程序名为SASTJFX4_13.SAS。

SAS程序第1步为建立数据集，style代表“摄食方式”，forage代表“饲料种类”，time代表“时间”，subject代表“小鼠个体”，y代表观测指标“体重”。第2、3、4、5、6步分别调用MIXED过程，采用VC、CS、UN、AR(1)、SP(POW)5种协方差结构模型对资料进行方差分析，其中model语句后的“style|forage|time”相当于“style forage time style*forage style*time forage*time style*forage*time”，即这三个因素的主效应及各级交互作用的效应。第7步为建立宏shuju，以实现对数据集中已有变量value的更名。第8、9步均用来实现对不同数据集的横向合并。第10、11步均用来将数据集中的内容输出到output窗口中。第2、3、4、5、6步分别调用MIXED过程，几个过程所用语句基本相同，仅在“type=”后的选项不同，5个过程分别指定了5种协方差结构模型。MIXED过程中的repeated语句用来规定个体的重复测量的协方差结构，“/”后的sub（也可写为subject）用来指定数据集中的个体，若不含有分组因素，直接在“sub=”后面给出受试对象个体变量名称即可，如程序SASTJFX4_12；若含有分组因素，则在“sub=”后面给出受试对象个体变量名称的同时，还需在后面加注“（）”，在括号内填入分组变量名称，如程序SASTJFX4_13。在调用MIXED过程进行方差分析时，使用了两个ods语句，分别用来将模型拟合的有关信息（fitstatistics）和模型维度有关参数（dimensions）输出。

对例4-14进行具有两个重复测量的两因素设计一元定量资料方差分析，SAS程序名为SASTJFX4_14.SAS。

SAS程序第1步为建立数据集，patient代表“病例号”，treat代表“处理”，type代表“神经节类型”，y代表观测指标“长出凸起的神经节的比例”。第2、3、4、5步分别调用MIXED过程，采用VC、CS、UN、AR(1)4种协方差结构模型对资料进行方差分析。第6、7步将不同数据集的内容进行合并查询，这与程序SASTJFX4_13中的宏shuju功能相同。

4.11.5 主要分析结果及解释

下面给出例4-12的分析结果，即程序SASTJFX4_12.SAS的部分输出结果。

以下是采用VC（方差分量型）协方差结构模型进行混合效应模型分析的输出结果。

以下是采用CS（复合对称型）协方差结构模型进行混合效应模型分析的输出结果。

以下是采用UN（无结构型）协方差结构模型进行混合效应模型分析的输出结果。

以下是采用AR(1)（一阶自回归型）协方差结构模型进行混合效应模型分析的输出结果。

这是采用SP（POW）（空间幂相关型）协方差结构模型进行混合效应模型分析的输出结果。

以上是采用MIXED过程进行混合效应模型分析的结果。这里给出了模型的拟合信息和固定效应假设检验的结果。通常关注的是固定效应假设检验的结果。但从上面的结果中可以看出，采用不同的协方差结构模型得出的固定效应假设检验的结果是不完全相同的，其中GLM过程计算出来的结果与CS协方差结构模型（采用REML即约束最大似然估计法）得出的结果是相同的，它是最简单的模型。在采用这5种协方差结构模型计算所得的结果中，应以哪个结果为准呢？这就是模型选择的问题了。

通常情况下，可以Akaike的信息准则AIC值或Schwarz的信息准则BIC（或称SBC）值来选择协方差结构模型。AIC值和BIC值越小，协方差结构模型拟合越好；若两个协方差结构模型拟合的AIC值和BIC值接近，还可参考−2 LogL（−2 Res Log Likelihood）的数值，小者为优。若两个协方差结构模型分别包含q+v和q个参数时，可用这两个模型的−2倍的对数似然函数值构造出似然比统计量，采用χ2检验进行推断，见式（4-1）。若似然比统计量对应的P值大于设定的临界值，则两个协方差结构模型对资料的拟合效果之间的差异无统计学意义，此时可选择参数个数较少的那个协方差结构模型。若似然比统计量对应的P值小于设定的临界值，则两个协方差结构模型对资料的拟合效果之间的差异有统计学意义，此时可选择参数个数较多的那个协方差结构模型。式中，χ2v服从自由度为v的χ2分布，−2 logLq和−2 logLq+v分别为含q和q+v个协方差参数模型的−2倍的对数似然函数值。

现将本资料输出结果的模型拟合信息部分的主要内容汇总在表4-17中，以利于比较。

由上面的结果可知：采用UN（无结构型）协方差结构模型拟合本资料，其迭代过程并没有收敛，拟合效果不好。对比其他4种协方差结构模型的拟合情况，可以看出采用CS协方差结构模型拟合本资料效果最好，因为在评价拟合效果的三个准则中，CS模型对应的值均为最小且其−2 Log L值较参数个数比其少的模型VC的值差异也比较大，所以此时应以CS（复合对称型）协方差结构模型进行混合模型分析的输出结果为准，其中F=1267.75，P<0.0001。说明不同时间点上测得的HbH患者血红蛋白含量均值之间的差异存在统计学意义。欲得出在各时间点上患者血红蛋白含量均值之间两两比较的结果，可将原程序中所有过程步删除，在data步后换上以下过程步：

        ods html;
        proc mixed;
            class person month;
            model y=month/ddfm=sat;
            repeated/type=CS sub=person;
            lsmeans month/pdiff;
        run;
        ods html close;

其中，“ddfm=sat”说明采用Satterthwaite近似的方法来计算分母的自由度，从而得到一个更加准确的近似程度更高的F值。但是，此法计算的结果有时会与GLM过程计算的结果有所不同。

程序运行后，得结果如下：

这是各时间点上患者血红蛋白含量均值之间两两比较的结果。首先，给出了各时间点上患者血红蛋白含量均值与0比较的检验结果，没有太大实际意义。然后，给出了两两比较的结果，month和_month列中的0、1、2、3分别代表时间的4个水平，即治疗前、服药1个月、服药2个月、服药3个月。读者可查看最后两列给出的两两比较的结果，可发现这4个时间点上患者血红蛋白含量均值之间的差异均有统计学意义，说明4个时间点上患者血红蛋白含量各不相同，根据上面给出的均值大小可判断，各时间点患者血红蛋白含量的大小关系为：服药1个月>服药3个月>服药2个月>治疗前。

以下是对例4-13的分析结果，即程序SASTJFX4_13.SAS的输出结果。

这是上述程序中ODS（Output Delivery System）输出的结果。首先，给出了5种协方差结构模型拟合本资料的有关情况，然后给出了协方差结构的有关信息（Covariance Parameters表示模型中待估计的协方差结构中的参数个数）。比较5种模型拟合资料情况的AIC、BIC数值，可发现UN与AR(1)两种协方差结构模型拟合资料较好。现比较这两种协方差结构模型拟合资料的效果之间的差异是否有统计学意义。

χv2=−2 logLq−(−2 logLq+v)=235.5−219.1=16.4

由ODS输出结果的第2部分可知，q=2，q+v=10，所以v=8。因χ0.05(8)2=15.51<16.4，故P<0.05。因为可认为不适合用AR(1)模型取代UN模型，所以最后的结论应按UN协方差结构模型计算出来的结果来下。其假设检验结果为：

由上述结果可知：摄食方式（style）、观测时间（time）、摄食方式与观测时间的交互作用（style*time）、饲料种类与观测时间的交互作用（forage*time）均有统计学意义。

欲得出各均数之间两两比较的结果，可将原程序中所有过程步删除，换上以下过程步：

        ods html;
        proc mixed data=sastjfx4_13;
          class style forage subject time;
          model y=style|forage|time/ddfm=sat;
          repeated/type=UN sub=subject(style forage);
          lsmeans style*forage*time/pdiff;
        run;
        ods html close;

此步运行结果大部分与上述TYPE=UN计算结果相同，仅多了多个均数两两比较的结果，因输出结果所占篇幅较多，故此处从略。

其实，对于含分组因素的重复测量设计定量资料，如果观测了各受试对象某指标的“基础”状态（即实验前状态），对资料进行分析时，最合适的做法是以“基础值”为定量的影响因素，即协变量，采用相应设计类型定量资料的协方差分析。例如本例，若第一个时间点观测的数据不是实验后第1周小鼠的体重，而是实验前小鼠的体重，则本资料最合适的分析方法是以实验前小鼠的体重为基础值（协变量），采用具有一个重复测量的三因素设计定量资料一元协方差分析。有关做法可参阅4.12节。

以下是对例4-14的分析结果，即程序SASTJFX4_14.SAS的输出结果。

这是上述程序中ODS输出的结果。首先给出了4种协方差结构模型拟合本资料的有关情况，然后给出了协方差结构的有关信息。比较4种模型拟合资料情况的AIC、BIC数值（越小越好），可发现VC协方差结构模型拟合资料较好且协方差结构中参数的个数最少。因此，最后的结论应按VC协方差结构模型计算出来的结果来下。其假设检验结果为：

由上述结果可知：处理（treat）和神经节类型（type）均有统计学意义。

欲得出各均数之间两两比较的结果，可将原程序中所有过程步删除，换上以下过程步：

        ods html;
        proc mixed data=sastjfx4_16;
          class patient treat type;
          model y=treat|type/ddfm=sat;
          repeated/type=VC sub=patient;
          lsmeans treat*type/pdiff;
        run;
        ods html close;

此步运行结果大部分与上述type=VC计算结果基本相同，仅多了多个均数两两比较的结果。因输出结果所占篇幅较多，故此处从略。