Interval-censored Cox model

Let's see it work

• Fitting a model for single-record-per-subject interval-censored data

• Stratification

• Checking proportional-hazards assumption

• Fitting a model for multiple-record-per-subject interval-censored data

• Graphing survivor functions

• Time-varying covariates

• Additional resources

Fitting a model for single-record-per-subject interval-censored data

We use data from Zeng et al. (2016) that studies time to HIV infection in a cohort study of injecting drug users in Thailand. The dataset contains 1,124 subjects. Those subjects initially tested negative for the HIV virus. They were then followed and assessed for HIV-1 seropositivity through blood tests approximately every four months. Because the subjects were tested periodically, the exact time of HIV infection was not observed, but the times are known to fall in intervals between blood tests with the lower and upper endpoints recorded in variables ltime and rtime.

The dataset also records several baseline factors, such as centered age at recruitment (age_mean), sex (male), history of needle sharing (needle), history of drug injection (inject), and whether a subject has been in jail at the time of recruitment (jail). The dataset contains one record per subject with the event-time information recorded as interval data and is called single-record-per-subject data. Here is a subset of the dataset for subjects 271–274:

. use https://www.stata-press.com/data/r18/idu
(Modified Bangkok IDU Preparatory Study)

. format ltime rtime age_mean %6.2f

. list ltime rtime age_mean male needle inject jail if _n >=271
     & _n<=274



      
 ltime   rtime   age_mean   male   needle   inject   jail 


 271. 
 22.00       .      -6.46    Yes      Yes       No     No 

 272. 
  3.80    9.41       8.54     No       No       No    Yes 

 273. 
 20.66       .     -11.46    Yes      Yes       No     No 

 274. 
  0.00    3.87      -4.46    Yes      Yes      Yes    Yes 

We can fit a Cox proportional hazards model in which time to HIV infection depends on above factors of interest. In a single-record-per-subject format, we must specify option interval() with stintcox.

. stintcox age_mean i.male i.needle i.inject i.jail, interval(ltime rtime)
note: using adaptive step size to compute derivatives.

Performing EM optimization (showing every 100 iterations):
Iteration 0:    Log likelihood = -1086.2564
Iteration 100:  Log likelihood = -597.65634
Iteration 200:  Log likelihood = -597.57555
Iteration 295:  Log likelihood = -597.56443

Computing standard errors: ............................ done

Interval-censored Cox regression                    Number of obs     =  1,124
Baseline hazard: Reduced intervals                         Uncensored =      0
                                                        Left-censored =     41
                                                       Right-censored =    991
                                                       Interval-cens. =     92

                                                    Wald chi2(5)      =  17.10
Log likelihood = -597.56443                         Prob > chi2       = 0.0043



             
 
                OPG                                             
             
 
 Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]

    age_mean 
 
   .9684341   .0126552    -2.45   0.014     .9439452    .9935582
             
 

        male 
 

        Yes  
 
   .6846949   .1855907    -1.40   0.162     .4025073    1.164717
             
 

      needle 
 

        Yes  
 
   1.275912   .2279038     1.36   0.173     .8990401    1.810768
             
 

      inject 
 

        Yes  
 
   1.250154   .2414221     1.16   0.248     .8562184    1.825334
             
 

        jail 
 

        Yes  
 
   1.567244   .3473972     2.03   0.043     1.014982    2.419998

Note: Standard-error estimates may be more variable for small datasets and
      datasets with low proportions of interval-censored observations.

We find that age at the recruitment is significantly associated with lower risk of HIV infection and being in jail at enrollment is associated with higher risk of HIV infection. Other factors are not statistically significant.

Stratification

If we assume that the baseline hazard function for female is different from that for male, we can fit a stratified Cox proportional hazards model by using the strata() option:

. stintcox age_mean i.needle i.inject i.jail, interval(ltime rtime) strata(male)
note: using adaptive step size to compute derivatives.

Performing EM optimization (showing every 100 iterations):
Iteration 0:   Log likelihood = -1087.0536
Iteration 100: Log likelihood = -585.59848
Iteration 200: Log likelihood = -585.53143
Iteration 282: Log likelihood =  -585.5222

