### 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)