临床随访资料的Markov模型构建与应用研究
熊林平 曹秀堂 孟岳良 郭祖超
摘 要 目的：对一般的临床随访资料构建Markov模型，预测各观察点的生存人数和死亡人数。方法：将患者在各观察点的状态表示为生存、死亡、失访三类，作为一类特殊的纵向观测分类数据，采用加权最小二乘法估计未知参数。结果：提出了随访资料的Markov模型，并结合实例，给出状态间的转移概率，以及各观察点生存人数、死亡人数的预测值。结论：采用Markov模型对临床随访资料进行预测分析可作为寿命表等分析方法的有效补充。
关键词：随访资料 生存分析 Markov模型
疾病的治疗情况，一方面看结局好坏，另一方面还要看出现这种结局所经历的时间长短，这类资料一般通过随访收集。常见的随访起点是确诊日期、手术日期等；最明确的阳性结局是死亡，此外还有复发、致残等。随访资料常因失访等原因造成某些数据观察不完全，要用专门的方法进行统计处理，这类方法称为生存分析，常用的有寿命表法、乘积极限法、Cox比例风险模型分析法等[1,2]。由于随访资料具有时间序列的属性，因此可以用纵向研究[3,4]的方法分析这类资料。本研究将临床随访资料作为一类特殊的纵向观测分类数据，建立随访资料生存分析的Markov模型，以加权最小二乘法估计未知参数，获得各状态间的转移概率，以及各观察点生存人数、死亡人数的预测值等。
, 百拇医药
1 方法和结果
1.1 有穷状态空间Markov过程 下面给出有穷状态空间连续时间Markov过程的有关结果[5]:状态空间为{1,2,…,k}的Markov过程{X(t):0ij(s,t)}(0≤s≤t)确定。pij(s,t)为转移概率:
pij(s,t)=Pr{X(t)=j｜X(s)=i} i,j=1,…,k
转移强度矩阵Q(t)={qij(t)}。若qij(t)=qij，即与时间t无关，称为齐次马尔柯夫过程，pij(s,t)仅是t-s的函数。记P(t)=P(0,t)，Chapman-Kolmogorov向前方程dp(t)/dt=P(t)Q有唯一解:
(1)
, 百拇医药
边界条件为P(0)=I(单位阵)。若Q的k个特征根为d1,…,dk,记D=diag(d1,…,dk)。以A表示相应的特征向量构成的矩阵，则有分解Q=ADA-1。代入(1)式可得：
P(t)=Adiag{exp(d1t),…,exp(dkt)}A-1 (2)
由(2)式可知，一旦确定了A和d1，…,dk,对任意t可确定P(t)。
在以下的讨论中假定齐次强度是参数向量θ=(θ1,…,θr)′的可微函数，即qij=qij(θ)，强度矩阵为Q={qij(θ)}。
, 百拇医药
1.2 随访资料观测数据 假设对一群相互独立的个体分别在时刻t01<…m进行观测，每个个体的状态构成一个时间齐次的Markov过程,其状态空间定义为{1,2,…,k}。设在时刻tυ处于状态j的个体数为Nj(tυ)(j=1,…,k;υ=0,1,…,m)，考虑封闭系统，即
，假定初始值Nj(t0)(j=1,…,k)已知。令观测时间间隔为uυ=tυ-tυ-1(υ=1,…,m)。记为：
Nυ={N1(tυ),…,Nk(tυ)}′ (3)
, 百拇医药
Mυ={N1(tυ),…,Nk-1(tυ)}′ (4)
其中υ=0,1,…,m。给定Nυ-1，Mυ的分布为渐近非齐次(k-1)维正态，有：
E(Mυ｜Nυ-1)=P′1(uυ)Nυ-1 (5)
Cov(Mυ｜Nυ-1)=V1v=diag{P′1(uυ)Nυ-1}-P′1(uυ)diag(Nυ-1)P1(uυ) (6)
, 百拇医药
其中P1(uυ)是去掉P(uυ)的最后一列构成的k×(k-1)阶矩阵，V1υ是Vυ的(k-1)阶主子阵。
1.3 加权最小二乘估计 利用式(5)和(6)可给出{qij(θ)}中参数θ的估计方法。记：
(7)
则使S达到极小的θ即为加权最小二乘估计，记作
。记Zυ=Mυ-P′1(uυ)Nυ-1，估计的迭代公式为[6]：
(8)
, http://www.100md.com
其中Cυ(θ)是(k-1)×r矩阵,其第j列为(
/θj)P′1(uυ-1)Nυ-1(j=1,…,r)。V1υ、Zυ均在θ=θ0处取值。的渐近协方差阵估计为：
(9)
1.4 随访资料的Markov模型 这里考虑的状态空间为{1，2，3}，即每个观察对象在任一观察点可处于三种状态之一:生存、死亡、失访。状态2(死亡)和状态3(失访)为吸收状态，进入到吸收状态的观察对象不可能转移到其他状态。转移强度矩阵为：
(10)
, 百拇医药
Q的特征根为，d1=-θ1-θ2,d2=d3=0，以A表示相应的特征向量构成的矩阵，则Q的分解表达式为：
Q=Adiag(d1,d2,d3)A-1
记x=exp{-(θ1+θ2)t}，则由(2)式可得如下转移概率矩阵：
(11)
而(6)式中的P1(t)可表示为：
(12)
, 百拇医药
(12)式分别关于θ1,θ2，求导，然后右乘向量Nυ-1获得相应于迭代公式(8)中的矩阵Cυ(θ)为：
其中λ=θ2(1-x)+θ1xuυ-1(θ1+θ)
μ=θ1xuυ-1(θ1+θ2)-θ1(1-x)
至此完成了随访资料的Markov模型构造。在此基础上我们编制了较为实用的计算机程序，下面以实例说明随访资料的Markov模型分析过程。
, http://www.100md.com
2 实例分析与讨论
表1为585例乳腺癌术后随访资料(括号内为理论值)。
表 1 585例乳腺癌术后随访结果
Tab 1 Follow-up results of 585 patients with breast cancer after operation
Yearsafter
operation
No.of
beginning
No.of
death
, 百拇医药
No.of
censor
[0,1)
585(585)
59(69)
59(99)
[1,2)
463(417)
69(50)
71(70)
[2,3)
323(297)
, http://www.100md.com
43(36)
55(50)
[3,4)
225(211)
30(25)
38(36)
[4,5)
157(151)
13(18)
31(25)
[5,6)
113(107)
, 百拇医药
7(13)
26(18)
[6,7)
80(76)
14(10)
21(13)
[7,8)
45(54)
4(7)
11(9)
[8,9)
30(38)
3(5)
, http://www.100md.com
15(7)
[9,10)
12(27)
0(4)
12(5)
应用前面叙述的方法，先将数据进行调整。分别以N1(ti)、N2(ti)、N3(ti)表示术后第ti年观测到的生存人数、死亡人数、失访人数。N1(ti)对应于表1中的期初观察人数，而N2(ti)、N3(ti)分别为表1中病死人数、失访人数按时间顺序的累加值，本例观测时间间隔μυ≡1(年)。
, 百拇医药
根据迭代公式(8)，确定(10)式中的参数θ1和θ2。选初始值θ1＝θ2＝0.5，得(θ1,θ2)的加权最小二乘估计为：1=0.1404,2=0.1990。
于是由矩阵(11)的第1行可得各观测点的转移概率，计算结果列于表2。
表 2 各观测点的转移概率
Tab 2 Transition probability at each observation point
, http://www.100md.com
t(year)
p11(t)
p12(t)
p13(t)
1
0.712 2
0.119 1
0.168 7
2
0.507 2
0.203 9
0.288 9
, 百拇医药
3
0.361 2
0.264 3
0.374 5
4
0.257 3
0.307 3
0.435 3
5
0.183 2
0.337 9
0.478 9
6
, http://www.100md.com
0.130 5
0.359 7
0.509 8
7
0.092 9
0.375 3
0.531 8
8
0.066 2
0.386 3
0.547 5
9
0.047 1
, 百拇医药
0.394 2
0.558 7
10
0.033 6
0.399 8
0.566 6
表中转移概率均有其实际意义，如p11(4)=0.2573,表示到术后第4年有25.73%的患者处于生存状态；p13(2)=0.2889，表示到术后第2年有28.89%的患者失访，其中包括生存和死亡。根据表2数据计算理论频数，将计算结果填入表1(括号内数据)，由表1可见,除去两端个别值外，理论值与实际值符合良好。
临床随访资料的分析是应用统计中的一个重要领域。我们应用Markov模型分析随访资料，将观察对象在各观察点的状态分为生存、死亡、失访三个状态，利用加权最小二乘法获得待估参数θ1和θ2，进而计算各观测点的转移概率、频数，收到了满意的效果。对于一个Markov过程，在给定过去状态和现在状态的条件下，将来时点状态的分布只依赖于当前状态，而与过去状态无关。因此在应用Markov模型进行实际分析时，可以不考虑初始时刻前的状态。在实际的医学研究中，常常由于某些原因使得研究者提前终止随访观测，在这种情况下，就可以采用本研究获得的结果，根据现有数据对观察对象的未来结局进行预测。然后利用统计软件，如SAS(statistics analysis system)[7]作分析处理。因此,本研究的结果可作为生存分析的有益补充。
, 百拇医药
作者简介：熊林平，男，1962年12月生，博士，讲师
作者单位 熊林平 孟岳良 第二军医大学卫生勤务学系卫生统计学教研室，上海，200433；
曹秀堂 郭祖超 第四军医大学卫生统计学教研室
参考文献
1 McCullagh P, Nelder JA. Generalized linear models. 2nd ed. London: Chapman and Hall, 1989.419～430
2 方积乾，徐勇勇，余松林，等.医学统计学与电脑实验.上海：上海科学技术出版社， 1997.332～359
3 Follman D. Modelling transitional and joint marginal distributions in re peated categorical data. Statist Med, 1994，13(5):467
, 百拇医药
4 Hendricks SA, Wassell JT, Collins JW, et al. Power determination for geo graphically clustered data using generalized estimating equations. Statist M ed, 1996, 15(9):1951
5 Ross SM. Stochastic processes. New York: John Wiley and Sons, 1983.100～140
6 Kalbfleisch JD, Lawless JF, Vollmer WM. Estimation in Markov models from aggregate data. Biometrics, 1983, 39(4):907
7 胡良平.现代统计学与SAS应用.北京：军事医学科学出版社，1996.295～315, http://www.100md.com
摘 要 目的：对一般的临床随访资料构建Markov模型，预测各观察点的生存人数和死亡人数。方法：将患者在各观察点的状态表示为生存、死亡、失访三类，作为一类特殊的纵向观测分类数据，采用加权最小二乘法估计未知参数。结果：提出了随访资料的Markov模型，并结合实例，给出状态间的转移概率，以及各观察点生存人数、死亡人数的预测值。结论：采用Markov模型对临床随访资料进行预测分析可作为寿命表等分析方法的有效补充。
关键词：随访资料 生存分析 Markov模型
疾病的治疗情况，一方面看结局好坏，另一方面还要看出现这种结局所经历的时间长短，这类资料一般通过随访收集。常见的随访起点是确诊日期、手术日期等；最明确的阳性结局是死亡，此外还有复发、致残等。随访资料常因失访等原因造成某些数据观察不完全，要用专门的方法进行统计处理，这类方法称为生存分析，常用的有寿命表法、乘积极限法、Cox比例风险模型分析法等[1,2]。由于随访资料具有时间序列的属性，因此可以用纵向研究[3,4]的方法分析这类资料。本研究将临床随访资料作为一类特殊的纵向观测分类数据，建立随访资料生存分析的Markov模型，以加权最小二乘法估计未知参数，获得各状态间的转移概率，以及各观察点生存人数、死亡人数的预测值等。
, 百拇医药
1 方法和结果
1.1 有穷状态空间Markov过程 下面给出有穷状态空间连续时间Markov过程的有关结果[5]:状态空间为{1,2,…,k}的Markov过程{X(t):0
pij(s,t)=Pr{X(t)=j｜X(s)=i} i,j=1,…,k
转移强度矩阵Q(t)={qij(t)}。若qij(t)=qij，即与时间t无关，称为齐次马尔柯夫过程，pij(s,t)仅是t-s的函数。记P(t)=P(0,t)，Chapman-Kolmogorov向前方程dp(t)/dt=P(t)Q有唯一解:
(1), 百拇医药
边界条件为P(0)=I(单位阵)。若Q的k个特征根为d1,…,dk,记D=diag(d1,…,dk)。以A表示相应的特征向量构成的矩阵，则有分解Q=ADA-1。代入(1)式可得：
P(t)=Adiag{exp(d1t),…,exp(dkt)}A-1 (2)
由(2)式可知，一旦确定了A和d1，…,dk,对任意t可确定P(t)。
在以下的讨论中假定齐次强度是参数向量θ=(θ1,…,θr)′的可微函数，即qij=qij(θ)，强度矩阵为Q={qij(θ)}。
, 百拇医药
1.2 随访资料观测数据 假设对一群相互独立的个体分别在时刻t0
，假定初始值Nj(t0)(j=1,…,k)已知。令观测时间间隔为uυ=tυ-tυ-1(υ=1,…,m)。记为：Nυ={N1(tυ),…,Nk(tυ)}′ (3)
, 百拇医药
Mυ={N1(tυ),…,Nk-1(tυ)}′ (4)
其中υ=0,1,…,m。给定Nυ-1，Mυ的分布为渐近非齐次(k-1)维正态，有：
E(Mυ｜Nυ-1)=P′1(uυ)Nυ-1 (5)
Cov(Mυ｜Nυ-1)=V1v=diag{P′1(uυ)Nυ-1}-P′1(uυ)diag(Nυ-1)P1(uυ) (6)
, 百拇医药
其中P1(uυ)是去掉P(uυ)的最后一列构成的k×(k-1)阶矩阵，V1υ是Vυ的(k-1)阶主子阵。
1.3 加权最小二乘估计 利用式(5)和(6)可给出{qij(θ)}中参数θ的估计方法。记：
(7)则使S达到极小的θ即为加权最小二乘估计，记作
。记Zυ=Mυ-P′1(uυ)Nυ-1，估计的迭代公式为[6]： (8), http://www.100md.com
其中Cυ(θ)是(k-1)×r矩阵,其第j列为(
/θj)P′1(uυ-1)Nυ-1(j=1,…,r)。V1υ、Zυ均在θ=θ0处取值。的渐近协方差阵估计为： (9)1.4 随访资料的Markov模型 这里考虑的状态空间为{1，2，3}，即每个观察对象在任一观察点可处于三种状态之一:生存、死亡、失访。状态2(死亡)和状态3(失访)为吸收状态，进入到吸收状态的观察对象不可能转移到其他状态。转移强度矩阵为：
(10), 百拇医药
Q的特征根为，d1=-θ1-θ2,d2=d3=0，以A表示相应的特征向量构成的矩阵，则Q的分解表达式为：
Q=Adiag(d1,d2,d3)A-1
记x=exp{-(θ1+θ2)t}，则由(2)式可得如下转移概率矩阵：
(11)而(6)式中的P1(t)可表示为：
(12), 百拇医药
(12)式分别关于θ1,θ2，求导，然后右乘向量Nυ-1获得相应于迭代公式(8)中的矩阵Cυ(θ)为：
其中λ=θ2(1-x)+θ1xuυ-1(θ1+θ)
μ=θ1xuυ-1(θ1+θ2)-θ1(1-x)
至此完成了随访资料的Markov模型构造。在此基础上我们编制了较为实用的计算机程序，下面以实例说明随访资料的Markov模型分析过程。
, http://www.100md.com
2 实例分析与讨论
表1为585例乳腺癌术后随访资料(括号内为理论值)。
表 1 585例乳腺癌术后随访结果
Tab 1 Follow-up results of 585 patients with breast cancer after operation
Yearsafter
operation
No.of
beginning
No.of
death
, 百拇医药
No.of
censor
[0,1)
585(585)
59(69)
59(99)
[1,2)
463(417)
69(50)
71(70)
[2,3)
323(297)
, http://www.100md.com
43(36)
55(50)
[3,4)
225(211)
30(25)
38(36)
[4,5)
157(151)
13(18)
31(25)
[5,6)
113(107)
, 百拇医药
7(13)
26(18)
[6,7)
80(76)
14(10)
21(13)
[7,8)
45(54)
4(7)
11(9)
[8,9)
30(38)
3(5)
, http://www.100md.com
15(7)
[9,10)
12(27)
0(4)
12(5)
应用前面叙述的方法，先将数据进行调整。分别以N1(ti)、N2(ti)、N3(ti)表示术后第ti年观测到的生存人数、死亡人数、失访人数。N1(ti)对应于表1中的期初观察人数，而N2(ti)、N3(ti)分别为表1中病死人数、失访人数按时间顺序的累加值，本例观测时间间隔μυ≡1(年)。
, 百拇医药
根据迭代公式(8)，确定(10)式中的参数θ1和θ2。选初始值θ1＝θ2＝0.5，得(θ1,θ2)的加权最小二乘估计为：
1=0.1404,2=0.1990。于是由矩阵(11)的第1行可得各观测点的转移概率，计算结果列于表2。
表 2 各观测点的转移概率
Tab 2 Transition probability at each observation point
, http://www.100md.com
t(year)
p11(t)
p12(t)
p13(t)
1
0.712 2
0.119 1
0.168 7
2
0.507 2
0.203 9
0.288 9
, 百拇医药
3
0.361 2
0.264 3
0.374 5
4
0.257 3
0.307 3
0.435 3
5
0.183 2
0.337 9
0.478 9
6
, http://www.100md.com
0.130 5
0.359 7
0.509 8
7
0.092 9
0.375 3
0.531 8
8
0.066 2
0.386 3
0.547 5
9
0.047 1
, 百拇医药
0.394 2
0.558 7
10
0.033 6
0.399 8
0.566 6
表中转移概率均有其实际意义，如p11(4)=0.2573,表示到术后第4年有25.73%的患者处于生存状态；p13(2)=0.2889，表示到术后第2年有28.89%的患者失访，其中包括生存和死亡。根据表2数据计算理论频数，将计算结果填入表1(括号内数据)，由表1可见,除去两端个别值外，理论值与实际值符合良好。
临床随访资料的分析是应用统计中的一个重要领域。我们应用Markov模型分析随访资料，将观察对象在各观察点的状态分为生存、死亡、失访三个状态，利用加权最小二乘法获得待估参数θ1和θ2，进而计算各观测点的转移概率、频数，收到了满意的效果。对于一个Markov过程，在给定过去状态和现在状态的条件下，将来时点状态的分布只依赖于当前状态，而与过去状态无关。因此在应用Markov模型进行实际分析时，可以不考虑初始时刻前的状态。在实际的医学研究中，常常由于某些原因使得研究者提前终止随访观测，在这种情况下，就可以采用本研究获得的结果，根据现有数据对观察对象的未来结局进行预测。然后利用统计软件，如SAS(statistics analysis system)[7]作分析处理。因此,本研究的结果可作为生存分析的有益补充。
, 百拇医药
作者简介：熊林平，男，1962年12月生，博士，讲师
作者单位 熊林平 孟岳良 第二军医大学卫生勤务学系卫生统计学教研室，上海，200433；
曹秀堂 郭祖超 第四军医大学卫生统计学教研室
参考文献
1 McCullagh P, Nelder JA. Generalized linear models. 2nd ed. London: Chapman and Hall, 1989.419～430
2 方积乾，徐勇勇，余松林，等.医学统计学与电脑实验.上海：上海科学技术出版社， 1997.332～359
3 Follman D. Modelling transitional and joint marginal distributions in re peated categorical data. Statist Med, 1994，13(5):467
, 百拇医药
4 Hendricks SA, Wassell JT, Collins JW, et al. Power determination for geo graphically clustered data using generalized estimating equations. Statist M ed, 1996, 15(9):1951
5 Ross SM. Stochastic processes. New York: John Wiley and Sons, 1983.100～140
6 Kalbfleisch JD, Lawless JF, Vollmer WM. Estimation in Markov models from aggregate data. Biometrics, 1983, 39(4):907
7 胡良平.现代统计学与SAS应用.北京：军事医学科学出版社，1996.295～315, http://www.100md.com