Simulating a continuous time-varying covariate, properly

In this post we’ll take a look at how to simulate survival data with a continuous, time-varying covariate. The aim is to simulate from a data-generating mechanism appropriate for evaluating a joint longitudinal-survival model. We’ll use avalon, the merlin family’s simulator, to simulate the survival times, and the merlin command to fit the corresponding true model.

The model we are simulating from

Let’s assume a proportional hazards survival model, with the current value parameterisation. So for the $i$ th patient, we have the observed longitudinal outcome,

y_{i} (t) = m_{i} (t) + ϵ_{i} (t)

where

m_{i} (t) = X_{2 i} (t) β_{2} + Z_{i} (t) b_{i}

and $ϵ_{i} (t)$ is our normally distributed residual variability. We call $m_{i} (t)$ our trajectory function, representing the true underlying value of the continuous outcome at time t, which is a function of fixed and random effects, with associated design matrices and coefficients. We assume normally distributed random effects,

b_{i} ~ N (0, Σ)

Our survival model can be defined in terms of the hazard function,

h_{i} (t) = h_{0} (t) \exp (X_{1 i} β_{1} + α m_{i} (t))

where $h_{0} (t)$ is the baseline hazard function, $X_{1 i}$ is a vector of baseline covariates with associated log hazard ratios $β_{1}$ . We then link the current value of the biomarker directly to survival, where $α$ is a log hazard ratio for a one-unit increase in the biomarker, at time t.

Let’s assume a simple random intercept and random linear trend for the biomarker trajectory, i.e.

m_{i} (t) = (β_{20} + b_{0 i}) + (β_{21} + b_{1 i}) t

The challenge with simulating survival times from such a model is that the cumulative hazard function, $H_{i} (t) = \int_{0}^{t} h_{i} (u) d u$ , must integrate the hazard over the biomarker’s continuous trajectory $m_{i} (u)$ ; it has no closed form, and so cannot be inverted analytically. This is one of the situations avalon’s user subcommand was designed for: it integrates the hazard numerically, and finds each survival time by root finding (more details on the algorithm can be found in Crowther and Lambert (2013)). Assuming a step function for the biomarker would be simpler, but would not reflect the true data-generating mechanism, nor the benefit of the joint model, which can model time continuously.

Simulating the survival times

We’ll simulate a dataset of 500 patients, generating an id variable.

Stata

. clear

. set seed 249587

. set obs 500
Number of observations (_N) was 0, now 500.

. gen id = _n

Next we simulate the things we need at the patient level. This includes a binary treatment group variable, trt, and the random effects for later use in the longitudinal data. We’ll simulate two independent random effects, with standard deviations 1 and 0.1 (note you can use drawnorm to specify a covariance structure), representing the subject-specific deviations of the intercept and slope.

Stata

. gen trt = runiform()>0.5

. gen b0 = rnormal()

. gen b1 = rnormal(0,0.1)

Now we can simulate our survival times, still at the patient level. The key is being explicit in your definition of the hazard function. From above, substituting in the longitudinal trajectory, our hazard function becomes,

h_{i} (t) = h_{0} (t) \exp (X_{1 i} β_{1} + α [X_{2 i} (t) β_{2} + Z_{i} (t) b_{i}])

avalon user allows you to specify a user-defined hazard function, meaning you have complete generality. The hazard function needs to be written in Mata code (Stata’s matrix programming language), but is fairly self-explanatory. Key things to note:

You refer to time using {t}
You can directly include Stata variables in the definition of your hazard function, and they will be read in automatically
You must use the colon operator, which makes Mata do element-by-element operations to allow the simulation calculations to be vectorised, and hence much faster

Stata

. avalon user stime died,                                 /// new variables
>       hazard(                                           /// hazard function
>              0.1:*1.2:*{t}:^0.2 :*                      /// Weibull hazard
>              exp(0.2 :* (b0 :+ (0.1 :+ b1) :* {t}))     /// current value
>              )                                          ///
>       covariates(trt -0.5)                              /// treatment effect
>       maxtime(5)                                        //  admin. censoring

which assumes a scale and shape of $λ = 0.1$ and $γ = 1.2$ for the baseline Weibull hazard function, and $β_{20} = 0$ and $β_{21} = 0.1$ for the mean intercept and slope of the longitudinal trajectory. The association parameter is set to $α = 0.2$ . We also have a constant treatment effect with a log hazard ratio of −0.5.

Generating the observed biomarker measurements

Now we can move on to generating the observed longitudinal data. We know the true trajectory function, so we can generate whatever observation scheme we like. We’ll simulate up to 5 observations per patient, measured at baseline, 1, 2, 3, 4 “years”. This is easiest to do using expand, and then generating our time variable to represent the time of observation, using the observation number (_n) minus 1.

Stata

. expand 5
(2,000 observations created)

. bys id : gen time = _n-1

. drop if time>stime
(410 observations deleted)

