rm(list=ls())
source("http://www.stat.washington.edu/tsr/s566/labs/dependence-functions.R")     
table.freq <- c(11038,38988,8962,41012) # Example from lecture notes

dim(table.freq) <- c(2,2)
dimnames(table.freq) <- list(c("Exposed","Unexposed"),c("Disease","No Disease"))

table.freq

risk.unexposed <- table.freq["Unexposed","Disease"]/(
                 table.freq["Unexposed","Disease"] + table.freq["Unexposed","No Disease"]
							)
risk.unexposed

risk.exposed   <- table.freq["Exposed","Disease"]/(
                    table.freq["Exposed","Disease"] + table.freq["Exposed","No Disease"])

risk.exposed 

risk.unexposed <- table.freq["Unexposed","Disease"]/sum(table.freq["Unexposed",])
risk.unexposed 

### (1) RISK DIFFERENCE

par(mfrow=c(1,1))

## Plot the risk in exposed and unexposed
do.prob.plot("Risk Difference")
points(risk.unexposed,risk.exposed,pch="+")

risk.difference <- risk.exposed - risk.unexposed
do.rd.line(risk.difference,color="green")


## (1a) Redo with different values in the table, or with a different RD


## (2) Relative Risk

relative.risk <- risk.exposed / risk.unexposed


do.prob.plot("Relative Risk")

points(risk.unexposed,risk.exposed,pch="+")
do.rr.line(relative.risk,color="red")

## (2a) Redo with different values in the table, or with a different RR



###(3) ODDS RATIO

# find the odds ratio
odds.ratio <- (risk.exposed / (1-risk.exposed))/
                         (risk.unexposed / (1-risk.unexposed))
odds.ratio <- getor(risk.exposed,risk.unexposed)

# for plot
do.prob.plot("Odds Ratio")
points(risk.unexposed,risk.exposed,pch="+")
do.or.curve(odds.ratio,col="blue")


###(3a) Re do with other values

##(4) Put all on one plot

do.prob.plot("RD,RR and OR")

points(risk.unexposed,risk.exposed,pch="+")
do.rd.line(risk.difference,color="green",lwd=1)
do.rr.line(relative.risk,color="red",lwd=1)
do.or.curve(odds.ratio,col="blue",lwd=1)


############

### (5) Stratified analysis

table.freq.strat <- c(9948,38388,52,1612,1090,600,8910,39400) # Example from lecture notes

dim(table.freq.strat) <- c(2,2,2)
dimnames(table.freq.strat) <- list(c("Exposed","Unexposed"),
                             c("Disease","No Disease"),
                             c("Female","Male"))
table.freq.strat

risk.unexposed.m <- table.freq.strat["Unexposed","Disease","Male"]/sum(table.freq.strat["Unexposed",,"Male"])

risk.unexposed.f <- table.freq.strat["Unexposed","Disease","Female"]/sum(table.freq.strat["Unexposed",,"Female"])

risk.unexposed.m
risk.unexposed.f

risk.exposed.m <- table.freq.strat["Exposed","Disease","Male"]/sum(table.freq.strat["Exposed",,"Male"])
risk.exposed.f <- table.freq.strat["Exposed","Disease","Female"]/sum(table.freq.strat["Exposed",,"Female"])

risk.exposed.m
risk.exposed.f

## Risk Differences

rd.f <- risk.exposed.f - risk.unexposed.f
rd.m <- risk.exposed.m - risk.unexposed.m

rd.f
rd.m


do.prob.plot("RD for Females vs Males")
points(risk.unexposed.f,risk.exposed.f,pch="F",)
points(risk.unexposed.m,risk.exposed.m,pch="M")
do.rd.line(rd.f,color="green")
do.rd.line(rd.m,lwd=1,color="red")

## Add back in combined
points(risk.unexposed,risk.exposed,pch="C")

points(c(risk.unexposed.f,risk.unexposed.m),c(risk.exposed.f,risk.exposed.m),lty=2,col="blue",type="l")

do.rd.line(risk.difference,lwd=1,color="black")



### Relative Risk

rr.f <- risk.exposed.f / risk.unexposed.f
rr.m <- risk.exposed.m / risk.unexposed.m

rr.f
rr.m

do.prob.plot("RR for Females vs Males")
points(risk.unexposed.f,risk.exposed.f,pch="F",)
points(risk.unexposed.m,risk.exposed.m,pch="M")
do.rr.line(rr.f,color="green")
do.rr.line(rr.m,lwd=1,color="red")


## Add back in combined
points(risk.unexposed,risk.exposed,pch="C")
points(c(risk.unexposed.f,risk.unexposed.m),c(risk.exposed.f,risk.exposed.m),lty=2,col="blue",type="l")

do.rr.line(relative.risk,lwd=1,color="black")


### Odds Ratio

or.f <- getor(risk.exposed.f,risk.unexposed.f)

or.m <- getor(risk.exposed.m,risk.unexposed.m)

or.f
or.m


do.prob.plot("OR for Females vs Males")
points(risk.unexposed.f,risk.exposed.f,pch="F",)
points(risk.unexposed.m,risk.exposed.m,pch="M")
do.or.curve(or.f,color="green")
do.or.curve(or.m,lwd=1,color="red")

points(risk.unexposed,risk.exposed,pch="C")
points(c(risk.unexposed.f,risk.unexposed.m),c(risk.exposed.f,risk.exposed.m),lty=2,col="blue",type="l")

do.or.curve(odds.ratio,col="blue",lwd=1)