机器算法验证 - 这种在 Cox 回归中拟合时间相关系数的方法有什么问题？ - 吾爱随笔录

这种在 Cox 回归中拟合时间相关系数的方法有什么问题？

机器算法验证生存 cox模型假设

2022-03-24 01:28:40

我有一个 Cox 比例风险模型。从 Schoenfeld 残差与时间图和零斜率的相应测试来看，有几个变量明显违反了 PH 假设（太多而无法分层）。

如果我们为简单起见采用单协变量情况，原始模型如下所示：

\hat{λ} (t) = λ_{0} (t) * \exp (X_{1} {\hat{β}}_{X_{1}})

$\hat\lambda(t) = \lambda_0(t)*\exp(X_1\hat\beta_{X_1})$

鉴于违反 PH 假设，我想尝试为问题变量拟合时间相关系数。这似乎很直观：

\hat{λ} (t) = λ_{0} (t) * \exp (X_{1} {\hat{β}}_{X_{1}} + X_{1} t * {\hat{β}}_{X_{1} t})

$\hat\lambda(t) = \lambda_0(t)*\exp(X_1\hat\beta_{X_1}+X_1t*\hat\beta_{X_1t})$

这简化为

\hat{λ} (t) = λ_{0} (t) * \exp (X_{1} ({\hat{β}}_{X_{1}} + t * {\hat{β}}_{X_{1} t}))

$\hat\lambda(t) = \lambda_0(t)*\exp(X_1(\hat\beta_{X_1}+t*\hat\beta_{X_1t}))$

因此，有效地，这将参数化拟合系数，它是时间的线性函数： $\hat\beta^\prime_{X_1}$