The last line simply drops any time points that occur after a patient’s event time. Now we generate the observed longitudinal responses, at the observation time points, incorporating some residual variability from the true normal distribution, with standard deviation 0.5.

Stata

. gen xb = b0 + (0.1 + b1) * time

. gen y = rnormal(xb,0.5)

Finally, because we used expand, our survival times were replicated in all the extra rows. merlin, which we’re going to use to estimate our joint model, allows more than one event time per patient, so it would count each copy as a separate event. It’s crucial we remove the repeats, and only pass a single event time and event indicator per id to the estimation command.

Stata

. bys id (time) : replace stime = . if _n>1
(1590 real changes made, 1590 to missing)

. bys id (time) : replace died = . if _n>1
(1,590 real changes made, 1,590 to missing)

So our final dataset looks like:

Stata

. list id trt time y stime died if id==1 | id==13, sepby(id)

      +------------------------------------------------+
      | id   trt   time           y       stime   died |
      |------------------------------------------------|
   1. |  1     0      0    1.332419           5      0 |
   2. |  1     0      1     1.35186           .      . |
   3. |  1     0      2    1.907379           .      . |
   4. |  1     0      3    1.459253           .      . |
   5. |  1     0      4    .6210317           .      . |
      |------------------------------------------------|
  53. | 13     1      0   -.9602892   2.8579337      1 |
  54. | 13     1      1    .1079234           .      . |
  55. | 13     1      2   -1.303469           .      . |
      +------------------------------------------------+

Where we have our patient identifier, id, patient treatment group, trt, the observation times of the longitudinal outcome, time, the observed values of the longitudinal outcome, y, the survival time, stime, and associated event indicator, died.

Fitting the true model back

We can fit the true data-generating joint model very simply using merlin, as follows

Stata

. merlin (stime                           /// survival time
>                 trt                     /// baseline treatment (PH)
>                 EV[y] ,                 /// Expected Value of y
>          family(weibull, failure(died)) /// Weibull distribution & event ind.
>                 timevar(stime))         /// time-dependent, because of EV[y]
>        (y                               /// response
>                 time                    /// fixed effect of time
>                 time#M2[id]@1           /// random effect on time
>                 M1[id]@1 ,              /// random intercept
>                 family(gaussian)        /// distribution
>                 timevar(time))          //  time-dependent

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood = -3601.7803
Iteration 1:  Log likelihood = -2917.0355
Iteration 2:  Log likelihood = -2899.6704
Iteration 3:  Log likelihood = -2899.5605
Iteration 4:  Log likelihood = -2899.5604

Mixed effects regression model                           Number of obs = 2,090
Log likelihood = -2899.5604
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
         trt |  -.5093395   .1406999    -3.62   0.000    -.7851063   -.2335727
        EV[] |   .1965032   .0753041     2.61   0.009     .0489099    .3440966
       _cons |  -2.383013   .1511786   -15.76   0.000    -2.679318   -2.086708
  log(gamma) |   .1812595   .0668036     2.71   0.007     .0503268    .3121922
-------------+----------------------------------------------------------------
y:           |
        time |   .1058906   .0096587    10.96   0.000     .0869598    .1248214
 time#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .0550409   .0461364     1.19   0.233    -.0353849    .1454666
  sd(resid.) |   .4920555   .0098604                      .4731041    .5117662
-------------+----------------------------------------------------------------
id:          |
      sd(M1) |    .952584   .0332475                      .8895989    1.020029
      sd(M2) |   .1013752   .0127413                      .0792408    .1296923
------------------------------------------------------------------------------

which links the expected value (current value) of the longitudinal outcome directly to survival by using the EV[] element type. Note the use of the timevar() options in both sub-models, which makes sure merlin knows which variables represent time in each, and allows the appropriate calculation of the likelihood. The model converges nicely, with parameter estimates around the true values. For the survival sub-model, the log hazard ratio for treatment is −0.509 (we simulated −0.5), the association with the current value of the biomarker is 0.197 (0.2), and the log scale and log shape of the baseline hazard are −2.383 and 0.181 (log 0.1 = −2.303 and log 1.2 = 0.182). For the longitudinal sub-model, the mean intercept and slope are 0.055 and 0.106 (0 and 0.1), the standard deviations of the random intercept and slope are 0.95 and 0.10 (1 and 0.1), and the residual standard deviation is 0.49 (0.5). Every estimate is within 1.5 standard errors of its true value. To convince ourselves that everything is working, we can increase the sample size to check the estimates converge on the true values, or repeat the simulation many times and average the estimates to check for bias.

This example should hopefully provide a base case on which to expand for your own work. Check out the avalon and merlin pages for more.

MethodSimulation studies

Simulating survival data with a continuous time-varying covariate…the right way

The model we are simulating from

Simulating the survival times

Generating the observed biomarker measurements

Fitting the true model back

Need this applied to your own data?