Mutiple variables violate Cox proportional-hazards assumption

[yaacovlawrence]

I am analyzing a large cancer registry (n=5000), focussing on cervical cancer. Median follow up about 9 years. Stata 11.
I performed a multivariate cox regression with overall survival as the endpoint. Variables include tumor size, age, ethnicity, histological grade, surgery type, use of radiotherapy.

When testing for the proportional-hazards assumption using Therneau and Grambsch method (with estat phtest) and lg-lg plots (with stphplot) I was surprised to discover that almost all the variables violated the assumptions.

Questions:
Is this a truly surprising result, or a known problem with large data-sets and long follow-up?
Would it be legitimate to use another regression model, e.g. Weibull distribution in this situation? or should I make almost every variable time-dependant within the Cox model?
Any other solutions? I tried limiting follow-up to 4 years, after which a few more variables fulfilled the assumption.

thank you!

Yaacov Lawrence MRCP
Dep. Radiation Oncology
Sheba Medical Center

[Maarten Buis]

A Weibull model won't help, as that model also imposes the proportial hazards assumption.

Other suggestions are harder to make, as it depends on the extend of the model violation and the purpose of the analysis. A Small violation may not be a problem, as a model is supposed to simply reality, i.e. be a bit wrong. Larger violations could be acceptable, as long as they don't have a big impact on the parameter of interest.

[yaacovlawrence]

1) Dear Maarten - thank you for your invaluable comments!
2) your comments prompted me to begin to understand the problem:
our data set is looking at non-metastatic disease, which will either be cured at time of surgery or recur. In the case of cervical cancer most recurrences (and hence cancer deaths) are early and occur within about 2 years.
Hence all the factors that play a role in determining cure rate (size of tumor, histology, use of radiotherapy) have all their effect in quite a small time window. Therefore it is entirely to be expected that the proportional hazard assumptions will vary over time, and the assumption therefore violated.

[Maarten Buis]

Depending on the exact purpose of the study non-proportional hazards could actually be of substantive interest rather than a "methodological inconvenience": You can test whether the effects are gone after 2 years by adding the variables in the tvc() option and in the main variable list and add the texp(_t > 730) (assuming analysis time is measured in days). Consider the example below:

webuse drugtr2
stset time, failure(cured)
stcox age drug1 drug2, tvc(drug1 drug2) texp(exp(_t>10))

         failure _d:  cured
   analysis time _t:  time

Iteration 0:   log likelihood = -116.54385
Iteration 1:   log likelihood = -97.718202
Iteration 2:   log likelihood = -96.258724
Iteration 3:   log likelihood = -96.231448
Iteration 4:   log likelihood =  -96.23143
Refining estimates:
Iteration 0:   log likelihood =  -96.23143

Cox regression -- Breslow method for ties

No. of subjects =           45                     Number of obs   =        45
No. of failures =           36
Time at risk    =  677.9000034
                                                   LR chi2(5)      =     40.62
Log likelihood  =    -96.23143                     Prob > chi2     =    0.0000

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
main         |
         age |   .8544492   .0286611    -4.69   0.000     .8000813    .9125117
       drug1 |   1.045584   .0130673     3.57   0.000     1.020283    1.071511
       drug2 |   1.035401   .0138725     2.60   0.009     1.008565     1.06295
-------------+----------------------------------------------------------------
tvc          |
       drug1 |   .9820064   .0056074    -3.18   0.001     .9710773    .9930586
       drug2 |   .9839315   .0059858    -2.66   0.008     .9722692    .9957337
------------------------------------------------------------------------------
Note: variables in tvc equation interacted with exp(_t>10)

In this case the hazard ratio before 10 days(?) for drug 1 is 1.05 and after 10 days this is multiplied by a factor .98, i.e. is decreased by 2%. This decrease is significant (the test reported after drug1 in the tvc equation. You can also test whether the hazard ratio after 10 days is significantly different from 1 (= no effect). Remember that the coefficients are in the form of log hazard ratios, so if you add the coefficient of drug 1 in the main equation to the coefficient of drug 1 in the tvc equation you get the log hazard ratio after 10 days, and to test if that hazard ratio equals 1 you can test whether that log hazard ratio equals 0, which in my example is not the case:

. test _b[main:drug1]+_b[tvc:drug1] = 0

 ( 1)  [main]drug1 + [tvc]drug1 = 0

           chi2(  1) =   12.36
         Prob > chi2 =    0.0004

Having said all that, if this is a situation where respondents are "cured" you may look at so called cure models:

Fitting and modeling cure in population-based cancer studies within the framework of flexible parametric survival models
T. M.-L. Andersson and P. C. Lambert. 2012.
Stata Journal Volume 12 Number 4.
http://www.stata-journal.com/article...ticle=st0165_1

Modeling of the cure fraction in survival studies
P. C. Lambert. 2007.
Stata Journal Volume 7 Number 3.
http://www.stata-journal.com/article...article=st0131