# Function "ttestsize"  
# plots rejection regions for one-sided (nominal size alpha=.05)
# t-tests with v=0, 0.5, 1, and the Smirnov test
# calculating actual sizes 
# for 2x3 contingency table 
#  x1,x2,x3
#  y1,y2,y3
# with ordered outcome levels
# Author: Anastasia Ivanova <aivanova@bios.unc.edu>

ttestsize_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 (col.m[1] == 0)  d1_1 else d1_dbinom(Gamma1[,1],Gamma1[,1]+Gamma1[,4],.5)
 if (col.m[2] == 0)  d2_1 else d2_dbinom(Gamma1[,2],Gamma1[,2]+Gamma1[,5],.5)
 if (col.m[3] == 0)  d3_1 else d3_dbinom(Gamma1[,3],Gamma1[,3]+Gamma1[,6],.5)
nprob_d1*d2*d3/dbinom(row.m[1],sum(row.m),.5)

#bellow are test statistics for different tests
scores1.test_Gamma1[,3]
scores2.test_Gamma1[,2]+Gamma1[,3]
scores3.test_Gamma1[,1]+2*Gamma1[,2]+3*Gamma1[,3]
Smirnov_pmax(Gamma1[,1]/row.m[1]-Gamma1[,4]/row.m[2],
             (Gamma1[,1]+Gamma1[,2])/row.m[1]-(Gamma1[,4]+Gamma1[,5])/row.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 rejection regions for a size <=.05 test (not randomized)

size_0
teststat_min(scores1.test)
while (size <= 0.05) 
{ size_size+sum(nprob[scores1.test == teststat])
  teststat_teststat+1
}
scores1.rejreg_Gamma1[scores1.test < teststat-1,]
if (is.array(scores1.rejreg)) scores1.size_sum(scores1.rejreg[ ,7]) else if (length(scores1.rejreg)>0) scores1.size_sum(scores1.rejreg[7]) else scores1.size_0


size_0
porog_min(scores2.test)
while (size <= 0.05) 
{ size_size+sum(nprob[scores2.test == porog])
  porog_porog+1
}

scores2.rejreg_Gamma1[scores2.test < porog-1,]
if (is.array(scores2.rejreg)) scores2.size_sum(scores2.rejreg[ ,7]) else if (length(scores2.rejreg)>0) scores2.size_sum(scores2.rejreg[7]) else scores2.size_0


size_0
porog_min(scores3.test)
while (size <= 0.05) 
{ size_size+sum(nprob[scores3.test == porog])
  porog_porog+1
}

scores3.rejreg_Gamma1[scores3.test < porog-1,]
if (is.array(scores3.rejreg)) scores3.size_sum(scores3.rejreg[ ,7]) else if (length(scores3.rejreg)>0) scores3.size_sum(scores3.rejreg[7]) else scores3.size_0

size_0
porogvector_sort(Smirnov)
porog_max(porogvector)
while (size <= 0.05) 
{ size_size+sum(nprob[Smirnov == porog])
  porogvector_porogvector[porogvector<porog]
  porogold_porog
  porog_max(porogvector)
}
Smirnov.rejreg_Gamma1[Smirnov > porogold,]
if (is.array(Smirnov.rejreg)) Smirnov.size_sum(Smirnov.rejreg[ ,7]) else if (length(scores1.rejreg)>0) scores1.size_sum(Smirnov.rejreg[7]) else Smirnov.size_0

# 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
if (is.array(scores1.rejreg)) points(scores1.rejreg[ ,1],scores1.rejreg[ ,2], type="p",pch=16,cex=1.2) else
    if (length(scores1.rejreg)>0) points(scores1.rejreg[1],scores1.rejreg[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 1", "(size=",round(scores1.size,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
if (is.array(scores2.rejreg)) points(scores2.rejreg[ ,1],scores2.rejreg[ ,2], type="p",pch=16,cex=1.2) else
    if (length(scores2.rejreg)>0) points(scores2.rejreg[1],scores2.rejreg[2], type="p",pch=16,cex=1.2) 
title(main=paste("Linear Rank Test 0 1 1","(size=",round(scores2.size,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
if (is.array(scores3.rejreg)) points(scores3.rejreg[ ,1],scores3.rejreg[ ,2], type="p",pch=16,cex=1.2) else
    if (length(scores3.rejreg)>0) points(scores3.rejreg[1],scores3.rejreg[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", "(size=",round(scores3.size,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
if (is.array(Smirnov.rejreg)) points(Smirnov.rejreg[ ,1],Smirnov.rejreg[ ,2], type="p",pch=16,cex=1.2) else
    if (length(Smirnov.rejreg)>0) points(Smirnov.rejreg[1],Smirnov.rejreg[2], type="p",pch=16,cex=1.2) 
points(x1,x2, type="p", pch=1, cex=2)
title(main=paste("Smirnov Test", "(size=",round(Smirnov.size,digits=4),")"))
mtext(outer=T,paste("Data=",x1,",",x2,",",x3,",",y1,",",y2,",",y3),side=1,cex=1.5)

c(scores1.size,scores2.size,scores3.size,Smirnov.size)
}
ttestsize(12,3,7,3,7,12)

 11,15,1,25)