## Abstract

We describe how to estimate progression-free survival while dealing with interval-censored data in the setting of clinical trials in oncology. Three procedures with SAS and R statistical software are described: one allowing for a nonparametric maximum likelihood estimation of the survival curve using the EM-ICM (Expectation and Maximization-Iterative Convex Minorant) algorithm as described by Wellner and Zhan in 1997; a sensitivity analysis procedure in which the progression time is assigned (i) at the midpoint, (ii) at the upper limit (reflecting the standard analysis when the progression time is assigned at the first radiologic exam showing progressive disease), or (iii) at the lower limit of the censoring interval; and finally, two multiple imputations are described considering a uniform or the nonparametric maximum likelihood estimation (NPMLE) distribution. *Clin Cancer Res; 22(23); 5629–35. ©2016 AACR*.

## Introduction

In oncology, progression-free survival (PFS) is often used as the primary endpoint in clinical trials (1), as using the response rate to evaluate targeted therapies has been strongly criticized. PFS appears as a relevant alternative, although it does not represent a direct clinical benefit, such as overall survival. Although this endpoint was recommended by a consortium of American experts (2), it remains controversial. Fleming and colleagues showed that this endpoint was interesting to document efficacy of a treatment, but many biases could occur (3). One of them is related to the imaging schedule. In fact, for patients with solid tumors, the disease progression is based on the radiologic RECIST (Response Evaluation Criteria in Solid Tumors Version 1.1) criteria. Thus, imaging visits follows a previously planned schedule, but the precise date of progression is then unknown but occurs between two assessment dates. It might be thought that the solution is to suggest a more compact imaging schedule, but that is not in the interest of patients and implies too much financial expenditure (4). Moreover, in the context of data sharing or retrospective analysis of “real-life” data, we have to deal with interval censoring while assessing PFS over time.

Since the 1950s, several methods have emerged to deal with interval-censored data, which may be parametric or not, as described precisely by Kim (5). Beyond those methods, multiple imputation procedures may also be considered, because the issue of interval censoring might be defined as a missing data problem. Any fixed-point imputation method, such as the usual PFS estimation, does not take the whole uncertainty into account in its variance estimate, and multiple imputation methods may, therefore, be very useful to estimate the PFS without bias and its variance with all the variability due to the interval censoring.

Some of these methods are implemented in statistical software (6): the Expectation–Maximization (EM) algorithm (7), the Iterative Convex Minorant (ICM) algorithm (8), and the EM-ICM algorithm (9). Although tutorials for them exist, some may be obsolete or not only for medical research purposes (10, 11).

Although these methods are now widely implemented and recommended in oncology by the European Medicines Agency and the U.S. Food and Drug Administration as part of a sensitivity analysis (12–14), only a few clinical trials have assessed PFS with methods that take interval-censored data into account (15–19).

To make these methods as accessible as possible, we provide a step-by-step guide from data formatting to output interpretation of interval-censoring methods using SAS and R software. To do so, each step was conducted in the database of a real phase II clinical trial assessing the efficacy of a first-line therapeutic regimen for patients with metastatic colorectal cancer.

## Preamble

### Dataset context

As an example of interval-censored data in oncology, 58 patients with unresectable, metastatic colorectal cancer in a nonrandomized phase II trial were reanalyzed. They received biweekly a 90-minute intravenous infusion of bevacizumab (5 mg/kg), followed by a 90-minute intravenous infusion of irinotecan (180 mg/m^{2}), followed by a simplified LV5FU2 regimen [leucovorin (400 mg/m^{2}) and bolus fluorouracil (400 mg/m^{2}) on day 1 and a 46-hour infusion of fluorouracil (2,400 mg/m^{2})] (20).

### Data formatting

PFS was defined as the time from study entry to the occurrence of the “event,” namely progression or death. To compute any type of survival analyses, data have to be properly formatted with regard to the following:

– Complete observations dataset: For PFS, this means that all patients have progressed or died, and the only variable needed is the time elapsed until the occurrence of the event. This is really rare in clinical studies.

