Fitting a survival model to see if there is interaction between gender and age

Question

hi I am new to survival analysis. I am trying to fit a cox regression model with age, sex, type of case(local vs imported cases), and regions(urban vs rural). I tries to fit two model

ModelA

> modelcox1 <- coxph(Surv(LOS1,censored)~ Gender + Age.Years. + 
                     Case.Type,data=covid, method="efron")
> summary(modelcox1)
Call:
coxph(formula = Surv(LOS1, censored) ~ Gender + Age.Years. + 
    Case.Type, data = covid, method = "efron")

  n= 524, number of events= 493 

                        coef exp(coef) se(coef)      z Pr(>|z|)    
GenderMale           0.06298   1.06501  0.09514  0.662 0.507965    
Age.Years.Over 40   -0.33155   0.71781  0.09422 -3.519 0.000434 ***
Case.TypeLocal Case  0.70353   2.02088  0.09333  7.538 4.76e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                    exp(coef) exp(-coef) lower .95 upper .95
GenderMale             1.0650     0.9390    0.8838    1.2833
Age.Years.Over 40      0.7178     1.3931    0.5968    0.8634
Case.TypeLocal Case    2.0209     0.4948    1.6831    2.4265

Concordance= 0.677  (se = 0.018 )
Likelihood ratio test= 62.92  on 3 df,   p=1e-13
Wald test            = 62.73  on 3 df,   p=2e-13
Score (logrank) test = 64.38  on 3 df,   p=7e-14

> 
> cox.zph(modelcox1, transform="km", global=TRUE)
            chisq df       p
Gender      0.713  1   0.399
Age.Years.  4.622  1   0.032
Case.Type  27.438  1 1.6e-07
GLOBAL     30.431  3 1.1e-06

Model B

> #Interaction between gender & age
> modelcox2 <- coxph(Surv(LOS1,censored)~ Gender*Age.Years. + 
                     Case.Type,data=covid, method="efron")
> summary(modelcox2)
Call:
coxph(formula = Surv(LOS1, censored) ~ Gender * Age.Years. + 
    Case.Type, data = covid, method = "efron")

  n= 524, number of events= 493 

                                 coef exp(coef) se(coef)      z Pr(>|z|)    
GenderMale                    0.16898   1.18409  0.12344  1.369    0.171    
Age.Years.Over 40            -0.16260   0.84993  0.15371 -1.058    0.290    
Case.TypeLocal Case           0.70727   2.02844  0.09346  7.568  3.8e-14 ***
GenderMale:Age.Years.Over 40 -0.26501   0.76719  0.19225 -1.378    0.168    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                             exp(coef) exp(-coef) lower .95 upper .95
GenderMale                      1.1841     0.8445    0.9296     1.508
Age.Years.Over 40               0.8499     1.1766    0.6288     1.149
Case.TypeLocal Case             2.0284     0.4930    1.6889     2.436
GenderMale:Age.Years.Over 40    0.7672     1.3035    0.5263     1.118

Concordance= 0.678  (se = 0.018 )
Likelihood ratio test= 64.81  on 4 df,   p=3e-13
Wald test            = 64.42  on 4 df,   p=3e-13
Score (logrank) test = 66.09  on 4 df,   p=2e-13
>
> test.ph<-cox.zph(modelcox2, transform="km", global=TRUE)
> test.ph
                   chisq df       p
Gender             0.499  1   0.480
Age.Years.         4.289  1   0.038
Case.Type          27.854  1 1.3e-07
Gender:Age.Years.  3.453  1   0.063
GLOBAL            30.524  4 3.8e-06

.

> 
> anova(modelcox1,modelcox2)
Analysis of Deviance Table
 Cox model: response is  Surv(LOS1, censored)
 Model 1: ~ Gender + Age.Years. + Case.Type
 Model 2: ~ Gender * Age.Years. + Case.Type
   loglik Chisq Df P(>|Chi|)
1 -2570.4                   
2 -2569.4 1.893  1    0.1689
> 

