# Function "adaptive.test"  
# plots extreme and rejection regions (one-sided, nominal alpha=0.05)
# for the adaptive test
# calculating p-value and the actual size
# for 2x3 contingency table 
#  x1,x2,x3
#  y1,y2,y3
# with ordered outcome levels
# other parameters are 
# vx: the direction of the treatment effect 
# degree: the level of confidence in vx
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>


adaptive.test_function(x1,x2,x3,y1,y2,y3,vx,degree)
{
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 (!is.matrix(Gamma1)) Gamma1_t(Gamma1)
 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)
dlina_length(Gamma1[,1])

aa_matrix(rep(Gamma1[,1],dlina),dlina,dlina)
bb_matrix(rep(Gamma1[,2],dlina),dlina,dlina)
vv_(bb-t(bb))/(aa-t(aa))
vv_ifelse(vv == Inf,10000,vv)
vv_ifelse(vv == -Inf,10000,vv)
diag(vv)_0
k1_vv				#coeff with zeros
k2_vv-.001			#perturb coeff
k3_ifelse(vv == 0,10000,vv)	#coeff where zeros are changed for 10000

pv1_rep(0,dlina)
pv2_rep(0,dlina)
pv5_rep(0,dlina)

kkk_-1/(1-vx)

for (i in (1:dlina))
{
koef1_ifelse(t(aa-Gamma1[i,1])-(t(bb-Gamma1[i,2]))/k3[i,]>=0,1,0)
pv1[i]_min(apply(rep(nprob,dlina)*t(koef1),2,sum)*((1+abs(vx-1-1/k3[i,]))^degree))
koef2_ifelse(t(aa-Gamma1[i,1])-(t(bb-Gamma1[i,2]))/k2[i,]>=0,1,0)
kv_k1[i,]
kv2_kv
for (j in (1:dlina)) 	kv2[j]_max(ifelse(kv >= kv[j], -1000000, kv))
kv2_ifelse((kv2-kkk)*(kv-kkk) < 0, kkk, kv2)
ko2_pmin((1+abs(vx-1-1/kv))^degree,(1+abs(vx-1-1/kv2))^degree)
pv2[i]_min(apply(rep(nprob,dlina)*t(koef2),2,sum)*ko2)
koef1_ifelse(t(aa-Gamma1[i,1])>=0,1,0)
pv5[i]_min(apply(rep(nprob,dlina)*t(koef1),2,sum))*((1+abs(vx-1))^degree)
}

testpv_pmin(pv5,pv1,pv2)
ranktestpv_rank(-testpv)

observed_as.numeric(Gamma1[,1] == x1 & Gamma1[,2] == x2)

Gamma1_cbind(Gamma1, nprob, testpv, observed,ranktestpv) 
observedpv_Gamma1[observed == 1][8]
obsrank_Gamma1[observed == 1][10]
adap.region_Gamma1[testpv <= observedpv,]
if (is.matrix(adap.region)) adap.p_sum(adap.region[,7]) else if (length(adap.region)>0) adap.p_sum(adap.region[7]) else adap.p_0

size_0
porog_max(ranktestpv)
while (size <= 0.05) 
{ size_size+sum(nprob[ranktestpv == porog])
  porog_porog-1
}
adapsize.region_Gamma1[ranktestpv > porog+1,]
noadapsize.region_Gamma1[ranktestpv <= porog+1,]
if (is.array(adapsize.region)) adap.size_sum(adapsize.region[ ,7]) else { if (length(adapsize.region)>0) 
								adap.size_sum(adapsize.region[7]) else adap.size_0 }

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

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

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="n", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,7))
title(main=paste("p=",round(adap.p,digits=4)))
if (is.array(adap.region)) text(adap.region[ ,1],adap.region[ ,2], 
	adap.region[,10]-maxobpo-1) else text(adap.region[1],adap.region[2], adap.region[10]-maxobpo-1)
	noadap.region_Gamma1[testpv > observedpv,]
if ( is.array(noadap.region) ) points(noadap.region[,1],noadap.region[,2], type="p",pch=1) else
	if (length(noadap.region)>0) points(noadap.region[1],noadap.region[2], type="p",pch=1)
points(x1,x2, type="p", pch=1, cex=2)

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="n", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,7))
title(main=paste("size=",round(adap.size,digits=4)))
if (is.array(adapsize.region)) text(adapsize.region[ ,1],adapsize.region[ ,2], adapsize.region[,10]-maxobpo-1) else
	 if (length(adapsize.region)>0)   text(adapsize.region[1],adapsize.region[2], adapsize.region[10]-maxobpo-1)
if ( is.array(noadapsize.region) ) points(noadapsize.region[,1],noadapsize.region[,2], type="p",pch=1) else
	 if (length(noadapsize.region)>0) points(noadapsize.region[1],noadapsize.region[2], type="p",pch=1)
points(x1,x2, type="p", pch=1, cex=2)

mtext(outer=T,paste("Adaptive p(v)[1+|v-",vx,"|]^",degree," Test for (",x1,x2,x3,";",y1,y2,y3,")"),side=3,cex=1.5)

c(adap.p,adap.size)
}