### Bootstrap: Relative Risk

###

### Consider a prospective cohort study to investigate the effect of apsirin on heart attacks.

### A group of patients who are at risk of heart attack are randomly assigned to take either

### a placebo or an aspirin. At the end of one year the number of patients suffering a heart

### attack is recorded. The results were

### Aspirin: Heart Attack 104 No 10933 No. Subjects 11037

### Placebo: Heart Attack 189 No 10845 No. Subjects 11034

### Get a 95% CI for the Relative Risk = p1/p2 of aspirin to placebo using the traditional

### Mantel-Haenszel test. Using SAS run Proc Freq.

### Use the Bootstrap to plot the Bootstrap distribution of the RR and give a 95% Bootstrap

### Empirical Confidence Interval for RR.

# create the data

asp.ht <- 104

n1 <- 11037

x.aspirin <- numeric(n1)

x.aspirin[1:asp.ht] <- 1

sum(x.aspirin)

mean(x.aspirin)

pla.ht <- 189

n2 <- 11034

x.placebo <- numeric(n2)

x.placebo[1:pla.ht] <- 1

sum(x.placebo)

mean(x.placebo)

mean(x.aspirin)/mean(x.placebo) # RR

### Implement the bootstrap directly.

B <- 1000

theta.star <- numeric(B)

for(i in 1:B){

x.aspirin.sample <- sample(x.aspirin, size = n1, replace = T)

x.placebo.sample <- sample(x.placebo, size = n2, replace = T)

theta.star[i] <- mean(x.aspirin.sample)/mean(x.placebo.sample)

}

mean(theta.star)

hist(theta.star) # histogram of the bootstrap samples

plot(density(theta.star),type="l") # density plot

quantile(theta.star, c(0.025,0.975)) # empirical bootstrap CI.

### Implement the bootstrap using bootstrap()

HeartAttack <- c(x.aspirin,x.placebo) # create column of data

Treatment <- c(rep(1,n1),rep(0,n2)) # create column indicating treatment

X <- cbind(HeartAttack,Treatment) # combine data into a matrix with 2 columns

mean(X[X[,2]==1,1])/mean(X[X[,2]==0,1]) # compute the RR for the data

RR <- function(X){

# This function computes the RR for data created with the data in the

# first column and the indicator in the second column.

mean(X[X[,2]==1,1])/mean(X[X[,2]==0,1])

}

RR(X) # check that the RR function works

boot.RR <- bootstrap(X, RR) # use the bootstrap() with the RR function

# options(object.size = 1000000000)

summary(boot.RR)

limits.emp(boot.RR)

limits.bca(boot.RR)

plot(boot.RR)