# Function "ttestpv"  
# plots extreme regions for t-tests with v=0, 0.5, 1, and the Smirnov test
# calculating corresponding p-values
# for 2x3 contingency table 
#  x1,x2,x3
#  y1,y2,y3
# with ordered outcome levels
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>

ttestpv_function(x1,x2,x3,y1,y2,y3)
{
# constracting the reference set, enumerating all the tables with  
# given col.m, column marging, and row.m, row margin
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)
# calculating the null probability of each table
 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)

scores1.test_Gamma1[,3]
scores2.test_Gamma1[,2]+Gamma1[,3]
scores3.test_Gamma1[,1]+2*Gamma1[,2]+3*Gamma1[,3]
Smirnov_pmax(Gamma1[,1],(Gamma1[,1]+Gamma1[,2])-row.m[1]/sum(row.m)*col.m[2],0)
observed_as.numeric(Gamma1[,1] == x1 & Gamma1[,2] == x2)

Gamma1_cbind(Gamma1, nprob, Smirnov, scores1.test, scores2.test,scores3.test,observed)

# calculating p-value
scores1.region_Gamma1[scores1.test <= x3, ]
if (!is.matrix(scores1.region)) scores1.region_matrix(scores1.region,nrow=1)
scores1.p_sum(scores1.region[ ,7])

scores2.region_Gamma1[scores2.test <= x2+x3, ]
if (!is.matrix(scores2.region)) scores2.region_matrix(scores2.region,nrow=1)
scores2.p_sum(scores2.region[ ,7])  

scores3.region_Gamma1[scores3.test <= x1+2*x2+3*x3, ]
if (!is.matrix(scores3.region)) scores3.region_matrix(scores3.region,nrow=1)
scores3.p_sum(scores3.region[ ,7])  

Smirnov.region_Gamma1[Smirnov >= pmax(x1,x1+x2-row.m[1]/sum(row.m)*col.m[2],0),]
if (!is.matrix(Smirnov.region)) Smirnov.region_matrix(Smirnov.region,nrow=1)
Smirnov.p_sum(Smirnov.region[,7])

# plotting the rejection regions for t-test with v=0, .5, and 1
# and Smirnov test
minx_min(Gamma1[,1])
maxx_max(Gamma1[,1])
miny_min(Gamma1[,2])
maxy_max(Gamma1[,2])

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

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="p", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,0))
text(c(minx:maxx), miny-.1*(maxy-miny), c(minx:maxx)) # label x
text(minx-.1*(maxx-minx),c(miny:maxy), c(miny:maxy))  # label y
points(scores1.region[ ,1],scores1.region[ ,2], type="p",pch=16,cex=1.2) #mark points in extreme region
points(x1,x2, type="p", pch=1, cex=2) #mark the observed point
title(main=paste("Linear Rank Test 0 0 1","(p=",round(scores1.p,digits=4),")"))  

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="p", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,0))
text(c(minx:maxx), miny-.1*(maxy-miny), c(minx:maxx)) # label x
text(minx-.1*(maxx-minx),c(miny:maxy), c(miny:maxy))  # label y
points(scores2.region[ ,1],scores2.region[ ,2], type="p",pch=16,cex=1.2)
points(x1,x2, type="p", pch=1, cex=2)
title(main=paste("Linear Rank Test 0 1 1","(p=",round(scores2.p,digits=4),")"))

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="p", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,0))
text(c(minx:maxx), miny-.1*(maxy-miny), c(minx:maxx)) # label x
text(minx-.1*(maxx-minx),c(miny:maxy), c(miny:maxy))  # label y
points(scores3.region[ ,1],scores3.region[ ,2], type="p",pch=16,cex=1.2)
points(x1,x2, type="p", pch=1, cex=2)
title(main=paste("Linear Rank Test 0 0.5 1","(p=",round(scores3.p,digits=4),")"))

plot(Gamma1[,1],Gamma1[,2],xlab="C11",type="p", ylab="C12",lab=c(maxx-minx+1,maxy-miny+1,0))
text(c(minx:maxx), miny-.1*(maxy-miny), c(minx:maxx)) # label x
text(minx-.1*(maxx-minx),c(miny:maxy), c(miny:maxy))  # label y
points(Smirnov.region[ ,1],Smirnov.region[ ,2], type="p",pch=16,cex=1.2)
points(x1,x2, type="p", pch=1, cex=2)
title(main=paste("Smirnov Test","(p=",round(Smirnov.p,digits=4),")"))

mtext(outer=T,paste("Data=",x1,",",x2,",",x3,",",y1,",",y2,",",y3),side=1,cex=1.5)
c(scores1.p,scores2.p,scores3.p,Smirnov.p)
}

ttestpv(12,3,7,3,7,12)