Non-linear covariate effects in survival analysis

Let’s begin. There will be a single continuous covariate, representing age, with a non-linear effect influencing survival. We’ll simulate survival times under a data-generating model that incorporates a non-linear effect of age. We’ll then fit some models accounting for the non-linear effect of age, and finally make predictions for specified values of age. Sounds simple, but it’s often not.

Simulating survival data with a non-linear covariate effect

To simulate survival data we will use avalon, the merlin family’s simulator, and its standard subcommand, which simulates from the standard parametric distributions. Our data-generating model is,

h (t) = λ γ t^{γ - 1} \exp (X β_{1} + X^{2} β_{2})

i.e. a Weibull baseline, and both age (X) and its square impacting survival, under proportional hazards. We choose sensible values for our baseline parameters, $λ = 0.01$ and $γ = 1.5$ , and for the coefficients for age and age^2, $β_{1} = 0.01$ and $β_{2} = 0.001$ . Coding it up, generating a dataset of 1000 observations, simulating age from N(30,5^2), and applying administrative censoring at 10 years, we have:

Stata

. clear

. set seed 98798

. set obs 1000
Number of observations (_N) was 0, now 1,000.

. gen age = rnormal(30,5)

. gen age2 = age^2

. avalon standard                         /// command and setting
>         stime died,                     /// new variables to store survival
>                                         ///   time & event indicator
>         distribution(weibull)           /// survival time distribution
>         lambda(0.01) gamma(1.5)         /// scale and shape of the baseline
>         maxtime(10)                     /// max. follow-up time (admin.
>                                         ///   censoring)
>         covariates(age 0.01 age2 0.001) /// linear predictor - age and age^2
>                                         ///   with coefficients (log hazard
>                                         //    ratios)

Fitting a model with a non-linear covariate effect

We’re working with survival data, but there’s no need to stset it: merlin takes the survival time as the response in its model specification, and the event indicator through the failure() suboption of family().

We can now start fitting models. We’re going to start by fitting the true model using the merlin command, with a Weibull baseline hazard, entering age and its square as age and age#age. One of the differences between merlin and other estimation commands is in how non-linear effects can be modelled. But we’ll get to that.

Stata

. merlin (stime age age#age , family(weibull, failure(died)))

Fitting full model:

Iteration 0:  Log likelihood = -6756.6356
Iteration 1:  Log likelihood = -2159.8185
Iteration 2:  Log likelihood = -2153.1273
Iteration 3:  Log likelihood = -2078.2946
Iteration 4:  Log likelihood = -2073.9919
Iteration 5:  Log likelihood = -2073.8839
Iteration 6:  Log likelihood = -2073.8839

Fixed effects regression model                           Number of obs = 1,000
Log likelihood = -2073.8839
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
         age |  -.0263253   .0707263    -0.37   0.710    -.1649463    .1122957
     age#age |   .0018451   .0011354     1.63   0.104    -.0003802    .0040704
       _cons |   -4.31365   1.092658    -3.95   0.000     -6.45522    -2.17208
  log(gamma) |   .4132446   .0345778    11.95   0.000     .3454734    .4810158
------------------------------------------------------------------------------

All nice and simple. Because we simulated the data, we know what this model should find. The log scale and log shape of the baseline are estimated as −4.31 and 0.413, against log 0.01 = −4.61 and log 1.5 = 0.405. The coefficients for age and age^2 are −0.026 and 0.0018, against the 0.01 and 0.001 we simulated from. That looks worse than it is: over the ages we simulated, age and its square are almost perfectly correlated, so neither coefficient is well determined on its own (standard errors 0.071 and 0.0011), and each is within one standard error of its true value. What the data do pin down is their combined effect, which is what we care about when we predict.

Now merlin has a few tricks in its syntax. It has element functions, notably fp(), rcs() and bs() to directly create fractional polynomial, restricted cubic spline, and B-spline basis functions for you. Unlike other commands in Stata, we don’t have to create such basis functions prior to model estimation; we can let merlin do it. Let’s take a look.

The above model can be specified by using the fp() element function:

Stata

. merlin (stime fp(age, powers(1 2)) , family(weibull, failure(died)))
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_2

Fitting full model:

Iteration 0:  Log likelihood = -6756.6356
Iteration 1:  Log likelihood = -2159.8185
Iteration 2:  Log likelihood = -2153.1273
Iteration 3:  Log likelihood = -2078.2946
Iteration 4:  Log likelihood = -2073.9919
Iteration 5:  Log likelihood = -2073.8839
Iteration 6:  Log likelihood = -2073.8839

Fixed effects regression model                           Number of obs = 1,000
Log likelihood = -2073.8839
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
      fp():1 |  -.0263253   .0707263    -0.37   0.710    -.1649463    .1122957
      fp():2 |   .0018451   .0011354     1.63   0.104    -.0003802    .0040704
       _cons |   -4.31365   1.092658    -3.95   0.000     -6.45522    -2.17208
  log(gamma) |   .4132446   .0345778    11.95   0.000     .3454734    .4810158
------------------------------------------------------------------------------

which reproduces the previous results. To explain further, it will create fractional polynomial basis functions of age, with the powers of 1 and 2, so age and age^2. If we favoured restricted cubic splines to model continuous covariates, we could use:

Stata

. merlin (stime rcs(age, df(3) orthog) , family(weibull, failure(died)))
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_3

Fitting full model:

Iteration 0:  Log likelihood = -6756.6356
Iteration 1:  Log likelihood = -2159.8723
Iteration 2:  Log likelihood = -2152.5812
Iteration 3:  Log likelihood = -2077.5861
Iteration 4:  Log likelihood = -2073.2997
Iteration 5:  Log likelihood = -2073.1925
Iteration 6:  Log likelihood = -2073.1925

Fixed effects regression model                           Number of obs = 1,000
Log likelihood = -2073.1925
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
     rcs():1 |   .4146441   .0408094    10.16   0.000     .3346591     .494629
     rcs():2 |  -.0780577   .0408809    -1.91   0.056    -.1581828    .0020673
     rcs():3 |   .0452994   .0395197     1.15   0.252    -.0321577    .1227565
       _cons |  -3.402843   .1165609   -29.19   0.000    -3.631298   -3.174388
  log(gamma) |   .4129331   .0345708    11.94   0.000     .3451755    .4806906
------------------------------------------------------------------------------

which would create a spline function with three degrees of freedom, and orthogonalise the basis functions (which often improves convergence).

merlin creates the basis functions for you, tells you what they are called, and leaves them behind in the dataset. Importantly, it labels them in the output so you can see what’s what. The elegant thing is that we only have to specify the variable age to be expanded in splines, and as such changing the degrees of freedom etc. is extremely simple to do. No (re-)creating basis functions prior to model fitting.

Predicting from a model with a non-linear covariate effect

But! The benefit of this kind of syntax really comes in when obtaining predictions from a model. Since merlin creates the spline variables internally, when we want predictions at different values of age, the spline gets automatically recreated, and as such, getting a survival prediction for someone aged 40 is as simple as…

Stata

. predict s1, survival at(age 40)

Which is rather simple, we hope you agree. The prediction call wouldn’t change if you used fractional polynomials, or more degrees of freedom to model the effect of age.

Because we simulated the data, we can also check such predictions against the truth. Let’s predict survival from the spline model over ten years of follow-up for people aged 20, 30 and 40 (the mean age, and two standard deviations either side), on a grid of 101 time points using the timevar() option, with 95% confidence intervals from the ci option. Alongside, we calculate the true survival function we simulated from,

S (t) = \exp {- λ t^{γ} \exp (X β_{1} + X^{2} β_{2})}

Stata

. range tvar 0 10 101
(899 missing values generated)

. foreach a in 20 30 40 {
  2.     predict s`a', survival at(age `a') timevar(tvar) ci
  3.     gen true`a' = exp(-0.01 * tvar^1.5 * exp(0.01*`a' + 0.001*`a'^2))
  4. }
