重复测量资料的统计分析方法简介.doc

重复测量资料的统计分析方法简介在医学研究中，一些干预研究和纵向研究都需要对研究对象进行随访，每次随访进行观测或测量一些效应指标，考察同一研究对象同一指标的变化情况。同一个对象的多次观察或测量所获得的资料称为重复测量的资料。由于同一对象同一指标的相邻两个时间点的效应指标观测值往往是相关的，也就是重复测量的资料存在不独立的问题，然而大多数的医学统计方法都要求资料是独立，所以这些资料的统计分析需要用比较特殊的统计方法进行分析。重复测量资料的统计分析方法可以用重复测量的方差分析，也可以用混合回归模型（Mixed regression Model），由于重复测量的方差分析要求资料满足球形对称性（可以理解为相关资料情况下的方差齐性），而Mixed回归模型并不要求资料满足球形对称，并可以借助计算机统计软件对未知参数进行限制的最大似然估计，其他统计分析的思想都是类似的。本节将主要介绍如何借助统计软件应用Mixed回归模型对重复测量资料进行统计分析。为了帮助读者对重复测量资料分析有一个简单的了解，本节将举一个非常简单的例子初步说明重复测量资料的统计分析概况。例1 为了比较A药和B药在疗程为6个月中的持续减肥的疗效，现有10个身高为160cm的女性肥胖者志愿参加这项研究。随机分成2组，每组各5人。分别考察这2组肥胖者在服药前、服药3个月和服药6个月的体重变化。这2组肥胖者在服用该药前、服药3个月和服药6个月的体重测量值(kg)见表1。表 1 2组肥胖者在服用该药前、服药3个月和6个月的体重组别和肥胖者编号服药前（）3个月（）6个月（）A药组1号524942A药组2号515046A药组3号504941A药组4号514944A药组5号494740B药组1号515453B药组2号494746B药组3号504744B药组4号494841B药组5号525048这是两组观察对象的多个测量时间点的重复观察测量资料，同样对于同一对象的不同观察时间点的观察资料是相关的，但由于需要比较两个药的减肥疗效，所以采用两因素方差分析，随机区组设计的方差分析或Friedman秩检验的统计方法都不适用于本例的数据统计分析，但可用Mixed模型对本例资料进行统计分析。设A药组对象在服药前体重总体均数为；服药3个月后的总体体重改变量为，故服药3个月时的体重总体均数为；服药6个月时，体重比服药前的总体改变量为，即服药6个月时的体重总体均数为；同理，设B药组对象在服药前体重总体均数为；服药3个月后的总体体重改变量为，故服药3个月时的体重总体均数为；服药6个月时，体重比服药前的总体改变量为，即服药6个月时的体重总体均数为；为了便于两组比较，引入两组比较的差异参数如下：记服药前的两组差异为b0m1m0，即：B组服药前的总体均数可以表示为；在服药3个月时A药和B药的体重总体改变量分别为b1和b11 ，记服用B药和A药3个月时的体重总体改变量的差异为，即在服药3个月时B药的体重总体改变量可以表示为；在服药6个月时A药和B药的体重总体改变量分别为b2和b12 ，记服用B药和A药6个月时的体重总体改变量的差异为，即在服药6个月时B药的体重总体改变量可以表示为，把两组差异的参数代入上述表达式，得到下列2组肥胖者在服用该药前、服药3个月和6个月的体重总体均数表达式如表2所示：表2 两组3个时点总体均数表达式服药前(t1=0，t2=0)服药3个月(t1=1，t2=0)服药6个月(t1=0，t2=1)A组总体均数(g=0)B组总体均数(g=1)不难验证表2的总体均数表达式可以用下列总体回归方程(式1)表示： (1)由于不同对象之间存在个体差异,同一对象不同时点之间也存在随机差异,因此第g组第i个对象第t时刻的体重观察值可以用式（2）表示为 （2）并且假定，称式（2）为混合线性模型（Mixed Model）。若和不全为0，则称两种药物与服药时间对疗效有交互作用。两组在3个时间点的总体均数差异分别为，和，因此只需检验H0：、H0：和 H0：就可以推断两组总体均数差异。反之若和全为0，则称两种药物与服药时间对疗效无交互作用，并且两组各个时间点的总体均数差异均为，因此只需检验H0：b30就可以推断两组的总体均数差异。我们同样借助Stata软件对上述资料用混合模型进行统计分析，相应的Stata软件的数据格式如下ygnot1t2520100490110420101510200500210460201521100050110104811001Stata操作命令如下：gen gt1=g*t1 产生交互作用项变量gt1gen gt2=g*t2 产生交互作用项变量gt2xtreg y t1 t2 g gt1 gt2,i(no)Random-effects GLS regression Number of obs = 30Group variable (i) : no Number of groups = 10R-sq: within = 0.8288 Obs per group: min = 3 between = 0.0973 avg = 3.0 overall = 0.5927 max = 3Random effects u_i Gaussian Wald chi2(5) = 78.32corr(u_i, X) = 0 (assumed) Prob chi2 = 0.0000- y | Coef. Std. Err. z P|z| 95% Conf. Interval-+- t1 | -1.8 1.053565 -1.71 0.088 -3.86495 .2649502 t2 | -8 1.053565 -7.59 0.000 -10.06495 -5.93505 g | -.4 1.612452 -0.25 0.804 -3.560347 2.760347 gt1 | .8 1.489966 0.54 0.591 -2.120281 3.720281 gt2 | 4.2 1.489966 2.82 0.005 1.279719 7.120281 _cons | 50.6 1.140175 44.38 0.000 48.3653 52.8347-+- sigma_u | 1.9300259 sigma_e | 1.6658331 rho | .57307692 (fraction of variance due to u_i)-由此得到下列统计结论：在服药前，A组的体重总体均数的估计值为50.6(kg)服药3个月时，A组的体重改变量的总体均数为b1的估计值为-1.8(kg)，P值0.0880.05，因此没有足够证据推断在服用A药3个月时，服用A药的人群平均体重发生改变。服A药6个月时，A组的体重改变量的总体均数为b2的估计值为-8(kg)，P值0.05，差别无统计学意义。服药3个月时，服药A药和服拥B药的两个人群平均体重下降的差异的估计值为0.8，P值=0.591，差异无统计学意义。服药6个月时，服药A药和服拥B药的两个人群平均体重下降的差异的估计值为4.2，P值=0.005 chi2 = 0.8041相应的P值0.8041a，差异无统计学意义，故无证据显示两组总体均数不等。服药6个月时的两组平均体重比较的Stata命令和输出结果如下：test ggt20 （H0：0）( 1) g + gt2 = 0.0 chi2( 1) = 5.55 Prob chi2 = 0.0184相应的P值0.01840，说明logit(PA) logit（PB），相应的A组转阴率PA高于B组转阴率PB；若b10，说明logit(PA) chi2 = 0.1974Log likelihood = -61.856052 Pseudo R2 = 0.0133- y | Coef. Std. Err. z P|z| 95% Conf. Interval-+- g | .5564203 .4345529 1.28 0.200 -.2952877 1.408128 _cons | -1.045969 .322408 -3.24 0.001 -1.677877 -.4140606-得到b0的估计值为-1.045969，b1的估计值为0.5564203，b1的P值=0.200，因此对于检验水准a0.05而言，两组疗效的差异无统计学意义。故没有充分证据说明两组疗效有差异。由于HP经常呈现反复，这次检查可能是阴性，但下次检查可能又呈阳性了，所以往往要考察几个疗程才能评价两种治疗方案的疗效。例3 收集100名HP阳性患者，随机分成两组，分别用两种治疗方案进行治疗（A方案g=1，B方案g=0），2个月为一个疗程并检查一次HP，连续治疗3个疗程，考察两种治疗HP方案的疗效（阴性y=1,仍为阳性y=0），比较两组疗效的差别。表5 两种治疗方案3个疗程的HP转阴情况一个疗程t1=0，t2=0二个疗程t1=1，t2=0三个疗程t1=0，t2=1A组（g=0）人数B组(g=1)人数阳性(y=0)阳性(y=0)阳性(y=0)11阳性(y=0)阳性(y=0)阴性(y=1)34阳性(y=0)阴性(y=1)阳性(y=0)71阳性(y=0)阴性(y=1)阴性(y=1)2015阴性(y=1)阳性(y=0)阳性(y=0)11阴性(y=1)阳性(y=0)阴性(y=1)53阴性(y=1)阴性(y=1)阳性(y=0)54阴性(y=1)阴性(y=1)阴性(y=1)821由于患者转阴后还有一定的概率转阳，因此不能用生存分析的方法进行比较，但本例资料可以认为是重复测量的二分类资料，可以用下列随机效应的Logistic模型：.（6）其中uN(0,s2)该模型的特点是常数项含有随机误差项以及随机效应与协变量无关以及与回归系数无关。也可以写成log(Odds)=（7）b0，b1，bm为参数，一般采用限制的最大似然估计，由于这种估计方法非常复杂，本书不作介绍，我们仅介绍借助Stata软件对参数进行估计和检验。本例模型可以表示为（8）式，其中变量定义见表12.9。（8）其中u服从N(0,s2)，称为Logit(P)的确定性部分，我们主要关心的是Logit(P)的确定性部分，其与各个变量之间的对应列表如下。表6 两种治疗方案与多个疗程的Logit(P)的确定性部分关系一个疗程（t1=0,t2=0）二个疗程（t1=1,t2=0）三个疗程（t1=0,t2=1）方案A（g=1）b0b3b0b1b3b4b0b2b3b5方案B（g=0）b0b0b1b0b2两种方案Log(Odds)的差异b3b3b4b3b5Odds Ratioexp(b3)exp(b3b4)exp(b3b5)因此两种治疗方案在一个疗程治疗的转阴率比较就是检验H0：b30；在二个疗程治疗的转阴率比较就是检验H0：b3b40；在三个疗程治疗的转阴率比较就是检验H0：b3b50。如果两种治疗方案与随访时间无交互作用，即b40且b50，则上述模型可以简化为 （9）式，相应的两种治疗方案与多个疗程的Logit(P)的确定性部分关系见表7。（9）表7 两种治疗方案与多个疗程的Logit(P)的确定性部分关系一个疗程（t1=0,t2=0）二个疗程（t1=1,t2=0）三个疗程（t1=0,t2=1）方案A（g=0）b0b0b1b0b2方案B（g=1）b0b3b0b1b3b0b2b3两种方案的差异b3b3b3第二疗程与第一疗程比较：B方案：（b0b1b3）（b0b3）b1 A方案：b0b1b0=b1第三疗程与第一疗程比较：B方案：（b0b2b3）（b0b3）b2 A方案：b0b2b0=b2第三疗程与第二疗程比较：B方案：（b0b2b3）（b0b1b3）=b2-b1 A：方案：b0b2(b0+b1)=b2-b1因此各个疗程的两种治疗方案转阴率比较就是检验H0：b30。本例数据和Stata数据结构如下t1t2nogy001001010001100002001020001201001001110100110110011本例的Stata操作命令gen gt1=g*t1 产生gt1项gen gt1=g*t2 产生gt2项xtlogit y t1 t2 g gt1 gt2,i(no) 即：xtlogit y 自变量,i(观察对象编号)结果如下Random-effects logit Number of obs = 300Group variable (i) : no Number of groups = 100Random effects u_i Gaussian Obs per group: min = 3 avg = 3.0 max = 3 Wald chi2(5) = 35.10Log likelihood = -165.70344 Prob chi2 = 0.0000- y | Coef. Std. Err. z P|z| 95% Conf. Interval-+- t1 | 1.875843 .4581372 4.09 0.000 .9779101 2.773775 t2 | 1.43401 .4290642 3.34 0.001 .5930595 2.27496 g | .8123216 .408646 1.99 0.047 .0113902 1.613253 gt1 | -.6822685 .6538294 -1.04 0.297 -1.96375 .5992135 gt2 | .0585067 .6575034 0.09 0.929 -1.230176 1.34719 _cons | -.4895482 .2913583 -1.68 0.093 -1.0606 .0815036-+- /lnsig2u | -14 672.3872 -1331.855 1303.855-+- sigma_u | .0009119 .3065689 6.2e-290 1.3e+283 rho | 8.32e-07 .0005591 0 .-Likelihood ratio test of rho=0: chibar2(01) = 0.00 Prob = chibar2 = 1.000交互项b4的估计值为-0.6822685，相应的P值为0.297，b5的估计值为0.0585067，相应的P值为0.929，因此我们需要检验H0：b4b50 vs H1：b4和b5不全为0，a=0.05,相应的Stata操作命令如下test gt1=0test gt2=0,a结果如下 ( 1) ygt1 = 0.0 chi2( 1) = 1.09 Prob chi2 = 0.2967 ( 1) ygt1 = 0.0 ( 2) ygt2 = 0.0 chi2( 2) = 1.38 Prob chi2 = 0.5024相应的P值为0.5024a，因此交互项没有统计学意义，因此可以选用（12.9）式的无交互作用logistic模型，相应的Stata命令为xtlogit y t1 t2 g,i(no)Random-effects logit Number of obs = 300Group variable (i) : no Number of groups = 100Random effects u_i Gaussian Obs per group: min = 3 avg = 3.0 max = 3 Wald chi2(3) = 33.89Log likelihood = -166.38816 Prob chi2 = 0.0000- y | Coef. Std. Err. z P|z| 95% Conf. Interval-+- t1 | 1.563917