pdf("ex0504.pdf") # Read in data from ASCII formatted (Fixed With File) ex0504dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0504.dat", width=c(8,8,8,8), col.names=c("trt","improve","trt_r","improve_r")) # Create a dataset (frame) from input data ex0504 <- data.frame(ex0504dat) # Attach the dataset for analysis attach(ex0504) # Obtain crosstabulation of trt and improve freqti <- table(trt, improve) freqti # Obtain the probabilities of improve for each trt group # Note that improve is labelled "0" and "1", and we are modelling probability of "1" # Note that trt is labelled "0" and "1", and we are modelling P(improve=1|trt=1)-p(improve=1|trt=0) probrow1 <- freqti[1,2]/(freqti[1,1]+freqti[1,2]) probrow2 <- freqti[2,2]/(freqti[2,1]+freqti[2,2]) absrisk <- probrow2-probrow1 selnar <- sqrt(((probrow1*(1-probrow1))/(freqti[1,1]+freqti[1,2]))+((probrow2*(1-probrow2))/(freqti[2,1]+freqti[2,2]))) absrisklo <- absrisk-(1.96*selnar) absriskhi <- absrisk+(1.96*selnar) probrow1 probrow2 absrisk absrisklo absriskhi dev.off()