###################################################################################
# R code used to estimate the ACE with bootstrap 95% CIs, and to create a plot
#
# 2-10-2006
###################################################################################
#
# Read in the data
datamat1 <- matrix(scan('/home/pgilbert/teaching/vaccine2006/Data/Vax004data578.12.05.part1.dat'),ncol=26,byrow=T)

sidnumb1 <- datamat1[,1]
vacc <- datamat1[,2]
infected <- datamat1[,18]
startart <- datamat1[,22]
daysstartart <- datamat1[,23]

datamat2 <- matrix(scan('/home/pgilbert/teaching/vaccine2006/Data/Vax004data578.12.05.part2.dat'),ncol=4,byrow=T)

sidnumb2 <- datamat2[,1]
grptime <- datamat2[,2]
vl <- datamat2[,3]

# Define the set-point viral loads:
n <- length(vacc)
setptvl <- rep(NA,n)

for (i in 1:n) {
if (infected[i]==1 & !is.na(startart[i]) & (startart[i]==0 | (startart[i]==1 & daysstartart[i] > 30))) {
# condition on being infected and not starting ART by day 30 post infection diagnosis

tempvls <- vl[sidnumb2==sidnumb1[i] & grptime==30]
# get all the viral loads at day 30

if (length(tempvls[!is.na(tempvls)])>0) {
# set the set-point value to be the mean of the available values

setptvl[i] <- mean(tempvls[!is.na(tempvls)]) }}}

source('Viral.Load.Sensitivity.Analysis.Program.Meanonly.r')

Nn1 <- length(vacc[vacc==1])
Nn0 <- length(vacc[vacc==0])
y1.samp <- setptvl[!is.na(setptvl) & vacc==1]
y0.samp <- setptvl[!is.na(setptvl) & vacc==0]

# Estimate the ACE with 95% CIs over a range of beta values
# Use 500 bootstrap replications
b <- 500
betavect <- seq(-3,3,length=30)

for (i in 1:length(betavect)) {
beta <- betavect[i]
viralloadprogram(b,Nn1,Nn0,y1.samp,y0.samp,beta) }

# The answers are written to the file output.dat

# Now plot the estimated ACE (with 95% pointwise bootstrap CIs) versus beta

# Write a plotting function:

Plotsensitivity.ACE <- function(datafile,filename="ACEplot.ps") {

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

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

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

plot(beta,ACE,xlab="beta",ylab="Estimated ACE(beta)",xlim=c(-3.2,3.2),type='n',cex=1.1)

lines(beta,ACE,lwd=3,lty=1)
lines(beta,lowCI,lwd=2,lty=3)
lines(beta,upCI,lwd=2,lty=3)

abline(h=0)

title("Estimated ACE and 95% pointwise CIs as a function of beta")
dev.off()
}

# Call the plotting function

Plotsensitivity.ACE("output.dat")

# The graph is in the file ACEplot.ps