摘　要　目的：对一般的临床随访资料构建Markov模型，预测各观察点的生存人数和死亡人数。方法：将患者在各观察点的状态表示为生存、死亡、失访三类，作为一类特殊的纵向观测分类数据，采用加权最小二乘法估计未知参数。结果：提出了随访资料的Markov模型，并结合实例，给出状态间的转移概率，以及各观察点生存人数、死亡人数的预测值。结论：采用Markov模型对临床随访资料进行预测分析可作为寿命表等分析方法的有效补充。

　　关键词：随访资料　生存分析　Markov模型

　　疾病的治疗情况，一方面看结局好坏，另一方面还要看出现这种结局所经历的时间长短，这类资料一般通过随访收集。常见的随访起点是确诊日期、手术日期等；最明确的阳性结局是死亡，此外还有复发、致残等。随访资料常因失访等原因造成某些数据观察不完全，要用专门的方法进行统计处理，这类方法称为生存分析，常用的有寿命表法、乘积极限法、Cox比例风险模型分析法等^［1,2］。由于随访资料具有时间序列的属性，因此可以用纵向研究^［3,4］的方法分析这类资料。本研究将临床随访资料作为一类特殊的纵向观测分类数据，建立随访资料生存分析的Markov模型，以加权最小二乘法估计未知参数，获得各状态间的转移概率，以及各观察点生存人数、死亡人数的预测值等。

　　1　方法和结果

　　1.1　有穷状态空间Markov过程　下面给出有穷状态空间连续时间Markov过程的有关结果^［5］:状态空间为｛1,2,…,k｝的Markov过程｛X(t):0＜t＜∞｝可由k阶转移概率矩阵列P(s,t)=｛p_ij(s,t)｝(0≤s≤t)确定。p_ij(s,t)为转移概率:

　　p_ij(s,t)=Pr｛X(t)=j｜X(s)=i｝　　　　i,j=1,…,k

　　转移强度矩阵Q(t)=｛q_ij(t)｝。若q_ij(t)=q_ij，即与时间t无关，称为齐次马尔柯夫过程，p_ij(s,t)仅是t-s的函数。记P(t)=P(0,t)，Chapman-Kolmogorov向前方程dp(t)/dt=P(t)Q有唯一解:

　　　　(1)

　　边界条件为P(0)=I（单位阵）。若Q的k个特征根为d₁,…,d_k,记D=diag(d₁,…,d_k)。以A表示相应的特征向量构成的矩阵，则有分解Q=ADA^-1。代入（1）式可得：

　　P(t)=Adiag｛exp(d₁t),…,exp(d_kt)｝A^-1　　　　（2）

　　由（2）式可知，一旦确定了A和d₁，…,d_k,对任意t可确定P(t)。

　　在以下的讨论中假定齐次强度是参数向量θ=(θ₁,…,θ_r)^′的可微函数，即q_ij=q_ij(θ)，强度矩阵为Q=｛q_ij(θ)｝。

　　1.2　随访资料观测数据　假设对一群相互独立的个体分别在时刻t₀＜t₁＜…＜t_m进行观测，每个个体的状态构成一个时间齐次的Markov过程,其状态空间定义为｛1,2,…,k｝。设在时刻t_υ处于状态j的个体数为N_j(t_υ)(j=1,…,k;υ=0,1,…,m)，考虑封闭系统，即 ，假定初始值N_j(t0)(j=1,…,k)已知。令观测时间间隔为u_υ=t_υ-t_υ-1(υ=1,…,m)。记为：

　　N_υ=｛N₁(t_υ),…,N_k(t_υ)｝′　　　　(3)

　　M_υ=｛N₁(t_υ),…,N_k-1(t_υ)｝′　　　　(4)

　　其中υ=0,1,…,m。给定N_υ-1，M_υ的分布为渐近非齐次(k-1)维正态，有：

　　E(M_υ｜N_υ-1)=P^′₁(u_υ)N_υ-1　　　　(5)