– Right-censored dataset: This is when the follow-up was not long enough to allow all events to occur during the study. In such cases, two variables are needed in the dataset: a status variable to know whether the event occurred (status = 1) or not (status = 0) and a time variable describing the time elapsed until the last news (if status = 0) or until the occurrence of the event (if status = 1).

– Interval-censored dataset: This means that the exact moment at which the event occurred is unknown, although there has been a time interval during which the event occurred. This is the focus of this article. In such cases, two variables are needed in the dataset: T1 and T2. The time period between T1 and T2 is the interval during which the event occurred. Data have to be formatted as follows:

For patients who died without progression (uncensored data):

T1 = T2 = survival time

For patients who have not progressed or died during the study (right-censored data):

T1 = time elapsed to last visit

T2 = infinite time (missing data for SAS or R)

For patients who have progressed (interval-censored data):

T1 = time elapsed to last visit showing no progression

T2 = time elapsed to first visit showing disease progression

If only a few or no patients are interval censored in the dataset, it is pointless to use methods that deal with interval-censored data.

### Statistical analysis

In this article, we estimate PFS over time using three different procedures to compare advanced and more standard methods. The pros and cons of these procedures are presented in Table 1.

– “Interval-censored procedure”: To assess PFS precisely over time, we first computed a nonparametric maximum likelihood estimation (NPMLE) of the survival curve using the EM-ICM algorithm (9).

– “Sensitivity analysis procedure”: For the sake of simplicity, we also suggested a method proposed by Panageas and colleagues using standard PFS estimations (Kaplan–Meier estimator) that assumes the progression time is (i) the middle of the interval, (ii) the upper limit of the interval (which corresponds to the usual procedure while considering the date of progression as being the date of the first radiographic scan showing progression), and (iii) the lower limit of the interval (21).

– “Multiple imputations”: In parallel, we decided to present two multiple imputation methods to estimate PFS and confidence interval (CI), as the issue of interval censoring may be regarded as a missing data problem. The first one consisted of imputed data according to a uniform distribution in the censoring interval, and the second one consisted of imputed data according to the NPMLE distribution. Those methods might be less time-consuming than the interval-censored procedure presented above, but they are not straightforward when implemented in SAS and R. Thus, we have decided not to present them in details, but to provide only the results to compare them with the two other methods we proposed.

For the sensitivity procedure, CIs of PFS estimation were computed on the basis of the log hazard, namely log(−log(PFS)).

Statistical analyses were done using the SAS software 9.3 (Copyright 2002–2010 SAS Institute Inc.) and the R-3.0.1 software (Copyright 2014 The R Foundation for Statistical Computing) with interval, survival, splines, perm, Icens, and MLEcens packages (22–28).

The dataset used for this article is available in Supplementary Material S1. SAS and R input and output are presented in Supplementary Materials S2 and S3. The source code of multiple imputations methods is available in Supplementary Material S4.

## Data Visualization and Description

After importing the data, it is necessary to first describe them in terms of number of patients with complete, right-censored, and interval-censored data for PFS.

### In SAS

**DATA WORK.MyData;**

** set WORK.MyData;**

** if T2>T1 then censoring=‘interval’;**

** if T2=. AND T1 NE. then censoring =‘right’;**

** if T2=T1 then censoring=‘none’;**

**run;**

**PROC FREQ DATA=WORK.MyData;**

** tables censoring;**

**run;**

### In R

**MyData <- within(MyData, censoring <-ifelse(is.na(T2) & !is.na(T1), ‘right’,**

** ifelse(T2!=T1, ‘interval’,**

** ifelse(T2==T1, ‘none’, NA))))**

**table(MyData$censoring)**

### Interpretation

The dataset is composed of three variables: the subject identifier, T1, and T2. Our dataset contains 53 (91%) patients with interval-censored data, 3 (5%) patients with right-censored data, and 2 (3%) patients with complete data for PFS estimation. In such a case, it is recommended to estimate PFS with methods dealing with interval-censored data, as described below.

