pdf("ex0502.pdf") # Read in data from ASCII formatted (Fixed With File) ex0502dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0502.dat", width=c(8,8,8,8), col.names=c("pipesmk","lipcncr","pipesmk_r","lipcncr_r")) # Create a dataset (frame) from input data ex0502 <- data.frame(ex0502dat) # Attach the dataset for analysis attach(ex0502) # Obtain crosstabulation of pipesmk and lipcncr freqpl <- table(pipesmk, lipcncr) freqpl # Obtain the probabilities of lipcncr for each pipesmk group # Note that lipcncr is labelled "0" and "1", and we are modelling odds of "1" # Note that pipesmk is labelled "0" and "1", and we are modelling # odds(lipcncr=1|pipesmk=1)/odds(lipcncr=1|pipesmk=0) oddsrow1 <- freqpl[1,2]/freqpl[1,1] oddsrow2 <- freqpl[2,2]/freqpl[2,1] oddsratio <- oddsrow2/oddsrow1 selnor <- sqrt((1/freqpl[1,1])+(1/freqpl[1,2])+(1/freqpl[2,1])+(1/freqpl[2,2])) oddsratiolo <- exp(log(oddsratio)-(1.96*selnor)) oddsratiohi <- exp(log(oddsratio)+(1.96*selnor)) oddsrow1 oddsrow2 oddsratio oddsratiolo oddsratiohi dev.off()