#This S-Plus program will produce Figure 5 which contains the relative times of period 2, 3, and 4 to period 1 for the
#best fitting two-parameter distribution in the GG family; 90% confidence bands are also shown.
#Boxplots of relative times for p=0.1, 0.2, and 0.3 based on 500 bootstrap samples for period2,3,4 to 1 are also shown.

#3 data files (P2 vs. P1, P3 vs. P1, P4 vs. P1) with five variables: period, percentile, relative time, 
#  lower limit of confidence band, upper limit of confidence band
p2p1_matrix(scan("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/RelTime2vs1FIG5.dat"), ncol=5, byrow=T)
p3p1_matrix(scan("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/RelTime3vs1FIG5.dat"), ncol=5, byrow=T)
p4p1_matrix(scan("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/RelTime4vs1FIG5.dat"), ncol=5, byrow=T)

#Change path below for postscript file location
postscript("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/GGfig5.ps", horizontal=T)

period_p2p1[,1]
cumprob_100 * p2p1[,2]
#relative times for period 2 to period 1 along with lower and upper limits of 90% confidence bands for different values of p
reltimep2p1_p2p1[,3]
reltimep2p1lb_p2p1[,4]
reltimep2p1ub_p2p1[,5]

period_p3p1[,1]
cumprob_100 * p3p1[,2]
#relative times for period 3 to period 1 along with lower and upper limits of 90% confidence bands for different values of p
reltimep3p1_p3p1[,3]
reltimep3p1lb_p3p1[,4]
reltimep3p1ub_p3p1[,5]

period_p4p1[,1]
cumprob_100 * p4p1[,2]
#relative times for period 4 to period 1 along with lower and upper limits of 90% confidence bands for different values of p
reltimep4p1_p4p1[,3]
reltimep4p1lb_p4p1[,4]
reltimep4p1ub_p4p1[,5]

par(mfrow=c(1,3))
par(mgp=c(1.8, 0.8, 0), mar=c(7, 5, 4, 2))

############################################################################################################
############################################################################################################
######################################### P2 versus P1 #####################################################
############################################################################################################
############################################################################################################

#Gamma Period 2 and Ammag in Period 1
#draw three boxplots representing the bootstrap distribution of relative times for p=0.1, 0.2, and 0.3 respectively
#p=0.1
x95to75gg10 <- c( 9.8, 9.8)
y95to75gg10 <- c(1.7789, 1.5952)
x75gg10 <- c(5.5, 14.5)
y75gg10 <- c(1.5952, 1.5952)
xleftgg10 <- c(5.5, 5.5)
yleftgg10 <- c(1.5952,1.3705)
x25gg10 <- c(5.5, 14.5)
y25gg10 <- c(1.3705, 1.3705)
xrightgg10 <- c(14.5, 14.5)
yrightgg10 <- c(1.5952, 1.3705)
x25to5gg10 <- c(9.8, 9.8)
y25to5gg10 <- c(1.3705, 1.2599)
x50gg10 <- c(5.5, 14.5)
y50gg10 <- c(1.4792, 1.4792)

plot(x95to75gg10,y95to75gg10,log='y',xlim=c(5,35),ylim=c(0.8,15),xlab="",ylab="",axes=F, type="n")
title(main="\n Period 2 to Period 1", cex=0.8)	
title(ylab="Relative Time", cex=1.1)	
axis(2,at=seq(5,15,5),label=c("5","10","15"),cex=0.80,tck=-0.01)
axis(2,at=seq(1,4,1),label=c("1","2","3","4"),cex=0.80,tck=-0.01)
axis(1,at=seq(10,30,10),label=c("10", "20", "30"),cex=0.8)
abline(h=1,lty=2)
lines(x95to75gg10,y95to75gg10, lty=1, lwd=1.5)
lines(x75gg10, y75gg10, lty=1, lwd=1.5)
lines(xleftgg10, yleftgg10, lty=1, lwd=1.5)
lines(x25gg10, y25gg10, lty=1, lwd=1.5)
lines(xrightgg10, yrightgg10, lty=1, lwd=1.5)
lines(x25to5gg10, y25to5gg10, lty=1, lwd=1.5)
lines(x50gg10, y50gg10, lty=1, lwd=1.5)