　　Cov(M_υ｜N_υ-1)=V_1v=diag｛P^′₁(u_υ)N_υ-1｝-P^′₁(u_υ)diag(N_υ-1)P₁(u_υ)　　　　(6)

　　其中P₁(u_υ)是去掉P（u_υ）的最后一列构成的k×(k-1)阶矩阵，V_1υ是V_υ的(k-1)阶主子阵。

　　1.3　加权最小二乘估计　利用式（5）和（6）可给出｛q_ij(θ)｝中参数θ的估计方法。记：

　　　　(7)

　　则使S达到极小的θ即为加权最小二乘估计，记作 。记Z_υ=M_υ-P^′₁(u_υ)N_υ-1，估计 的迭代公式为^［6］：

　　　　(8)

　　其中C_υ(θ)是(k-1)×r矩阵,其第j列为( / θ_j)P′₁(u_υ-1)N_υ-1(j=1,…,r)。V_1υ、Z_υ均在θ=θ₀处取值。 的渐近协方差阵估计为：

　　　　(9)

　　1.4　随访资料的Markov模型　这里考虑的状态空间为｛1，2，3｝，即每个观察对象在任一观察点可处于三种状态之一:生存、死亡、失访。状态2(死亡)和状态3(失访)为吸收状态，进入到吸收状态的观察对象不可能转移到其他状态。转移强度矩阵为：

　　　　(10)

　　Q的特征根为，d₁=-θ₁-θ₂,d₂=d₃=0，以A表示相应的特征向量构成的矩阵，则Q的分解表达式为：

Q=Adiag(d₁,d₂,d₃)A^-1

　　记x=exp｛-(θ₁+θ₂)t｝，则由（2）式可得如下转移概率矩阵：

　　　　(11)

　　而（6）式中的P₁(t)可表示为：

　　　　(12)

　　（12）式分别关于θ₁,θ₂，求导，然后右乘向量N_υ-1获得相应于迭代公式（8）中的矩阵C_υ(θ)为：

　　其中λ=θ₂(1-x)+θ₁xu_υ-1(θ₁+θ)

　　μ=θ₁xu_υ-1(θ₁+θ₂)-θ₁(1-x)

　　至此完成了随访资料的Markov模型构造。在此基础上我们编制了较为实用的计算机程序，下面以实例说明随访资料的Markov模型分析过程。

　　2　实例分析与讨论

　　表1为585例乳腺癌术后随访资料(括号内为理论值)。

表 1　585例乳腺癌术后随访结果

　　Tab 1　Follow-up results of 585 patientswith breast cancer after operation

　　应用前面叙述的方法，先将数据进行调整。分别以N₁(t_i)、N₂(t_i)、N₃(t_i)表示术后第t_i年观测到的生存人数、死亡人数、失访人数。N₁(t_i)对应于表1中的期初观察人数，而N₂(t_i)、N₃(t_i)分别为表1中病死人数、失访人数按时间顺序的累加值，本例观测时间间隔μ_υ≡1（年）。

　　根据迭代公式（8），确定（10）式中的参数θ₁和θ₂。选初始值θ₁＝θ₂＝0.5，得（θ₁,θ₂）的加权最小二乘估计为： ₁=0.1404, ₂=0.1990。

　　于是由矩阵（11）的第1行可得各观测点的转移概率，计算结果列于表2。

表 2　各观测点的转移概率

　　Tab 2　Transition probability at each observation point

　　表中转移概率均有其实际意义，如p₁₁(4)=0.2573,表示到术后第4年有25.73%的患者处于生存状态；p₁₃(2)=0.2889，表示到术后第2年有28.89%的患者失访，其中包括生存和死亡。根据表2数据计算理论频数，将计算结果填入表1（括号内数据），由表1可见,除去两端个别值外，理论值与实际值符合良好。

　　临床随访资料的分析是应用统计中的一个重要领域。我们应用Markov模型分析随访资料，将观察对象在各观察点的状态分为生存、死亡、失访三个状态，利用加权最小二乘法获得待估参数θ₁和θ₂，进而计算各观测点的转移概率、频数，收到了满意的效果。对于一个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