library(postscriptfonts) ps.options(setfont=ps.setfont.latin1) options(mixed.text.screen=T) ps.options( reset=T, colors = ps.colors.rgb[c("black", "green", "blue"),] ) 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 } 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[timet0)*((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 } postscript("C:/mikefolder/mikeFolder/survivalmethods/programsforweb/Figures_revisedpaper/GGfig2.ps", horizontal=T, height=6.5, width=8, color.p=T) xx_seq(.001, 5.0, .01) par(mfrow=c(2,2), cex=1.0, tck=-0.02, mgp=c(2, 1.0, 0), mar=c(4,4,1,2)+.1, lwd=2) ibeta=0 plot(xx, hggamma(ibeta, 4/3, 1.0, xx), lwd=2, type="l", xlim=c(0, 5.0), ylim=c(0, 1.19), ylab="Hazard", xlab="Time", cex=1.2, xaxs='s', xaxt='n', yaxs='s', yaxt="n") axis(1, at=seq(0, 5.0, 1.0), labels=F, line=0) axis(1, at=seq(0, 5.0, 1.0), labels=T, ticks=F, line=-0.5) axis(2, at=seq(0, 1.2, 0.2), labels=T, line=0) #mixed.text(2.5, 1.6, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 4/3, ~f13~.l~f1~.= 1", slope=0, cex=0.8) #mixed.text(2.5, 1.45, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 8/3, ~f13~.l~f1~.= 1", slope=0, cex=0.8) #key(x=1.7, y=1.67, lines=list(lty=1, type="l", size=2), adj=0, cex=0.8) #key(x=1.7, y=1.52, lines=list(lty=2, type="l", size=2), adj=0, cex=0.8) lines(xx, hggamma(ibeta, 1.15, 0.9, xx), lty=2) lines(xx, hggamma(ibeta, 2, 0.6, xx), lty=2) lines(xx, hggamma(ibeta, 1.1, 1.05, xx), lty=2) plot(xx, hggamma(ibeta, 3/4, 1.0, xx), lwd=2, type="l", xlim=c(0, 5.0), ylim=c(0, 5), ylab="Hazard",xlab="Time", cex=1.2, xaxs='s', xaxt='n', yaxs='s', yaxt="n") axis(1, at=seq(0, 5.0, 0.5), labels=F, line=0) axis(1, at=seq(0, 5.0, 1.0), labels=T, ticks=F, line=-0.5) axis(2, at=seq(0, 5.0, 1.0), labels=T, line=0) #mixed.text(1.8, 4.5, "~f13~.b~f1~.= 0, ~f13~.s~f1~.=3/4, ~f13~.l~f1~.= 1", slope=0, cex=0.8) #mixed.text(1.8, 4.18, "~f13~.b~f1~.= 0, ~f13~.s~f1~.=0.6, ~f13~.l~f1~.= 1", slope=0, cex=0.8) #key(x=1.0, y=4.65, lines=list(lty=1, type="l", size=2), adj=0, cex=0.8) #key(x=1.0, y=4.33, lines=list(lty=2, type="l", size=2), adj=0, cex=0.8) lines(xx, hggamma(ibeta, 0.4, 0.41, xx), lty=2) lines(xx, hggamma(ibeta, 0.4, 2.0, xx), lty=2) lines(xx, hggamma(ibeta, 0.9, 0.95, xx), lty=2) plot(xx, hggamma(ibeta, 1.00, 4, xx), lwd=2, type="l", xlim=c(0, 5), ylim=c(0, 5), ylab="Hazard",xlab="Time", cex=1.2, xaxs='s', xaxt='n', yaxs='s', yaxt="n") axis(1, at=seq(0, 5.0, 1.0), labels=F, line=0) axis(1, at=seq(0, 5.0, 1.0), labels=T, ticks=F, line=-0.5) axis(2, at=seq(0, 5.0, 1.0), labels=T, line=0) #mixed.text(1.8, 7, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 1, ~f13~.l~f1~.= 3.0", slope=0, cex=0.8) #mixed.text(1.8, 6.4, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 1, ~f13~.l~f1~.= 2.0", slope=0, cex=0.8) #key(x=1.0, y=7.3, lines=list(lty=1, type="l", size=2), adj=0, cex=0.8) #key(x=1.0, y=6.7, lines=list(lty=2, type="l", size=2), adj=0, cex=0.8) lines(xx, hggamma(ibeta, 1.25, 2.5, xx), lty=2) lines(xx, hggamma(ibeta, 1.0, 2.50, xx), lty=2) lines(xx, hggamma(ibeta, 2, 4, xx), lty=2) plot(xx, dlnorm(xx, 0, 0.75)/(1-plnorm(xx, 0, 0.75)), lwd=2, type="l", xlim=c(0, 5.0), ylim=c(0, 1.19), ylab="Hazard",xlab="Time", cex=1.2, xaxs='s', xaxt='n', yaxs='s', yaxt="n") axis(1, at=seq(0, 5.0, 1.0), labels=F, line=0) axis(1, at=seq(0, 5.0, 1.0), labels=T, ticks=F, line=-0.5) axis(2, at=seq(0, 1.2, 0.2), labels=T, line=0) #mixed.text(2.3, 0.9, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 1.0, ~f13~.l~f1~.= 1/4", slope=0, cex=0.8) #mixed.text(2.3, 0.83, "~f13~.b~f1~.= 0, ~f13~.s~f1~.= 1.4, ~f13~.l~f1~.= 1/4", slope=0, cex=0.8) #key(x=1.5, y=0.94, lines=list(lty=1, type="l", size=2), adj=0, cex=0.8) #key(x=1.5, y=0.87, lines=list(lty=2, type="l", size=2), adj=0, cex=0.8) lines(xx, hggamma(ibeta, 1.40, 0.25, xx), lty=2) lines(xx, hggamma(ibeta, 1.00, 0.50, xx), lty=2) lines(xx, hggamma(ibeta, 1, -3, xx), lty=2) dev.off()