## jh Nov 8, 2010

## See also the book
## R.A. Fisher, An Appreciation. by Stepehen E Fienberg and DV Hinkley.
## Preview via Googe Books  http://books.google.ca/books

## relevent portion starts on page 182 in the chapter Smoking and Lung Cancer.


data=scan()
 7  17  141  67  133  63  96  78  21  24
17  21  162  80  157  44  74  44  16   7 
 3   3    7   8    7   5   5   3   0   0
 2  10    2   7    6   0   0   0   1   0

#
A=array(data,dim=c(2,5,2,2))

a = A[1, ,1,1] ; b = A[2, ,1,1] ; 
c = A[1, ,2,1] ; d = A[2, ,2,1] ;
inhale = a+c ; inhale.not= b+d
cases  = a+b ; size.denom.series =c+d ; n = cases+size.denom.series
mh.num = a*d/n ; mh.den=b*c/n
a.null = inhale*cases/n
v.a.null = cases*size.denom.series*inhale*inhale.not/(n*n*n)

a.null.fisher = cases * (c / size.denom.series)

reduced.deficit.fisher =  - (a*d - b*c ) / n 

v.reduced.deficit.fisher = cases*size.denom.series*inhale*inhale.not/(n*n*n)

def = a - a.null.fisher
v.def = (1/size.denom.series^2) * cases*size.denom.series*inhale*inhale.not / n

# for Woolf

or = a*d/(b*c) ; 
log.or = log(or);
v.log.or = 1/a + 1/b + 1/c + 1/d
weight  = (1/v.log.or) / sum(1/v.log.or)

# for RGB variance of log(or.mh)

P = (a+d)/n; Q = (b+c)/n ; R = a*d/n ; S = b*c/n ;

cbind(a,b,c,d,inhale,inhale.not,cases, size.denom.series,
  n, mh.num, mh.den, a.null, v.a.null, 
  a.null.fisher, reduced.deficit.fisher, v.reduced.deficit.fisher, 
  def, v.def, or, log.or, v.log.or, weight, P, Q, R, S)

or.mh = sum(mh.num) / sum(mh.den) ; round(or.mh,2)

X.sq = ( sum(a - a.null) )^2 / sum(v.a.null) ;

c(X.sq, sqrt(X.sq) )

 20.299 / 8.274
 
sum(reduced.deficit.fisher)
sum(v.reduced.deficit.fisher)
 
X.Fisher = sum(reduced.deficit.fisher)  / sqrt( sum(v.reduced.deficit.fisher) ) ;

X.Fisher

sum(def)
sum(v.def)
sum(def)/sqrt(sum(v.def))

test.based.ci = or.mh^(1 + c(1.96,-1.96) / sqrt(X.sq) ) ; round(test.based.ci,2)

or.woolf = exp( sum(weight*log.or) ) ; round(or.woolf,2)

v.log.or.woolf = 1 / sum(1/v.log.or) ; v.log.or.woolf

woolf.ci = or.woolf * ( exp( c(-1.96,1.96) * sqrt(v.log.or.woolf ) ) ) 

round(woolf.ci,2)

# rgb

rgb.var.log.or.mh = sum(P*R)      / ( 2 *   ( (sum(R) )^2 ) ) 

                  + sum(P*S + Q*R)/ ( 2 *   sum(R)*sum(S)   )
              
                  + sum(Q*S)      / ( 2 *   ( (sum(S) )^2 ) ) 
              
rgb.var.log.or.mh

 
rgb.ci = or.mh * ( exp( c(-1.96,1.96) * sqrt(rgb.var.log.or.mh ) ) ) ;

round(rgb.ci,2)