### 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 <- 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[i] <- mean(x.aspirin.sample)/mean(x.placebo.sample)
}
hist(theta) # histogram of the bootstrap
samples
plot(density(theta),type="l") #
density plot
quantile(theta, c(0.025,0.975)) #
empirical bootstrap CI.