0

I'm interested in testing sex differences in the outcome Y (e.g., mortality), and if there are sex differences in the effect of BMI on mortality.

I planned to use a Cox regression model

primary<-coxph(Surv(time,status) ~ sex, data=data) second<-coxph(Surv(time,status) ~ sex*bmi, data=data)

And my first question is 1) if I can test proportionality assumption for an interaction. I use cox.zph function for the test

cox.zph(second)

and I'm wondering if it makes sense to test an assumption for an interaction term, and if this is the right way to do it.

Also, when I tested the assumption of proportionality for the first model (no interaction term), it was violated for sex. I read that the nonproportionality could be handled by including the time*covariate interaction term, so I did that...

coxph(Surv(start, stop, status.time) ~ sex+sex:stop, data=longdata)

And I think I can report that sex differences decreased over time (the interaction was significant), but I'm wondering what this means for my research question about the bmi*sex interaction ('second' model)...

2) should I include a three-way interaction term for that?

coxph(Surv(start, stop, status.time) ~ sex*bmi+sex:stop+sex:bmi:stop, data=longdata)

Thanks in advance for your help!!

llbia
  • 331

1 Answers1

1

It's OK to use the cox.zph() function (and inspection of the plots it can produce) to examine non-proportionality of hazard for an interaction term. You perhaps shouldn't worry too much about the non-proportionality in the first model using only sex as a predictor, as Cox models typically are improved by including as many predictors reasonably related to outcome as possible (without going so far as to overfit the model). For example, it's possible that omitting bmi or the bmi:sex interaction helps make sex appear to have a non-proportional hazard when used as a single predictor.

If you do find non-proportional hazards, an interaction of time with a covariate can be a way to deal with that problem. The particular way you propose to code it, however, doesn't work.

The stop time represents the actual last observation time (death or loss to follow-up), which isn't really a predictor: you didn't know it at the beginning of the study. It's a (potentially censored) outcome of the study. So using stop time in the regression is essentially using the outcome to help predict the outcome, which is circular reasoning. See this answer for a bit more detail and a link to a vignette that shows how to use a correct time-transform function to accomplish what you wish.

If you aren't interested in the association of outcome with sex per se but rather in the interaction of bmi with sex, you might be able to accomplish your goal more simply by using a predictor bmi*strata(sex). Analyzing by sex strata could also simplify any modeling of bmi with splines or other transformations that might be needed to meet the linearity assumption of Cox regression.

One final thought: although not necessary for statistical analysis, it might help to start with bmi values centered around some typical value like the median. Otherwise the coefficient value reported for sex in a bmi*sex interaction term will be the hazard ratio for sex at a bmi value of 0, which of course can't ever happen. Starting by centering such continuous predictors included in interaction terms can make it easier to explain to others what the hazard ratios mean in practice.

EdM
  • 92,183
  • 10
  • 92
  • 267
  • Thank you so much this is great help! I tested the assumption in the model with the interaction included and none of the main or interaction effects are significant there, but the global model is. Should I still say that the assumption is violated (and go down the time interaction road...?) or could this be because of the large sample size I have? – llbia Jun 16 '19 at 18:35
  • @llbia look carefully at the curves from cox.zph() to get an idea about what might be going on. Based on your initial finding that PH is violated by sex as a single predictor, I'd suggest starting with stratification by sex. Also, not modeling bmi properly to meet the linearity assumption of the regression could lead to loss of PH; see this question for a worked-out example. With a large data set use regression splines for bmi to help with any linearity problem, provided by rcs() terms from the R rms package. – EdM Jun 16 '19 at 19:42
  • The first issue is that I actually haven't been able to plot things cz I keep getting an error "Error in qr.default(t(const)) : NA/NaN/Inf in foreign function call (arg 1)" whenever I'm trying to use plot() or ggsurvplot()... Do you have any thoughts on this? – llbia Jun 16 '19 at 21:21
  • Second, thanks a lot for the suggestion. When I test a model, coxph(Surv(time,status)~ strata(sex) + rcs(bmi) +//there are other covariates//, data=data), there is no significant effect with the 4 bmi terms. Would that suggest linearity of the variable? When I compare the models Harrell's way as the answer here suggests, https://stackoverflow.com/questions/47937065/how-can-i-interpret-the-p-value-of-nonlinear-term-under-rms-anova, I see that the model is significantly different (p = .01). – llbia Jun 16 '19 at 21:21
  • @llbia questions about plotting problems and error messages are best asked on R help or other sites (with a reproducible example provided); they are off topic here. A significant difference between nested models with and without rcs(bmi) means that bmi is important overall as a predictor. The spline fits the nonlinear shape of the relation of bmi to outcome; with a large data set it makes sense to include the spline to improve linearity of the model as a whole even if the higher-order spline terms for bmi aren't individually "significant." – EdM Jun 16 '19 at 21:48