Computing standard errors: ........................... done

Stratified interval-censored Cox regression         Number of obs     =  1,124
Baseline hazard: Reduced intervals                         Uncensored =      0
Strata variable: male                                   Left-censored =     41
                                                       Right-censored =    991
Event-time interval:                                   Interval-cens. =     92
  Lower endpoint: ltime
  Upper endpoint: rtime
                                                    Wald chi2(4)      =  14.84
Log likelihood = -585.5222                          Prob > chi2       = 0.0051



             
 
                OPG                                             
             
 
 Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]

    age_mean 
 
   .9682508   .0126326    -2.47   0.013     .9438052    .9933295
             
 

      needle 
 

        Yes  
 
   1.276222   .2270302     1.37   0.170     .9005422    1.808625
             
 

      inject 
 

        Yes  
 
   1.245357   .2393768     1.14   0.254     .8544367    1.815131
             
 

        jail 
 

        Yes  
 
    1.57314   .3490687     2.04   0.041     1.018337    2.430205

Note: Standard error estimates may be more variable for small datasets and
      datasets with low proportions of interval-censored observations.

Our conclusion here is similar to the case without stratification.

Checking proportional hazards assumption

It is important to evaluate the validity of the underlying assumption for the Cox proportional hazards model, which is that the hazard ratio is constant over time. One way of testing the proportional-hazards assumption for a covariate is to test whether the coefficient associated with that covariate is time invariant. This can be accomplished by including an interaction between this covariate and a function of time in the model and testing whether the corresponding coefficient equals zero.

For example, to check the proportional-hazards assumption for inject, we can include an interaction term between inject and time using the tvc() option:

. stintcox age_mean i.male i.needle i.inject i.jail, interval(ltime rtime)
     tvc(i.inject) nohr
note: using adaptive step size to compute derivatives.

Performing EM optimization (showing every 100 iterations):
Iteration 0:   Log likelihood = -1086.2564
Iteration 100: Log likelihood = -597.58166
Iteration 200: Log likelihood = -597.50007
Iteration 295: Log likelihood = -597.48895

Computing standard errors: ......................................... done

Interval-censored Cox regression                    Number of obs     =  1,124
Baseline hazard: Reduced intervals                         Uncensored =      0
                                                        Left-censored =     41
Event-time interval:                                   Right-censored =    991
  Lower endpoint: ltime                                Interval-cens. =     92
  Upper endpoint: rtime
                                                    Wald chi2(6)      =  17.36
Log likelihood = -597.48895                         Prob > chi2       = 0.0080



             
 
                OPG                                             
             
 
 Coefficient  std. err.      z    P>|z|     [95% conf. interval]

main         
 
                                                                
    age_mean 
 
  -.0321889   .0131026    -2.46   0.014    -.0578696   -.0065082
             
 

        male 
 

        Yes  
 
  -.3817124   .2714232    -1.41   0.160    -.9136921    .1502674
             
 

      needle 
 

        Yes  
 
   .2434583   .1791214     1.36   0.174    -.1076131    .5945298
             
 

      inject 
 

        Yes  
 
   .3194498   .3209978     1.00   0.320    -.3096943    .9485939
             
 

        jail 
 

        Yes  
 
   .4498814   .2230002     2.02   0.044      .012809    .8869538

tvc          
 
                                                                
      inject 
 

        Yes  
 
  -.0080255   .0208446    -0.39   0.700    -.0488801    .0328292

Notes: Standard error estimates may be more variable for small datasets and
       datasets with low proportions of interval-censored observations.
       Variables in tvc equation interacted with _t.

Wald test that [tvc] = 0: chi2(1) = 0.1482                Prob > chi2 = 0.7002

With the nohr option, the first equation, main, reports the coefficients for the covariates that do not vary over time; the second equation, tvc, reports the estimated coefficients for covariates interacted with time. There is no evidence that the coefficient of inject interacted with time is different from 0. Hence the proprtional-hazards assumption appears to hold for covariate inject.

