#######################################
# function "trio" calls "trio_ftns.R" #
#######################################

# N2 vector with length 2 (n20,n21)
# N4 vector with length 3 (n40,n41,n42)
# N5 vector with length 2 (n51,n52)
# N_bootstrap is the number of simulations to obtain p-values of MAX3 and CLRT
# output (ZREC, pREC) for the test statistic and p-value using REC model
# output (ZADD, pADD) for the test statistic and p-value using ADD model
# output (ZDOM, pDOM) for the test statistic and p-value using DOM model
# output (ZMERT, pMERT) for the test statistic and p-value using MERT
# output (ZMAX3, pMAX3) for the test statistic and p-value using MAX3
# output (ZCLRT, pCLRT) for the test statistic and p-value using CLRT



trio = function(N2,N4,N5,N_bootstrap){

	source("trio_ftns.R")
	
	ZCLRT = ftn_ZCLRT(N2,N4,N5)

	n20 = N2[1]
	n21 = N2[2]

	n40 = N4[1]
	n41 = N4[2]
	n42 = N4[3]

	n51 = N5[1]
	n52 = N5[2]

	n2 = sum(N2)
	n4 = sum(N4)
	n5 = sum(N5)
	
	ZR = ((n42-n4/4)+(n52-n5/2))/sqrt(3*n4/16+n5/4)
	ZA = ((n21-n20)+2*(n42-n40)+(n52-n51))/sqrt(n2+2*n4+n5)
	ZD = ((n21-n2/2)+(n42+n41-3*n4/4))/sqrt(n2/4+3*n4/16)

	corr_ZRA = (2*n4+2*n5)/sqrt(4*n5+3*n4)/sqrt(n2+2*n4+n5)
	corr_ZAR = corr_ZRA

	corr_ZRD = (n4)/sqrt(4*n2+3*n4)/sqrt(3*n4+4*n5)
	corr_ZDR = corr_ZRD

	corr_ZAD = (2*n2+2*n4)/sqrt(4*n2+3*n4)/sqrt(n2+2*n4+n5)
	corr_ZDA = corr_ZAD

	ZMAX = max(ZR,ZA,ZD)
	ZMERT = (ZR+ZD)/sqrt(2*(1+corr_ZRD))

	pR = 1-pnorm(ZR)
	pA = 1-pnorm(ZA)
	pD = 1-pnorm(ZD)
	pMERT = 1-pnorm(ZMERT)

	pMAX = ftn_pMAX(N_bootstrap,ZMAX,corr_ZRA,corr_ZRD,corr_ZAD)
	pCLRT = ftn_pCLRT(N_bootstrap,ZCLRT,corr_ZRD)

	out = list(REC = c(ZR,pR),ADD = c(ZA,pA), DOM = c(ZD,pD), MERT = c(ZMERT, pMERT), MAX3 = c(ZMAX, pMAX), CLRT = c(ZCLRT,pCLRT))

	return(out)

} # end function ftn_TDT

###########
# Example #
###########

n2 = c(286,248)
n4 = c(32,40,37)
n5 = c(6,13)

trio(n2,n4,n5,100000)