# Function "chull2x2.test"  
# plots extreme and rejection regions (one-sided, nominal alpha=0.05)
# for the convex hull 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>

chull2x2.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

probb_0
flag_0
level_1
chulllev_rep(1000,length(v1))
nprob1_nprob

while (flag == 0 || probb < .05)
{
che_sum(abs((x-x[1])*(y[length(y)]-y[1])-(y-y[1])*(x[length(x)]-x[1])))
ind_chull(x,y)

	ymax_max(y[ind])
	bmax_match(y[ind],ymax,nomatch="0")
	xmax_max(x[ind]*bmax)
	xmin_max(x[ind])
	bmin_match(x[ind],xmin,nomatch="0")
	ymin_max(y[ind]*bmin) 
	check_(x[ind]-xmin)*(ymax-ymin)-(y[ind]-ymin)*(xmax-xmin)
	if (sum(check^2) == 0) dirind_ifelse((y[ind]-ymin)^2-(x[ind]-xmin)^2 == 0,ind,0) else	dirind_ifelse(check >= 0,ind,0)

for (k in (1:length(dirind))) if (dirind[k]!=0) chulllev_ifelse((Gamma1[,1]-x[dirind[k]])^2+(Gamma1[,2]-y[dirind[k]])^2==0,level,chulllev)
probb_probb+sum(nprob1[dirind])
flag_max(sum(observed[dirind]),flag)
level_level+1
x_x[-dirind]
y_y[-dirind]
observed_observed[-dirind]
nprob1_nprob1[-dirind]
}

chulllev_-chulllev+level-1
Gamma1_cbind(Gamma1,chulllev)
chullobs_sum(chulllev*Gamma1[,4])
chull.region_Gamma1[chulllev >= chullobs ,]
nochull.region_Gamma1[chulllev < chullobs ,]
if (!is.array(chull.region)) chull.region_t(chull.region)
chull.p_sum(chull.region[,3])

minx_min(Gamma1[,1])
maxx_max(Gamma1[,1])
miny_min(Gamma1[,2])
maxy_max(Gamma1[,2])

par(oma=c(0, 0, 3, 0),mfcol=c(1,2))
plot(Gamma1[,1],Gamma1[,2],xlab="X11",type="n", ylab="Y11",lab=c(maxx-minx+1,maxy-miny+1,2))
title(main=paste("p=",round(chull.p,digits=4)))
text(chull.region[,1],chull.region[,2],chull.region[,5])
if (is.array(nochull.region)) points(nochull.region[ ,1],nochull.region[ ,2], type="p",pch=1) else
    if (length(nochull.region)>0) points(nochull.region[1],nochull.region[2], type="p",pch=1) 
points(x1,y1, type="p", pch=1, cex=2)

size_0
porog_max(chulllev)
while (size <= 0.05) 
{ size_size+sum(nprob[chulllev == porog])
  porog_porog-1
}
chull.rejreg_Gamma1[chulllev > porog+1,]
nochull.rejreg_Gamma1[chulllev <= porog+1,]
if (is.array(chull.rejreg)) chull.size_sum(chull.rejreg[ ,3]) else if (length(chull.rejreg)>0) scores1.size_sum(Smirnov.rejreg[3]) else chull.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(chull.size,digits=4)))
if (is.array(chull.rejreg)) text(chull.rejreg[,1],chull.rejreg[,2],chull.rejreg[,5]) else
    if (length(chull.rejreg)>0) text(chull.rejreg[1],chull.rejreg[2],chull.rejreg[5])
if (is.array(nochull.rejreg)) points(nochull.rejreg[ ,1],nochull.rejreg[ ,2], type="p",pch=1) else
    if (length(nochull.rejreg)>0) points(nochull.rejreg[1],nochull.rejreg[2], type="p",pch=1) 
points(x1,y1, type="p", pch=1, cex=2)

mtext(outer=T,paste("Convex Hull test for (",x1,x2,x3,x4,"),(",y1,y2,y3,y4,")"),side=3,cex=1.3)
c(chull.p,chull.size)
}

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