# Function "chull.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 2x3 contingency table 
#  x1,x2,x3
#  y1,y2,y3
# with ordered outcome levels
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>

chull.test_function(x1,x2,x3,y1,y2,y3)
{
col.m_c(x1+y1,x2+y2,x3+y3)
row.m_c(x1+x2+x3,y1+y2+y3)
v1_rep(c(0:col.m[1]),col.m[2]+1)
v2_rep(c(0:col.m[2]),rep(col.m[1]+1,col.m[2]+1))
v3_row.m[1]-v1-v2
w1_col.m[1]-v1
w2_col.m[2]-v2
w3_col.m[3]-v3
Gamma1_cbind(v1,v2,v3,w1,w2,w3)
Gamma1_Gamma1[v3 >=0 & w3 >=0, ]
 if ((x1+y1) == 0)  d1_1 else d1_dbinom(Gamma1[,1],Gamma1[,1]+Gamma1[,4],.5)
 if ((x2+y2) == 0)  d2_1 else d2_dbinom(Gamma1[,2],Gamma1[,2]+Gamma1[,5],.5)
 if ((x3+y3) == 0)  d3_1 else d3_dbinom(Gamma1[,3],Gamma1[,3]+Gamma1[,6],.5)
nprob_d1*d2*d3/dbinom(x1+x2+x3,x1+x2+x3+y1+y2+y3,.5)
observed_as.numeric(Gamma1[,1] == x1 & Gamma1[,2] == x2)
Gamma1_cbind(Gamma1,nprob,observed)
x_Gamma1[,1]
y_Gamma1[,2]

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])))
if (che != 0) ind_chull(x,y) else if (y[1]-y[length(y)] == 0) ind_c(length(x)) else ind_c(1,length(x))  

	ymax_max(y[ind])
	bmax_match(y[ind],ymax,nomatch="0")
	xmax_max(x[ind]*bmax)
	ymin_min(y[ind])
	bmin_match(y[ind],min(y[ind]),nomatch="0")
	xmin_max(x[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[,8])
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[ ,7]) 

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))
par(mfcol=c(1,2))
plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="n", ylab="C12",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[,9])
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,x2, 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[ ,7]) else if (length(chull.rejreg)>0) scores1.size_sum(Smirnov.rejreg[7]) else chull.size_0

plot(Gamma1[,1],Gamma1[,2],xlab="C11",ylab="C12",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[,9]) else
    if (length(chull.rejreg)>0) text(chull.rejreg[1],chull.rejreg[2],chull.rejreg[9])
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,x2, type="p", pch=1, cex=2)

mtext(outer=T,paste("Convex Hull Test for (",x1,x2,x3,";",y1,y2,y3,")"),side=3,cex=1.5)
c(chull.p,chull.size)
}
chull.test(12,3,7,3,7,12)