I have a general question regarding which p-value for interaction to report when modelling the interaction between a continuous predictor with restricted cubic splines and a categorical predictor.

mRCS_k3_int <- cph(Surv(pyears, event) ~ rcs(bmi, 3)*mets + age + sex + smoker, data = data)

Which gives this output:

##                 Wald Statistics          Response: Surv(pyears, event) 
##  Factor                                          Chi-Square d.f. P     
##  bmi     (Factor+Higher Order Factors)             53.87    4    <.0001
##   All Interactions                                  6.91    2    0.0317
##   Nonlinear (Factor+Higher Order Factors)          46.59    2    <.0001
##  mets     (Factor+Higher Order Factors)            50.67    3    <.0001
##   All Interactions                                  6.91    2    0.0317
##  age                                             2630.58    1    <.0001
##  sex                                              126.99    1    <.0001
##  smoker                                           194.27    2    <.0001
##  bmi    * mets     (Factor+Higher Order Factors)    6.91    2    0.0317
##   Nonlinear                                         3.87    1    0.0492
##   Nonlinear Interaction : f(A,B) vs. AB             3.87    1    0.0492
##  TOTAL NONLINEAR                                   46.59    2    <.0001
##  TOTAL NONLINEAR + INTERACTION                     51.41    3    <.0001
##  TOTAL                                           2935.04    9    <.0001

Which p-value do I report for the interaction between bmi and mets?

The way I interpret it, the first p-value for the interaction term, p=0.0317, is for the overall interaction (does bmi vary over categories of mets). The second and the third I am not so sure about.



The tests reported by anova() for a cph object are Wald tests on specific combinations of coefficients. They take into account the variances and covariances of the coefficients in question. Your model is relatively simple so several of the tests end up evaluating the same coefficient combinations. The print() function for anova.rms() allows for much more complicated models in general.

You can display the specific coefficient combinations corresponding to each test with the call:


In your model you have one linear and one nonlinear coefficient related to "bmi", called bmi and bmi' respectively. What you will find when you examine the details will be:

The "All Interactions" tests under both "bmi" and "mets", and the highest-level test for "bmi * mets", include interaction coefficients for both bmi and bmi' with mets.

The "Nonlinear" test under "bmi" includes both nonlinear coefficients involving "bmi", bmi' itself and the interaction bmi'*mets.

The "Nonlinear" and "Nonlinear Interaction" tests reported under "bmi * mets" represent the one nonlinear coefficient involving "bmi", the interaction bmi'*mets.

So your initial sense is correct. All 3 tests reporting p = 0.0317 in the table are evaluating both the linear and nonlinear interactions of "bmi" with "mets"; that would presumably be the comparison of most interest to yourself and your audience.

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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