Cox regression baseline risk and interpretation of coefficients

Question

I recently encountered a problem with Cox regression that I couldn't find an answer to (please direct me to the appropriate resources if you know of any). So the problem goes as follows:

1. I have a gene expression dataset for two different tumour grades (A and B).
2. I want to check how gene expression affects HR depending on tumour grade and if the difference is statistically significant.

The thing is that I want to keep the gene expression as numerical values, I don't want to do median splits or any other transformation of these values.

I'm working in R and I'm using the survival package for the analysis.

I found some very useful information in this post (stratification in cox model) but I still feel a bit lost... mainly about the baseline hazard in some of the models and the correct way to compare the models. I have looked at several models:

• Stratified models

fit.strat <- summary(coxph(Surv(time, death == 1) ~ gene + strata(grade), data = test))

.            coef        exp(coef)       se(coef)      z      Pr(>|z|)
gene        0.2636      1.3016          0.1026        2.569     0.0102 *

The model assumes that the baseline hazard is different for the two grades (we don't know the difference) and that a unit increase in gene expression increases HR by 1.3016 for both grades.

• Two models

fit.grpA <- summary(coxph(Surv(time, death == 1) ~ gene, subset=(grade== "A"), data = test))

.         coef         exp(coef)        se(coef)       z      Pr(>|z|)
gene     0.7343       2.0840            0.1787         4.11     3.96e-05 ***

fit.grpB <- summary(coxph(Surv(time, death == 1) ~ gene, subset=(grade== "B"), data = test))

.          coef         exp(coef)        se(coef)     z  Pr(>|z|)
gene      0.0837      1.0873           0.12020       0.69      0.49

fit.grpA and fit.grpB are two different models with different baseline hazards and different effects of change in gene expression.

These models could be compared using the ANOVA test.

• No stratification

fit <- summary(coxph(Surv(time, death == 1) ~ grade:gene, data = test))

.            coef      exp(coef)        se(coef)          z            Pr(>|z|)
gradeA:gene  0.28458    1.32920         0.11894           2.393       0.01673 *
gradeB:gene  0.27266    1.31346         0.09799           2.783       0.00539 **

Here we have information about the change in HR from baseline, but what is the baseline? Some baseline level of gene expression without differentiation by grade?

Is it possible to check if the influence of the change in gene expression is statistically different between grades (if grade2:gene and grade3:gene results are significantly different)?

Answer 1

First, the way that you are modeling gene assumes that its values, in the scale that you are using, has a strictly linear association with log-hazard. That's probably not a good assumption. If it doesn't hold, you might end up with a proportional hazards (PH) violation. If you have enough events, you should consider fitting gene more flexibly, for example with a spline. That's a recommended type of transformation, unlike dichotomization.

Second, your first two approaches don't let you "check how gene expression affects HR depending on tumour grade and if the difference is statistically significant." The first assumes that the association of gene with outcome is independent of grade. The two models aren't nested in the second approach, so ANOVA can't be used to compare them directly.

Third, it's not clear how you have coded grade. Your first two approaches suggest that there are only 2 levels. The third shows coefficients involving grade2 and grade3, suggesting that there is also a grade1 being used as a reference.

Fourth, your third approach has a couple of problems. Without stratification, it assumes that the baseline hazard is the same for all levels of grade. The use of : instead of * for the interaction term means that the coefficients might just represent a hazard difference due to grade rather than the hazard associated with gene differing as a function of grade.

Once you have straightened out your coding of grade, modify your third approach to use * instead of :. If the PH assumption is met well enough, then the interaction term will cleanly document whether "gene expression affects HR depending on tumour grade." You can use survfit.coxph() to display estimated survival curves for any combination of gene and grade, including for the reference level of grade that doesn't show a coefficient in the model.

If PH is violated with respect to grade in the suggested interaction model, then combine stratification of grade with an interaction: gene*strata(grade). Then you have separate baselines for each grade and allow for different associations of gene with outcome as a function of grade.

Some final warnings. First, it doesn't look like your model contains any clinical information besides grade. In a survival model it's best to include as many outcome-associated clinical variables as you can (with precautions against overfitting), as omitting them can lead to bias or confounding. For example, what if your gene expression is related to age, and it's age that's directly associated with the hazard? Second, if this is a large-scale gene-expression data set so that you are doing this analysis for many genes, be very careful about the multiple comparisons problem.