# Function "MH2x2.test"  
# plots extreme and rejection regions (one-sided, nominal alpha=0.05)
# for the Mantel-Haenszel test
# calculating p-value and the actual size
# for a pair of 2x2 contingency tables 
#  x1,x2	y1,y2
#  x3,x4	y3,y4
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>

MH2x2.test_function(x1,x2,x3,x4,y1,y2,y3,y4)
{
xx1_min(x1+x2,x1+x3)
yy1_min(y1+y2,y1+y3)
v1_rep(c(0:xx1),yy1+1)
v2_as.vector(t(matrix(rep(c(0:yy1),xx1+1),ncol=(xx1+1))))
nprob_dhyper(v1,x1+x2,x3+x4,x1+x3)*dhyper(v2,y1+y2,y3+y4,y1+y3)
Gamma1_cbind(v1,v2,nprob)
observed_as.numeric(Gamma1[,1] == x1 & Gamma1[,2] == y1)
Gamma1_cbind(Gamma1,observed)
x_v1
y_v2

par(oma=c(0, 0, 3, 0),mfcol=c(1,2))

minx_min(Gamma1[,1])
maxx_max(Gamma1[,1])
miny_min(Gamma1[,2])
maxy_max(Gamma1[,2])
Manhen.region_Gamma1[v1+v2 >= x1+y1,]
noManhen.region_Gamma1[v1+v2 < x1+y1,]
Manhen.p_sum(Manhen.region[,3])
plot(Gamma1[,1],Gamma1[,2],xlab="X11",ylab="Y11",type="n",lab=c(maxx-minx+1,maxy-miny+1,2))
title(main=paste("p=",round(Manhen.p,digits=4)))
text(Manhen.region[,1],Manhen.region[,2],Manhen.region[,1]+Manhen.region[,2])
if (is.array(noManhen.region)) points(noManhen.region[ ,1],noManhen.region[ ,2], type="p",pch=1) else
    if (length(noManhen.region)>0) points(noManhen.region[1],noManhen.region[2], type="p",pch=1) 
points(x1,y1, type="p", pch=1, cex=2)

size_0
porog_max(v1+v2)
while (size <= 0.05) 
{ size_size+sum(nprob[v1+v2 == porog])
  porog_porog-1
}
Manhen.rejreg_Gamma1[v1+v2 > porog+1,]
noManhen.rejreg_Gamma1[v1+v2 <= porog+1,]
if (is.array(Manhen.rejreg)) Manhen.size_sum(Manhen.rejreg[ ,3]) else if (length(Manhen.rejreg)>0) Manhen.size_sum(Manhen.rejreg[3]) else Manhen.size_0


plot(Gamma1[,1],Gamma1[,2],xlab="X11",ylab="Y11",type="n",lab=c(maxx-minx+1,maxy-miny+1,2))
title(main=paste("size=",round(Manhen.size,digits=4)))
if (is.array(Manhen.rejreg)) text(Manhen.rejreg[,1],Manhen.rejreg[,2],Manhen.rejreg[,1]+Manhen.rejreg[,2]) else
    if (length(Manhen.rejreg)>0) text(Manhen.rejreg[1],Manhen.rejreg[2],Manhen.rejreg[1]+Manhen.rejreg[2])
if (is.array(noManhen.rejreg)) points(noManhen.rejreg[ ,1],noManhen.rejreg[ ,2], type="p",pch=1) else
    if (length(noManhen.rejreg)>0) points(noManhen.rejreg[1],noManhen.rejreg[2], type="p",pch=1) 
points(x1,y1, type="p", pch=1, cex=2)

mtext(outer=T,paste("Mantel-Haentszel Test for (",x1,",",x2,",",x3,",",x4,");(",y1,",",y2,",",y3,",",y4,")"),side=3,cex=1.3)
c(Manhen.p,Manhen.size)
}

MH2x2.test(4,14,0,19,11,15,1,25)