Interval censoring regression: predicted survival probabilities and CIs

[Innocent Karangwa]

Hi everyone
With interval censoring regression models, stintreg is used. To predict survival analysis, the following is done:
. predict surv_l surv_u, surv

. list surv_l surv_u in 1/5

	surv_l	surv_u
1.	1	.95814
2.	1	.948338
3.	1	.9754614
4.	.9828176	.9151379
5.	.9754614	.9029849

Can anyone explain to me how the above probabilities are combined to plot a survival curve? What is the command to calculate the probability at a particular time (just a single value)? Eg. Pr(5) = 0.80 instead of 0.97 and 0.90 as in the table above.
What is the command for confidence intervals? In other survival models, they use stci but with stintreg I don't know. How can I predict percentiles?
Please help!

[Isabel Canette (StataCorp)]

The behavior of predict after stintreg is similar to its behavior after streg. The only difference is that for stintreg we have two values of time to evaluate the prediction at: the two limits of the censoring interval.

Predictions are performed at the covariate pattern present in the data for each observation. Different covariate patterns will produce different survival estimates. Let's look at the cosmesis example:

. use http://www.stata-press.com/data/r15/cosmesis
. stintreg i.treat, interval(ltime rtime) distribution(weibull)
. predict S1 S2, surv
. list treat ltime rtime S1 S2 if _n<3 | _n>91


     +---------------------------------------------------+
     |       treat   ltime   rtime         S1         S2 |
     |---------------------------------------------------|
  1. |       Radio       0       7          1     .95814 |
  2. |       Radio       0       8          1    .948338 |
 92. | Radio+Chemo      35      39   .2382198   .1811756 |
 93. | Radio+Chemo      44      48   .1255049   .0917869 |
 94. | Radio+Chemo      48       .   .0917869          . |
     +---------------------------------------------------+

In observation 1, S1 and S2 show the value of the survivor function evaluated at treatment = Radio (i.e., corresponding to patients treated with radiotherapy only). S1 is evaluated at time 0 (ltime), and S2 is evaluated at time 7 (rtime). For observation 94, the predictions are computed for a patient who is treated with radiotherapy and chemotherapy. Predictions at different covariate patterns are not comparable, but you can plot survivor functions for a certain covariate pattern, for example:

. twoway line S1 ltime if treat == 0, sort
. twoway line S2 rtime if treat == 0, sort

The two lines above produce the same survivor function, corresponding to treatment == 0.

However, the easiest way to plot a survivor function is to use the stcurve command. To plot the survivor curve evaluated at the mean of the covariates, type

. stcurve, survival

You are also interested in quantiles and their confidence intervals. If you want to compute the median, you can use margins with the predict(median) option, for example:

.  margins, predict(median) at(treat=0)

Adjusted predictions                            Number of obs     =         94
Model VCE    : OIM

Expression   : Predicted median for (ltime,rtime], predict(median)
at           : treat           =           0

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |    39.3324   5.342493     7.36   0.000      28.8613    49.80349
------------------------------------------------------------------------------

For other quantiles, you can use the expression() option for margins. For example, let's assume I want the quantile 75 for treat==0.

If F(t) is the distribution function for treat=0, the quantile q_75 is the value that satisfies:

F(q_75) = .75;

F is the Weibull survivor distribution

F(t) = 1- exp(-lamda . t^p)

where p is the ancillary parameter, and lambda is exp(xb) with xb being the linear prediction for treat=0.

In other words,

q_75 = F^(-1)(.75)


Because we fit a Weibull model in the PH parameterization, the inverse distribution function can be obtained using the invweibullph() function with arguments lambda=exp(xb), p, and .75.


The ancillary parameter p can be obtained from the output from stintreg. We can use the coeflegend option with any estimation command to see how parameters are named.

. stintreg, coeflegend

Weibull PH regression                           Number of obs     =         94
                                                   Uncensored     =          0
                                                   Left-censored  =          5
                                                   Right-censored =         38
                                                   Interval-cens. =         51

                                                LR chi2(1)        =      10.93
Log likelihood = -143.19228                     Prob > chi2       =     0.0009

------------------------------------------------------------------------------
             |      Coef.  Legend
-------------+----------------------------------------------------------------
       treat |
Radio+Chemo  |   .9157008  _b[1.treat]
       _cons |  -6.292388  _b[_cons]
-------------+----------------------------------------------------------------
       /ln_p |   .4785787  _b[/ln_p]
-------------+----------------------------------------------------------------
           p |   1.613779
         1/p |   .6196635
------------------------------------------------------------------------------
Note: Estimates are transformed only in the first equation.
Note: _cons estimates baseline hazard.

The ancillary parameter p is not estimated directly. stintreg estimates ln(p), which is stored as _b[/ln_p]. This transformation is made to ensure that the parameter will be in the right domain (in this case, that it will be positive).

I will use exp(predict(xb)) as lambda and (exp(_b[/ln_p]) as the value of p in the invweibullph() function.

I use margins and its at() option to obtain a confidence interval for the quantile 75 for treat=0. I specify the expression for the quantile in the expression() option.

. margins, expression(invweibullph(exp(_b[/ln_p]), exp(predict(xb)), .75)) at(treat=0) force

Notice that if you use .5 instead of .75 in the above, you will obtain the same results as we obtained earlier for the median.

This is just the tip of the iceberg. There are many other options. For example, if you don't use the at() option, predictions will be averaged over your dataset. If you are looking for a specific result not covered in this answer, please contact our Technical Services.

-- Isabel

[Innocent Karangwa]

Dear Isabel
Thank you very much for the response. By confidence intervals, I meant the confidence intervals for the survival probabilities S1 and S2. I would also like to ask if one can add confidence bands on stcurve, survival. How can we plot cumulative survival (cumsurv) and add confidence bands on the plot?

[Isabel Canette (StataCorp)]

Currently, we don't have an option to add confidence bands to the plots produced by stcurve. The user-written command survci can be used to plot confidence bands for the predicted survivor function and the cumulative hazard function after fitting a model with streg, but it doesn't support models for interval-censored data.

If you want to obtain confidence intervals for S1 and S2, you can pass the formula for your prediction to the predictnl command; when you use the option ci(), confidence intervals will be computed. predictnl obtains confidence intervals using linearization, which might produce a lower limit that is negative. To avoid this, you can use a transformation, like the log function. Here is an example:

. clear
. webuse cosmesis
. stintreg i.treat, interval(ltime rtime) distribution(weibull)
. predict S1 S2, surv


. predict xb, xb
. predictnl logS1b = log(1-weibullph( exp(_b[/ln_p])),   ///
               exp(_b[1.treat]*treat + _b[_cons]), ltime) , ci(logl1 logl2) force

. gen S1b = exp(logS1b)
. gen l1 = exp(logl1)
. gen l2 = exp(logl2)

. list l1 S1 S1b l1 l2

The weibullph() function computes the cumulative Weibull distribution function, and the parameters p and lambda are entered the same way as in invweibullph(), shown in my previous post. The third parameter is the time at which I want to evaluate the function, in this case variable ltime. Because I want the survivor function, I compute one minus the distribution function. Once log(S1) and the logarithm of the limits of the CI are computed, I need to transform them back to the original scale by using the exponential function.

You will see that variables S1 and S2 contain the same values, and l1 and l2 contain the limits of the confidence intervals. By default, predictnl produces 95% confidence intervals, but you can customize the level with the level() option.

--Isabel

[Innocent Karangwa]

Dear Isabel
Thanks a lot for your answers