# datafile is output.dat, the output of Viral.Load.Sensitivity.Analysis.Program.r
# filename1 is the name of the postscript file that shows the estimated AVE(beta) with confidence intervals
# filename2 is the name of the postscript file that shows the 1-sided p-values for the three test statistics
#           versus beta
# smoothed = option for lowess smoothed confidence limits and 1-sided p-values are not smoothed

Viral.Load.plots <- function(datafile="output.dat",filename1="figure1.ps",filename2="figure2.ps",smoothed=T) {

postscript(filename1,horizontal=T)
par(mar=c(4.1,4.1,4.1,4.1),las=1)

mat <- matrix(scan(datafile),ncol=11,byrow=T)

beta <- mat[,1]
n <- length(beta)
ACE <- mat[,9]
lowCI <- mat[,10]
upCI <- mat[,11]

plot(beta,ACE,xlab="Sensitivity parameter beta",ylab="Estimated ACE(beta)",type='n',ylim=c(min(lowCI),max(upCI)),cex=1.1)

if (smoothed) {
lines(lowess(beta,ACE),lwd=3,lty=1)
lines(lowess(beta,lowCI),lwd=2,lty=3)
lines(lowess(beta,upCI),lwd=2,lty=3) }
else {
lines(beta,ACE,lwd=3,lty=1)
lines(beta,lowCI,lwd=2,lty=3)
lines(beta,upCI,lwd=2,lty=3) }

abline(h=0,lty=2)

title("Estimated ACE(beta) with 95% Pointwise Confidence Intervals")

dev.off()

postscript(filename2,horizontal=T)
par(mar=c(4.1,4.1,4.1,4.1),las=1)

y1 <- mat[,6]
y2 <- mat[,7]
y3 <- mat[,8]

plot(beta,y1,xlab="Sensitivity parameter beta",ylab="1-sided p-value",type='n',cex=1.1,ylim=c(0,1))

if (smoothed) {
lines(lowess(beta,y1),lty=1)
lines(lowess(beta,y2),lty=2)
lines(lowess(beta,y3),lty=3) }
else { 
lines(beta,y1,lty=1)
lines(beta,y2,lty=2)
lines(beta,y3,lty=3) }

legend(x=0,y=.5,legend=c("Semiparametric Means","Kolmogorov-Smirnov","Anderson-Darling"),lty=c(1:3))

title("1-sided P-values for 3 Semiparametric Test Statistics") 

dev.off()

}