pdf("ex0501.pdf")

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

ex0501dat <- read.fwf("http://www.stat.ufl.edu/~winner/data/biostat/ex0501.dat", width=c(8,8,8,8), 
        col.names=c("pamid","skev","pamid_r","skev_r"))


# Create a dataset (frame) from input data

ex0501 <- data.frame(ex0501dat)

# Attach the dataset for analysis 

attach(ex0501)

# Obtain crosstabulation of pamid and skev

freqps <- table(pamid, skev)

freqps

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

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

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

relrisk <- probrow2/probrow1

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

relrisklo <- exp(log(relrisk)-(1.96*selnrr))

relriskhi <- exp(log(relrisk)+(1.96*selnrr))

probrow1

probrow2

relrisk 

relrisklo 

relriskhi




dev.off()