Biostat CDA - Codes for Chapter 03

## smoking data table
t33
## CMH test
n = rowSums(t33)
u11 = (t33[,1]+t33[,2])*(t33[,1]+t33[,3])/n
v11 = (t33[,1]+t33[,2])/n*(t33[,1]+t33[,3])/n*(t33[,3]+t33[,4])/(n-1)*(t33[,2]+t33[,4])
cmh = (sum(u11)-sum(t33[,1]))^2/sum(v11)
cmh
1-pchisq(cmh,1)
## MH common odds ratio estimation
mh = sum(t33[,1]*t33[,4]/n)/sum(t33[,2]*t33[,3]/n)
mh
## Breslow-Day test
bd = apply(t33, 1, function(x){
  n1p = x[1]+x[2]; np1 = x[1]+x[3]
  np2 = x[2]+x[4]
  b2 = mh-1; b1 = -(n1p+np1)*mh-(np2-n1p); b0 = n1p*np1*mh
  u11 = (-b1+sqrt(b1^2-4*b2*b0))/2/b2
  v11 = (-b1-sqrt(b1^2-4*b2*b0))/2/b2
  if( all(np1-v11>0,n1p-v11>0,np2-n1p+v11>0) ){
    u11 = v11
  }  
  return((u11-x[1])^2*(np1/u11/(np1-u11)+np2/(n1p-u11)/(np2-n1p+u11)))
})
sum(bd)
1-pchisq(sum(bd),7)

## Exact inference: artificial example
t1 = c(1,8,5,17)
t2 = c(1,8,5,14)
t3 = c(1,9,3,14)
n123 = expand.grid(0:3,0:3,0:3)
a = n123[rowSums(n123)<=3,]
## Prob
p1 = dhyper(a[,1], t1[1]+t1[2],t1[3]+t1[4],t1[1]+t1[3]) 
p2 = dhyper(a[,2], t2[1]+t2[2],t2[3]+t2[4],t2[1]+t2[3]) 
p3 = dhyper(a[,3], t3[1]+t3[2],t3[3]+t3[4],t3[1]+t3[3]) 
## p-value
sum(p1*p2*p3)