For single-record interval-censored data, Stata also provides two graphical commands to visually assess the proportional-hazards assumption. These two commands can be used without running stintcox first.

For a single categorical covariate in a Cox model, you can use stintphplot, which plots -ln{-ln(survival)} curves for each category versus ln(analysis time). If the plots are parallel, the proportional-hazards assumption has not been violated. Alternatively, you can use stintcoxnp to plot the nonparametric maximum-likelihood estimation survival curve versus the Cox predicted survival curve for each category. When the two curves are close together, the proportional-hazards assumption is valid.

When a Cox model contains multiple covariates, as mentioned in the above example, only stintphplot is appropriate for graphically checking the proportional-hazards assumption. In that case, we need to use the adjustfor() option.

For example, to graphically check the proportional-hazards assumption for inject, we include all covariates except inject in the adjustfor() option:

. stintphplot, interval(ltime rtime) by(inject) adjustfor(age_mean i.male i.needle i.jail)

Fitting Cox model with covariates from option adjustfor()
for inject = No ...

Fitting Cox model with covariates from option adjustfor()
for inject = Yes ...

A separate Cox model, which contains all covariates from the adjustfor() option, is fit for each level of inject. And the two plots are almost parallel, which indicates that the proportional-hazards assumption has not been violated for the categorical variable inject.

Fitting a model for multiple-record-per-subject interval-censored data

Interval-censored data can also be recorded as multiple-record-per-subject format. In this format, each subject may contain multiple records with multiple examination times, the event status at each examination time, and potential time-varying covariates at each examination time. Here we use an extended version of previous single-record-per-subject dataset, which contains all the baseline covariates as well as a time-varying imprisonment indicator variable (jail_vary) that indicates whether the subject has been imprisioned since the last clinic visit. It also records a subject identifier (id), the examination time of the blood test (time), and whether the blood test is positive at the examination time (is_seropos). Here is a subset of this dataset for subject 271-274:

. use https://www.stata-press.com/data/r18/idu2
(Modified Bangkok IDU Preparatory Study with time-varying variable jail_vary)

. format time age_mean %6.2f

. list id time is_seropos age_mean male needle inject jail_vary
     if id >=271 & id <=274, sepby(id) noobs abbreviate(10)


    

    
  id    time   is_seropos   age_mean   male   needle   inject   jail_vary 

    

    
 271    4.89           No      -6.46    Yes      Yes       No          No 

    
 271    9.31           No      -6.46    Yes      Yes       No          No 

    
 271   13.38           No      -6.46    Yes      Yes       No         Yes 

    
 271   17.97           No      -6.46    Yes      Yes       No         Yes 

    
 271   22.00           No      -6.46    Yes      Yes       No          No 

    

    
 272    3.80           No       8.54     No       No       No         Yes 

    
 272    9.41          Yes       8.54     No       No       No          No 

    

    
 273    3.93           No     -11.46    Yes      Yes       No          No 

    
 273    8.00           No     -11.46    Yes      Yes       No          No 

    
 273   12.07           No     -11.46    Yes      Yes       No         Yes 

    
 273   15.97           No     -11.46    Yes      Yes       No         Yes 

    
 273   20.66           No     -11.46    Yes      Yes       No         Yes 

    

    
 274    3.87          Yes      -4.46    Yes      Yes      Yes         Yes 

    

We use stintcox to fit a Cox proportional-hazards model in which the time to HIV infection depends on baseline covariates age_mean, male, needle, inject, and time-varying covariate jail_vary. In a multiple-record-per-subject format, we must specify options id(), time(), and status() with stintcox.

. stintcox age_mean i.male i.needle i.inject i.jail_vary, id(id) time(time)
     status(is_seropos)
note: time-varying covariates detected in the data; using method nearleft to
      impute their values between examination times.
note: using adaptive step size to compute derivatives.

Performing EM optimization (showing every 100 iterations):
Iteration 0:   Log likelihood = -1086.2564
Iteration 100: Log likelihood = -598.45375
Iteration 200: Log likelihood = -598.35872
Iteration 285: Log likelihood = -598.34887

