重复测量资料的统计分析方法简介.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号
52
49
42
A药组2号
51
50
46
A药组3号
50
49
41
A药组4号
51
49
44
A药组5号
49
47
40
B药组1号
51
54
53
B药组2号
49
47
46
B药组3号
50
47
44
B药组4号
49
48
41
B药组5号
52
50
48
这是两组观察对象的多个测量时间点的重复观察测量资料，同样对于同一对象的不同观察时间点的观察资料是相关的，但由于需要比较两个药的减肥疗效，所以采用两因素方差分析，随机区组设计的方差分析或Friedman秩检验的统计方法都不适用于本例的数据统计分析，但可用Mixed模型对本例资料进行统计分析。
设A药组对象在服药前体重总体均数为；服药3个月后的总体体重改变量为，故服药3个月时的体重总体均数为；服药6个月时，体重比服药前的总体改变量为，即服药6个月时的体重总体均数为；
同理，设B药组对象在服药前体重总体均数为；服药3个月后的总体体重改变量为，故服药3个月时的体重总体均数为；服药6个月时，体重比服药前的总体改变量为，即服药6个月时的体重总体均数为；
为了便于两组比较，引入两组比较的差异参数如下：记服药前的两组差异为b0＝m1－m0，即：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：b3＝0就可以推断两组的总体均数差异。我们同样借助Stata软件对上述资料用混合模型进行统计分析，相应的Stata软件的数据格式如下
y
g
no
t1
t2
52
0
1
0
0
49
0
1
1
0
42
0
1
0
1
51
0
2
0
0
50
0
2
1
0
46
0
2
0
1
…
…
…
…
…
52
1
10
0
0
50
1
10
1
0
48
1
10
0
1
Stata操作命令如下：
gen gt1=g*t1 产生交互作用项变量g´t1
gen gt2=g*t2 产生交互作用项变量g´t2
xtreg y t1 t2 g gt1 gt2,i(no)
Random-effects GLS regression Number of obs = 30
Group variable (i) : no Number of groups = 10
R-sq: within = 0.8288 Obs per group: min = 3
between = 0.0973 avg = 3.0
overall = 0.5927 max = 3
Random effects u_i ~ Gaussian Wald chi2(5) = 78.32
corr(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.088>0.05，因此没有足够证据推断在服用A药3个月时，服用A药的人群平均体重发生改变。
服A药6个月时，A组的体重改变量的总体均数为b2的估计值为-8(kg)，P值<0.001，因此可以认为在服用A药6个月时，服用A药人群的平均体重已有下降。
服药前两组平均体重的差异的估计值为－0.4（kg）,相应的P值为=0.804>0.05，差别无统计学意义。
服药3个月时，服药A药和服拥B药的两个人群平均体重下降的差异的估计值为0.8，P值=0.591，差异无统计学意义。
服药6个月时，服药A药和服拥B药的两个人群平均体重下降的差异的估计值为4.2，P值=0.005<0.05，差异有统计学意义。
服用A药和B药在各个时间的总体均数估计见表3。

表3 在各个时间点的两组总体均数及其估计值

A组（g=0）
B组（g=1）
B组－A组
两组差异检验
服药前
(t1=0,t2=0)
服药3个月时
(t1=1,t2=0)
服药6个月时
(t1=0,t2=1)
总体
均数
总体均数
估计值
总体
均数
总体均数
估计值
总体
均数之差
总体均数
估计值
P值

50.6

50.2

-0.4
0.804

48.8

49.2

0.4
0.804

42.6

46.4

3.8
0.018
注：表中均数估计值是参数估计值和总体均数参数表达式计算所得。如：服药3个月时A组的总体均数估计值＝50.6-1.8＝48.8。

服药后的2个时间点的两组平均体重比较的Stata统计检验命令和输出结果如下
设a=0.05
服药3个月时的两组平均体重比较的Stata命令和输出结果如下：
test g＋gt1＝0 （H0：＝0 即服药3个月时的两组体重总体均数相等）

( 1) g + gt1 = 0.0
chi2( 1) = 0.06
Prob > chi2 = 0.8041
相应的P值＝0.8041>a，差异无统计学意义，故无证据显示两组总体均数不等。
服药6个月时的两组平均体重比较的Stata命令和输出结果如下：
test g＋gt2＝0 （H0：＝0）
( 1) g + gt2 = 0.0
chi2( 1) = 5.55
Prob > chi2 = 0.0184
相应的P值＝0.0184<a，差异有统计学意义，故可以认为服A药6个月时的人群平均体重低于服B药6个月的人群平均体重。

（四）随机效应的Logistic模型
在第十章中介绍了常见的两分类logistic模型为
（3）
或者可以写为
（4）
并由此可知：P越大，logit(P)也越大；P越小，logit(P)也越小，并且两者一一对应。logistic模型的主要适应范畴为：应变量为两分类变量，每个观察记录是独立的。为了使读者能较连贯地理解随机效应的logistic模型，故通过下面一个简单例子复习相关内容。
例2 收集100名HP阳性患者，随机分成两组，分别用两种治疗方案进行治疗（A方案g=1，B方案g=0），2个月后考察HP治疗的疗效（阴性Y=1,仍为阳性y=0），比较两组疗效的差别。
表4 两组100名HP阳性患者治疗结果

阴性(Y=1)
阳性(Y=0)
合计
A组(g=1)
19(38%)
31(62%)
50
B组(g=0)
13(26%)
37(74%)
50
由于本例的应变量为二分类变量，所以可以用Logistic模型对本例资料进行统计分析。本例的Logistic模型为

（5）
对于A组，g=1，logit(PA)＝b0＋b1；对于B组，g=0，logit（PB）＝b0，因此若b1>0，说明logit(PA)> logit（PB），相应的A组转阴率PA高于B组转阴率PB；若b1<0，说明logit(PA)< logit（PB），相应的A组转阴率PA低于B组转阴率PB；若b1＝0，说明两组转阴率相同。我们仍用Stata软件进行参数估计和统计检验。Stata数据格式为（数据文件为ex12-5.dta）
g
y
w
1
1
19
1
0
31
0
1
13
0
0
37
Stata命令： logit y g [freq=w]
结果如下
Logit estimates Number of obs = 100
LR chi2(1) = 1.66
Prob > chi2 = 0.1974
Log 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，因此对于检验水准a＝0.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=1
A组（g=0）
人数
B组(g=1)
人数
阳性(y=0)
阳性(y=0)
阳性(y=0)
1
1
阳性(y=0)
阳性(y=0)
阴性(y=1)
3
4
阳性(y=0)
阴性(y=1)
阳性(y=0)
7
1
阳性(y=0)
阴性(y=1)
阴性(y=1)
20
15
阴性(y=1)
阳性(y=0)
阳性(y=0)
1
1
阴性(y=1)
阳性(y=0)
阴性(y=1)
5
3
阴性(y=1)
阴性(y=1)
阳性(y=0)
5
4
阴性(y=1)
阴性(y=1)
阴性(y=1)
8
21

由于患者转阴后还有一定的概率转阳，因此不能用生存分析的方法进行比较，但本例资料可以认为是重复测量的二分类资料，可以用下列随机效应的Logistic模型：.

（6）
其中u~N(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）
b0＋b3
b0＋b1＋b3＋b4
b0＋b2＋b3＋b5
方案B（g=0）
b0
b0＋b1
b0＋b2
两种方案Log(Odds)的差异
b3
b3＋b4
b3＋b5
Odds Ratio
exp(b3)
exp(b3＋b4)
exp(b3＋b5)
因此两种治疗方案在一个疗程治疗的转阴率比较就是检验H0：b3＝0；在二个疗程治疗的转阴率比较就是检验H0：b3＋b4＝0；在三个疗程治疗的转阴率比较就是检验H0：b3＋b5＝0。
如果两种治疗方案与随访时间无交互作用，即b4＝0且b5＝0，则上述模型可以简化为 （9）式，相应的两种治疗方案与多个疗程的Logit(P)的确定性部分关系见表7。

（9）
表7 两种治疗方案与多个疗程的Logit(P)的确定性部分关系

一个疗程（t1=0,t2=0）
二个疗程（t1=1,t2=0）
三个疗程（t1=0,t2=1）
方案A（g=0）
b0
b0＋b1
b0＋b2
方案B（g=1）
b0＋b3
b0＋b1＋b3
b0＋b2＋b3
两种方案的差异
b3
b3
b3
第二疗程与第一疗程比较：B方案：（b0＋b1＋b3）－（b0＋b3）＝b1 A方案：b0＋b1－b0=b1
第三疗程与第一疗程比较：B方案：（b0＋b2＋b3）－（b0＋b3）＝b2 A方案：b0＋b2－b0=b2
第三疗程与第二疗程比较：B方案：（b0＋b2＋b3）－（b0＋b1＋b3）=b2-b1 A：方案：b0＋b2－(b0+b1)=b2-b1
因此各个疗程的两种治疗方案转阴率比较就是检验H0：b3＝0。
本例数据和Stata数据结构如下
t1
t2
no
g
y
0
0
1
0
0
1
0
1
0
0
0
1
1
0
0
0
0
2
0
0
1
0
2
0
0
0
1
2
0
1
…
…
…
…
…
0
0
100
1
1
1
0
100
1
1
0
1
100
1
1

本例的Stata操作命令
gen gt1=g*t1 产生g×t1项
gen gt1=g*t2 产生g×t2项
xtlogit y t1 t2 g gt1 gt2,i(no) 即：xtlogit y 自变量,i(观察对象编号)
结果如下
Random-effects logit Number of obs = 300
Group variable (i) : no Number of groups = 100
Random effects u_i ~ Gaussian Obs per group: min = 3
avg = 3.0
max = 3
Wald chi2(5) = 35.10
Log 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：b4＝b5＝0 vs H1：b4和b5不全为0，a=0.05,相应的Stata操作命令如下
test gt1=0
test gt2=0,a
结果如下
( 1) [y]gt1 = 0.0
chi2( 1) = 1.09
Prob > chi2 = 0.2967
( 1) [y]gt1 = 0.0
( 2) [y]gt2 = 0.0
chi2( 2) = 1.38
Prob > chi2 = 0.5024
相应的P值为0.5024>a，因此交互项没有统计学意义，因此可以选用（12.9）式的无交互作用logistic模型，相应的Stata命令为
xtlogit y t1 t2 g,i(no)
Random-effects logit Number of obs = 300
Group variable (i) : no Number of groups = 100
Random effects u_i ~ Gaussian Obs per group: min = 3
avg = 3.0
max = 3
Wald chi2(3) = 33.89
Log likelihood = -166.38816 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
t1 | 1.563917 .328479 4.76 0.000 .9201103 2.207724
t2 | 1.436818 .321107 4.47 0.000 .8074599 2.066176
g | .6409473 .2705639 2.37 0.018 .1106518 1.171243
_cons | -.402588 .2447068 -1.65 0.100 -.8822044 .0770284
-------------+----------------------------------------------------------------
/lnsig2u | -14 650.5396 -1289.034 1261.034
-------------+----------------------------------------------------------------
sigma_u | .0009119 .2966076 1.2e-280 6.8e+273
rho | 8.32e-07 .0005409 0 .
------------------------------------------------------------------------------
Likelihood ratio test of rho=0: chibar2(01) = 0.00 Prob >= chibar2 = 1.000
公式（12.9）的参数b3的估计值为0.6409473，相应的，P值为0.018，对于统计检验水准为0.05而言，两组疗效的差别有统计学意义，可以认为B方案的疗效优于A方案；参数b1的估计值为1.563917，相应的，P值<0.001，因此可以认为二个疗程的HP转阴率高于一个疗程，并且差别有统计学意义；参数b2的估计值为1.436818，相应的，相应的P值<0.001，因此可以认为三个疗程的HP转阴率高于一个疗程，并且差别有统计学意义;三个疗程与二个疗程的HP转阴率比较需要作下列检验。
H0：b1＝b2 vs H1：b1¹b2 a＝0.05
用Stata检验命令 test t1=t2
结果如下
( 1) [y]t1 - [y]t2 = 0.0
chi2( 1) = 0.13
Prob > chi2 = 0.7216
P值＝0.7216，二个疗程与三个疗程的转阴率的差异无统计学意义，因此没有足够证据可以推断可以认为三个疗程的HP转阴率与二个疗程的疗效有差别。
结论：可以认为B方案的疗效优于A方案，并且可以认为二个疗程治疗的疗效高于一个疗程治疗的疗效，但没有足够的证据显示有必要进行三个疗程的治疗。