Allow for non-linear relationship between continuous independent variable and binary outcome in Cox regression

Question

One thing that is not clear to me if how one can model non-linear relationship between continuous independent variables and binary outcome (i.e. dependent variable) in a Cox regression model. Suppose that I want to explore the association of the age with death in a cohort of unselected patients. A linear regression may not provide the best fit, while a natural spline may does. The questions are:

1. I suppose that modelling non-linear relationship is not really possible in software like SPSS, but it may be in R. If so, how?
2. How can I report the results of the analysis coming from non-linear regression in a Cox model? If I use a spline regression, how can I report the coefficients so that the beta coefficients (and therefore the HR) are easily interpretable?

Answer 1

First, although you are proposing to model with a non-linear transformation of a continuous predictor, that is not technically a "non-linear regression," as the model will still be linear in the coefficients. With that terminological point out of the way:

1. Modeling continuous predictors with splines is straightforward in R (and I suspect with SPSS, too, though I don't use it). For survival models, you can model penalized splines (lots of regression coefficients, but with penalized magnitudes) with the pspline() function from the survival package, or restricted cubic splines (smooth curves with a few "knots" at specified positions along the range of the predictor, linear beyond the outer knots, with coefficients not penalized) with the rcs() function in the rms package. Frank Harrell's course notes discusses those starting in Sections 2.4.

2. If you fit a continuous predictor with a spline, there is no longer a single HR associated with it. The best way to display your results is to plot the relative hazard associated with the predictor (and its standard error) as a continuous function of its value, based on the spline fit. Examples with survival analysis are in Chapter 20 of Harrell's notes.

Answer 2

An answer to your second question: One thing that is very clear is that the effect of age on morbidity/mortality is highly non-linear. The very young and old are at higher risk, and the lowest risk is for those in the "prime of their lives."

Despite the up-then-down-then-up curvature, the quadratic function is not sufficiently flexible to model this effect, so cubic functions and splines as you mention are preferred to quadratics.

Because there are extra parameters, determining an "age risk effect" becomes complicated with such models. In addition, modelling interaction effects involve age also becomes more cumbersome when using quadratics, cubics, and splines.

The following paper was written to address these problems:

Moreau AR, Westfall PH, Cancio LC, Mason AD Jr, "Development and validation of an age-risk score for mortality predication after thermal injury," The Journal of Trauma, 30 Apr 2005, 58(5):967-972.

In it, the authors develop and validate an "Age Risk Score" which can be used to convert age to "age risk", which can be used in regression models as an ordinary linear term. The score and also be used to assess interaction effects in the usual way, by using a simple multiplicative term.

The specific "Age Risk" score is given by $$ AGESCORE = –5*Age + 14*(Age^2/100) – 7*(Age^3/10,000)$$

Lower values of the score correspond to lower risk, and higher values to higher risk.

Then your regression model (Cox, logistic, ordinary, etc) models the link of $Y$ in the usual simple regression way as $\beta_0 + \beta_1 AGESCORE + \beta_2 X_2 + \cdots$.

Hope this is helpful!