Computing standard errors: ........................................ done

Interval-censored Cox regression                   Number of obs      =  6,453
Baseline hazard: Reduced intervals                 Number of subjects =  1,124
                                                           Uncensored =      0
ID variable: id                                         Left-censored =     41
Examination time: time                                 Right-censored =    991
Status indicator: is_seropos                           Interval-cens. =     92

                                                   Wald chi2(5)       =  17.03
Log likelihood = -598.34887                        Prob > chi2        = 0.0044



             
 
                OPG                                             
        time 
 
 Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]

    age_mean 
 
   .9714605    .012757    -2.20   0.027     .9467762    .9967884
             
 

        male 
 

        Yes  
 
   .6678044   .1816576    -1.48   0.138     .3918353    1.138138
             
 

      needle 
 

        Yes  
 
   1.271409   .2275426     1.34   0.180     .8952546    1.805609
             
 

      inject 
 

        Yes  
 
   1.370672   .2575405     1.68   0.093     .9484142    1.980928
             
 

   jail_vary 
 

        Yes  
 
   1.440966   .2916178     1.81   0.071     .9691488    2.142481

Time varying: jail_vary
Note: Standard error estimates may be more variable for small datasets and
      datasets with low proportions of interval-censored observations.

Compared with the previous example, after we account for time-varying imprisonment, the hazard ratio for inject increases from 1.25 to 1.37, but the effect of imprisonment decreases from 1.57 for baseline jail to 1.44 for time-varying jail_vary.

Graphing survivor functions

After fitting the model, we can use stcurve to plot the estimated survivor, failure, hazard, or cumulative hazard function. By default, stcurve evaluates the functions at the overall means of covariates.

. stcurve, survival
note: function evaluated at overall means of covariates.

With multiple-record data with time-varying covariates, we may want to incorporate the time-varying nature of those covariates when plotting the functions. In this case, we can use the attmeans option to evaluate the function at time-specific means.

. stcurve, survival attmeans
note: function evaluated at time-specific means of covariates.

We can also use the atframe() option to specify our own time-varying covariate values to be used to evaluate the survivor function. Suppose we want to plot the survivor curve for an individual with the same covariate pattern as subject 1 in our dataset. We can create a new frame called id1 and use frame put to copy the relevant information for subject 1 in this new frame. Then we list the data we just saved in frame id1.

. frame put time age_mean male needle inject jail_vary
     if id==1, into(id1)

. frame id1: list



      
  time   age_mean   male   needle   inject   jail_v~y 


  1. 
  4.13       5.54    Yes       No      Yes        Yes 

  2. 
  8.26       5.54    Yes       No      Yes         No 

  3. 
 12.33       5.54    Yes       No      Yes         No 

  4. 
 16.10       5.54    Yes       No      Yes         No 

  5. 
 20.10       5.54    Yes       No      Yes         No 


  6. 
 24.23       5.54    Yes       No      Yes         No 

  7. 
 28.62       5.54    Yes       No      Yes         No 

  8. 
 32.59       5.54    Yes       No      Yes         No 

  9. 
 36.20       5.54    Yes       No      Yes         No 

 10. 
 40.56       5.54    Yes       No      Yes        Yes 

We can graph the survivor curve for this particular profile by typing

. stcurve, survival atframe(id1)
note: function evaluated at specified values of selected covariates and overall
      means of other covariates (if any).
note: covariate values from frame id1 used to evaluate function.

References

Zeng, D., L. Mao, and D. Lin. 2016. Maximum likelihood estimation for semiparametric transformation models with interval-censored data. Biometrika 103: 253–271.S

Zeng, D., F. Gao, and D. Lin. 2017. Maximum likelihood estimation for semiparametric regression models with multivariate interval-censored data. Biometrika 104: 505–525.

Turnbull, B. W. 1976. The empirical distribution function with arbitrarily grouped, censored and truncated data. Journal of the Royal Statistical Society, Series B 38: 290–295.

Farrington, C. P. 2000. Residuals for proportional hazards models with interval-censored survival data. Biometrics 56: 473–482.