#p=0.2
x95to75gg20 <- c(19.8, 19.8)
y95to75gg20 <- c(1.4984, 1.3743)
x75gg20 <- c(15.5, 24.5)
y75gg20 <- c(1.3743, 1.3743)
xleftgg20 <- c(15.5, 15.5)
yleftgg20 <- c(1.3743, 1.2165)
x25gg20 <- c(15.5, 24.5)
y25gg20 <- c(1.2165, 1.2165)
xrightgg20 <- c(24.5, 24.5)
yrightgg20 <- c(1.3743, 1.2165)
x25to5gg20 <- c(19.8, 19.8)
y25to5gg20 <- c(1.2165, 1.1279)
x50gg20 <- c(15.5, 24.5)
y50gg20 <- c(1.2912, 1.2912)

lines(x95to75gg20,y95to75gg20, lty=1, lwd=1.5)
lines(x75gg20, y75gg20, lty=1, lwd=1.5)
lines(xleftgg20, yleftgg20, lty=1, lwd=1.5)
lines(x25gg20, y25gg20, lty=1, lwd=1.5)
lines(xrightgg20, yrightgg20, lty=1, lwd=1.5)
lines(x25to5gg20, y25to5gg20, lty=1, lwd=1.5)
lines(x50gg20, y50gg20, lty=1, lwd=1.5)

#p=0.3
x95to75gg30 <- c(29.8, 29.8)
y95to75gg30 <- c(1.3697, 1.2756)
x75gg30 <- c(25.5, 34.5)
y75gg30 <- c(1.2756, 1.2756)
xleftgg30 <- c(25.5, 25.5)
yleftgg30 <- c(1.2756, 1.1545)
x25gg30 <- c(25.5, 34.5)
y25gg30 <- c(1.1545, 1.1545)
xrightgg30 <- c(34.5, 34.5)
yrightgg30 <- c(1.2756, 1.1545)
x25to5gg30 <- c(29.8, 29.8)
y25to5gg30 <- c(1.1545,1.0816)
x50gg30 <- c(25.5, 34.5)
y50gg30 <- c(1.2112,1.2112)

lines(x95to75gg30,y95to75gg30, lty=1, lwd=1.5)
lines(x75gg30, y75gg30, lty=1, lwd=1.5)
lines(xleftgg30, yleftgg30, lty=1, lwd=1.5)
lines(x25gg30, y25gg30, lty=1, lwd=1.5)
lines(xrightgg30, yrightgg30, lty=1, lwd=1.5)
lines(x25to5gg30, y25to5gg30, lty=1, lwd=1.5)
lines(x50gg30, y50gg30, lty=1, lwd=1.5)