I dropped region because it was not significant in the univariate cox analysis. Gender also was not significant but I kept it because it is an important prognostic factor in my research. I have to analyse if there is an interaction between gender and age keeping gender as the main effects. I dont know how to proceed. Can anyone pls help to interprete the r output and also how to proceed further

After adding age as binary

 modelcox1 <- coxph(Surv(LOS1,censored)~ Gender + Age + Case.Type + Urban.Rural,data=covid, method="efron")
> summary(modelcox1)
Call:
coxph(formula = Surv(LOS1, censored) ~ Gender + Age + Case.Type + 
    Urban.Rural, data = covid, method = "efron")

  n= 524, number of events= 493 

                         coef exp(coef)  se(coef)      z Pr(>|z|)    
GenderMale           0.070114  1.072630  0.095402  0.735  0.46238    
Age                 -0.007457  0.992571  0.002881 -2.588  0.00965 ** 
Case.TypeLocal Case  0.694778  2.003265  0.097012  7.162 7.97e-13 ***
Urban.RuralVillage   0.069353  1.071815  0.160996  0.431  0.66663    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                    exp(coef) exp(-coef) lower .95 upper .95
GenderMale             1.0726     0.9323    0.8897    1.2932
Age                    0.9926     1.0075    0.9870    0.9982
Case.TypeLocal Case    2.0033     0.4992    1.6564    2.4228
Urban.RuralVillage     1.0718     0.9330    0.7818    1.4695

Concordance= 0.674  (se = 0.019 )
Likelihood ratio test= 57.84  on 4 df,   p=8e-12
Wald test            = 58.32  on 4 df,   p=7e-12
Score (logrank) test = 60  on 4 df,   p=3e-12

> cox.zph(modelcox1, transform="km", global=TRUE)
             chisq df       p
Gender       0.833  1  0.3614
Age         10.600  1  0.0011
Case.Type   28.058  1 1.2e-07
Urban.Rural  0.363  1  0.5469
GLOBAL      38.270  4 9.9e-08


modelcox2 <- coxph(Surv(LOS1,censored)~ Gender*Age + Case.Type +Urban.Rural ,data=covid, method="efron")
> summary(modelcox2)
Call:
coxph(formula = Surv(LOS1, censored) ~ Gender * Age + Case.Type + 
    Urban.Rural, data = covid, method = "efron")

  n= 524, number of events= 493 

                         coef exp(coef)  se(coef)      z Pr(>|z|)    
GenderMale           0.504504  1.656164  0.239259  2.109   0.0350 *  
Age                 -0.001406  0.998595  0.004179 -0.336   0.7366    
Case.TypeLocal Case  0.693814  2.001334  0.097065  7.148 8.81e-13 ***
Urban.RuralVillage   0.047343  1.048481  0.161239  0.294   0.7690    
GenderMale:Age      -0.011282  0.988782  0.005659 -1.993   0.0462 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                    exp(coef) exp(-coef) lower .95 upper .95
GenderMale             1.6562     0.6038    1.0362    2.6470
Age                    0.9986     1.0014    0.9905    1.0068
Case.TypeLocal Case    2.0013     0.4997    1.6546    2.4207
Urban.RuralVillage     1.0485     0.9538    0.7644    1.4382
GenderMale:Age         0.9888     1.0113    0.9779    0.9998

Concordance= 0.68  (se = 0.019 )
Likelihood ratio test= 61.82  on 5 df,   p=5e-12
Wald test            = 62.3  on 5 df,   p=4e-12
Score (logrank) test = 63.85  on 5 df,   p=2e-12

> test.ph<-cox.zph(modelcox2, transform="km", global=TRUE)
> test.ph
             chisq df       p
Gender       0.486  1  0.4856
Age         10.115  1  0.0015
Case.Type   28.706  1 8.4e-08
Urban.Rural  0.348  1  0.5555
Gender:Age   3.985  1  0.0459
GLOBAL      38.689  5 2.7e-07

the global test is insignificant for both model. I have to fit a survival model to see if there is interaction. How to do this, can anyone tell me how to proceed further?

Answer 1

I dropped region because it was not significant in the univariate cox analysis.

