a<-1 # set prior parameters a0<-a1<-b0<-b1<-a #---------------- #set up the data #---------------- y <- c( 25, 30) n <- c(300,900) ex2_1.risk <- data.frame(y=y,n=n) # # y1<- y[1] y0<- y[2] n1<- n[1] n0<- n[2] # generate 1000 observations from the posterior of p1 and p0 p1 <- rbeta(1000, y1+a1, n1+b1) p0 <- rbeta(1000, y0+a0, n0+b0) # obtain values for risk measures AR <- p1-p0 RR <- p1/p0 OR <- p1*(1-p0)/(p0*(1-p1)) # mean(AR);mean(RR);mean(OR)# calculate the posterior means # calculate the posterior standard deviations sd(AR);sd(RR);sd(OR) # calculate a 95% credible interval postscript('chap1_plot_sens1.ps', width = 10.0, height = 6.0, horizontal=FALSE) par( xaxs='i', yaxs='r', bty='l' , cex=1.5, mfrow=c(3,2)) # plots of posterior densities hist(AR, main='', xlab='Attributable Risk', ylab='', probability=TRUE) plot(density(AR), main='', xlab='Attributable Risk', ylab='Posterior Density') hist(RR, main='', xlab='Relative Risk', ylab='', probability=TRUE) plot(density(RR), main='', xlab='Relative Risk', ylab='Posterior Density') hist(OR, main='', xlab='Odds Ratio', ylab='', probability=TRUE) plot(density(OR), main='', xlab='Odds Ratio', ylab='Posterior Density') graphics.off()