重复测量数据分析及结果详解（之三）——多⽔平模型
上⼀篇⽂章继续介绍了基于重复测量数据的⼴义估计⽅程，对⼴义估计⽅程的结果进⾏了详细解释。本⽂承接上⽂，继续介绍医学中重复测量数据的常⽤⽅法：多⽔平模型（multilevel model）。
三、多⽔平模型
要了解多⽔平模型，得先了解多⽔平数据，或者叫多层次数据，简单来说就是具有⼀定层次结构的数据。⽐如学校-学⽣、村-村民、单位-职⼯，这些都是⼆层结构，同样，重复测量数据也是⼀个⼆层次数据，即：个⼈-时间点。每个⼈包含多个时间点，所以是个⼆层次结构数据。
多⽔平模型在不同领域有不同叫法，⽐如随机系数模型、随机效应模型、混合效应模型等等，其实说的都是差不多⼀回事。它主要是⼀种处理⾮独⽴数据的⽅法（重复测量数据就是⾮独⽴数据，因为每次测量结果可能是相关的，⽽不是独⽴的）。
多⽔平模型应⽤场景很多，这⾥只介绍它在重复测量数据中的应⽤情况。对于重复测量数据，它是将个体与时间点分别看做⽔平2和⽔平1两个层次上的数据，⽔平2包含⽔平1（即每个⼈包含多次测量）。
（⼀）基于重复测量数据的多⽔平模型思想二手骊威
对于重复测量数据，传统⽅法（如⼀般线性模型）的做法是将数据合并，估计⾃变量（如组别因素、时间）的效应，这种做法忽略了个体间的变异，因为它假定个体之间是没有差异的。⽽现实情况是，数据的差异不仅体现在个体内（不同时间点），往往个体之间也是有差异的，因此分析时需要考虑到个体间的差异。
多⽔平模型正是基于这种思想，将个体间（⽔平2）的变异分解出来，并估计其⼤⼩（通常⽤⽅差表⽰）。然后在此基础上，再估计⾃变量（如组别）的效应。
（⼆）多⽔平模型中的随机效应
多⽔平模型中，随机效应是⼀个⾮常关键的词。多⽔平模型不同于传统⽅法的特点就是可以估计随机效应。
在重复测量数据中，随机效应可以理解为个体（⽔平2）间的随机差异。常见的随机效应有两类：
⼀是个体间只有截距存在差异，斜率相同（图2左，表现为⼏条平⾏曲线或直线），即每⼀个体⾼低不同，但时间变化趋势相同，这时候如何体现个体间（⽔平2）的差异呢？只要估计出截距的变异（⽅差）就⾏了，所以这种情况下叫做随机截距模型；
⼆是个体间不仅截距存在差异，斜率也不相同（图2右，表现为⼏条不平⾏的曲线或直线），即每⼀
个体不仅⾼低不同，⽽且每个⼈随时间变化的趋势也各不相同。这种情况下，要体现出个体间的差异，仅估计截距的变异是不够，还得估计斜率的变异，即随机系数模型。
重复测量数据分析中，通常先考虑随机截距模型，如果出于专业需要，或经图形发现每个⼈的时间变化趋势差别较⼤，也可以考虑随机系数模型。
（三）多⽔平模型中的⽤途
多⽔平模型⽤途⼗分⼴泛，它可以⽤于具有层次结构的数据分析，如“学校-学⽣”“乡镇-村民”“患者-肿瘤部位”“个体-时间”等数据，重复测量数据可以看做是多⽔平数据的特例。
（四）多⽔平模型的SAS软件实现
仍以前⾯⽂章的例1作为⽰例。在多⽔平模型和⼴义估计⽅程中，是以模型的形式，主要是明确⾃变量和因变量。这⾥组别是⾃变量（试验组和对照组分别赋值为1和0），因变量⽤y表⽰。下⾯列出数据的图⽰和每个时间点的均值情况。
多⽔平模型的分析通常需要进⾏⼀定的初步探索过程，然后才可考虑最终确定的模型。
重复测量数据的多⽔平模型分析，建议第⼀步先绘制每⼀个体随时间的变化趋势图，根据图形为进⼀步的分析提供参考。例1中20名个体随时间变化的趋势图见图3，可以看出，每个⼈的Y值⾼低不同，差异较⼤，因此⾄少可以考虑随机截距模型。从时间变化的斜率来看，尽管在某些个体中起伏较⼤，但总的⼤都是呈增加趋势，⽽且可以发现，试验组增加的趋势似乎幅度更⼤（换句话说，试验组和对照组的变化趋势似乎有差异，当然，有没有的话，需要后⾯统计分析来确定，但起码可以给我们⼀些初步提⽰）。
确定，但起码可以给我们⼀些初步提⽰）。
如果每个⼈随时间变化趋势基本⼀致，拟合随机截距模型即可，否则需要考虑随机系数模型。当然图形只是给我们⼀个初步的主观判断，是否需要采⽤随机系数模型，可结合统计分析结果来确定。
根据图 3提⽰，我们先考虑拟合随机截距模型，以确定个体间（⽔平2）是否变异较⼤。如果变异（⽅差）较⼤，那就不能忽视这种变异，需要采⽤多⽔平模型将其估计出来；如果变异（⽅差）较⼩，那可以忽略，不⽤多⽔平模型也没事，直接采⽤传统的线性模型，把数据合并分析即可。
实际中通常先拟合⼀个随机截距的空模型（即模型中不包含任何变量），以观察个体间变异情况。SAS程序如下：
广州本田锋范报价data ex2；
input id group time y；
cards；……；新奥迪a7
proc mixed covtest data=ex2；
iconaclass id；
model y=/solution；
random int/subject=id type=un solution；
/
*random指定随机效应，此处指定 int （截距），拟合随机截距模型；
type 指定⽔平2变量，重复测量数据中即个体id；type指定协⽅差结构，常见为vc和un；solution指定输出随机效应估计结果*/
run；
表 6 显⽰了随机效应估计结果，其中 id 对应的估计值反映了个体间的变异⼤⼩（即⽔平2的⽅差），Z值和P值反映了⽅差与0相⽐的统计学检验结果，本例分析结果显⽰，个体间差异较⼤，且这种变异有统计学意义；residual 则显⽰了⽔平1的误差，即个体间变异外的随机误差。
这⾥需要提醒⼀点，由于此处估计的是⽅差，⽅差不可能⼩于0，因此此处的P值是单侧P值。有的统计软件给出的是双侧P值，是有问题的。
通过对空模型的分析，可以看出，个体间的确存在较⼤的变异，且有统计学意义，因此需要考虑多⽔平模型。如果个体间变异很⼩，⽆统计学意义，也可以不⽤多⽔平模型，直接采⽤普通的⼀般线性模型即可。
基于表 6 结果，进⼀步拟合随机截距模型。此时的分析与⼴义估计⽅程类似，如果要考虑组别和时间的主效应，可先不加⼊交互项。SAS程序如下：
proc mixed covtest data=ex2；
class id time（ref=first）；
model y=time group/solution；
10月新规random int/subject=id type=un solution；
run；
结果显⽰，时间和组别的主效应均有统计学意义（表 7）。与前⾯⽂章的表4的⼴义估计⽅程结构相⽐，两种⽅法的参数估计值相同，但标准误不同，从⽽导致不同的统计量和 P值。但总的来说，结论是⼀致的。⼆者的解释也相同，这⾥不再赘述。可参考前⾯⽂章中⼴义估计⽅程的解释结果，这两种⽅法都是从模型的⾓度解释的。
此时拟合的是基于随机截距模型的结果，我们进⼀步考虑是否有必要采⽤随机系数模型，即每个⼈随时间变化的趋势（每个⼈的斜率）是否有差异。SAS程序如下：
data ex2；
input id group time y；
input id group time y；
timec=time；
…… ；东风huv报价
proc mixed covtest data=ex2；
class id time（ref=first）；
model y=group time/solution；
random int timec/subject=id type=un solution；
/*random指定int和timec作为随机效应，拟合随机系数模型 */
run；
注意程序中新产⽣了 1个新变量 timec，其值完全等同于 time，即 0~4。产⽣的⽬的主要是因为 random 语句指定的随机效应通常为连续变量⽽不是分类变量（否则容易导致估计不出结果，在SAS中是如此，不同软件规定不同，不是所有软件都需要这么做），然⽽固定效应估计中我们期望将其作为分类变量，这样才能发现每⼀时间点与参照时间点的差异。因此分析时将time作为分类变量，放在
固定效应的估计中，⽽timec作为连续变量放在随机效应的估计中。当然，如果固定效应分析本⾝就将time作为连续变量，就⽆需这⼀步。所以，这取决于你把time作为连续变量纳⼊模型，还是作为分类变量纳⼊模型，作为连续变量，那就是默认时间变化是直线趋势，如果不能确定这⼀点，最好作为分类变量纳⼊。
表 8 给出了随机系数模型的随机效应估计结果，与表 6 相⽐，尽管只增加了 1个随机斜率，但结果却增加了 2个，除了时间的随机效应外，还有1个是时间随机效应与截距随机效应的协⽅差。
表8结果中，第⼀⾏反映的是截距的⽅差，提⽰个体间差异较⼤，且有统计学意义（P=0.020）；第三⾏反映的是（时间变量的）斜率的⽅差，提⽰每⼀个体随时间变化的斜率差异并不⼤，且⽆统计学意
义（P=0.257）；第⼆⾏是截距和斜率的协⽅差，这⼀结果很有实际意义，它反映了截距与斜率的相关性，表8结果中，该估计值为负数，提⽰如果个体初始时间点的值较⼤，其斜率较⼩，反之，如果某个⼈初始时间点的值较低，其斜率变化则较⼤。通俗地说，如果某⼈⼀开始Y值很低，那么后期很可能会增长（或降低）很快，⽽如果⼀开始Y 值较⾼，后期增长（或降低）速度会较慢。当然由于该值并⽆统计学意义（P=0.188），也可能并⽆这种关联。
表 8结果提⽰，time 作为随机效应并⽆统计学意义，因此实际中考虑随机截距模型即可。
进⼀步在模型中加⼊交互项，以考量两组随时间变化的趋势是否相等。SAS程序如下：
proc mixed covtest data=ex2；
class id time（ref=first）；
model y=time group time*group/solution；