|
重复测量数据的分析方法已经有多篇笔记进行过演示： 多个分组因素的重复测量资料的方差分析【SPSS】 因此本次笔记不再做理论上的介绍，直接演示实例，尽可能详细地对过程和结果进行解读。示例操作参照[多层统计分析模型：方法与应用/王济川,谢海义,姜宝法等,-北京:高等教育出版社，2008.1]。原文采用的软件是SAS，将干预前的因变量值按照重复测量的一个点来处理。我们采用SPSS来完成，也做了一些相应的调整。因为研究的目的是考察心理治疗的作用，可以将开始治疗后的所有随访结果作为一个整体来对待，我们关注的是这个“整体”的变化，而将干预前的因变量值作为协变量来处理。另外在操作的过程中也发现了一个SPSS26中的水平2协方差结构设置的一个bug，当然我更觉得是由于使用的是破解试用版导致的，还好问题可以简单通过程序命令的修改来补救。 示例来自来自美国New Hampshire State的一个关于心理健康和毒品滥用的研究项目。该研究对患有严重精神疾病的毒品滥用者进行了为期3年（1998-2001年）的对比观察试验。在203名参与者中，98名被分配到了治疗组，105位分配到了对照组中，不同组的受试者接受了不同的心理治疗。在基线调查后，每6个月对受试者进行一个随访调查。 结局变量： QOL（general quality of life）：总体生活质量/总体生活满意度，1-7分，1表示很糟糕，7表示很快乐； Hospdum：二分类结局变量，0表示过去6个月中未曾在精神病院住院，1表示过去6个月中曾在精神病院住院。以后会利用这个结局变量示例二分类资料的多层混合效应模型。 水平1解释变量，也称低水平变量、个体内变量，随时间变化而变化的受试者个体内测量（within-subject measures），随测量时间的变化可在个体内和个体间取不同的值。 time：随访时间。0表示基线，1-6表示基线后每6个月随访一次； SATS：Substance Abuse Treatment Scale，毒品滥用治疗程度，随时间变化而变化的个体内协变量（时依协变量），取值1-8分，分值越高代表越好的解毒效果。 水平2解释变量，也称高水平变量、个体间变量，是时间恒定变量，不随时间变化而变化的受试者个体间测量（between-subject measures）。 Trt：基线调查时的心理治疗，0=Control，1=Treatnment； gender：性别，0=女，1=男； age：基线调查时的年龄； Edu：教育水平，0=初中及以下，1=高中及以上； Diagnosis：最初的精神疾病诊断，0=躁狂抑郁症，1=精神分裂及压抑症。 变量的处理与数据按需调整： 正如开始所言，我们将基线值作为协变量来处理，因此在数据集中去掉了time=0的记录，并增加了一个变量Qt0，其值为QOL基线值，Qt0是时间恒定变量，属于水平2变量； 为了使模型的截距和交互作用的解释具有实际意义，对连续变量age、Qt0、SATS进行总均数中心化处理，处理后的变量分别对应agec、Qt0c、SATSc；本例将time也作为连续变量处理（最后我们也会按分类变量纳入做下对比），对time进行重新编码，原1-6分别编码为-5，-4，-3，-2，-1，0，新编码变量名称为time36，这样新编码的变量便有了有实际意义的0值。 【0】数据初判：结局测量随时间的变化趋势 在进行建模分析前，我们先作图（Graphs>>Chart Builder…：line）看一下数据的趋势。每个受试者QOL随时间的变化趋势如下图（图中仅选择了其中20例）。研究对象的QOL在研究开始时（基线）各不相同，每个对象的QOL变化率也不同。在整个随访期间，有的研究对象QOL呈现上升趋势，有些则出现下降。初始QOL比较高的受试者在研究结束的时候QOL也倾向高一些，表现出数据的相关性。另外我们也发现并非所有的受试者都完成了3年的随访研究。
本次笔记我们考虑使用的是多层混合效应模型，混合效应模型也在线性模型的范畴。我们可以先考察一下QOL是否随时间变化呈线性趋势。总体而言，所有受试者的平均QOL在3年期间呈上升态势，但看起来似乎并不是直线上升趋势。
统计学上的直线趋势未必就是视觉上的直线。本例去除基线后有6个时间点，理论上最多可以拟合5次方的时间多项式模型。我们关注的是去掉time=0基线的数据，去除基线后，曲线剩余1个弯曲，我们可以考虑拟合时间的2次方模型（当然也可考虑从更高阶开始拟合）。我们可以构建时间的2次方和3次方变量time2和time3，然后将time、time2和time3与QOL进行直线回归，当高次项系数没有统计学意义时表明模型不满足该次项的模型，然后再去掉该次项后再进行回归。注意线性、平方和高次方项之间可能存在共线性问题，可以将time先进行中心化再进行检验。结果表明2次方项和3次方项的系数均无统计学意义。 【1】随机截距模型 （方差成分模型） 1.1 空模型（最简单的随机截距模型） QOLij=γ00+μ0j+eij 水平1（个体内）模型：QOLij=β0j+eij 水平2（个体间）模型：β0j=γ00+μ0j 空模型，也称零模型、截距模型、无条件均值模型。建立该模型的目的就是计算组内相关系数（ICC），评估组内同质性或者组间异质性，确定是否有必要进行多水平模型分析。 QOLij表示第j个研究对象（水平2，个体间）第i次测量（水平1，个体内）的QOL值，空模型中只有常数项γ00，表示总截距，是所有研究对象QOL的总均数；μ0j表示水平2残差，是研究对象j的QOL均值与QOL总均值间的差异，其方差为σμ02；eij表示水平1残差，是个体j内各次测量间的差异，其方差为σ2。在多层混合效应模型中，用表示数据离散程度的方差来作为随机效应的衡量指标的。β0j表示研究对象j的QOL测量均数。 Analyze>>Mixed Models>>Linear… 个体变量选入id，重复变量选入time36，重复协方差结构选择Scaled Identity；
主体变量对话框中指定的指示变量，取值相同的个体，个体内的数据可能不独立（个体内存在集聚性）。在纵向数据发展模型中，表示对该变量代表的个体进行重复测量； 在建立空模型时，模型不含任何解释变量。重复测量指示变量选入或者不选入并不影响空模型的结果，但是选入状态下可以获得残差的协方差矩阵。选入原始的变量time还是选入重新编码的time36结果也都是一致的。 重复测量的协方差常用结构介绍如下：
协方差结构详情可参见与本文同时发送的“Covariance Structures”一文，IBM SPSS出品。 将QOL选入因变量框； 【随机】按钮对话框中选中“包含截距”的复选框来指定截距的随机项μ0j；将id从主体框中选入联合框，用以设定μ0j在不同的id个体中是不同的。模型的截距（固定效应γ00）会在在【固定】按钮中默认的“包含截距”指定。随机效应对话框中也存在一个协变量类型，这里指的是水平2残差的方差/协方差结构。当只有一个随机效应时，比如随机截距模型只有截距是随机的，当然就不会存在与其他随机效应变量的关系，此时不论选择什么结构，都会按照identity执行。这里面有一个前面个体内协方差结构中没有的结构：方差成分（Variance Components）， This structure assigns a scaled identity (ID) structure to each of the specified random effects。 【统计量】中选中模型的一些参数。
空模型主要结果如下： 模型基本信息
当前模型中只有一个常数项γ00，随机效应也只有μ0j。QOL以个人（id）为集聚水平（每个个体不同的测量间可能相关），个体间协方差结构默认方差成分，如在预定义对话框中指定time36为重复测量指示变量，表中将显示个体内协方差结构，本例指定为Scale Identity。有三个参数需要估计：γ00、μ0j、eij。 信息准则
信息准则本身不能说明单一模型拟合的优劣，可用于不同模型的比较，在比较时，一般都是遵循“小即好”的原则。-2LL即-2倍的对数似然值，在作为拟合优度统计量考察模型的拟合优度时，实际上是考察的设定模型与饱和模型之间的自然对数似然值之差乘以负2。饱和模型完美拟合数据，似然值为1，对数似然值为0，设定模型的-2LL实际上就是偏差统计量，反映的是设定模型与饱和模型拟合数据的差别。-2LL也用于嵌套模型间的比较，当一个模型是另外一个模型的子模型时，如前面说的设定模型与饱和模型，可以采用似然比检验（LR test）进行检验，似然比检验的检验统计量就是两个模型-2LL之差（偏差），偏差统计量之差近似满足卡方分布，自由度为两模型的参数数目之差，LR检验可以很好地检验新加入的变量是否具有统计学意义。信息准则一般都是对-2LL进行参数数目的调整，甚至考虑样本量后的校正（AICC、BIC）。信息准则的优势在于不仅可以比较嵌套模型也可以比较非嵌套模型。 另外，需要特别说明的是混合效应模型的参数估计方法，具体可在【估计】按钮中进行设置，默认是REML。当要比较的两个模型随机效应相同、固定效应不同时，模型估计需采用最大似然法（ML）；当随机效应（方差/协方差成分）不同、固定效应相同时，模型估计需采用限制性最大似然法（REML）；当随机效应和固定效应都不同，但两模型是嵌套模型时，也可采用ML。模型的这些比较原则不仅适用于LR检验，也适用于信息准则。当模型完成比较后，仍需要REML对最终模型进行估计，毕竟REML比ML估计更精确。 固定效应结果
固定效应的检验和参数估计：本例为空模型，只有常数项γ00，估计值为4.42，即所有研究对象测量值的总均数为4.42；原假设为常数项为0，结果有统计学意义，即所有个体的常数项不全为0。固定效应和参数估计的统计学检验是一致的，F=t2。 随机效应结果（随机效应使用协方差来表示的）
协方差参数估计：重复测量方差（或残差方差）即模型中水平1随机效应（eij）的方差，其估计值为0.954，且有统计学意义，说明研究对象不同测量之间存在差异（显而易见的结论）。截距方差即水平2截距的随机效应（μ0j）的方差，空模型中代表的研究对象的QOL的测量均值与总均值之间的差异，其估计值为0.968，具有统计学意义，表明QOL的变异在不同的研究对象（水平2）上存在集聚性，即在治疗后的QOL在不同的个体间存在差异。组内相关系数ICC=0.954/(0.954+0.968)=0.4964，QOL的总变异中有相当一部分（49.64%）是由不同研究对象引起的。个体间差异显著，应进行多水平模型的分析。 值得注意的是在SPSS中直接给出var(μ0j)是否为0的Wald Z的检验结果和估计参数的95%CI。一般情况下“置信区间是否包含无效假设”与“假设检验的结果”是一致的，但这种一致是在判断“是否等于”的双侧检验时。在随机效应（用方差表示）检验时应注意，方差不＜0，无效假设是方差=0，备择假设就是方差＞0，检验为单侧检验，所以此处用95%CI来大体判断随机效应是否有统计学意义是不合适的。 后面两个表可以查看G矩阵（水平2残差的方差/协方差矩阵）和R矩阵（水平1残差的方差/协方差矩阵）的关系。在多层混合效应模型中，要求水平1残差符合正态分布，水平2残差满足多元正态分布，水平1残差和水平2残差相互独立（G和R协方差为0，相关系数为0）。 1.2 带重复测量变量的随机截距模型 QOLij=γ00+γ10timeij+μ0j+eij 水平1（个体内）模型： QOLij=β0j+β1jtimeij+eij 水平2（个体间）模型： β0j=γ00+μ0j，β1j=γ10 QOLij表示第j个研究对象（水平2，个体间）第i次测量（水平1，个体内）的QOL值。在加入时间变量后，模型中的常数项γ00，已不再对应所有研究对象所有QOL的总均数，而是time取值为0时的研究对象的QOL均数；μ0j表示个体特异性，因其存在，QOL在time=0时便因人而异；eij是个体j每次测量间的差异；模型要求不同观察水平的随机残差项μ0j和eij是相互独立的。β0j和β1j表示研究对象j在time=0时的QOL测量值和对象j随时间的QOL变化率。当前模型将time按连续变量处理并设为固定效应，其变化率β1j不随时间发生变化。 个体特异性效应μ0j是随机的效应，假定为服从均数为0，方差为σμ02的正态分布；误差eij被假设服从均数为0方差为σ2的正态分布。 在当前模型中，γ00的实际意义不大也并不是很重要。因为实际的time=0是基线值，但我们的处理中并未将time=0作为重复测量的一个点来分析，所以当前模型中的γ00只是根据模型推断出来的一个值。为了使截距具有实际意义，可以对time重新进行编码，比如将原来的1-6编码为0-5，或者数据中心化，或者虑将1年半时（18个月）作为中心点将原编码1-6改为-2、-1、0、1、2、3。本例变量time36将重复测量的最后一个随访赋值为0，以其为主效应更能代表治疗后的效果，以此建立的模型常数项γ00表示time36=0时的QOL均值，此时time36=0已经具有了实际意义。如果将time按分类变量处理，在SPSS里面常将高赋值的作为默认参照水平，这样可以很好的体现治疗最后测量的效应，time=0表示的往往是最后一个测量点的QOL均值。甚至在有些软件中，该模型γ00依旧表示所有研究对象QOL的总均数，通过μ0j和eij的调整，依旧可以获得类似的模型。这些问题其实都不是问题，因为不同的编码改变的只是截距值，并不影响我们最为关注的系数和方差估计结果。 当前模型： QOLij=γ00+γ10time36ij+μ0j+eij Level 1 Model: QOLij=β0j+β1jtime36ij+eij Level 2 Model: β0j=γ00+μ0j，β1j=μ0j Analyze>>Mixed Models>>Linear… 主体及重复测量指定对话框同空模型； 将QOL选入因变量框；本例随访时间按连续变量处理，time36选入协变量框。 【固定】按钮中将time选入模型，【随机】、【统计量】同1.1空模型。
主要结果 基本信息：当前模型的主要信息如上表，详情可参照1.1空模型，需要估计4个参数：γ00、γ10、μ0j、eij。
信息准则：-2LL和信息准则介绍可参照空模型部分。
与空模型相比，带重复测量因素的随机截距模型的-2LL变小了，但信息准则却变大了。但正如空模型部分所言，当要用-2LL和信息准则比较随机效应相同、固定效应不同的两个模型时，模型估计需采用ML法，SPSS中默认的是REML法，我们可以重新运行空模型和带时间变量的随机截距模型，采用ML估计法（【估计】按钮中修改即可），获得空模型和带时间因素的随机截距模型的-2LL和信息准则如下：
结果表明-2LL和信息准则均有下降，但幅度并不大。我们可以用似然比检验（LR test）来检验一下新加入的变量（模型的改变）是否有统计学意义： Chi2≈-2LL0-(-2LL)=3564.674-3553.619=11.055，df=4-3=1，只需要在excel的一个单元格中输入“=CHISQ.DIST.RT(11.055,1)”或者“=1-CHISQ.DIST(11.055,1,TRUE)”，结果会直接输出P值，本例P=0.00088<0.05，表明加入变量time，对模型的改善是有统计学意义的。 固定效应 详情可参照空模型部分。常数项γ00估计值为为4.56，具有统计学意义，其实际意义表示所有观测对象在随访结束的时（3年）的QOL均值；不同时间点的QOL是不同的，平均而言，每过半年QOL约升高0.057分，且具有统计学意义（P=0.001）。
随机效应 个体内水平（水平1）的值既有组内（个体内）变异，也有个体间（水平2）变异，因此个体内水平的解释变量既可在个体内水平上影响结局的变异，也可以在个体间水平影响结局的变异。而个体间水平（水平2）解释变量仅影响个体间水平的结局变异，并不影响低水平结局的变异，因为个体间水平变量对同一个组内的个体值相同。
与空模型相比，在加入重复测量的时间变量后，个体内残差和个体间残差均有变化。个体内残差方差由原来的0.9539变成0.9438，依旧具有统计学意义（Wald Z=8.546，P<0.001）,说明有一部分变异虽然被新引入的变量解释掉了，但仍然存在显著的个体内变异。 个体间变异也会受到个体内水平解释变量的影响。与空模型相比，当前模型的个体间变异在加入时间变量后变大反而变大（0.9682→0.96962），这只有在新加入个体间方差为负值才讲得通，但理论上方差不为负。这种负方差的问题在重复测量数据的对水平模型分析时较为常见，主要是由于结果的重复测量部分（水平1）的个体变异通常比抽样模型假设要大，这导致了只含截距的模型中的测量间变异被高估，而个体间方差被低估。组间方差比例由49.64%变为50.67%，仍有超过一半的总变异是由于个体间异质性引起的。在王等的《多层统计分析模型：方法与应用》中并不是把空模型作为计算ICC的模型，而是从带重复测量变量的随机截距模型中计算。 【2】建立含重复测量变量的随机截距-斜率模型 QOLij=γ00+γ10timeij+(μ0j+μ1jtimeij+eij) 水平1（个体内）模型：QOLij=β0j+β1jtimeij+eij 水平2（个体间）模型：β0j=γ00+μ0j ,β1j=γ10+μ1j 模型解读同1.2含重复测量变量的随机截距模型, 不同的是在1.2模型中将time设为固定效应，其变化率β1j=γ10不随时间发生变化，当前模型将time设为随机效应，γ10表示所有个体的平均变化率（总斜率），μ1j表示个体j的变化率与总的平均变化率（总斜率）的差异。 当前模型随机部分包括三部分：个体间水平上的μ0j和μ1j，个体内水平上的eij。个体内水平上的随机效应eij涉及多个重复测量点，其方差协方差可以形成一个重复测量的协方差矩阵R，个体间存在2个及2个以上随机效应，其方差协方差可以形成一个组件协方差矩阵G。 采用重新编码后的time36作为重复测量因素建模： QOLij=γ00+γ10time36ij+(μ0j+μ1jtime36ij+eij) Level 1 Model: QOLij=β0j+β1jtime36ij+eij Level 2 Model: β0j=γ00+μ0j ,β1j=γ10+μ1j Analyze>>Mixed Models>>Linear… 主体及重复测量指定对话框同截距模型； 将QOL选入因变量框；本例随访时间按连续变量处理，time36选入协变量框。 【随机】按钮中协方差结构选择“无结构”，选中“包含截距”复选框，将time36选入模型，id选入联合框；【固定】【统计量】同1.2带重复测量变量的截距模型。此处需要特别说明一下将模型斜率也设为随机效应时，水平2就有2个随机效应：μ0j和μ1j，需要指定协方差结构。本例选择无结构（UN），选择无结构协方差结构时建议在【估计】选择采用Kenward-Roger法来计算固定效应检验的分母自由度。另外在我选用UN结构时，软件运行的时候依旧用的ID（Identity），不知道是SPSS26的一个Bug还是试用版破解的问题，这个问题可以通过[Paste]按钮修改程序来进行。
主要结果如下：
当前模型需要估计6个参数。当前模型有两个固定效应：γ00=4.561(t=54.580,P<0.001)和γ01=0.057(t=2.945,P=0.004)；有3个随机效应：随机截距方差估计σμ02的=0.909(Wald Z=6.414,P<0.001)，表示QOL在time36=0时（最后1次随访）在研究对象之间存在显著性差异；重复测量变量的随机斜率方差估计σμ12的=0.021(Wald Z=2.595,P=0.009)，表示QOL随时间的变化的斜率在不同的研究对象之间也存在显著差异；随机截距和斜率间的协方差估计σμ012的=0.110(Wald Z=0.423,P=0.672)，协方差为正表明是正相关关系，表示QOL在最后1次随访时得分越高，其值随时间的变化率就越大。模型在加入了随机斜率后，个体内残差方差由0.9438进一步变小为0.8716。-2LL与信息准则与原来相比有进一步下降，大家可以参照1.2模型中采用LR进行检验新加入的随机效应是否有统计学意义，不一样的是当随机效应（方差/协方差成分）不同、固定效应相同时，模型估计需采用REML法[Chi2=12.962,df=2,P=0.0015]。 【3】治疗效应评估（在随机截距-斜率模型中加入个体间的治疗效应） QOLij=γ00+γ01Trtj+γ10timeij+γ11Trtj*timeij+(μ0j+μ1jtimeij+eij) 水平1（个体内）模型： QOLij=β0j+β1jtimeij+eij 水平2（个体间）模型： β0j=γ00+γ01Trtj+μ0j ； β1j=γ10+γ11Trtj+μ1j 通过方程来可以会对系数的实际意义有更深刻地理解。γ00是Trt=0的研究对象在time=0时的结局测量（QOL）。此处要特别注意，Trt=0和time=0指的都是参照水平，并不一定是实际赋值，在SPSS中会将分类变量高水平为参照，连续变量较小的值作为参照水平。γ01则是Trt=1与Trt=0的研究对象的QOL在time=0时的平均水平的差值。对斜率而言，γ10是Trt=0的QOL的平均变化率，γ11则是Trt=1与Trt=0的QOL的随时间的平均变化率差异。eij是水平1残差，表示个体内变异，假设其满足均数为0方差为σ2的正态分布，μ0j和μ1j是水平2间或跨组的变异，假设其满足均数为0方差分别为σμ02和σμ12的多元正态分布。
回归到本例的模型： QOLij=γ00+γ01Trtj+γ10time36ij+γ11Trtj*time36ij+(μ0j+μ1jtime36ij+eij) Level 1 Model: QOLij=β0j+β1jtime36ij+eij Level 2 Model: β0j=γ00+γ01Trtj+μ0j ； β1j=γ10+γ11Trtj+μ1j 本例Trt是分类变量，高水平（治疗组，实际赋值为1）为参照水平，time36=0是最后一个随访点，γ00是治疗组的研究对象的最后一个随访点的QOL。γ01就是对照组与治疗组研究对象的QOL在最后一个随访点的平均差值。γ10是治疗组的QOL随时间的平均变化率，γ11则对照组与治疗组的QOL的随时间的平均变化率差异。 Analyze>>Mixed Models>>Linear… 主体及重复测量指定对话框同空模型； 将QOL选入因变量框；Trt选入因子变量框，本例随访时间按连续变量处理，time36选入协变量框。 【固定】按钮中将Trt、time36及两者交互项Trt*time36选入模型。 【随机】、【统计量】、【估计】、【统计量】同2带重复测量变量的截距模型。
主要结果如下：
结果显示，Trt的主效应不显著，γ01=-0.044，表示在最后一次随访时时，对照组比治疗组的QOL要低，但统计不显著（P=0.791）。时间主效应γ10=0.069，表明治疗组的QOL随时间的变化率为0.069（P=0.012）。Trt与时间变量time的交互作用γ11=-0.025（P=0.527），表明相对治疗组而言，对照组的QOL增长的更少一些，只是这个额外的增长不具有显著的统计学意义。与带时间变量的随机截距-斜率模型相比，模型拟合程度并未得到显著改善（Chi2=3540.485-3540.893，df=2，P=0.815，注意随机效应相同、固定效应不同时模型估计需采用ML法）。从模型拟合角度，无需将Trt放入模型。但我们的研究目的是研究心理治疗作用的效果，因此我们需要将其保留在模型中，但交互项无统计学意义可以去掉。 个体间解释变量的加入只对水平2残差产生影响。从随机效应上看，随机截距和随机斜率的方差变化也不大，仍然具有统计学意义。表明Trt的加入对模型QOL的离散程度的解释力度也不大。截距随机效应和斜率随机效应都具有统计学意义，而两者的协方差不具有统计学意义（P=0.672）,当前随机效应的方差协方差结可能满足Diagonal结构。 模型估计方法为REML，几种不同的误差/残差方差协方差结构对应的-2LL和信息准则如下：
笔者最终AIC为标准选定DIAG-AR1的方差协方差结构，该结构下固定效应和随机效应的结果就不再放上来了。当然SPSS里面提供了非常多的协方差结构，笔者也只是选择了有限的几种。 这里也涉及到以AIC和BIC为标准结果发生冲突的时候该以哪个未标准的问题。-2LL是AIC和BIC的基础，BIC的惩罚项比AIC大，考虑了样本个数，可以防止模型精度过高造成的模型复杂度过高。AICC是对AIC的校正，CIAC是对BIC的校正，随着样本量的增加AICC收敛于AIC、CIAC收敛于BIC。来自网络的观点是“AIC是从预测角度选择一个好的模型用来预测。而BIC是从拟合角度，选择一个对现有数据拟合最好的模型，从贝叶斯因子的解释来讲，就是边际似然最大的那个模型”。 【4】控制其他影响因素 QOLij=γ00+γ01Trtj+γ02genderj+γ03ageij+γ04Eduj+γ05diagnosisj+γ06Qt0j+γ10timeij+γ11Trtj*timeij+γ20SATSij+(μ0j+μ1jtimeij+eij) 水平1（个体内）模型： QOLij=β0j+β1jtimeij+β2jSATSij+eij 水平2（个体间）模型： β0j=γ00+γ01Trtj+γ02genderj+γ03ageij+γ04Eduj+γ05diagnosisj+γ06Qt0j+μ0j, β1j=γ10+γ11Trtj+μ1j β2j=γ20 这里面需要正确判断哪些测量属于水平1（个体内）变量，哪些测量水域水平2（个体间）变量。在重复测量模型中，个体水平是水平2，时间恒定的变量（如性别、种族等）只有个体间变异，属于水平2变量；随测量时间的变化可在个体内和个体间取不同的值的个体内变量，属于水平1变量。本例中个体水平上的性别、初始年龄、教育水平、疾病诊断及QOL的基线水平都是时间恒定变量， 属于个体间变量；SATS（毒品滥用治疗程度）是随时间变化而变化的个体内协变量（时依协变量），跟time一样是个体内变量。对时依协变量的处理，在王等《多层统计分析模型：方法与应用》中提供了两种解决策略，一是直接按个体内变量纳入，二是进一步分解成两个变量：个体均值和均值的离均差，按两个变量纳入建模与直接纳入相比较。 同样的，为使解释更有实际意义，我们对时间变量进行了重新编码，其他连续变量最终都进行了中心化处理： QOLij=γ00+γ01Trtj+γ02genderj+γ03agecij+γ04Eduj+γ05diagnosisj+γ06Qt0cj+γ10time36ij+γ11Trtj*time36ij+γ20SATScij+(μ0j+μ1jtime36ij+eij) Analyze>>Mixed Models>>Linear… 主体框选入id，重复测量指定对话框选入time36，重复协方差类型选择AR(1)； 将QOL选入因变量框；Trt、gender、Edu、diagnosis选入因子变量框，本例随访时间按连续变量处理，time36选入协变量框，其他选入协变量框的还有agec、Qt0c； 【固定】将Trt、time36及两者交互项Trt*time36选入模型； 【随机】随机效应协方差类型选入Diagonal，选择“包含截距”复选框，将time36选入模型，id选入联合框； 【统计量】选择需要的统计量； 【估计】采用默认。 如结果不能按设定的协方差结构运行，在主对话框点击【Paste】修改一下命令即可。
主要结果如下：
在加入控制变量后，模型的信息准则出现了较大的下降，感兴趣的可以参考前面做下模型的LR检验。我们发现在加入控制变量后，变量Qt0c(γ06=0.425,P<0.001)、agec(γ03=-0.020,P=0.010)及SATSc(γ20=0.136,P<0.001)对结局QOL有影响外，其他变量均无统计学意义。总体来说，QOL基线值越高、年龄越小、毒品滥用治疗越好，患者的的生活满意度（QOL）越高。 另外一个比较有意思的结果是在加入变量SATSc之前，时间（time36）和诊断结果（diagnosis）都是具有统计学意义的，而在加入SATSc之后，这两个变量就都没有统计学意义了，这说明吸毒者的QOL随时间而改善的效果是他们接受了毒品滥用，治疗效果随时间变化而得到改善的结果，且不论是躁狂抑郁症还是精神分裂压抑症（diagnosis），随着毒品滥用治疗效果的出现，都可以有类似的生活满意度。加入变量SATSc之前的结果如下：
结果中我们也发现，交互项Trt*time是没有统计学意义的。为简化模型，我们可以考虑去掉交互项。但我们的目的并不是寻找生活满意度的影响因素，而是考察心理治疗作用对生活满意度的效应，因此Trt无论有没有统计学意义，都应该作为主要研究因素保留在模型中，其他因素作为校正因素也保留在模型中。最终主要结果如下：
最后，我们也可以将时间变量按分类变量纳入。当时间按分类变量纳入时，相当于同时纳入了K-1个哑变量。比如本例时间因素按分类变量纳入、随机变量的方差协方差类型设置为CS-AR1，固定效应的分析结果可以参考如下（未收敛）：
当然我们建立的这些模型有一个前提假设，就是这些控制变量与Trt、time之间并没有交互作用，大家感兴趣的可以去检验一下到底有没有。另外SATS会不会也是随机变量，大家也可以取尝试一下看。篇幅已经很长了，就不再啰嗦了。 总算是结尾了！ |
|
|
来自： Memo_Cleon > 《待分类》