(899 missing values generated)
(899 missing values generated)
(899 missing values generated)

Comparing them at ten years,

Stata

. foreach a in 20 30 40 {
  2.     display "Age `a': fitted " %5.3f s`a'[101]                ///
>             " (95% CI " %5.3f s`a'_lci[101]                   ///
>             " to " %5.3f s`a'_uci[101]                        ///
>             "), true " %5.3f true`a'[101]
  3. }
Age 20: fitted 0.554 (95% CI 0.445 to 0.649), true 0.562
Age 30: fitted 0.359 (95% CI 0.322 to 0.396), true 0.350
Age 40: fitted 0.062 (95% CI 0.030 to 0.111), true 0.097

and plotting them over the whole of follow-up,

Stata

. twoway (rarea s20_lci s20_uci tvar, color(navy%20) lwidth(none))         ///
>        (rarea s30_lci s30_uci tvar, color(maroon%20) lwidth(none))       ///
>        (rarea s40_lci s40_uci tvar, color(forest_green%20) lwidth(none)) ///
>        (line s20 s30 s40 tvar, lcolor(navy maroon forest_green))         ///
>        (line true20 true30 true40 tvar, lcolor(black black black)        ///
>                 lpattern(dash dash dash))                                ///
>        , xtitle("Follow-up time (years)") ytitle("Survival")             ///
>        ylabel(0(0.2)1, angle(h) format(%2.1f)) xlabel(0(2)10)            ///
>        legend(pos(6) rows(1) order(4 "Age 20" 5 "Age 30"                 ///
>                 6 "Age 40" 7 "True survival"))

Predicted survival over ten years from the spline model for people aged 20, 30 and 40, each with a shaded 95% confidence band, and the true survival curves as black dashed lines. At ten years the predicted survival is 0.55, 0.36 and 0.06 at ages 20, 30 and 40, against true values of 0.56, 0.35 and 0.10, and each true curve stays inside its band throughout. — Survival predicted from the spline model at ages 20, 30 and 40 (solid lines, with 95% confidence intervals), against the true survival functions the data were simulated from (dashed).

The spline model tracks the truth closely at ages 20 and 30. At age 40 the predicted survival at ten years, 0.062, sits below the true 0.097, but the truth lies inside the 95% confidence interval (0.030 to 0.111), as it does at every time point for all three ages. The bands widen towards either end of the age distribution, where there are fewer people to learn from, which is worth remembering before predicting for anyone further out.

Now just imagine if you had a second continuous covariate…and then a spline-spline interaction…and then time-dependent effects on them… Convinced merlin makes life a bit easier…?

MethodsSimulation studies Prediction and risk modelling

Simulation, modelling and prediction with a non-linear covariate effect

Simulating survival data with a non-linear covariate effect

Fitting a model with a non-linear covariate effect

Predicting from a model with a non-linear covariate effect

Need this applied to your own data?