# MLE-conditional of common OddsRatio from k 2 x 2 tables

# using as example, the data in Table 5.2 B&D Vol II (BrCa estrogen)

ds=read.table("BD-II-5-2-BrCa-estrogen.txt",header=T)
ds
attach(ds)

cases=c.H1+c.H0
nc.H0=nc.H0+nc.H2
non.cases=nc.H1+nc.H0

# M-H summary OR

  sum(c.H1*nc.H0/n) / sum(c.H0*nc.H1/n)
  
## Newton-Raphson iterations to OR-ml-condn'l

# for explanation of the notation, see the article 
# An 'Unconditional-like' Structure for the Conditional
# Estimator of Odds Ratio from 2 x 2 Tables
# Biometrical Journal 48 (2006) 1, 23—34 DOI: 10.1002/bimj.200510167,
# available on http://www.epi.mcgill.ca/hanley/reprints

# function to get (non-null) expectation and variance
# of a r.v. with a non-central hypergeometric distr'n

# uses the (diagonal) cell beside the row or column
# with the smallest marginal frequency (m  for 'min'). 
# This way, the range of this cell frequency is from 0 
# to m and so the highest power of the normalizing polynomial
# is also m.

# examples

# 1 2 | 3    the a cell is opposite the column with
# 1 3 | 4    the smallest margin(2) .. 
# -------      min.margin.opposite.f   = 2 [col 1] (f is 1)
# 2 5 | 7      other.margin.opposite.f = 3 [row 1] ...  n = 7

# 4 4 | 8    the d cell is opposite the row with
# 1 2 | 3    the smallest margin(2) .. 
# -------       min.margin.opposite.f   = 3 [row 2] (f is 2) 
# 5 6 | 11      other.margin.opposite.f = 6 [col 2] ... n = 11  


condnl.mean.var=
 function(min.margin.opposite.f,other.margin.opposite.f,n,OR){
 	min.margins.neighbour = n - min.margin.opposite.f;
 	other.margins.neighbour = n - other.margin.opposite.f;
 	b=c();
 	for (i in 0:min.margin.opposite.f){
 		rest = min.margin.opposite.f-i ;
 		b=c(b,choose(other.margin.opposite.f,i)*
 		      choose(other.margins.neighbour,rest))
 		}
 	values=seq(0,min.margin.opposite.f);
 	powers=OR^values;
 	mean=sum(values*b*powers)/sum(b*powers);
 	var=sum(values*values*b*powers)/sum(b*powers) - mean^2;
 	return(c(mean,var))
 	}
 	
# test condnl.mean.var(2,3,7,1); condnl.mean.var(2,3,7,sqrt(2))

# function to iterate via Newton-Raphson...

# arguments 4 vectors, each of size k, for the traditional a b c d

ml.condnl=function(a,b,c,d){
 	k=length(a);
 	target=c(); n=c(); # target ca be different in each table
 	# for each table, find the 2 key margins...
 	min.opposite.target=c();other.opposite.target=c();
 	for (i in 1:k){
 		r1=a[i]+b[i]; r2=c[i]+d[i]; c1=a[i]+c[i];c2=b[i]+d[i];
 		n=c(n,r1+r2);
 		           which="row1"; m=r1; 
 		if (r2 < m) {which="row2"; m=r2};
 		if (c1 < m) {which="col1"; m=c1};
 		if (c2 < m) {which="col2"; m=c2};
 		if (which=="row1") {
 			 target=c(target,a);
 			 min.opposite.target=c(min.opposite.target,r1);
 			 other.opposite.target=c(other.opposite.target,c1) 
 			                 };
 		if (which=="row2") {
 			 target=c(target,d);
 			 min.opposite.target=c(min.opposite.target,r2);
 			 other.opposite.target=c(other.opposite.target,c2) 
 			                 };
 		if (which=="col1") {
 			 target=c(target,a);
 			 min.opposite.target=c(min.opposite.target,c1);
 			 other.opposite.target=c(other.opposite.target,r1) 
 			                 };
 		if (which=="col2") {
 			 target=c(target,d);
 			 min.opposite.target=c(min.opposite.target,c2);
 			 other.opposite.target=c(other.opposite.target,r2) 
 			                 }
   };
 	sum.target=sum(target); # estimating eqn: sum(targets) = sum(fitted targets)
 	
 	logOR.prev=0; OR.prev=exp(logOR.prev) ;
 	
 	for (iteration in 1:10){
 		U = 0; I=0;
 		for (i in 1:k){
 			m.v=condnl.mean.var(min.opposite.target[i],
 			                 other.opposite.target[i],n[i],OR.prev);
 			U=U+(target[i]-m.v[1]) ; 
 			I = I + m.v[2]
 			};
 		V=1/I;
 		print(round(c(logOR.prev,OR.prev,U,I,V),6));
 		logOR.prev = logOR.prev + U/I; 
 		OR.prev=exp(logOR.prev)
 	}
}	
 
# 1 stratum

ml.condnl(1,2,1,3)

# the strata in the estrogen study above

a=c.H1; b=cases-a; c= nc.H1; d=non.cases-c; 

ml.condnl(a,b,c,d)