1 General Information

This reduced piecewise exponential survival software implements the likelihood ratio test and backward elimination procedure in Han, Schell, and Kim (2012¹, 2014²), and Han et al. (2016³). Inputs to the program can be either times when events/censoring occur or the vectors of total time on test and the number of events. Outputs of the programs are times and the corresponding p-values in the backward elimination. Details about the model and implementation are given in Han et al. 2014. Adelson (2016⁴) also mentioned the application of the method. This program can run in R version 3.2.2 and above.

2 Inputs and Outputs

2.1 Inputs

This software has one driver files RPEXEv1_2.R. Inputs to RPEXEv1_2.R include

‘EventTime’ = time, a vector having the times of the events occurred and of censoring. Note that ’EventTime’ is a required input.
‘Censor’ = censor, a vector with 0 or 1 indicating whether an observation is an event or censored. We use 0 to denote censoring and 1 to denote event. The length of time and censor are identical. Note that ’Censor’ is a required input.
‘CutTimes’ = cuttimes, a vector of the potential times where the piecewise exponential model is divided. Note that ‘CutTimes’ is an optional input. By default, cuttimes is a vector of times when the events occur.
‘Monotone’ = monotone, an indicator with values and meanings
- 0: no monotonic assumption;
- 1: assuming that the hazard is decreasing over time;
- 2: assuming that the hazard is increasing over time;
- 3: assuming that the hazard is monotonic;
- 4: assuming that the hazard is increasing and then decreasing;
- 5: assuming that the hazard is decreasing and then decreasing;
- 6: assuming that the hazard is increasing and then decreasing with the peak removed first;
- 7: assuming that the hazard is decreasing and then increasing with the peak removed first;
Note that ‘Monotone’ is an optional input. The default value of ‘Monotone’ is 0.
‘Criticalp’ = The critical (naive) p-value cutoff where all p-values in the backward elimination that are lower than this will be regarded as being significant. The prediction of the survival probability will be made on 100 equally spaced time points within the range of the event times based on the piecewise exponential estimate determined by all the changepoints. Default == -1 (equivalent to NA).

2.2 Outputs

times: times to make the cuts
pvalues: pvalues correspond to the times
times_c: critical times to make the cuts
pvalues_c: critical p-values that are smaller than the
trend: trend information
struct: structure information for multiple order restrictions
changet: change point in time for umbrella alternatives.

3 Example with dataset data2

3.1 Load data and extract variables

library(RPEXE.RPEXT)
data(data2)
times = data2[,1]
censor = data2[,2]
group = data2[,3]

ID_nan = which(is.na(times))
times = times[-ID_nan] 
censor = censor[-ID_nan]
group = group[-ID_nan]
armsA_ID = which(group == 1)
armsB_ID = which(group == 2)

3.2 Plot the data

# figure(1): Kaplan Meier curve of Arm A without indicating censored points 
km(times[armsA_ID], censor[armsA_ID], 0)

# figure(2): Kaplan Meier curve of armA with censored points indicated
km_red(times[armsA_ID], censor[armsA_ID], 1)

# figure(3): Kaplan Meier curve of armB without indicating censored points 
km(times[armsB_ID], censor[armsB_ID], 0)

# figure(4): Kaplan Meier curve of Arm B with censored points indicated
km_red(times[armsB_ID], censor[armsB_ID], 0)

3.3 Plot the KME in the paper

# figure(5) : Combined plot of both armA and armB 
x1 = cbind(times[armsA_ID], censor[armsA_ID])
x2 = cbind(times[armsB_ID], censor[armsB_ID])
km_combine(x1,x2)

3.4 The reduced piecewise exponential analysis

# Fit the rpexe with monotonic order restriction;
pexeoutA     =  RPEXEv1_2(times[armsA_ID],censor[armsA_ID], monotone = 1,criticalp = 0.05)

## [1] "last"
## [1] 3

pexeoutB     =  RPEXEv1_2(times[armsB_ID],censor[armsB_ID],monotone = 1,criticalp = 0.05)

## [1] "last"
## [1] 4

# combined
pexeout = RPEXEv1_2(times,censor,monotone = 1,criticalp = 0.05)

## [1] "last"
## [1] 4

3.5 Additional analysis

Given the RPEXE estimates, using the total time on test and number of events to compare the two arms where the hazard rates are costant.

# Calculate the ttot and n from a), 0-2.777, b), 2.777-8.959, c), 8,959-end;

