pdf("ex0503.pdf")

# Read in data from ASCII formatted (Fixed With File)

ex0503dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0503.dat", width=c(8,8,8), 
        col.names=c("drgcls","cancer","cancer_r"))


# Create a dataset (frame) from input data

ex0503 <- data.frame(ex0503dat)

# Attach the dataset for analysis 

attach(ex0503)

# Obtain crosstabulation of drgcls and cancer

freqps <- table(drgcls, cancer)

freqps

# Obtain the probabilities of cancer for each drgcls group
# Note that cancer is labelled "0" and "1", and we are modelling probability of "1"
# Note that drgcls is labelled "1", "2", "3", and we are modelling P(cancer=1|drgcls=2)/p(cancer=1|drgcls=1)
# and P(cancer=1|drgcls=3)/p(cancer=1|drgcls=1)

probrow1 <- freqps[1,2]/(freqps[1,1]+freqps[1,2])

probrow2 <- freqps[2,2]/(freqps[2,1]+freqps[2,2])

probrow3 <- freqps[3,2]/(freqps[3,1]+freqps[3,2])

relrisk2 <- probrow2/probrow1

selnrr2 <- sqrt(((1-probrow1)/freqps[1,2])+((1-probrow2)/freqps[2,2]))

relrisklo2 <- exp(log(relrisk2)-(1.96*selnrr2))

relriskhi2 <- exp(log(relrisk2)+(1.96*selnrr2))

relrisk3 <- probrow3/probrow1

selnrr3 <- sqrt(((1-probrow1)/freqps[1,2])+((1-probrow3)/freqps[3,2]))

relrisklo3 <- exp(log(relrisk3)-(1.96*selnrr3))

relriskhi3 <- exp(log(relrisk3)+(1.96*selnrr3))

probrow1

probrow2

probrow3

relrisk2 

relrisklo2 

relriskhi2

relrisk3 

relrisklo3 

relriskhi3


dev.off()