Some background: For my dissertation I'm working on a survival analysis on the hazard of exiting a firm for men vs women. After doing the basic analysis I want to investigate if the control variables in my model have different effects on hazard for men and women (e.g. age has more impact on men). To do this I made a full interaction model by interacting all of my control variables with the gender variable using R.

Coxinteractfull <- coxph(Surv(Duration,Censoring) ~ Control1 + Control2 + ... + Control1*Gender + Control2*Gender + ...)  

Then I checked if each of the interaction terms were significant by making a restricted model without the interaction term for each resprective control variable and comparing this with the full interaction model.

RestrictedControl1 <- coxph(Surv(Duration,Censoring) ~ Control1 + Control2 + ... + Control2*Gender + ...)  
anova(RestrictedControl1, Coxinteractfull)

I am planning on reporting this in a table, with the coefficients/hazard ratios seperately for men and women. Is this the correct way to go about this?

The current issue: I want to check if an interaction term is different for men and women, which essentially means doing a three-way interaction. The interaction is a continuous x nominal interaction. I'm not sure 1. how to model/code this and 2. how to interpret the results.

Any way of coding will lead to something like this:

variable coef exp(coef)
Male:Group1:Continuous -0.401 0.667 *
Male:Group2:Continuous -0.600 0.549
Female:Group1:Continuous 0.053 1.055
Female:Group2:Continuous -3.449 0.032

So based on this, I would say the continuous variable only has significant impact on the hazard for men in Group 1, correct? How do I now check if this interaction is significantly different for men and women?

Thank you!


1 Answer 1


This isn't really different from the issue in any regression model with interaction terms. The only difference is whether baseline predictor values are included in the baseline hazard (Cox) or contribute to an intercept term (other regression models). When you write in R with its default treatment coding of predictors:

outcome ~ Control * Gender

where Control is continuous, you get a set of 3 coefficients: a slope coefficient for Control at the baseline value of Gender (assuming a linear association of Control with log-hazard in a Cox model), a coefficient for the Gender difference at a value of Control = 0, and an interaction coefficient that can be thought of as the difference between the Gender values in the slope coefficient for Control.*

If you have an additional 2-level predictor Group then you can simply include that in your model:

outcome ~ Control * Gender * Group

and examine the significance of the single reported 3-way interaction term. That will tell you directly whether the Control:Group interaction term differs as a function of Gender.

If you have a large enough data set you can include 3-way interactions for all of your Continuous control variables:

outcome ~ Continuous1 * Gender * Group + Continuous2 * Gender * Group + ...

The Gender coefficient then represents the Gender-associated difference in the baseline Group with all Continuous variables at 0. With 2 levels of Group, its coefficient represents the corresponding difference for cases in the baseline Gender with all Continuous variables at 0.

Each Continuous:Gender interaction coefficient is as for Control:Gender in the first model, but only for cases in the baseline Group. Each Control:Group interaction coefficient is the Group difference in the Control coefficient at the baseline Gender. The significance of each 3-way interaction coefficient will answer your question about Group:Control interaction differences between the Genders.

Be careful, as the table of coefficients that you show wouldn't typically be reported that way by default. If Male and Group1 are the reference categories, then those categories won't show up in the default output. There would be no Male:Group1:Continuous interaction term; the coefficient for Continuous will represent the Continuous effect at that set of reference values.

What you show as Male:Group2:Continuous would be reported as the Group2:Continuous interaction term. Female:Group1:Continuous would be the Female:Continuous interaction term. The only 3-way interaction reported would be Female:Group2:Continuous. Remember that under usual treatment coding each of those interactions represents differences from lower levels in the interaction hierarchy.

You can use those coefficient estimates and their covariances to construct survival curves for any particular combination of predictor values. It's very easy to make mistakes in putting all that information together. With interactions in survival models I use the survplot(), Predict(), and contrast() functions in the R rms package.

*It's not clear that you have to do all of the restricted-model comparisons that you propose, as the significance of the Control:Gender interaction terms should be reported directly from the full model.

  • $\begingroup$ Thank you for this extensive reply, it's very helpfull. I have a more indepth question about the reporting of the tables that you talked about: I saw some papers showing a table where they show the HR of variables separately for two groups (in my case this would be men and women). As per your comment, I'm unsure if this is a fiable way to do this AND if so, how I should interpret the output. For example with the first coding I put down I would get something like: Variable 1 : Coef -0.125, exp(Coef) 0.883 * Female*Variable1: Coef 0.420, exp(Coef) 1.522 * $\endgroup$
    – KHT
    Dec 9, 2021 at 15:57
  • $\begingroup$ So if I would want to see if the HR of variable 1 is different for men and women, could I do that based on these results? If so how? $\endgroup$
    – KHT
    Dec 9, 2021 at 16:02
  • $\begingroup$ @KHT the details depend on the complexity of your model. For a simple model with a Gender:variable1 interaction term but no higher-order interactions like Gender:variable1:Group, the coefficient for the Gender:variable1 interaction gives you the extra log-hazard for females versus males with respect to variable1. The corresponding p-value indicates the significance. With higher-order interactions, you have to think carefully about just what you mean by "if the HR of variable 1 is different for men and women," as that difference might depend on the Group membership. $\endgroup$
    – EdM
    Dec 9, 2021 at 16:15
  • $\begingroup$ Right, so if I wanted to see if the variables have different effects on the probability of exit rates for men compared to women (e.g. aging raises the probability of exit for men but not for women), would I then be able to say (with the outcomes I just described) that for men the effect of variable 1 (e.g. age) is Coef -0.125, exp(Coef) 0.883, while for women it's Coef (-0.125+0.420), exp(Coef) 0.883*1.522. Or is this incorrect? $\endgroup$
    – KHT
    Dec 9, 2021 at 16:29
  • 1
    $\begingroup$ Side note: The sample size needed to support interaction analysis can be amazingly high. Basically, the number of subjects in the least frequent three-way group (binning age into say 4 intervals temporarily) needs to be in the neighborhood of 100. Other wise any non-astounding interaction effect (e.g., overall $P$-value for 3-way terms = 0.6) gives you no information because of power issues. $\endgroup$ May 2, 2023 at 12:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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