多状态Cox模型的参数估计和检验
作者：陈冠民 张选群 陈华
单位：陈冠民(湖北医科大学医学统计教研室 武汉430071)；张选群(湖北医科大学数学教研室)；陈华(湖北医科大学医学统计教研室 武汉430071)
关键词：多状态Cox模型；参数估计；生存分析
数理医药学杂志000324
摘 要 根据风险函数的定义，提出了多状态特定失效的时间的偏似然函数，给出了多状态Cox模型的参数和公共效应参数的估计方法，并对其假设检验进行了讨论。最后介绍了实现多状态Cox模型参数估计和假设检验的计算程序。
中图分类号：R 195.1 文献标识码：A
文章编号：1004-4337(2000)03-0228-02
, 百拇医药
一般的Cox模型只分析两种状态、一个生存时间的情形。如果事物存在三种以上的状态、两个以上的生存时间的情况晚，就需要用多状态的Cox模型进行分析。如糖耐量异常(IGT)、糖尿病(NIDDM)、糖尿病合并并发症共有三种状态、两个生存时间，本文拟对多状态Cox模型的参数估计和检验进行讨论。
1 参数估计
假定有K+1种失效事件，tki表示第i个个体发生第K种失效事件的时间，δki表示第k种失效事件的截尾指示变量。xki=[x1ki,…,xpki]表示第i个样本在生存时间t上、第K种失效事件所对应的p×1维协变量。在条件Xki上，假定生存时间、截尾指示变量和协变量是独立的，则对于第i个个体的第k种失效事件，其风险函数hki有如下的形式[1，2]：
, 百拇医药
hki(t)=hk0(t)exp[β｀kxki(t)] (1)
其中hk0(t)是非特定的基础风险函数，而βk=(β1k,…,βpk)是待估的失效回归参数。设Rk(t)={l:tkl≥t}是在K种失效情况下生存时间t以前的个体组成的危险集。因此，K种特定失效的偏似然函数为：
(2)
βk的最大似然估计可以通过解
而得到。当以上比例风险的假定成立时，则
k是特定的βk的一致估计。其估计方法采用Wei[3]提出的边缘分布估计方法。对于具有相同的截尾生存时间，采用Breslow渐进法进行计算[4]，其参数估计可采用MULCOX实现[5]。
, http://www.100md.com
如果协变量Xki是P维向量，有K+1种失效事件，则求出的回归系数应为一维的PK个元素的向量：
其中，
若比较第一个危险因素在K个阶段生存时间上的效应有无差异，可通过检验β11、β12、…、β1k之间有无差异来推断。若比较不同协变量在第一阶段生存时间上的效应有无差异，可通过比较β11、β21、…、βp1之间的有无差异来推断。
2 K个阶段回归系数之间的检验
, http://www.100md.com
对于大样本的情况，
趋近于均数为
联合协方差矩阵Q的正态分布。其中
的估计方法可参见文献[3]，同时，Wei还证明了Q矩阵的稳健性[6]。
为检验βk参数间(k=1,…,K)是否有显著性差异，其检验假设通常可以写成如下线性组合的形式：
H0:CβT=0
其中C是设定的r×PK阶的矩阵，称为比较矩阵(contrast matrix)。例如，如果要检验所有的多重生存时间均与任何协变量无关，则C矩阵是一个PK×PK阶的单位矩阵。其无效假设的Weld统计量为：
(3)
, 百拇医药
W服从于自由度为r的χ2分布。
假定要检验某个协变量在K+1种事件上的效应是否有统计学意义，其相应K个参数用ηk来表示，其中k=1,…,K。ηk可以通过设定的比较矩阵C与βT相乘而得到CβT=(η1,…,ηk)。
的协方差矩阵用
来表示，
,其无效假设Hk：ηk=0,(k=1,…,K)。则其二次项
(4)
, http://www.100md.com
可用来检验无效假设。可以看出式(4)是式(3)的特例，他们是完全等价的。
假定有4种状态、三个阶段生存时间、有两个协变量，则其回归系数为：
如果要比较第一个协变量在三个阶段生存时间上的效应是否有统计学意义，其比较的矩阵为：
不同的比较矩阵，可以比较同一个协变量在不同状态之间的参数是否有显著性差异，也可以比较同一生存时间阶段上不同协变量之间的参数是否有显著性差异。
在上述的假定中，如果要比较三个阶段生存时间上的回归系数是否相等，其比较矩阵为：
, http://www.100md.com
其H0为CβT=0，矩阵相乘后，
，即β11=β21=β31，其统计量可以根据式(3)得到。
3 K个阶段生存时间公共效应参数的估计
K个阶段生存时间公共参数用η来表示，则可以假定：η1=…=ηk=η。自然η的估计可以通过ηk的线性组合而得到。即：
(5)
其中，
，估计
所用的权重矩阵可用下式求得： (6)
, http://www.100md.com
其中，e=(1,…,1),
,用式(6)估计的权重系数在所有的线性组合中具有较小的渐近方差[7]。
的方差可以通过
求得。
4 计算程序的数据格式和控制参数
4.1 MULCOX规定的数据格式见表1。表中q是数据库中所包含的协变量的参数，它与模型中所分析的协变量不同，模型中所分析的协变量个数为p。当分类变量转为哑变量时，变量数可能增加，当没有协变量转化为哑变量时，则q=p。该程序配有一个子程序可以将分类变量转化为哑变量。数据可以Fortran的形式读入，也可以Free的形式读入。以Free的形式读入时，数据之间用空格或逗号隔开。
, 百拇医药
表1 多状态Cox比例风险回归的数据结构 生存时间和截尾指示变量
协变量
t11
δ11…
tk1
δk1
x11
x21
xq1
t12
, http://www.100md.com
δ12…
tk2
δk2
x12
x22
xq2
…
……
…
…
…
…
, http://www.100md.com
…
t1n
δ1n…
tkn
δkn
x1n
x2n
xqn
4.2 计算程序的控制参数见表2。表2 在拟合模型时程序要求输入的控制参数 参 数
类型
, http://www.100md.com 研究问题的题目
特征变量
数据文件名
特征变量
输出结果文件名
特征变量
样本含量
整数
事件状态数
整数
第一个生存时间变量
特征变量
…
, 百拇医药
…
第k个生存时间变量
特征变量
数据库中协变量数
整数
模型中要分析的协变量数
整数
第一个协变量名
特征变量
…
…
第p个协变量名
特征变量
, http://www.100md.com
输入数据的格式
特征变量
进行多变量检验的个数
整数
第一个检验C矩阵行数
整数
…
…
最后一个检验C矩阵行数
整数
要检验的公共参数的个数
整数
第一个公共参数检验时C矩阵的行数
, http://www.100md.com
整数
…
…
最后一个公共参数检验时C矩阵的行数
整数
第一个多变量检验的C矩阵
实际矩阵
最后一个多变量检验的C矩阵
实际矩阵
第一个公共参数检验的C矩阵
实际矩阵
…
, 百拇医药
…
最后一个公共参数检验的C矩阵
实际矩阵
湖北省自然科学研究基金课题和湖北省重大科研项目资助 参考文献
1，Prentice RL, et al. On the regression analysis of multivariate failure time data. Biometrika, 1981,68:373～379.
2，Clayton D, et al. Multivariate generalization of the proportional hazards model (with discussion). J R Statist Soc A, 1985,148:82.
, http://www.100md.com
3，Wei LJ, et al. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of American Statistics Association, 1989,84:1065～1073.
4，Breslow NE. Covariance analysis of censored survival data. Biometrics, 1974,30:89.
5，Lin DY. MULCOX: a computer program for the Cox regression analysis of multiple failure time variable. Computer Methods and Programs in Biomedicine, 1990,32:125～135.
6，Lin DY, Wei LJ. The robust inference for the Cox proportional hazard model. Journal of American Statistics Association, 1989,84:1074～1078.
7，Wei LJ, Johson WE. Combining dependent test with incomplete repeated measurements. Biometrika,1985,72:359～364.
收稿日期：1999-12-27, 百拇医药
单位：陈冠民(湖北医科大学医学统计教研室 武汉430071)；张选群(湖北医科大学数学教研室)；陈华(湖北医科大学医学统计教研室 武汉430071)
关键词：多状态Cox模型；参数估计；生存分析
数理医药学杂志000324
摘 要 根据风险函数的定义，提出了多状态特定失效的时间的偏似然函数，给出了多状态Cox模型的参数和公共效应参数的估计方法，并对其假设检验进行了讨论。最后介绍了实现多状态Cox模型参数估计和假设检验的计算程序。
中图分类号：R 195.1 文献标识码：A
文章编号：1004-4337(2000)03-0228-02
, 百拇医药
一般的Cox模型只分析两种状态、一个生存时间的情形。如果事物存在三种以上的状态、两个以上的生存时间的情况晚，就需要用多状态的Cox模型进行分析。如糖耐量异常(IGT)、糖尿病(NIDDM)、糖尿病合并并发症共有三种状态、两个生存时间，本文拟对多状态Cox模型的参数估计和检验进行讨论。
1 参数估计
假定有K+1种失效事件，tki表示第i个个体发生第K种失效事件的时间，δki表示第k种失效事件的截尾指示变量。xki=[x1ki,…,xpki]表示第i个样本在生存时间t上、第K种失效事件所对应的p×1维协变量。在条件Xki上，假定生存时间、截尾指示变量和协变量是独立的，则对于第i个个体的第k种失效事件，其风险函数hki有如下的形式[1，2]：
, 百拇医药
hki(t)=hk0(t)exp[β｀kxki(t)] (1)
其中hk0(t)是非特定的基础风险函数，而βk=(β1k,…,βpk)是待估的失效回归参数。设Rk(t)={l:tkl≥t}是在K种失效情况下生存时间t以前的个体组成的危险集。因此，K种特定失效的偏似然函数为：
(2)βk的最大似然估计可以通过解
而得到。当以上比例风险的假定成立时，则k是特定的βk的一致估计。其估计方法采用Wei[3]提出的边缘分布估计方法。对于具有相同的截尾生存时间，采用Breslow渐进法进行计算[4]，其参数估计可采用MULCOX实现[5]。, http://www.100md.com
如果协变量Xki是P维向量，有K+1种失效事件，则求出的回归系数应为一维的PK个元素的向量：
其中，
若比较第一个危险因素在K个阶段生存时间上的效应有无差异，可通过检验β11、β12、…、β1k之间有无差异来推断。若比较不同协变量在第一阶段生存时间上的效应有无差异，可通过比较β11、β21、…、βp1之间的有无差异来推断。
2 K个阶段回归系数之间的检验
, http://www.100md.com
对于大样本的情况，
趋近于均数为联合协方差矩阵Q的正态分布。其中的估计方法可参见文献[3]，同时，Wei还证明了Q矩阵的稳健性[6]。为检验βk参数间(k=1,…,K)是否有显著性差异，其检验假设通常可以写成如下线性组合的形式：
H0:CβT=0
其中C是设定的r×PK阶的矩阵，称为比较矩阵(contrast matrix)。例如，如果要检验所有的多重生存时间均与任何协变量无关，则C矩阵是一个PK×PK阶的单位矩阵。其无效假设的Weld统计量为：
(3), 百拇医药
W服从于自由度为r的χ2分布。
假定要检验某个协变量在K+1种事件上的效应是否有统计学意义，其相应K个参数用ηk来表示，其中k=1,…,K。ηk可以通过设定的比较矩阵C与βT相乘而得到CβT=(η1,…,ηk)。
的协方差矩阵用来表示，,其无效假设Hk：ηk=0,(k=1,…,K)。则其二次项 (4), http://www.100md.com
可用来检验无效假设。可以看出式(4)是式(3)的特例，他们是完全等价的。
假定有4种状态、三个阶段生存时间、有两个协变量，则其回归系数为：
如果要比较第一个协变量在三个阶段生存时间上的效应是否有统计学意义，其比较的矩阵为：不同的比较矩阵，可以比较同一个协变量在不同状态之间的参数是否有显著性差异，也可以比较同一生存时间阶段上不同协变量之间的参数是否有显著性差异。
在上述的假定中，如果要比较三个阶段生存时间上的回归系数是否相等，其比较矩阵为：
, http://www.100md.com
其H0为CβT=0，矩阵相乘后，
，即β11=β21=β31，其统计量可以根据式(3)得到。3 K个阶段生存时间公共效应参数的估计
K个阶段生存时间公共参数用η来表示，则可以假定：η1=…=ηk=η。自然η的估计可以通过ηk的线性组合而得到。即：
(5)其中，
，估计所用的权重矩阵可用下式求得： (6), http://www.100md.com
其中，e=(1,…,1),
,用式(6)估计的权重系数在所有的线性组合中具有较小的渐近方差[7]。的方差可以通过求得。4 计算程序的数据格式和控制参数
4.1 MULCOX规定的数据格式见表1。表中q是数据库中所包含的协变量的参数，它与模型中所分析的协变量不同，模型中所分析的协变量个数为p。当分类变量转为哑变量时，变量数可能增加，当没有协变量转化为哑变量时，则q=p。该程序配有一个子程序可以将分类变量转化为哑变量。数据可以Fortran的形式读入，也可以Free的形式读入。以Free的形式读入时，数据之间用空格或逗号隔开。
, 百拇医药
表1 多状态Cox比例风险回归的数据结构 生存时间和截尾指示变量
协变量
t11
δ11…
tk1
δk1
x11
x21
xq1
t12
, http://www.100md.com
δ12…
tk2
δk2
x12
x22
xq2
…
……
…
…
…
…
, http://www.100md.com
…
t1n
δ1n…
tkn
δkn
x1n
x2n
xqn
4.2 计算程序的控制参数见表2。表2 在拟合模型时程序要求输入的控制参数 参 数
类型
, http://www.100md.com 研究问题的题目
特征变量
数据文件名
特征变量
输出结果文件名
特征变量
样本含量
整数
事件状态数
整数
第一个生存时间变量
特征变量
…
, 百拇医药
…
第k个生存时间变量
特征变量
数据库中协变量数
整数
模型中要分析的协变量数
整数
第一个协变量名
特征变量
…
…
第p个协变量名
特征变量
, http://www.100md.com
输入数据的格式
特征变量
进行多变量检验的个数
整数
第一个检验C矩阵行数
整数
…
…
最后一个检验C矩阵行数
整数
要检验的公共参数的个数
整数
第一个公共参数检验时C矩阵的行数
, http://www.100md.com
整数
…
…
最后一个公共参数检验时C矩阵的行数
整数
第一个多变量检验的C矩阵
实际矩阵
最后一个多变量检验的C矩阵
实际矩阵
第一个公共参数检验的C矩阵
实际矩阵
…
, 百拇医药
…
最后一个公共参数检验的C矩阵
实际矩阵
湖北省自然科学研究基金课题和湖北省重大科研项目资助 参考文献
1，Prentice RL, et al. On the regression analysis of multivariate failure time data. Biometrika, 1981,68:373～379.
2，Clayton D, et al. Multivariate generalization of the proportional hazards model (with discussion). J R Statist Soc A, 1985,148:82.
, http://www.100md.com
3，Wei LJ, et al. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of American Statistics Association, 1989,84:1065～1073.
4，Breslow NE. Covariance analysis of censored survival data. Biometrics, 1974,30:89.
5，Lin DY. MULCOX: a computer program for the Cox regression analysis of multiple failure time variable. Computer Methods and Programs in Biomedicine, 1990,32:125～135.
6，Lin DY, Wei LJ. The robust inference for the Cox proportional hazard model. Journal of American Statistics Association, 1989,84:1074～1078.
7，Wei LJ, Johson WE. Combining dependent test with incomplete repeated measurements. Biometrika,1985,72:359～364.
收稿日期：1999-12-27, 百拇医药