---
title: "RPEXE.RPEXT"
author: "Gang Han, Yu Zhang"
date: "`r Sys.Date()`"
output: html_document
vignette: >
  %\VignetteIndexEntry{RPEXE.RPEXT}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

# 1 General Information

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.

# 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;
    
    <br/>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
```{r}
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

```{r}
# figure(1): Kaplan Meier curve of Arm A without indicating censored points 
km(times[armsA_ID], censor[armsA_ID], 0)
```
```{r}
# figure(2): Kaplan Meier curve of armA with censored points indicated
km_red(times[armsA_ID], censor[armsA_ID], 1)
```
```{r}
# figure(3): Kaplan Meier curve of armB without indicating censored points 
km(times[armsB_ID], censor[armsB_ID], 0)
```
```{r}
# 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
```{r}
# 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
```{r}
# 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)
```

## 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.
```{r}
# 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 
```{r}
# 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
```

# References