pdf("ex0507.pdf") # Read in data from ASCII formatted (Fixed With File) ex0507dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0507.dat", width=c(8,8,8,8,8,8), col.names=c("age","smoke","death","numcases","smoke_r","death_r")) # Create a dataset (frame) from input data ex0507 <- data.frame(ex0507dat) # Attach the dataset for analysis attach(ex0507) numages <- max(age) # Obtain crosstabulation of smoke and death separately by age with numcases being the count in each cell freqasd <- xtabs(numcases ~ smoke+death+age) freqasd # Obtain R=AD/n and S=BC/n for each age group and sum them r_sum <- 0.0 s_sum <- 0.0 r_mh <- c(rep(0,numages)) s_mh <- c(rep(0,numages)) for (i in 1:numages) { r_mh[i] <- freqasd[1,1,i]*freqasd[2,2,i]/(freqasd[1,1,i]+freqasd[1,2,i]+freqasd[2,1,i]+freqasd[2,2,i]) s_mh[i] <- freqasd[1,2,i]*freqasd[2,1,i]/(freqasd[1,1,i]+freqasd[1,2,i]+freqasd[2,1,i]+freqasd[2,2,i]) r_sum <- r_sum+r_mh[i] s_sum <- s_sum+s_mh[i] } # Obtain the Mantel-Haenszel estimate of odds ratio and print results or_mh <- r_sum/s_sum r_mh s_mh r_sum s_sum or_mh dev.off()