#This S-Plus program will produce Figure 4; left panel contains the graphical representation of the estimated 
#scale paramater (closed symbols) for the best fitting two-parameter model in the GG family for each of the four periods, also
#shown are the estimated 5th and 95th percentiles of the bootstrap distribution of the estimate.  The open symbols
#represent the parameter estimates for the saturated GG model.  The right panel contains the hazard function for the
#best fitting two-parameter distribution in each of the four periods.
library(postscriptfonts)
ps.options(setfont=ps.setfont.latin1)
ps.options( reset=T, colors = ps.colors.rgb[c("black", "green", "blue"),] )
options(mixed.text.screen=T)

#########################################################################################################
##########################################  Left Panel   ################################################
#########################################################################################################

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

par(mgp=c(2, 1.5, 0),mfrow=c(1,2),mar=c(5,6,4,2))

x<-c(0)
y<-c(0)

plot(x,y,xlim=c(0,3),ylim=c(0,2),xlab="",ylab="",axes=F,type="n",yaxs='s', xaxs='s',cex=1.2)
axis(2,at=seq(0,2,0.5),label=c("0.0","0.5","1.0","1.5","2.0"),cex=1.0,tck=-0.01,lwd=3.5)
par(mgp=c(2, 0.7, 0))
axis(1,at=seq(0,3,0.5),label=c("0.0","0.5","1.0","1.5","2.0","2.5","3.0"),cex=1.0,lwd=3.5)
par(col=1, cex=1.3)
mixed.text.vector(1.2, -0.17, c(1, 13, 1), c("Scale (", "\s", ")") )
mixed.text.vector(-0.45, 0.8, c(1, 13, 1), c("Shape (", "\l", ")"), srt=90)
mixed.text(x=1.6, y=1.7, "Gamma", slope=1 )
mixed.text(x=0.6, y=1.95, "Ammag", slope=-2.1 )

#Period 1: Ammag
ttx_seq(0,3,.01)
tty_ 1 / ttx[ttx>0]
#0.686 and 0.911 are 5th and 95th percentiles of estimated sigma in period 1
lines(ttx[ttx >= 0.686 & ttx < 0.911],tty[ttx >= 0.686 & ttx < 0.911],lwd=6,lty=1)
points(0.798,1/0.798,pch=17,cex=1.7)
points(0.773,1.369,pch=2,cex=1.7)
text(0.4,1.3,"Period 1",cex=1.2)
text(0.4,1.21,"Q1=0.59",cex=1.05)

#Period 2: Gamma
ciper2x <- c(0,0.806)
ciper2y <- c(0,0.806)
lines(ciper2x,ciper2y,lwd=1,lty=1)
llper2x <- c(0.806,0.890)
llper2y <- c(0.806,0.890)
lines(llper2x,llper2y,lwd=6,lty=1)
ulper2x <- c(0.890,1.98)
ulper2y <- c(0.890,1.98)
#0.806 and 0.890 are 5th and 95th percentiles of estimated sigma in period 2
lines(ulper2x,ulper2y,lwd=1,lty=1)
points(0.847,0.851,pch=0,cex=1.7)
points(0.848,0.848,pch=15,cex=1.2)
#Closed box for period 2 is not plotted since it will cover estimate of sigma and lambda for saturated general gamma model
text(0.4,0.9,"Period 2",cex=1.2)
text(0.4,0.81,"Q1=0.73",cex=1.05)


#Period 3: Ammag
lines(ttx[ttx >= 0.500 & ttx < 1.534],tty[ttx >= 0.500 & ttx < 1.534],lwd=1,lty=1)
#1.534 and 2.602 are 5th and 95th percentiles of estimated sigma in period 3
lines(ttx[ttx >= 1.534 & ttx < 2.602],tty[ttx >= 1.534 & ttx < 2.602],lwd=6,lty=1)
lines(ttx[ttx >= 2.602 & ttx < 3.000],tty[ttx >= 2.602 & ttx < 3.000],lwd=1,lty=1)
points(2.082,1/2.082,pch=16,cex=1.7)
points(2.008,0.5490,pch=1,cex=1.7)
text(2.55,0.6,"Period 3",cex=1.2)
text(2.55,0.51,"Q1=1.83",cex=1.05)