lines(cumprob[cumprob>=5],reltimep2p1[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep2p1lb[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep2p1ub[cumprob>=5],lty=1,lwd=1.5)

box()

##########################################################################################################
##########################################################################################################
######################################### P3 versus P1 ###################################################
##########################################################################################################
##########################################################################################################

#Ammag Period 3 and Ammag Period 1
#three boxplots representing the bootstrap distribution of relative times for p=0.1, 0.2, and 0.3 respectively
#p=0.1
x95to75gg10 <- c( 9.8, 9.8)
y95to75gg10 <- c(2.6897, 2.0406)
x75gg10 <- c(5.5, 14.5)
y75gg10 <- c(2.0406, 2.0406)
xleftgg10 <- c(5.5, 5.5)
yleftgg10 <- c(2.0406, 1.3100)
x25gg10 <- c(5.5, 14.5)
y25gg10 <- c(1.3100, 1.3100)
xrightgg10 <- c(14.5, 14.5)
yrightgg10 <- c(2.0406, 1.3100)
x25to5gg10 <- c(9.8, 9.8)
y25to5gg10 <- c(1.3100, 0.9730)
x50gg10 <- c(5.5, 14.5)
y50gg10 <- c(1.6356, 1.6356)

plot(x95to75gg10,y95to75gg10,log='y',xlim=c(5,35),ylim=c(0.8,15),xlab="",ylab="",axes=F,type="n")
axis(2,at=seq(5,15,5),label=c("5","10","15"),cex=0.80,tck=-0.01)
axis(2,at=seq(1,4,1),label=c("1","2","3","4"),cex=0.80,tck=-0.01)
axis(1,at=seq(10,30,10),label=c("10", "20", "30"),cex=0.8)
title(main="\n Period 3 to Period 1", cex=0.8)	
key(2, exp(-0.5), text="Cummulative Percent Deseased", cex=1.0, adj=0)	
key(-6, exp(-0.65), text="(Deciles for Period 1: 0.23, 0.47, 0.72 years)", cex=1.0, adj=0)	
abline(h=1,lty=2)
lines(x95to75gg10,y95to75gg10, lty=1, lwd=1.5)
lines(x75gg10, y75gg10, lty=1, lwd=1.5)
lines(xleftgg10, yleftgg10, lty=1, lwd=1.5)
lines(x25gg10, y25gg10, lty=1, lwd=1.5)
lines(xrightgg10, yrightgg10, lty=1, lwd=1.5)
lines(x25to5gg10, y25to5gg10, lty=1, lwd=1.5)
lines(x50gg10, y50gg10, lty=1, lwd=1.5)


#p=0.2
x95to75gg20 <- c(19.8, 19.8)
y95to75gg20 <- c(3.5001, 2.8939)
x75gg20 <- c(15.5, 24.5)
y75gg20 <- c(2.8939, 2.8939)
xleftgg20 <- c(15.5, 15.5)
yleftgg20 <- c(2.8939, 2.1220)
x25gg20 <- c(15.5, 24.5)
y25gg20 <- c(2.1220, 2.1220)
xrightgg20 <- c(24.5, 24.5)
yrightgg20 <- c(2.8939, 2.1220)
x25to5gg20 <- c(19.8, 19.8)
y25to5gg20 <- c(2.1220, 1.7229)
x50gg20 <- c(15.5, 24.5)
y50gg20 <- c(2.4927, 2.4927)

lines(x95to75gg20,y95to75gg20, lty=1, lwd=1.5)
lines(x75gg20, y75gg20, lty=1, lwd=1.5)
lines(xleftgg20, yleftgg20, lty=1, lwd=1.5)
lines(x25gg20, y25gg20, lty=1, lwd=1.5)
lines(xrightgg20, yrightgg20, lty=1, lwd=1.5)
lines(x25to5gg20, y25to5gg20, lty=1, lwd=1.5)
lines(x50gg20, y50gg20, lty=1, lwd=1.5)


#p=0.3
x95to75gg30 <- c(29.8, 29.8)
y95to75gg30 <- c(4.6420, 3.9188)
x75gg30 <- c(25.5, 34.5)
y75gg30 <- c(3.9188, 3.9188)
xleftgg30 <- c(25.5, 25.5)
yleftgg30 <- c(3.9188, 3.0697)
x25gg30 <- c(25.5, 34.5)
y25gg30 <- c(3.0697, 3.0697)
xrightgg30 <- c(34.5, 34.5)
yrightgg30 <- c(3.9188, 3.0697)
x25to5gg30 <- c(29.8, 29.8)
y25to5gg30 <- c(3.0697, 2.5380)
x50gg30 <- c(25.5, 34.5)
y50gg30 <- c(3.4792, 3.4792)

lines(x95to75gg30,y95to75gg30, lty=1, lwd=1.5)
lines(x75gg30, y75gg30, lty=1, lwd=1.5)
lines(xleftgg30, yleftgg30, lty=1, lwd=1.5)
lines(x25gg30, y25gg30, lty=1, lwd=1.5)
lines(xrightgg30, yrightgg30, lty=1, lwd=1.5)
lines(x25to5gg30, y25to5gg30, lty=1, lwd=1.5)
lines(x50gg30, y50gg30, lty=1, lwd=1.5)

lines(cumprob[cumprob>=5],reltimep3p1[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep3p1lb[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep3p1ub[cumprob>=5],lty=1,lwd=1.5)

box()

############################################################################################################
############################################################################################################
######################################### P4 versus P1 #####################################################
############################################################################################################
############################################################################################################

#Weibull Period 4 and Ammag Period 1
#three boxplots representing the bootstrap distribution of relative times for p=0.1, 0.2, and 0.3 respectively
#p=0.1
x95to75gg10 <- c( 9.8, 9.8)
y95to75gg10 <- c(10.2052, 7.2077)
x75gg10 <- c(5.5, 14.5)
y75gg10 <- c(7.2077, 7.2077)
xleftgg10 <- c(5.5, 5.5)
yleftgg10 <- c(7.2077, 3.9641)
x25gg10 <- c(5.5, 14.5)
y25gg10 <- c(3.9641, 3.9641)
xrightgg10 <- c(14.5, 14.5)
yrightgg10 <- c(7.2077, 3.9641)
x25to5gg10 <- c(9.8, 9.8)
y25to5gg10 <- c(3.9641, 2.1391)
x50gg10 <- c(5.5, 14.5)
y50gg10 <- c(5.3373, 5.3373)

plot(x95to75gg10,y95to75gg10,log='y',xlim=c(5,35),ylim=c(0.8,15),xlab="",ylab="",axes=F,type="n")
axis(2,at=seq(5,15,5),label=c("5","10","15"),cex=0.80,tck=-0.01)
axis(2,at=seq(1,4,1),label=c("1","2","3","4"),cex=0.80,tck=-0.01)
axis(1,at=seq(10,30,10),label=c("10", "20", "30"),cex=0.8)
title(main="\n Period 4 to Period 1", cex=0.8)	
abline(h=1,lty=2)
lines(x95to75gg10,y95to75gg10, lty=1, lwd=1.5)
lines(x75gg10, y75gg10, lty=1, lwd=1.5)
lines(xleftgg10, yleftgg10, lty=1, lwd=1.5)
lines(x25gg10, y25gg10, lty=1, lwd=1.5)
lines(xrightgg10, yrightgg10, lty=1, lwd=1.5)
lines(x25to5gg10, y25to5gg10, lty=1, lwd=1.5)
lines(x50gg10, y50gg10, lty=1, lwd=1.5)


#p=0.2
x95to75gg20 <- c(19.8, 19.8)
y95to75gg20 <- c(10.4848, 8.3940)
x75gg20 <- c(15.5, 24.5)
y75gg20 <- c(8.3940, 8.3940)
xleftgg20 <- c(15.5, 15.5)
yleftgg20 <- c(8.3940, 5.8297)
x25gg20 <- c(15.5, 24.5)
y25gg20 <- c(5.8297, 5.8297)
xrightgg20 <- c(24.5, 24.5)
yrightgg20 <- c(8.3940, 5.8297)
x25to5gg20 <- c(19.8, 19.8)
y25to5gg20 <- c(5.8297, 3.8537)
x50gg20 <- c(15.5, 24.5)
y50gg20 <- c(6.9108, 6.9108)

lines(x95to75gg20,y95to75gg20, lty=1, lwd=1.5)
lines(x75gg20, y75gg20, lty=1, lwd=1.5)
lines(xleftgg20, yleftgg20, lty=1, lwd=1.5)
lines(x25gg20, y25gg20, lty=1, lwd=1.5)
lines(xrightgg20, yrightgg20, lty=1, lwd=1.5)
lines(x25to5gg20, y25to5gg20, lty=1, lwd=1.5)
lines(x50gg20, y50gg20, lty=1, lwd=1.5)

#p=0.3
x95to75gg30 <- c(29.8, 29.8)
y95to75gg30 <- c(11.1717, 9.3637)
x75gg30 <- c(25.5, 34.5)
y75gg30 <- c(9.3637, 9.3637)
xleftgg30 <- c(25.5, 25.5)
yleftgg30 <- c(9.3637, 7.3225)
x25gg30 <- c(25.5, 34.5)
y25gg30 <- c(7.3255, 7.3255)
xrightgg30 <- c(34.5, 34.5)
yrightgg30 <- c(9.3637, 7.3225)
x25to5gg30 <- c(29.8, 29.8)
y25to5gg30 <- c(7.3225, 5.5503)
x50gg30 <- c(25.5, 34.5)
y50gg30 <- c(8.2301, 8.2301)

lines(x95to75gg30,y95to75gg30, lty=1, lwd=1.5)
lines(x75gg30, y75gg30, lty=1, lwd=1.5)
lines(xleftgg30, yleftgg30, lty=1, lwd=1.5)
lines(x25gg30, y25gg30, lty=1, lwd=1.5)
lines(xrightgg30, yrightgg30, lty=1, lwd=1.5)
lines(x25to5gg30, y25to5gg30, lty=1, lwd=1.5)
lines(x50gg30, y50gg30, lty=1, lwd=1.5)

lines(cumprob[cumprob>=5],reltimep4p1[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep4p1lb[cumprob>=5],lty=1,lwd=1.5)
lines(cumprob[cumprob>=5],reltimep4p1ub[cumprob>=5],lty=1,lwd=1.5)

box()
dev.off()