纵向观测二分类数据的广义线性模型分析
熊林平 曹秀堂 徐勇勇 陆健
摘 要 目的：利用广义线性模型对纵向观测二分类数据进行分析，充分考虑纵向观测间的相关性，给出一般分析方法。方法：采用Zeger和提出的广义估计方程，拟合logistic广义线性模型，讨论3种协方差矩阵结构。结果：同时获得回归参数、相关参数的估计，完成了较为实用的运行程序，并进行了实例分析。结论：医学研究和临床试验中经常接触到纵向观测数据，对这类数据需采用特殊的方法进行分析处理，以解决重复观测间的相关性问题。
关键词：纵向观测 二分类数据 广义线性模型
医学研究和临床试验中经常碰到重复观测二分类数据。如:进行临床试验, 比较两种处理对呼吸道疾病的疗效。在两个研究中心将111名患者(甲中心56名，乙中心55名)随机分为两个处理组：积极治疗组(54名)和安慰剂治疗组(57名)，治疗期间连续4次检查患者的呼吸道情况，结果为二分类变量，以0表示差，1表示好。可能的影响因素有研究中心、处理、性别、基准呼吸状况以及开始治疗时患者的年龄。数据如表1示，表中处理栏A为积极治疗，P为安慰治疗。 这种资料是按时间顺序对个体进行重复观测而得到的，其观测结果分为两类。对于这种纵向观测资料，若采用通常意义下的t 检验或方差分析
, 百拇医药
表 1 甲中心56名患者及乙中心55名患者呼吸道状况观测结果
Tab 1 The observations of breathing condition of 56 patients in center A and 55 patients in center B
No.
Tr.
Sex
Age
Status(0=bad,1=good)
No.
Tr.
Sex
, http://www.100md.com
Age
Status(0=bad,1=good)
Baseline
1
2
3
4
Baseline
1
2
3
4
1
, http://www.100md.com
P
M
46
0
0
0
0
0
57
P
F
39
0
0
, 百拇医药
0
0
0
2
P
M
28
0
0
0
0
0
58
A
, 百拇医药
M
25
0
0
1
1
1
3
A
M
23
1
1
1
, http://www.100md.com
1
1
59
A
M
58
1
1
1
1
1………………………………………………
55
P
, 百拇医药 M
24
1
1
1
1
1
110
A
F
63
1
1
1
, 百拇医药
1
1
56
A
M
25
0
1
1
0
1
111
A
M
, 百拇医药
31
1
1
1
1
1
等横向分析方法，将造成信息损失；更为重要的是，由于忽略了重复数据间可能存在的相关性，会导致参数估计不准确,甚至得出错误的结论。
纵向研究由于充分考虑了数据间的相互依赖性，具有可提高功效以及对模型的选择有稳健性等优点，正日益受到统计学界的广泛关注。纵向数据的分析目前多采用拟合线性模型和广义线性模型的方法，研究的焦点集中在如何解决重复观测间的相关性问题上[1,2]。Zeger等在1986年提出的广义估计方程，较好地解决了似然函数中多余参数的估计问题[3]，它将横向研究中的线性模型以及广义线性模型的参数估计推广到了纵向数据的分析。本文讨论纵向观测二分类数据的logistic广义线性模型分析，利用广义估计方程解决模型的回归参数和相关参数的估计问题。
, 百拇医药
1 模型及参数估计
一般地，设有m个个体，个体i的观测序列为yi=(yi1,…,yini)′，相应的观测时间为ti=(ti1,…,tini)′，对应于yij(j=1,…,ni)的协变量向量为xij=(xij1,…,xijp)′，以Yij表示对应于yij的随机变量，记μij=E(Yij)，μi=(μi1,…,μini)′。为方便起见，下面的叙述中，去掉ni的下标i。
Yij为二分类变量，观测结果为0或1，μij=Pr{Yij=1}，Var(Yij)=μij(1-μij)，二分类数据的logistic模型为：
, 百拇医药
logit(μij)=x′ijβ (i=1,…,m;j=1,…,n) (1)
即:
log{μij/(1-μij)}=xij1β1+…+xijpβp (2)
这里β=(β1,…,βp)′为p维回归参数向量。
回归参数β的估计值
为下列广义估计方程的解[4]：
(3)
, 百拇医药
其中Di=
/
=(djk)n×p，Vi为近似协方差阵，定义为：
Vi=A1/2iRi(α)A1/2i
Ai=diag{μi1(1-μi1),…,μin(1-μin)}
, http://www.100md.com
Ri(α)为近似相关阵，其结构完全由相关参数α确定。
根据(2)式可得：
μij=exp(x′ijβ)/{1+exp(x′ijβ)}
则Di的元素djk为：
j=1,…,n;k=1,…,p
至于相关参数α的估计，根据样本相关系数
, http://www.100md.com
记其均值E(Rijk=δijk)，Rijk的方差为[5]：
wijk=Var(Rijk)=1+(1-2μij)(1-2μik).
{μij(1-μij)μik(1-μik)}1/2δijk-δ2ijk
j
将δijk看作Ri(α)的(j,k)元，它是α的函数，α满足如下方程：
(4)
, http://www.100md.com
其中Ri=(Ri12,…,Ri1n,…,Rin-1,n)′为n(n-1)/2维向量，δi=E(Ri)=(δi12,…,δi1n,…,δin-1,n)′
Ei=
Wi=diag{wi12,…,wi1n,…,win-1,n}
于是可得到估计α,β的迭代公式为：
βs+1=βs+(∑iD′iV-1iDi)-1∑iD′iV-1i(yi-μi)
, 百拇医药
αs+1=αs+(∑iE′iW-1iEi)-1∑iE′iW-1i(Ri-δi)
至此，只要知道相关阵Ri(α)，就可迭代得到参数的估计。以下就Ri(α)的几种常用结构进行讨论。
2 相关阵结构
2.1 独立结构(IND) R(α)=I(单位阵),此时每个个体的各次重复观测间相互独立，可用一般的logistic回归模型进行分析[6,7]。
2.2 一阶依赖结构(DE-1) R(α)为三对角阵，次对角线元素为：
, 百拇医药
δij,j+1=α(j=1,…,n-1)
此时 Ei=δi/α=(ei12,…,ei1n,…,ein-1,n)′=(1,0,…,1)′
即关于Ei的元素有：
2.3 可交换结构(EXC) δijk=α (ji=(1,…,1)′为n(n-1)/2维元素均为1的向量。Wi的对角元为：
wijk=1+(1-2μij)(1-2μik)*
, 百拇医药
{μij(1-μij)μik(1-μik)}α-α2 (j
若已经获得估计,
，则可据此得到的协方差阵的稳健估计。我们已将上述方法编制成计算机运行程序，下面以表1的资料进行分析讨论。
3 实例分析与讨论
对表1资料进行分析。以yij表示第i个患者第j次观测时的呼吸道状况，i=1,…,111,j=1,…,4,μij表示yij的均值，协变量为：研究中心xi1、处理xi2、性别xi3、年龄xi4，以及基准呼吸状况xi5。令：
, http://www.100md.com
利用logit连接函数
h(x)=log{x/(1-x)} 针对两种相关阵结构：一阶依赖结构(DE-1)，可交换结构(EXC)，拟合logistic模型。
logi t(μij)=β0+xi1β1+xi2β2+xi3β3+xi4β4+xi5β5
i=1,…,111;j=1,…,4
由估计结果(表2)可见，性别效应(β3)及年龄效应(β4)均无显著性意义。去掉这两个因素，拟合新的logistic模型。最后得到的模型含有参数β2(处理效应)和β5(基准效应)，均有显著性意义，结果见表3。
, http://www.100md.com
表 2 logistic模型Ⅰ拟合结果
Tab 2 Fitted result of logistic model 1
Parameter
Correlated
structure
Estimate
Standard
error
U
β0
DE-1
, http://www.100md.com
-0.906 2
0.456 6
1.984 7
EXC
-0.856 1
0.456 4
1.875 9
β1
DE-1
0.713 0
0.351 9
2.026 1
, 百拇医药
EXC
0.649 5
0.353 2
1.838 8
β2
DE-1
1.206 7
0.347 2
3.475 9
EXC
1.265 4
0.346 7
, http://www.100md.com
3.649 9
β3
DE-1
0.135 1
0.445 4
0.303 2
EXC
0.136 8
0.440 2
0.310 7
β4
DE-1
, http://www.100md.com
-0.017 7
0.012 9
1.374 2
EXC
-0.018 8
0.013 0
1.446 7
β5
DE-1
1.865 4
0.345 6
5.397 5
, http://www.100md.com
EXC
1.845 7
0.346 0
5.334 8
α
DE-1
0.380
EXC
0.329
由表3参数估计结果，可得结论如下：
(1) 基准呼吸状况比较
基准为差(xi5=0)
, http://www.100md.com
基准为好
5≈2，于是基准呼吸状况好的患者的概率比值Pr(yij=1)/Pr(yij=0)为基准差的患者的exp(2)=7.39倍，即指标优势比OR=7.39。表 3 logistic模型拟合结果
Tab 3 Fitted result of logistic model
Parameter
Correlated
structure
, http://www.100md.com
Estimate
Standard
error
U
β0
DE-1
-1.170 0
0.301
3.887 5
EXC
-1.182 7
0.297 5
, 百拇医药
3.976 1
β2
DE-1
1.183 9
0.329 9
3.588 3
EXC
1.246 5
0.328 7
3.792 9
β5
DE-1
, 百拇医药
2.019 2
0.322 1
6.268 2
EXC
1.989 8
0.322
6.179 7
α
DE-1
0.390
EXC
0.352
(2)处理组比较
, http://www.100md.com
安慰剂组(xi2=0)
积极治疗组(xi2=1)
2≈1.2,积极治疗组的Pr(yij=1)/Pr(yij=0)比值为安慰剂组的exp(1.2)=3.32倍，即指标优势比OR=3.32。
参 考 文 献
1 Hendricks SA, Wassell JT, Collins JW, et al. Power determination for geographically clustered data using generalized estima-ting equations. Statist Med, 1996, 15(9): 1951
, 百拇医药
2 Follman D. Modelling transitional and joint marginal distributions in repeated categorical data. Statist Med, 1994, 13(5):467
3 Zeger SL，Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics, 1986, 42(1):121
4 McCullagh P , Nelder JA. Generalized linear models. 2nd ed. London: Chapman and Hall, 1989.419～430
5 Prentice RL. Correlated binary regression with covariates specific to each binary observation.Biometrics, 1988,44(4):1033
6 方积乾，徐勇勇，余松林，等. 医学统计学与电脑实验.上海：上海科学技术出版社，1997.322～331
7 胡良平. 现代统计学与SAS应用.北京:军事医学科学出版社,1996.199～214, 百拇医药
摘 要 目的：利用广义线性模型对纵向观测二分类数据进行分析，充分考虑纵向观测间的相关性，给出一般分析方法。方法：采用Zeger和提出的广义估计方程，拟合logistic广义线性模型，讨论3种协方差矩阵结构。结果：同时获得回归参数、相关参数的估计，完成了较为实用的运行程序，并进行了实例分析。结论：医学研究和临床试验中经常接触到纵向观测数据，对这类数据需采用特殊的方法进行分析处理，以解决重复观测间的相关性问题。
关键词：纵向观测 二分类数据 广义线性模型
医学研究和临床试验中经常碰到重复观测二分类数据。如:进行临床试验, 比较两种处理对呼吸道疾病的疗效。在两个研究中心将111名患者(甲中心56名，乙中心55名)随机分为两个处理组：积极治疗组(54名)和安慰剂治疗组(57名)，治疗期间连续4次检查患者的呼吸道情况，结果为二分类变量，以0表示差，1表示好。可能的影响因素有研究中心、处理、性别、基准呼吸状况以及开始治疗时患者的年龄。数据如表1示，表中处理栏A为积极治疗，P为安慰治疗。 这种资料是按时间顺序对个体进行重复观测而得到的，其观测结果分为两类。对于这种纵向观测资料，若采用通常意义下的t 检验或方差分析
, 百拇医药
表 1 甲中心56名患者及乙中心55名患者呼吸道状况观测结果
Tab 1 The observations of breathing condition of 56 patients in center A and 55 patients in center B
No.
Tr.
Sex
Age
Status(0=bad,1=good)
No.
Tr.
Sex
, http://www.100md.com
Age
Status(0=bad,1=good)
Baseline
1
2
3
4
Baseline
1
2
3
4
1
, http://www.100md.com
P
M
46
0
0
0
0
0
57
P
F
39
0
0
, 百拇医药
0
0
0
2
P
M
28
0
0
0
0
0
58
A
, 百拇医药
M
25
0
0
1
1
1
3
A
M
23
1
1
1
, http://www.100md.com
1
1
59
A
M
58
1
1
1
1
1………………………………………………
55
P
, 百拇医药 M
24
1
1
1
1
1
110
A
F
63
1
1
1
, 百拇医药
1
1
56
A
M
25
0
1
1
0
1
111
A
M
, 百拇医药
31
1
1
1
1
1
等横向分析方法，将造成信息损失；更为重要的是，由于忽略了重复数据间可能存在的相关性，会导致参数估计不准确,甚至得出错误的结论。
纵向研究由于充分考虑了数据间的相互依赖性，具有可提高功效以及对模型的选择有稳健性等优点，正日益受到统计学界的广泛关注。纵向数据的分析目前多采用拟合线性模型和广义线性模型的方法，研究的焦点集中在如何解决重复观测间的相关性问题上[1,2]。Zeger等在1986年提出的广义估计方程，较好地解决了似然函数中多余参数的估计问题[3]，它将横向研究中的线性模型以及广义线性模型的参数估计推广到了纵向数据的分析。本文讨论纵向观测二分类数据的logistic广义线性模型分析，利用广义估计方程解决模型的回归参数和相关参数的估计问题。
, 百拇医药
1 模型及参数估计
一般地，设有m个个体，个体i的观测序列为yi=(yi1,…,yini)′，相应的观测时间为ti=(ti1,…,tini)′，对应于yij(j=1,…,ni)的协变量向量为xij=(xij1,…,xijp)′，以Yij表示对应于yij的随机变量，记μij=E(Yij)，μi=(μi1,…,μini)′。为方便起见，下面的叙述中，去掉ni的下标i。
Yij为二分类变量，观测结果为0或1，μij=Pr{Yij=1}，Var(Yij)=μij(1-μij)，二分类数据的logistic模型为：
, 百拇医药
logit(μij)=x′ijβ (i=1,…,m;j=1,…,n) (1)
即:
log{μij/(1-μij)}=xij1β1+…+xijpβp (2)
这里β=(β1,…,βp)′为p维回归参数向量。
回归参数β的估计值
为下列广义估计方程的解[4]： (3), 百拇医药
其中Di=
/=(djk)n×p，Vi为近似协方差阵，定义为：Vi=A1/2iRi(α)A1/2i
Ai=diag{μi1(1-μi1),…,μin(1-μin)}
, http://www.100md.com
Ri(α)为近似相关阵，其结构完全由相关参数α确定。
根据(2)式可得：
μij=exp(x′ijβ)/{1+exp(x′ijβ)}
则Di的元素djk为：
j=1,…,n;k=1,…,p
至于相关参数α的估计，根据样本相关系数
, http://www.100md.com
记其均值E(Rijk=δijk)，Rijk的方差为[5]：
wijk=Var(Rijk)=1+(1-2μij)(1-2μik).
{μij(1-μij)μik(1-μik)}1/2δijk-δ2ijk
j
将δijk看作Ri(α)的(j,k)元，它是α的函数，α满足如下方程：
(4), http://www.100md.com
其中Ri=(Ri12,…,Ri1n,…,Rin-1,n)′为n(n-1)/2维向量，δi=E(Ri)=(δi12,…,δi1n,…,δin-1,n)′
Ei=
Wi=diag{wi12,…,wi1n,…,win-1,n}
于是可得到估计α,β的迭代公式为：
βs+1=βs+(∑iD′iV-1iDi)-1∑iD′iV-1i(yi-μi)
, 百拇医药
αs+1=αs+(∑iE′iW-1iEi)-1∑iE′iW-1i(Ri-δi)
至此，只要知道相关阵Ri(α)，就可迭代得到参数的估计。以下就Ri(α)的几种常用结构进行讨论。
2 相关阵结构
2.1 独立结构(IND) R(α)=I(单位阵),此时每个个体的各次重复观测间相互独立，可用一般的logistic回归模型进行分析[6,7]。
2.2 一阶依赖结构(DE-1) R(α)为三对角阵，次对角线元素为：
, 百拇医药
δij,j+1=α(j=1,…,n-1)
此时 Ei=δi/α=(ei12,…,ei1n,…,ein-1,n)′=(1,0,…,1)′
即关于Ei的元素有：
2.3 可交换结构(EXC) δijk=α (j
wijk=1+(1-2μij)(1-2μik)*
, 百拇医药
{μij(1-μij)μik(1-μik)}α-α2 (j
若已经获得估计
,，则可据此得到的协方差阵的稳健估计。我们已将上述方法编制成计算机运行程序，下面以表1的资料进行分析讨论。3 实例分析与讨论
对表1资料进行分析。以yij表示第i个患者第j次观测时的呼吸道状况，i=1,…,111,j=1,…,4,μij表示yij的均值，协变量为：研究中心xi1、处理xi2、性别xi3、年龄xi4，以及基准呼吸状况xi5。令：
, http://www.100md.com
利用logit连接函数
h(x)=log{x/(1-x)} 针对两种相关阵结构：一阶依赖结构(DE-1)，可交换结构(EXC)，拟合logistic模型。
logi t(μij)=β0+xi1β1+xi2β2+xi3β3+xi4β4+xi5β5
i=1,…,111;j=1,…,4
由估计结果(表2)可见，性别效应(β3)及年龄效应(β4)均无显著性意义。去掉这两个因素，拟合新的logistic模型。最后得到的模型含有参数β2(处理效应)和β5(基准效应)，均有显著性意义，结果见表3。
, http://www.100md.com
表 2 logistic模型Ⅰ拟合结果
Tab 2 Fitted result of logistic model 1
Parameter
Correlated
structure
Estimate
Standard
error
U
β0
DE-1
, http://www.100md.com
-0.906 2
0.456 6
1.984 7
EXC
-0.856 1
0.456 4
1.875 9
β1
DE-1
0.713 0
0.351 9
2.026 1
, 百拇医药
EXC
0.649 5
0.353 2
1.838 8
β2
DE-1
1.206 7
0.347 2
3.475 9
EXC
1.265 4
0.346 7
, http://www.100md.com
3.649 9
β3
DE-1
0.135 1
0.445 4
0.303 2
EXC
0.136 8
0.440 2
0.310 7
β4
DE-1
, http://www.100md.com
-0.017 7
0.012 9
1.374 2
EXC
-0.018 8
0.013 0
1.446 7
β5
DE-1
1.865 4
0.345 6
5.397 5
, http://www.100md.com
EXC
1.845 7
0.346 0
5.334 8
α
DE-1
0.380
EXC
0.329
由表3参数估计结果，可得结论如下：
(1) 基准呼吸状况比较
基准为差(xi5=0)
, http://www.100md.com
基准为好
5≈2，于是基准呼吸状况好的患者的概率比值Pr(yij=1)/Pr(yij=0)为基准差的患者的exp(2)=7.39倍，即指标优势比OR=7.39。表 3 logistic模型拟合结果Tab 3 Fitted result of logistic model
Parameter
Correlated
structure
, http://www.100md.com
Estimate
Standard
error
U
β0
DE-1
-1.170 0
0.301
3.887 5
EXC
-1.182 7
0.297 5
, 百拇医药
3.976 1
β2
DE-1
1.183 9
0.329 9
3.588 3
EXC
1.246 5
0.328 7
3.792 9
β5
DE-1
, 百拇医药
2.019 2
0.322 1
6.268 2
EXC
1.989 8
0.322
6.179 7
α
DE-1
0.390
EXC
0.352
(2)处理组比较
, http://www.100md.com
安慰剂组(xi2=0)
积极治疗组(xi2=1)
2≈1.2,积极治疗组的Pr(yij=1)/Pr(yij=0)比值为安慰剂组的exp(1.2)=3.32倍，即指标优势比OR=3.32。参 考 文 献
1 Hendricks SA, Wassell JT, Collins JW, et al. Power determination for geographically clustered data using generalized estima-ting equations. Statist Med, 1996, 15(9): 1951
, 百拇医药
2 Follman D. Modelling transitional and joint marginal distributions in repeated categorical data. Statist Med, 1994, 13(5):467
3 Zeger SL，Liang KY. Longitudinal data analysis for discrete and continuous outcomes. Biometrics, 1986, 42(1):121
4 McCullagh P , Nelder JA. Generalized linear models. 2nd ed. London: Chapman and Hall, 1989.419～430
5 Prentice RL. Correlated binary regression with covariates specific to each binary observation.Biometrics, 1988,44(4):1033
6 方积乾，徐勇勇，余松林，等. 医学统计学与电脑实验.上海：上海科学技术出版社，1997.322～331
7 胡良平. 现代统计学与SAS应用.北京:军事医学科学出版社,1996.199～214, 百拇医药