Survival analysis with interval censoring

Interval censoring occurs when we don’t know the exact time an event occurred, only that it occurred within a particular time interval. Such data is common in ophthalmology and dentistry, where events are only picked up at scheduled appointments, but they actually occurred at some point since the previous visit. Arguably, we could say all survival data is subject to some degree of interval censoring – recording time of death to the nearest day, for example, sets a time interval of one day.

In this post, we will look at how to simulate interval censored survival times, and then show how to use merlin to fit an interval censored flexible parametric survival model. We simulate continuous time survival data, and then apply interval censoring.

First we’ll simulate a dataset of 1000 observations, and include a binary variable which we will assume is a treatment group indicator coded 0 for control, and 1 for active treatment. Observations will be assigned to treatment with 50% probability. And of course we set a seed for reproducibility.

. clear

. set seed 72549

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

. gen trt = runiform()>0.5

We will use survsim to simulate survival times from a Weibull distribution with shape and scale parameters, 1.2 and 0.1. We’ll assume a log hazard ratio due to treatment of -0.5. Let’s apply administrative right censoring at 5 years.

. survsim stime event ,   dist(weib) lambda(0.1) gamma(1.2) ///
>                         covariates(trt -0.5) maxtime(5)

We’ve simulated continuous, exactly observed survival times, which is what we would observe in an ideal world. We will now apply interval censoring. Let’s assume observations were only made at annual intervals, i.e. every year, so any event times will be rounded up to the nearest year to give the right side of the interval, and floored to give the left interval. This code is easily adapted to widen or shorten the intervals. The left interval will be stored in st1,

. gen st1 = floor(stime) if stime < 5
(590 missing values generated)

. replace stime = st1 + 1 if st1 < 5
(410 real changes made)

merlin identifies interval censored observations through the failure() indicator being coded as a 2, rather than a 1 for exactly observed events or a 0 for right censored observations, so we make that change (as survsim codes all events by default as 1). We’re assuming all events are interval censored, but merlin allows any combination of exactly observed, interval and right censored observations.

. replace event = 2 if event == 1
(410 real changes made)

Let’s have a quick look at the dataset (always recommended),

. list stime st1 event in 1/5

     +---------------------+
     | stime   st1   event |
     |---------------------|
  1. |     5     .       0 |
  2. |     5     .       0 |
  3. |     5     .       0 |
  4. |     5     4       2 |
  5. |     1     0       2 |
     +---------------------+

The first three observations are right censored at 5 years, and hence their st1 is missing. Observation 4 is an interval censored event which occurred between 4 and 5 years, and observation 5 is an event which occurred between 0 and 1 years.

Now we are all setup and ready to fit an interval censored survival model with merlin. We will keep it relatively simple and show how to fit a Royston-Parmar flexible parametric survival model which uses restricted cubic splines to model the baseline log cumulative hazard function, but importantly now allowing for interval censoring.

The important option to focus on is the linterval() option, which lets you pass the variable which contains the left interval time for those observations which are interval censored, which must have their associated event indicator specified as a 2. The right side of the interval is simply that in your outcome variable, stime in our case. We specify our family() to be rp and use 3 degrees of freedom for the baseline log cumulative hazard function (1 df is equivalent to a Weibull – the true model), and add in the failure() and linterval() options,

. merlin  (stime trt , family(rp, df(3) failure(event) linterval(st1))) 
variables created: _rcs1_1 to _rcs1_3

Fitting full model:

Iteration 0:   log likelihood = -1792.7911  
Iteration 1:   log likelihood = -1507.7269  
Iteration 2:   log likelihood = -1321.6384  
Iteration 3:   log likelihood = -1316.8481  
Iteration 4:   log likelihood = -1316.8434  
Iteration 5:   log likelihood = -1316.8434  

Fixed effects regression model                           Number of obs = 1,000
Log likelihood = -1316.8434
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |            
         trt |  -.6211561   .1019953    -6.09   0.000    -.8210632    -.421249
       _cons |  -.6726925   .0631678   -10.65   0.000    -.7964992   -.5488858
------------------------------------------------------------------------------
    Warning: Baseline spline coefficients not shown - use ml display

Our model fits smoothly and we estimate our log hazard ratio due to treatment. You can increase the sample size to something large to remove Monte Carlo error, to make sure there’s no bias – we leave that to you to check. Next steps would of course be to assess proportional hazards in our treatment effect by including an interaction between a function of time and trt (just add stime#trt, for example), and we can do any number of other things (non-linear effects, frailties, …).

Hopefully this simple example gives a primer for both data setup and syntax. Remember, once a feature is added to merlin, it tends to work with everything – such as any of the other inbuilt survival distributions, including the user-defined ones.