临床随访资料的Markov模型构建与应用研究

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

    Construction and application of Markov model for follow-up data
    Xiong Linping，Cao Xiutang，Meng Yueliang，Guo Zuchao (Department of Health Statistics，Faculty of Medical Service, Second Military Medical University，Shanghai，200433)

    ABSTRACT　Objective: To construct Markov model for clinical follow-up data and predict the number of survival and death. Methods: The state of patients was clas sified into survival, death and censor. Being a kind of special categorical data, weighted least square method was used for estimating parameters. Results: Markov model for survival analysis of follow-up data was presented. With an example , the transition probability matrices were obtained. The number of survival and death at each observation point could be predicted respectively. Conclusion: Markov model is constructed to analyze clinical follow-up data. It can be used as an effective supplement for life table analysis.
K    EY WORDS　follow-up data; survival analysis; Markov model

［Acad J Sec Mil Med Univ, 1999, 20(7): 359～361］

　　疾病的治疗情况，一方面看结局好坏，另一方面还要看出现这种结局所经历的时间长短，这类资料一般通过随访收集。常见的随访起点是确诊日期、手术日期等；最明确的阳性结局是死亡，此外还有复发、致残等。随访资料常因失访等原因造成某些数据观察不完全，要用专门的方法进行统计处理，这类方法称为生存分析，常用的有寿命表法、乘积极限法、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 patients with breast cancer after operation