## Interval-Censored Procedure

### In SAS

#### PFS estimation.

If version 13.1 of SAS/STAT in SAS 9.4 TS Level 1M1 is installed, the ICLIFETEST procedure may be used with the help of the well-described SAS support (29). Otherwise, you may use the EMICM macro, which may be downloaded from its help page (http://support.sas.com/kb/24/980.html). The macro has to be first compiled in SAS before running it with the data imported:

**%EMICM(data=WORK.MyData, left=T1, right=T2,**

**method=EMICM, out=intervalPFS, options=plot, mRS=100);**

Outputs should appear in the Results Viewer tab, and a new table “Intervalpfs” should be created in the work library, containing survival probabilities and SEs for each time interval defined with “Lower” and “Upper” columns.

#### PFS comparison.

The ICLIFETEST procedure provides a TEST statement to compare several groups with a weighted generalized log-rank test. Otherwise, the ICSTEST macro allows for comparison of survival curves using generalized log-rank tests as well (30, 31). The SAS support of ICLIFETEST procedure well described the difference between both (29).

### In R

Nonparametric estimation and comparison of interval-censored survival data are available in an add-on package called “interval.”

#### PFS estimation.

**intervalPFS <- icfit(Surv(T1, T2, type = "interval2") ∼ 1, data = MyData,**

** conf.int=TRUE, control = icfitControl(B=100,seed=1234))**

The CI is estimated by bootstrapping, which is a resampling method that can be very time-consuming. By default, the number of bootstrap samples is set at 200. To speed up the procedure, we reduced this number to 100 samples with the icfitControl function. Moreover, as bootstrapping relies on random sampling, the random seed is set to an arbitrary value so that the estimation of the CI does not change between two calls of the icfit function with the same entry values.

As the "intervalPFS" object contains only the event probability and not the cumulative survival probability, to see and plot the interval times with associated PFS estimation, call:

**intPFS <- (1-cumsum(intervalPFS$pf))**

**intPFS <-data.frame(time.interval=names(intPFS),PFS=intPFS,row.names=NULL)**

**intPFS**

**plot(intervalPFS)**

#### PFS comparison.

The package interval also provides a function to compare several groups, called "ictest" (22). Three different types of scores are available: the most commonly used log-rank scores of Sun (32), the scores associated with the grouped proportional hazards model of Finkelstein (33), and the generalized Wilcoxon–Mann–Whitney scores. The other options for scores only allow the permutation methods and follow cases where the error under the grouped continuous model is either normally distributed or distributed by some other distribution (34).

## Sensitivity Analysis Procedure

Two new variables (a status and a time variable) are computed to conduct this analysis: Right-censored patients are labeled as censored (status = 0), with the time of censoring being T1; deceased patients are labeled as not censored (status = 1), with the time of the event at the time of death; and progressive patients are labeled as not censored (status = 1), with the time of the event at (i) the midpoint between T1 and T2, (ii) the lower limit of the interval, or (iii) the upper limit of the interval for the three separate analyses.

### Midpoint of the interval

#### In SAS.

**DATA WORK.MyData;**

**set WORK.MyData;**

**if T2=. then do;**

** status=0;**

** midpoint=T1; end;**

**if T2 NE. then do;**

** status=1;**

** midpoint=(T2+T1)/2; end;**

**run;**

Now the LIFETEST procedure may be used to estimate the “midpoint” PFS over time:

** PROC LIFETEST DATA=WORK.MyData PLOTS=s OUTSURV=WORK.midpointPFS;**

** TIME midpoint*status(0);**

** run;**

Outputs of the PFS should appear in the Results Viewer tab and a new table “Midpointpfs” should be created in the work library, containing survival probabilities and a 95% CI for each “event” time defined in the “Midpoint” column.

#### In R.

**MyData <- within(MyData, status <- ifelse(is.na(T2), 0, 1))**

**MyData <- within(MyData, midpoint <- ifelse(is.na(T2), T1, (T2+T1)/2))**

Now the “survfit” function of the “survival” package may be used to estimate the “midpoint” PFS over time:

**midpointPFS <-survfit(Surv(midpoint, status)∼1, data = MyData, conf.type="log-log")**

Outputs of this “midpoint” PFS are now saved in the “midpointPFS” object. The estimation is summarized with:

**midpointPFS**

Then the PFS estimation may be plotted over time by calling:

**summary(midpointPFS)**

**plot(midpointPFS)**

### Upper limit of the interval

With this procedure considering progression at the date of first appearance, the status and event variables are the same as those already computed in the “midpoint procedure,” but the time variable is equal to T2 in the event of interval censoring:

#### In SAS.

**DATA WORK.MyData;**

**set WORK.MyData;**

** upper=midpoint;**

** if T2 NE. AND T1 NE. then upper=T2;**

**run;**

The LIFETEST procedure may then be used to estimate the upper PFS over time, as described in the midpoint procedure section, while replacing the “midpoint” variable by the “upper” variable.

#### In R.

**MyData <- within(MyData, upper <- ifelse(!is.na(T2) & !is.na(T1),**

** T2, midpoint))**

The “survfit” function of the “survival” package may then be used to estimate the upper PFS over time, as described in the midpoint procedure section, while replacing the “midpoint” variable by the “upper” variable.

### Lower limit of the interval

As detailed above, the time variable is now equal to T1 in the event of interval censoring.

## Output interpretation

### Interval-censored procedure

The shape of PFS assessed by an, NPMLE taking into account interval-censored data is similar to a Kaplan–Meier curve except for gray rectangles. In fact, as the exact PFS is unknown (no unique solution) in these periods of time, the rectangles display all the possibilities that PFS may be.

For each time interval, output tables returned by SAS and R display PFS probability. The main difference between them only relies on the display: In R output, there are less time intervals than in SAS output because consecutive time intervals with the same PFS probabilities are pooled together.

To estimate median PFS, the point where the survival probability reaches 0.5 must be established, namely during the time interval of 9.66 to 9.79 months. As the PFS is between 0.5329 and 0.3380 in this time interval, it is then simple to estimate that the median PFS is approximately reached at 9.68 months (Table 2).

*95% CI of median PFS*.

##### In SAS.

It is not easy to estimate a CI of the median PFS with the EMICM macro, but SEs of survival estimations are available for each time interval.

##### In R.

To find the 95% CI of this median, the same must be done in both columns of the 95% CI of survival over time (lower and upper bound).

– For the lower bound, one may call:

**lower<-with(intervalPFS$CI,{i<-which(lower==min(lower[lower>=0.5]));**

**i<-c(i,tail(i,1)+1); data.frame**

**(time=round(time[i],2), lower=lower[i])})**

The probability at 5.65 months falls from 0.53 to 0.48. So, in Fig. 1, the horizontal line representing a PFS of 0.5 (definition of median survival) cuts the lower line of the CI at this time.

– For the upper bound, the same code may be called while replacing lower by upper: The probability at 10.18 months is 0.55, and at 13.54 months it is 0.37. So, the probability 0.5 is reached within this time interval, at 11.07 months.

#### PFS estimation at a precise time.

One may also need to estimate PFS at a precise time. To do so, the interval to which this time belongs must be found. For example, at 12 months, which belongs to the time interval of 10.18 to 13.54 months, as seen in the R output, or to 11.79 to 12.65 months, as seen in the SAS output, the associated survival probability is 0.248. At 10 months, which is not within any interval, the associated survival probability is between 0.338 and 0.248.

### Sensitivity analysis procedure

In this procedure, PFS is assessed using the standard Kaplan–Meier estimator. The median PFS is thus directly estimated by SAS and R (Table 2). To estimate PFS at a precise time, the nearest time that is lower than or equal to the target time must be found in the output table. For example, at 12 months, midpoint PFS is the same as at 11.385 months, namely 0.259 [95% CI (0.155–0.375)].

As described above for the midpoint PFS procedure, the estimation of the median PFS while considering the upper limit of the censoring interval is 10.215 months (95% CI, 8.310–11.530; Table 2). At 12 months, PFS is 0.310 (95% CI, 0.197–0.430). The estimation of the median PFS while considering the lower limit of the censoring interval is 7.62 months, (95% CI, 5.75–9.23; Table 2). At 12 months, the PFS is 0.207 (95% CI, 0.114–0.318).

### Multiple imputation procedures

Both multiple imputation procedures gave close results (Table 2), so only the multiple imputation using the NPMLE distribution is presented in Fig. 1. We observed that the CI obtained with the multiple imputation procedure was narrower than the one obtained with the bootstrap resampling of the EM-ICM algorithm. In fact, as detailed by Hsu and colleagues, even with the multiple imputation based on the NPMLE, not all the uncertainty is taken into account, because those imputations are conditional on the data distribution and may not be viewed as a perfect multiple imputation procedure (35). Hsu and colleagues also propose adding a bootstrap stage before the imputation stage, but that makes it more time-consuming, thus losing this advantage.

## Discussion

In this article, we explain how to estimate PFS by using three different methods with interval-censored data using SAS and R statistical software and applied to a real dataset of patients with metastatic colorectal cancer enrolled in a phase II clinical trial.

As expected, the standard procedure (which is the procedure considering progression at the lower limit of the censoring interval) overestimates PFS in this example. In clinical trials with short time intervals between visits, this error may be negligible. The “midpoint procedure” is an easy-to-use method, but its scope of use is limited when used alone (36). We have, therefore, proposed a sensitivity analysis with the three systematic imputations: midpoint, upper limit, and lower limit of the interval, as recommended by Panageas and colleagues (21). The results obtained with the sensitivity analysis are similar to the results obtained with the other two procedures, however more complex.

The three procedures proposed here (the sensitivity analysis procedure, the interval-censored procedures, and the multiple imputation methods) require that the time of last visit before progression of each patient has to be recorded in the dataset, which is not always the case in the case report form. Other methods, such as the maximum bias, exist to analyze the worst-case scenario.

Beyond the wish to estimate PFS without bias, all the uncertainty in its variance estimate has to be taken into account. One solution is to add a bootstrap resampling stage, but it makes the procedures more time-consuming. Semiparametric maximum likelihood estimators or methods for grouped data (where all subjects follow the same assessment schedule) may be preferred in the case of large dataset (4, 37, 38).

Finally, caution is required when dealing with cancer types where progression is not only assessed by radiologic techniques. In settings such as gastrointestinal stromal tumor, progression may be established during a clinical examination, so the interval to which the exact time of progression belongs may be unknown. However, it would certainly be shorter than with radiologic exams. Thus, we do not recommend the use of methods dealing with interval-censored data when these intervals are unknown or when the censoring mechanism is not independent of the event process.

As the dataset used here is available from the Supplementary Material, readers may practice before using such methods on their own data.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**Conception and design:** A.E. Dugué

**Development of methodology:** A.E. Dugué

**Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.):** S. Chabaud

**Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis):** A.E. Dugué, J. Gal

**Writing, review, and/or revision of the manuscript:** A.E. Dugué, M. Pulido, S. Chabaud, L. Belin, J. Gal

**Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases):** A.E. Dugué

## Grant Support

This study was supported by François Baclesse Cancer Centre (Caen, France).

## Acknowledgments

We thank Dr. Yves Becouarn from the Bergonié Cancer Institute for providing us the dataset of his clinical trial and Ray Cooke for his help in language and content editing of this manuscript. We thank the anonymous referees and the senior editor for their thoughtful comments, which have helped improve the presentation of this article.

## Footnotes

**Note:**Supplementary data for this article are available at Clinical Cancer Research Online (http://clincancerres.aacrjournals.org/).

- Received April 21, 2016.
- Revision received August 4, 2016.
- Accepted August 12, 2016.

- ©2016 American Association for Cancer Research.