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 our `survsim`

command. Our data-generating model is,

i.e. a Weibull baseline, and both age (X) and it’s square impacting survival, under proportional hazards. We choose sensible values for our baseline parameters, and for the coefficients for age and age^2. 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:

```
. 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
. survsim /// command
> 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
> maxt(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, so let’s now `stset`

our simulated data.

```
. stset stime, failure(died)
Survival-time data settings
Failure event: died!=0 & died<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
1,000 total observations
0 exclusions
--------------------------------------------------------------------------
1,000 observations remaining, representing
653 failures in single-record/single-failure data
6,756.636 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 10
```

We can now start fitting models. We’re going to start by fitting the true model using the `stmerlin`

command, which is a wrapper function for the more powerful `merlin`

command. One of the differences of `merlin`

to other estimation commands, is in how non-linear effects can be modelled. But we’ll get to that.

```
. stmerlin age age#age , distribution(weibull)
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
Survival model Number of obs = 1,000
Log likelihood = -2073.8839
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
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. Now `merlin`

and `stmerlin`

have a few tricks in their syntax. They have 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, you can let `merlin`

do it. Let’s take a look.

The above `stmerlin`

model can be specified by using the `fp()`

element function:

```
. stmerlin fp(age, powers(1 2)) , distribution(weibull)
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
Survival model Number of obs = 1,000
Log likelihood = -2073.8839
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
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:

```
. stmerlin rcs(age, df(3) orthog) , distribution(weibull)
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.8722
Iteration 2: log likelihood = -2152.5807
Iteration 3: log likelihood = -2077.5848
Iteration 4: log likelihood = -2073.2984
Iteration 5: log likelihood = -2073.1912
Iteration 6: log likelihood = -2073.1912
Survival model Number of obs = 1,000
Log likelihood = -2073.1912
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
rcs():1 | .4146476 .0408083 10.16 0.000 .3346649 .4946303
rcs():2 | -.0781111 .0408836 -1.91 0.056 -.1582414 .0020193
rcs():3 | .045297 .0395261 1.15 0.252 -.0321728 .1227667
_cons | -3.402849 .116561 -29.19 0.000 -3.631304 -3.174393
log(gamma) | .4129337 .0345708 11.94 0.000 .3451761 .4806912
------------------------------------------------------------------------------
```

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

`stmerlin`

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 50 is as simple as…

`. predict s1, survival at(age 50)`

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. 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…?

Subscribe to get new post alerts directly in your inbox:

## One thought on “Simulation, modelling and prediction with a non-linear covariate effect in survival analysis”