This reduced piecewise exponential survival software implements the likelihood ratio test and backward elimination procedure in Han, Schell, and Kim (2012^[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.], 2014^[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.]), and Han et al. (2016^[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.]). 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^[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.]) also mentioned the application of the method. This program can run in R version 3.2.2 and above.
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).
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.
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)
# 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)
# 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)
# Fit the rpexe with monotonic order restriction; pexeoutA = RPEXEv1_2(times[armsA_ID],censor[armsA_ID], monotone = 1,criticalp = 0.05) pexeoutB = RPEXEv1_2(times[armsB_ID],censor[armsB_ID],monotone = 1,criticalp = 0.05) # combined pexeout = RPEXEv1_2(times,censor,monotone = 1,criticalp = 0.05)
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 # second piece result=exact_pvalue(ttotA2,ttotB2,dA2,dB2,0) a12 = result[1] p12 = result[2] p12 # third piece result=exact_pvalue(ttotA3,ttotB3,dA3,dB3,0) a13 = result[1] p13 = result[2] p13 # One-sided test # first piece result=exact_pvalue(ttotA1,ttotB1,dA1,dB1,1) a21 = result[1] p21 = result[2] p21 # second piece result=exact_pvalue(ttotA2,ttotB2,dA2,dB2,1) a22 = result[1] p22 = result[2] p22 # third piece result=exact_pvalue(ttotA3,ttotB3,dA3,dB3,1) a23 = result[1] p23 = result[2] p23
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.