#Will generate the graph of the log of the cumulative hazard of death for each of the four periods 
#Path below reads in data file that contains for each unique event time for each of the four period three variables: 
#exit time, probability of survival, Period 1-4 indicator for the non-parametric estimated cumulative mortality
nonpara_matrix(scan("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/nonparametricFig3.dat"),ncol=3,byrow=T)

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


#time to death Period 1
timeaid1_nonpara[nonpara[,3]==1,1]
#time to death Period 2
timeaid2_nonpara[nonpara[,3]==2,1]
#time to death Period 3
timeaid3_nonpara[nonpara[,3]==3,1]
#time to death Period 4
timeaid4_nonpara[nonpara[,3]==4,1]

#Period 1 mortality
survaid1_(1-nonpara[nonpara[,3]==1,2])*100
#Period 2 mortality
survaid2_(1-nonpara[nonpara[,3]==2,2])*100
#Period 3 mortality
survaid3_(1-nonpara[nonpara[,3]==3,2])*100
#Period 4 mortality
survaid4_(1-nonpara[nonpara[,3]==4,2])*100

surv1 _ log(-log(1-(survaid1/100)))
surv2 _ log(-log(1-(survaid2/100)))
surv3 _ log(-log(1-(survaid3/100)))

#maximum likelihood parameter estimates from saturated Generalized Gamma models (Table 3 of paper)
#Period 1: beta = 0.674, sigma = 0.773, lamda = 1.369
#Period 2: beta = 0.649, sigma = 0.847, lamda = 0.851
#Period 3: beta = 2.446, sigma = 2.008, lamda = 0.549
#Period 4: beta = 3.071, sigma = 1.376, lamda = 0.845

tta1_stepfun(timeaid1,surv1) 
  par(mgp=c(3, 1.5, 0))
  plot(tta1$x,tta1$y,xlab='Years from AIDS diagnosis', ylab='Log of Cumulative Hazard of Death after AIDS',xlim=c(0,8),ylim=c(-3,1),type='l',lwd=3, yaxs='s', xaxs='s', axes=F)
  lines(tta1$x,tta1$y,lty=1,lwd=3) 
  sgam1_ pgamma((exp(-0.674)*tta1$x)^(1.369/.773)/(1.369^2),1/(1.369^2),1,1)
  sgam1_ log(-log(1-sgam1))
  lines(tta1$x,sgam1,lty=4,lwd=3)
  axis(2, at=seq(-3, 1, 1), labels=T, ticks=T, line=0, outer=F, lwd=3)
  key(x=5, y=95, lines=list(lty=c(1, 4), cex=1.5), text=list(c("Nonparametric", "Generalized Gamma")), lwd=3)
  par(mgp=c(2, 1, 0))
  axis(1, at=seq(0, 8, 1), labels=T, ticks=T, line=0, outer=F, lwd=3)

tta2_stepfun(timeaid2,surv2)
 lines(tta2$x,tta2$y,lty=1,lwd=3) 
 sgam2_ pgamma((exp(-0.649)*tta2$x)^(0.851/0.847)/(0.851^2),1/(0.851^2),1,1)
 sgam2_ log(-log(1-sgam2))
 lines(tta2$x,sgam2,lty=4,lwd=3)

tta3_stepfun(timeaid3,surv3)
 lines(tta3$x,tta3$y,lty=1,lwd=3)
 sgam3 _ pgamma((exp(-2.446)*tta3$x)^(0.549/2.008)/(0.549^2),1/(0.549^2),1,1)
 sgam3_ log(-log(1-sgam3))
 lines(tta3$x,sgam3,lty=4,lwd=3)

tta4_stepfun(timeaid4,survaid4)
 correction4 <-(1-(pgamma((exp(-3.071)*0.5)^(0.845/1.376)/(0.845^2),1/(0.845^2),1,1)))
 newsurv _ 100-((100 - tta4$y) * correction4)
 newsurv _ log(-log(1-(newsurv/100)))
 print(tta4$y)
 lines(tta4$x[tta4$x>0],newsurv[tta4$x>0],lty=1,lwd=3)
 sgam4_ pgamma((exp(-3.071)*tta4$x)^(0.845/1.376)/(0.845^2),1/(0.845^2),1,1)
 sgam4_ log(-log(1-sgam4))
 lines(tta4$x,sgam4,lty=4,lwd=3)

text(1.7,0.55,"July 1984 - Dec 1989")
text(3.7,0.35,"Jan 1990 - Dec 1994")
text(3.7,-0.65,"Jan 1995 - June 1998")
text(6.7,-1.15,"July 1998 - Dec 2003")

dev.off()