4
$\begingroup$

I have a Cox proportional hazards model. Judging by Schoenfeld residual vs. time plots and corresponding tests for zero slope, there is clear violation of the PH assumption for several of the variables (too many to stratify).

If we take the single-covariate case for simplicity, the original model looks like:

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

Given the violation of the PH assumption, I want to try fitting time-dependent coefficients for the problem variables. It seems intuitive to fit this:

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

which simplifies to

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

So, effectively, this would parametrically fit a coefficient $\hat\beta^\prime_{X_1}$ that is a linear function of time: $$ \hat\beta^\prime_{X_1} = \hat\beta_{X_1}+t*\hat\beta_{X_1t} $$

Before I could happily dance off into the sunset with this model, I learned it's invalid (Terry Therneau #1 and #2). But WHY is it invalid? Why do we have to construct a dataset with start-stop times instead?

$\endgroup$
3
$\begingroup$

I am not very familiar with the mechanics of survival analysis data management in R, but I think I can explain why this happens and show an example using Stata.

The hazard of the risk changes the instant the variable changes as time flows (no delays or anticipation allowed), though it remains constant in the intervals that forms the rows in the data. You can achieve the right result by splitting the data at the observed failure times and manually generating the time-varying covariates that you include in the model. Splitting the data allows you to estimate a separate HR for each episode and get the right time-varying coefficient. The cost is data inflation.

Here's an example using the hip fracture study where some elderly folks were given an inflatable device to protect them from falls and an initial dosage of bone-fortifying drug. We will treat the initial dosage as continuous and interact it with time. The HR model is

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

This probably makes very little sense pharmacologically since the effect of the initial dosage should decay over time, but we will roll with it to make the example more similar to your question.

First we load the data and take a peek at it:

. 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 |
  +--------------------------------------+

The variable id indexes patients, _t0 is the entry date, _t is study time in months, and _d indicates failure. The initial dosage was either 50 or 100 mg.

Here's wrong way to do things (create the interaction between dosage and time and use the data as is):

. 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
-------------------------------------------------------------------------------------

Here's the right way using the automated option (so there's no need to split the data):

  . 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

Here's how you would do it by hand, after splitting the records:

. 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 |
  +----------------------------------------------------------------------------------+

Note how different the data looks. The key point is that we have added a record for every time that time increments while the patient is still alive, so the data is much bigger. Note that for patient 9, we previously assumed that the current dosage was 250, equal to what it was at the end of the fifth month. Now we have current dosage vary in the months 0 to 5, and we're no longer pretending that it was the same level as it was at the end of month 5. Since Cox regression is a series of comparisons of those subjects who fail to those subjects at risk of failing for periods where there was some failure, we are now comparing apples to apples when we include the extra data that reflects this variation.

When we include this data, we get the same estimates as with the automated version:

. 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
-------------------------------------------------------------------------------------
$\endgroup$
  • 1
    $\begingroup$ Crystal-clear, and the example is perfect. Thank you! $\endgroup$ – half-pass Aug 20 '14 at 23:42
1
$\begingroup$

As a side note, you can archive your initial model in a "non-invalid" way in R by using the tt function. See the last example in the Using Time Dependent Covariates and Time Dependent Coefficients in the Cox Model vignette of the Survival package (at least the last example in version 2.41-3).

$\endgroup$
  • $\begingroup$ Yes! I believe this is a newer feature. tt() can be quite slow for large datasets, though. $\endgroup$ – half-pass Oct 25 '17 at 11:41
  • 1
    $\begingroup$ It should be an equivalent model to the stcox from stata that @Dimitriy V. Masterov mentions. The key part is that ... we have added a record for every time that time increments while the patient is still alive, so the data is much bigger. It should be much slower as @Dimitriy V. Masterov points out. $\endgroup$ – Benjamin Christoffersen Oct 25 '17 at 15:38
  • 1
    $\begingroup$ Makes sense. I haven't done a head-to-head comparison vs. STATA, but I have found that tt() is sometimes slower than manually splitting the data in an efficient way and then fitting the Cox model. Nevertheless, it's definitely the easiest solution for smallish datasets. $\endgroup$ – half-pass Oct 25 '17 at 19:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.