{\hat{β}}_{X_{1}}^{'} = {\hat{β}}_{X_{1}} + t * {\hat{β}}_{X_{1} t}

$\hat\beta^\prime_{X_1} = \hat\beta_{X_1}+t*\hat\beta_{X_1t}$

在我可以用这个模型愉快地跳舞到日落之前，我知道它是无效的（Terry Therneau #1和#2）。但为什么它无效？为什么我们必须构建一个包含起止时间的数据集呢？

2个回答

我对 R 中生存分析数据管理的机制不是很熟悉，但我想我可以解释为什么会发生这种情况并展示一个使用 Stata 的示例。

风险的危害会在变量随着时间流逝而发生变化的瞬间发生变化（不允许延迟或预期），尽管它在形成数据行的间隔中保持不变。您可以通过在观察到的故障时间拆分数据并手动生成模型中包含的随时间变化的协变量来获得正确的结果。拆分数据可以让您为每一集估计单独的 HR，并获得正确的时变系数。代价是数据膨胀。

这是一个使用髋部骨折研究的例子，其中一些老年人被给予一个充气装置以保护他们免受跌倒和初始剂量的骨骼强化药物。我们将初始剂量视为连续的，并与时间相互作用。人力资源模型是

h (t | x) = h_{0} (t) \cdot \exp (β \cdot p r o t e c t + γ \cdot i n i t_{-} d o s a g e + η \cdot i n i t_{-} d o s a g e \times t)

$h(t \vert x) = h_0(t) \cdot \exp(\beta \cdot protect + \gamma \cdot init_-dosage + \eta \cdot init_-dosage \times t)$

这可能在药理学上几乎没有意义，因为初始剂量的效果会随着时间的推移而衰减，但我们将继续使用它以使示例更类似于您的问题。

首先，我们加载数据并查看它：

. set more off

. use http://www.stata-press.com/data/cggm3/hip4, clear
(hip fracture study)

. sort id _t

. list id _t0 _t _d init_drug_level if inlist(id,1,5,9), sepby(id) noobs ab(30)

  +--------------------------------------+
  | id   _t0   _t   _d   init_drug_level |
  |--------------------------------------|
  |  1     0    1    1                50 |
  |--------------------------------------|
  |  5     0    4    1               100 |
  |--------------------------------------|
  |  9     0    5    0                50 |
  |  9     5    8    1                50 |
  +--------------------------------------+

变量 id 索引患者，_t0 是进入日期，_t 是以月为单位的学习时间，_d 表示失败。初始剂量为 50 或 100 mg。

这是错误的做事方式（创建剂量和时间之间的交互并按原样使用数据）：

. gen current_drug_level1 = init_drug_level *_t

. stcox protect init_drug_level current_drug_level1, nolog 

         failure _d:  fracture
   analysis time _t:  time1
                 id:  id

Cox regression -- Breslow method for ties

No. of subjects =           48                     Number of obs   =       106
No. of failures =           31
Time at risk    =          714
                                                   LR chi2(3)      =     56.88
Log likelihood  =   -70.129892                     Prob > chi2     =    0.0000

-------------------------------------------------------------------------------------
                 _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
            protect |   .1002764   .0563649    -4.09   0.000     .0333229    .3017554
    init_drug_level |   1.034508   .0152702     2.30   0.022     1.005008    1.064874
current_drug_level1 |   .9954672   .0011483    -3.94   0.000     .9932191    .9977204
-------------------------------------------------------------------------------------

这是使用自动选项的正确方法（因此无需拆分数据）：

  . stcox protect init_drug_level, tvc(init_drug_level) texp(_t) nolog

         failure _d:  fracture
   analysis time _t:  time1
                 id:  id


Cox regression -- Breslow method for ties

No. of subjects =           48                     Number of obs   =       106
No. of failures =           31
Time at risk    =          714
                                                   LR chi2(3)      =     33.23
Log likelihood  =    -81.95591                     Prob > chi2     =    0.0000

---------------------------------------------------------------------------------
             _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------+----------------------------------------------------------------
main            |
        protect |   .0868497   .0417166    -5.09   0.000     .0338774    .2226521
init_drug_level |   .9770202   .0134021    -1.69   0.090     .9511026    1.003644
----------------+----------------------------------------------------------------
tvc             |
init_drug_level |   .9999956   .0009067    -0.00   0.996     .9982201    1.001774
---------------------------------------------------------------------------------
Note: variables in tvc equation interacted with _t

以下是拆分记录后手动执行的方法：

. stsplit, at(failures)
(21 failure times)
(452 observations (episodes) created)

. gen current_drug_level2 = init_drug_level *_t 
. sort id _t
. list id _t0 _t _d *_drug_level* if inlist(id,1,5,9), sepby(id) noobs ab(30)

  +----------------------------------------------------------------------------------+
  | id   _t0   _t   _d   init_drug_level   current_drug_level1   current_drug_level2 |
  |----------------------------------------------------------------------------------|
  |  1     0    1    1                50                    50                    50 |
  |----------------------------------------------------------------------------------|
  |  5     0    1    0               100                   400                   100 |
  |  5     1    2    0               100                   400                   200 |
  |  5     2    3    0               100                   400                   300 |
  |  5     3    4    1               100                   400                   400 |
  |----------------------------------------------------------------------------------|
  |  9     0    1    0                50                   250                    50 |
  |  9     1    2    0                50                   250                   100 |
  |  9     2    3    0                50                   250                   150 |
  |  9     3    4    0                50                   250                   200 |
  |  9     4    5    0                50                   250                   250 |
  |  9     5    6    0                50                   400                   300 |
  |  9     6    7    0                50                   400                   350 |
  |  9     7    8    1                50                   400                   400 |
  +----------------------------------------------------------------------------------+

注意数据看起来有多么不同。关键是我们在患者还活着的情况下，每次时间增量都添加了一条记录，因此数据要大得多。请注意，对于患者 9，我们之前假设当前剂量为 250，等于第五个月末的剂量。现在我们在第 0 个月到第 5 个月有当前剂量变化，我们不再假装它与第 5 个月末的水平相同。因为 Cox 回归是对那些未能达到的受试者的一系列比较那些在失败期间有失败风险的受试者，当我们包含反映这种变化的额外数据时，我们现在正在将苹果与苹果进行比较。

当我们包含这些数据时，我们会得到与自动版本相同的估计值：

. stcox protect init_drug_level current_drug_level2, nolog 

         failure _d:  fracture
   analysis time _t:  time1
                 id:  id

Cox regression -- Breslow method for ties

No. of subjects =           48                     Number of obs   =       558
No. of failures =           31
Time at risk    =          714
                                                   LR chi2(3)      =     33.23
Log likelihood  =    -81.95591                     Prob > chi2     =    0.0000

-------------------------------------------------------------------------------------
                 _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
            protect |   .0868497   .0417166    -5.09   0.000     .0338774    .2226521
    init_drug_level |   .9770202   .0134021    -1.69   0.090     .9511026    1.003644
current_drug_level2 |   .9999956   .0009067    -0.00   0.996     .9982201    1.001774
-------------------------------------------------------------------------------------

作为旁注，您可以使用该tt函数在 R 中以“非无效”方式存档您的初始模型。请参阅包的 Cox 模型小插图中的Using Time Dependent Covariates and Time Dependent Coefficients 中Survival的最后一个示例（至少是 version 中的最后一个示例2.41-3）。

其它你可能感兴趣的问题

上一篇任意应用引导程序是否存在问题？下一篇全批 vs 在线学习 vs 小批