Maybe not a good idea. A point of multiple regression is to account for other covariates that might be masking a true association of a predictor with outcome. If your interest is in prediction then you generally want to include as many predictors that might be associated with outcome as possible, provided that you aren't overfitting. With nearly 500 events, you can include many more predictors than these (maybe up to 25, including interactions, for a 20/1 event/predictor ratio) and not be at much overfitting risk.

I have to analyse if there is an interaction between gender and age keeping gender as the main effects.

You already have a start to an answer in what you have shown. The (apparently binary) age predictor and the gender:age interaction both have estimated hazard ratios that might be of practical interest, but neither passes a test of "statistical significance" at the usual p < 0.05 (p-values of 0.29 and 0.17, respectively in the coxph summary; anova test comparing with and without interaction, p-value also 0.17). Again, if your interest is in prediction that doesn't really matter. Without the interaction term the age predictor did have a "statistically significant" association with outcome, so you clearly want to keep age in the model in some form.

Instead of arbitrarily cutting age into 2 groups, you should model age continuously. Binning a continuous predictor is almost never a good idea. Your two-group handling of age might be harming your ability to determine its true association with outcome and its interaction with gender. At least, evaluate age as a linear continuous predictor. With this many events you could do some very flexible nonlinear fitting of age, for example with restricted cubic splines. (For simplicity of interpretation, you might want to limit the interaction with gender to the linear part of the association.) That's likely to provide substantial improvement to your model.

Added in response to edited question:

Your modelcox2 is moving in the right direction. Although you say that age is coded as "binary" for that model, it seems that it's really been incorporated as a continuous linear predictor. (The break of age into groups above and below age 40 in the earlier models is binary in the sense I meant.) The coefficient for age in that model represents the association of age with outcome for Female (the reference level), an insignificant (both statistically and practically) hazard ratio (exp(coef)) of 0.999 per year (if age is expressed in years). The hazard ratio for the GenderMale:Age interaction, 0.989, represents the extra hazard per year for Male beyond the Female value, so the net hazard per year for Male is $0.999*0.989=0.988$, about 1% less risk of an event per year of age. The interaction passes (barely) the usual test of statistical significance (p = 0.046). You must apply your knowledge of the subject matter to decide if that "statistically significant" association represents a practically significant association.

There are some warnings here, however. First, your most substantial association with outcome is Case.Type, which shows a statistically significant departure from a proportional hazard. You need to evaluate that more closely; examine plots of coefficient estimates over time as described in the time dependence vignette to see if this also represents a practically significant departure. You might need to model that with time-dependent coefficients, as described there, or stratify your model by Case.Type. Second, you thus far have only modeled Age linearly. I suspect that its association with outcome is not linear. Try a non-linear modeling of Age as suggested above, with the rcs() function in the rms package. With the cph() Cox modeling in that package, you could write terms for Age interacting with Gender for example as rcs(Age,5) %ia% Gender to provide very flexible modeling of Age along with a linear interaction between Age and Gender. The anova() function then can document the overall significance of the predictors, the significance of the non-linear terms, and the significance of the interaction. Finally, although it doesn't make a difference with respect to predictions from the model, you might want to express Age as the difference from the median or mean age (if you haven't done so already). With the Age:Gender interaction, the coefficient for Male in modelcox2 represents the association of Gender with outcome when age is 0. Interpretation of that coefficient's practical importance might be easier if the 0 reference point for Age is something like the mean or median age.

General thoughts:

As you are just getting started, take advantage of the resources provided by Harrell's rms site. The course notes cover the principles that you need to learn for many types of models including Cox and other survival models.

Also, start using the associated tools provided by the rms package in R. Although there can be something of a learning curve, once you've mastered them they greatly simplify analysis and presentation of regression models, including survival models. For example, the rcs() function easily allows cubic spline fits to continuous predictors. For predictors involved in interactions or spline modeling, the anova() function designed for its objects nicely displays both overall significance including all interactions and nonlinearities, and the significance of the interactions/nonlinearities. It provides straightforward ways to examine combinations of predictors and to display results.