#	Table 3.1 data

Seat.belt<-c("no","yes") 
Fatality<-c("yes","no") 
tbl <-expand.grid(seat.belt=Seat.belt,fatality=Fatality) 
Data<-c(1601,510,162527,412368)
tbl.2<-cbind(tbl,count=Data) 
( tmp <- tapply(tbl.2$count,tbl.2[,1:2], sum)  ) 

tmp[1,1] / sum(tmp[1,])
tmp[2,1] / sum(tmp[2,])

Wald.ci<-function(Table, aff.response, alpha=.05){ 
# Gives two-sided Wald CI's for odds ratio, difference in proportions and relative risk. 
# Table is a 2x2 table of counts with rows giving the treatment populations 
# aff.response is a string like "c(1,1)" giving the cell of the beneficial response and the 
# treatment category 
# alpha is significance level 
pow<-function(x, a=-1) x^a 
z.alpha<-qnorm(1-alpha/2) 
if(is.character(aff.response)) 
where<-eval(parse(text=aff.response)) 
else where<-aff.response 
Next<-as.numeric(where==1) + 1 
# OR 
odds.ratio<-Table[where[1],where[2]]*Table[Next[1],Next[2]]/(Table[where[1],Next[2]]*Table[Next[1],where[2]]) 
se.OR<-sqrt(sum(pow(Table))) 
ci.OR<-exp(log(odds.ratio) + c(-1,1)*z.alpha*se.OR) 
# difference of proportions 
p1<-Table[where[1],where[2]]/(n1<-Table[where[1],Next[2]] + Table[where[1],where[2]]) 
p2<-Table[Next[1],where[2]]/(n2<-Table[Next[1],where[2]]+Table[Next[1],Next[2]]) 
se.diff<-sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2) 
ci.diff<-(p1-p2) + c(-1,1)*z.alpha*se.diff 
# relative risk 
RR<-p1/p2 
se.RR<-sqrt((1-p1)/(p1*n1) + (1-p2)/(p2*n2)) 
ci.RR<-exp(log(RR) + c(-1,1)*z.alpha*se.RR) 
list(OR=list(odds.ratio=odds.ratio, CI=ci.OR), proportion.difference=list(diff=p1-p2, CI=ci.diff), relative.risk=list(relative.risk=RR,CI=ci.RR)) 
} 

Wald.ci(tmp,aff.response=c(1,1))




# Haashara: log-odds ratio variance in three sampling schemes


#	1. poisson sampling


mu.11 <- 60 ; mu.12 <- 40 ; mu.21 <- 150 ; mu.22 <- 50
or.pois <- rep(NA,10000)

for(i in 1:10000)
{

n.11 <- rpois(1,mu.11)
n.12 <- rpois(1,mu.12)
n.21 <- rpois(1,mu.21)
n.22 <- rpois(1,mu.22)

or.pois[i] <- n.11 * n.22 / (n.21 * n.12)
}

#	2. multinomial sampling

pi.11 <- 60  / 300 ; pi.12 <- 40  / 300 
pi.21 <- 150 / 300 ; pi.22 <- 50  / 300
or.mult <- rep(NA,10000)

for(i in 1:10000)
{

n.vec <- rmultinom(1,300,prob=c(pi.11,pi.12,pi.21,pi.22))
n.11 <- n.vec[1,1]
n.12 <- n.vec[2,1]
n.21 <- n.vec[3,1]
n.22 <- n.vec[4,1]

or.mult[i] <- n.11 * n.22 / (n.21 * n.12)
}

#	3. product multinomial sampling

pi.11 <- 60  / 100  ;  pi.12 <- 40  / 100
pi.21 <- 150 / 200  ;  pi.22 <- 50  / 200

or.pmlt <- rep(NA,10000)
for(i in 1:10000)
{


n.vec <- rmultinom(1,100,prob=c(pi.11,pi.12))
n.11 <- n.vec[1,1]
n.12 <- n.vec[2,1]

n.vec <- rmultinom(1,200,prob=c(pi.21,pi.22))
n.21 <- n.vec[1,1]
n.22 <- n.vec[2,1]

or.pmlt[i] <- n.11 * n.22 / (n.21 * n.12)
}

boxplot(list(or.pois,or.mult,or.pmlt))
abline(h = (60 * 50) / (150 * 40),col=2)

(60 * 50) / (150 * 40)

t.test(or.pois - (60 * 50) / (150 * 40))
t.test(or.mult - (60 * 50) / (150 * 40))
t.test(or.pmlt - (60 * 50) / (150 * 40))



boxplot(list(log(or.pois),log(or.mult),log(or.pmlt)))
abline(h = log((60 * 50) / (150 * 40)),col=2)

t.test(log(or.pois) - log((60 * 50) / (150 * 40)))
t.test(log(or.mult) - log((60 * 50) / (150 * 40)))
t.test(log(or.pmlt)  - log((60 * 50) / (150 * 40)))

var(log(or.pois))
var(log(or.mult))
var(log(or.pmlt))

1/60 + 1/50 + 1/150 + 1/40