returnvA=totaltest(times[armsA_ID],censor[armsA_ID]) 
m=dim(returnvA)[2]/3
time_dieA=returnvA[,1:m]
ttotA=returnvA[,(m+1):(2*m)]
deathsA=returnvA[,(2*m+1):3*m]

returnvB=totaltest(times[armsB_ID],censor[armsB_ID]) 

m=dim(returnvB)[2]/3

time_dieB=returnvB[,1:m]
ttotB=returnvB[,(m+1):(2*m)]
deathsB=returnvB[,(2*m+1):3*m]

ttotA1 = 0
ttotA2 = 0
ttotA3 = 0
dA1 = 0
dA2 = 0
dA3 = 0
for (i in 1:length(time_dieA))
{
  if ( time_dieA[i]<=2.777)
  {
    ttotA1 = ttotA1+ttotA[i]
    dA1    = dA1+deathsA[i]
  }else if (time_dieA[i]<=8.959)
  {
    ttotA2 = ttotA2+ttotA[i]
    dA2    = dA2+deathsA[i]
  } else 
  {
    ttotA3 = ttotA3+ttotA[i]
    dA3    = dA3+deathsA[i]
  }
}
      
ttotB1 = 0
ttotB2 = 0
ttotB3 = 0
dB1 = 0
dB2 = 0
dB3 = 0
for (i in 1:length(time_dieB))
{
  if ( time_dieB[i]<=2.777)
  {
    ttotB1 = ttotB1+ttotB[i]
    dB1    = dB1+deathsB[i]
  }else if (time_dieB[i]<=8.959)
  {
    ttotB2 = ttotB2+ttotB[i]
    dB2    = dB2+deathsB[i]
  } else 
  {
    ttotB3 = ttotB3+ttotB[i]
    dB3    = dB3+deathsB[i]
  }
}

Compute the test statistic and p-value. Show the p-values

# Test the two side hypothesis;

# Two-sided test
# first piece
result=exact_pvalue(ttotA1,ttotB1,dA1,dB1,0)
a11 = result[1]
p11 = result[2]
p11

## [1] 0.7960157

# second piece
result=exact_pvalue(ttotA2,ttotB2,dA2,dB2,0)
a12 = result[1]
p12 = result[2]
p12

## [1] 0.05638443

# third piece
result=exact_pvalue(ttotA3,ttotB3,dA3,dB3,0)
a13 = result[1]
p13 = result[2]
p13

## [1] 0.5973843

# One-sided test
# first piece
result=exact_pvalue(ttotA1,ttotB1,dA1,dB1,1)
a21 = result[1]
p21 = result[2]
p21

## [1] 0.7960157

# second piece
result=exact_pvalue(ttotA2,ttotB2,dA2,dB2,1)
a22 = result[1]
p22 = result[2]
p22

## [1] 0.05139599

# third piece
result=exact_pvalue(ttotA3,ttotB3,dA3,dB3,1)
a23 = result[1]
p23 = result[2]
p23

## [1] 1.371367

References

Han, G., Schell, M. J., and Kim, J. (2012) “Comparing Two Exponential Distributions Using the Exact Likelihood Ratio Test,” Statistics in Biopharmaceutical Research, 4(4), 348-356.↩
Han, G., Schell, M. J., and Kim, J. (2014) “Improved Survival Modeling in Cancer Research Using a Reduced Piecewise Exponential Approach,” Statistics in Medicine, 33(1), 59-73.↩
Han, G., Schell, M., Zhang, H., Zelterman, D., Pusztai, L., Adelson, K., and Hatzis, C. (2016) “Testing Violations of the Exponential Assumption in Cancer Clinical Trials with Survival Endpoints,” Biometrics, DOI: 10.1111/biom.12590; PMID: 27669414.↩
Adelson, K. B., Ramaswamy, B., Sparano, J. A., Christos, P. J., Wright, J. J., Raptis, G., Han, G., Villalona-Calero, M., Ma, C., Hershman, D., Baar, J., Klein, P., Cigler, T., Budd, T., Novik, Y., Tan, A.R., Tannenbaum, S., Goel, A., Levine, E., Shapiro, C. L., Andreopoulou, E., Naughton, M., Kalinsky, K., Waxman, S., Germain, D. (2016) “Randomized Phase II Trial of Fulvestrant Alone or in Combination with Bortezomib in Hormone Receptor-Positive Metastatic Breast Cancer Resistant to Aromatase Inhibitors: A New York Cancer Consortium Trial,” Nature Partner Journals Breast Cancer, Volume 2, Article ID 16037, DOI: 10.1038/npjbcancer.2016.37.↩

RPEXE.RPEXT

Gang Han, Yu Zhang

2020-05-20