#Period 4: Weibull
llper4x <- c(0.950,1.803)
llper4y <- c(1.000,1.000)
#0.950 and 1.803 are 5th and 95th percentiles of estimated sigma in period 4
lines(llper4x,llper4y,lwd=6,lty=1)
points(1.292,1,pch=18,cex=1.7)
points(1.376,0.845,pch=5,cex=1.7)
text(2.2,1.05,"Period 4",cex=1.2)
text(2.2,0.96,"Q1=4.46",cex=1.05)

#########################################################################################################
##########################################  Right Panel  ################################################
#########################################################################################################

#computes the hazard function for the Generalized Gamma distribution
hggamma <- function(beta, sigma, lambda, timet)
{

s <- lambda^(-2)
r <- s * exp(( - beta * lambda)/sigma)

# consider timet below exp(-500) at exp(-500) to avoid division by zero
timet[timet<exp(-500)] <- exp(-500)

b <- lambda/sigma
haz <- (lambda>0)*((b*(timet^(b-1))*dgamma(timet^b,shape=s,rate=r))/(1-pgamma(timet^b,shape=s,rate=r)))
haz <- haz+(lambda<0)*((-b*(timet^(b-1))*dgamma(timet^b,shape=s,rate=r))/pgamma(timet^b,shape=s,rate=r))
haz
}


#computes the quantile function for the Generalized Gamma distribution
qggamma <- function(beta, sigma, lambda, p)
{
	s <- lambda^(-2)
	r <- s * exp(( - beta * lambda)/sigma)
	# consider p below 0.001 at 0.001 to avoid numerical problems
	p[p<0.001] <- 0.001
	# consider p above 0.999 at 0.999 to avoid numerical problems
	p[p>0.999] <- 0.999
	b <- lambda/sigma
	quan <- (lambda>0)*((qgamma(p,shape=s, rate=r))^(1/b))
	quan <- quan+(lambda<0)*((qgamma(1-p,shape=s, rate=r))^(1/b))
	quan
}

#Plot the estimated hazard function for best fitting two-parameter distribution in each period
#Period 1: Ammag; beta=0.637, sigma = 0.798 estimates from final model shown in last column of Table 3
ttx_seq(.01,8,.001)
ttyag_ hggamma(0.637,0.798,1/0.798,ttx)
tt95_qggamma(0.674,0.773,1.369,.95)

par(mgp=c(2, 1.5, 0))
plot(ttx[ttx<tt95],ttyag[ttx<tt95],xlab="Years from AIDS",ylab="",
	xlim=c(0,8),ylim=c(0,100),type='n', lty=1, lwd=1 ,axes=F,yaxs='s', xaxs='s',cex=1.2)
axis(2,at=seq(0,100,20),label=c("0","20","40","60","80","100"),cex=1.0,tck=-0.01,lwd=3.5)
par(mgp=c(2, 0.7, 0))
mixed.text.vector(-0.45*8/3, 0.8*50, c(1), c("Hazard"), srt=90)

axis(1,at=seq(0,8,2),label=c("0","2","4","6","8"), cex=1.0,lwd=3.5)
lines(ttx[ttx<tt95],100*ttyag[ttx<tt95],lwd=3)
text(5.2,94,"Period 1",cex=1.2)

#Period 2: Gamma; beta = 0.648, sigma = 0.848 estimates from final model shown in last column of Table 3
ttyga_ hggamma(0.648,0.848,0.848,ttx)
tt95_qggamma(0.649,0.847,0.851,.95)
lines(ttx[ttx<tt95],100*ttyga[ttx<tt95],lwd=3)
text(6.1,67,"Period 2",cex=1.2)


#Period 3: Ammag; beta = 2.432, sigma = 2.082 estimates from final model shown in last column of Table 3
ttyam_ hggamma(2.432,2.082,1/2.082,ttx)
tt95_min(c(8,qggamma(2.446,2.008,0.549,.95)))
lines(ttx[ttx<tt95],100*ttyam[ttx<tt95],lwd=3)
text(4.5,11,"Period 3",cex=1.2)


#Period 4: Weibull; beta = 3.106, sigma = 1.293 estimates from final model shown in last column of Table 3
ttywb_ hggamma(3.106,1.293,1,ttx)
tt95_min(c(8,qggamma(3.071,1.376,0.845,.95)))
lines(ttx[ttx<tt95],100*ttywb[ttx<tt95],lwd=3)
text(2,3.6,"Period 4",cex=1.2)


for (i in seq(0,100,20)) lines(c(0,8),c(i,i),lty=2,lwd=0.8)

dev.off()