simulate.study	<- function(n.trtmnt.1 = 100, n.trtmnt.2 = 100, n.placebo = 150,
 p.DD = 0.4, p.succ.placebo = 0.2, p.succ.trt.not.DD = 0.25, p.succ.trt.DD = 0.4)
{
	p.D	<- sqrt(p.DD)
	p.DI	<- 2 * p.D * (1 - p.D)
	p.II	<- (1 - p.D)^2
	
	result.table		<- array(dim=c(3,3,2))
	dimnames(result.table)	<- list(c("Treatment 1","Treatment 2","Placebo"),c("DD","DI","II"),c("Success","Failure"))
	
	n.trtmnt.1.genome	<- rmultinom(1,n.trtmnt.1,c(p.DD,p.DI,p.II))
	n.trtmnt.2.genome	<- rmultinom(1,n.trtmnt.2,c(p.DD,p.DI,p.II))
	n.placebo.genome	<- rmultinom(1,n.placebo,c(p.DD,p.DI,p.II))

#	Treatment 1 table entries

	result.table[1,1,1]	<- rbinom(1,n.trtmnt.1.genome[1],p.succ.trt.DD)
	result.table[1,1,2]	<- n.trtmnt.1.genome[1] - result.table[1,1,1]

	result.table[1,2,1]	<- rbinom(1,n.trtmnt.1.genome[2],p.succ.trt.not.DD)
	result.table[1,2,2]	<- n.trtmnt.1.genome[2] - result.table[1,2,1]

	result.table[1,3,1]	<- rbinom(1,n.trtmnt.1.genome[3],p.succ.trt.not.DD)
	result.table[1,3,2]	<- n.trtmnt.1.genome[3] - result.table[1,3,1]

#	Treatment 2 table entries

	result.table[2,1,1]	<- rbinom(1,n.trtmnt.2.genome[1],p.succ.trt.DD)
	result.table[2,1,2]	<- n.trtmnt.2.genome[1] - result.table[2,1,1]

	result.table[2,2,1]	<- rbinom(1,n.trtmnt.2.genome[2],p.succ.trt.not.DD)
	result.table[2,2,2]	<- n.trtmnt.2.genome[2] - result.table[2,2,1]

	result.table[2,3,1]	<- rbinom(1,n.trtmnt.2.genome[3],p.succ.trt.not.DD)
	result.table[2,3,2]	<- n.trtmnt.2.genome[3] - result.table[2,3,1]


#	Placebo table entries

	result.table[3,1,1]	<- rbinom(1,n.placebo.genome[1],p.succ.placebo)
	result.table[3,1,2]	<- n.placebo.genome[1] - result.table[3,1,1]

	result.table[3,2,1]	<- rbinom(1,n.placebo.genome[2],p.succ.placebo)
	result.table[3,2,2]	<- n.placebo.genome[2] - result.table[3,2,1]

	result.table[3,3,1]	<- rbinom(1,n.placebo.genome[3],p.succ.placebo)
	result.table[3,3,2]	<- n.placebo.genome[3] - result.table[3,3,1]

	return(result.table)
}


(result.table <- simulate.study())

sum(result.table)
apply(result.table,c(1,3),sum)
apply(result.table,c(1,2),sum)

result.table[,,1]
result.table[,,1] / apply(result.table,c(1,2),sum)


genome.table		<- array(dim=c(2,2,3))
dimnames(genome.table)	<- list(c("Treatment","Placebo"),c("Success","Failure"),c("DD","DI","II"))
genome.table[,,1]		<- rbind(result.table[1,1,] + result.table[2,1,],result.table[3,1,])
genome.table[,,2]		<- rbind(result.table[1,2,] + result.table[2,2,],result.table[3,2,])
genome.table[,,3]		<- rbind(result.table[1,3,] + result.table[2,3,],result.table[3,3,])

apply(genome.table,1,sum)

fisher.test(genome.table[,,1])	
fisher.test(genome.table[,,2])	
fisher.test(genome.table[,,3])	

apply(genome.table,c(1,2),sum)
fisher.test(apply(genome.table,c(1,2),sum))
mantelhaen.test(genome.table)

treatment.table			<- array(dim=c(2,2,3))
dimnames(treatment.table)	<- list(c("DD","not.DD"),c("Success","Failure"),c("Treatment 1","Treatment 2","Placebo"))
treatment.table[,,1]		<- rbind(result.table[1,1,], result.table[1,2,] + result.table[1,3,])
treatment.table[,,2]		<- rbind(result.table[2,1,], result.table[2,2,] + result.table[2,3,])
treatment.table[,,3]		<- rbind(result.table[3,1,], result.table[3,2,] + result.table[3,3,])

fisher.test(treatment.table[,,1])	
fisher.test(treatment.table[,,2])	
fisher.test(treatment.table[,,3])	
mantelhaen.test(treatment.table)


analyze.data	<- function(result.table)
{
	genome.table		<- array(dim=c(2,2,3))
	dimnames(genome.table)	<- list(c("Treatment","Placebo"),c("Success","Failure"),c("DD","DI","II"))
	genome.table[,,1]		<- rbind(result.table[1,1,] + result.table[2,1,],result.table[3,1,])
	genome.table[,,2]		<- rbind(result.table[1,2,] + result.table[2,2,],result.table[3,2,])
	genome.table[,,3]		<- rbind(result.table[1,3,] + result.table[2,3,],result.table[3,3,])

	mh.treatment.p	<- mantelhaen.test(genome.table)$p.value

	treatment.table			<- array(dim=c(2,2,3))
	dimnames(treatment.table)	<- list(c("DD","not.DD"),c("Success","Failure"),c("Treatment 1","Treatment 2","Placebo"))
	treatment.table[,,1]		<- rbind(result.table[1,1,], result.table[1,2,] + result.table[1,3,])
	treatment.table[,,2]		<- rbind(result.table[2,1,], result.table[2,2,] + result.table[2,3,])
	treatment.table[,,3]		<- rbind(result.table[3,1,], result.table[3,2,] + result.table[3,3,])

	mh.genome.p	<- 	mantelhaen.test(treatment.table)$p.value

	return(list(mh.treatment.p,mh.genome.p))
}

analyze.data(result.table)



########################################################################################

# Power simulation .....


trtmnt.p	<- rep(NA,1000)
genome.p	<- rep(NA,1000)

for(i in 1:1000)
{
	tbl	<-  simulate.study(n.trtmnt.1 = 100, n.trtmnt.2 = 100, n.placebo = 150, p.DD = 0.4, p.succ.placebo = 0.2, p.succ.trt.not.DD = 0.25, p.succ.trt.DD = 0.4)
	aa	<- analyze.data(tbl)

	trtmnt.p[i]	<- aa[[1]]
	genome.p[i]	<- aa[[1]]
}

mean(genome.p < 0.02 & trtmnt.p < 0.02)












#########################################################################






#	Study distribution of statistic testing diff in proportions and Fisher exact test statistic 


study.dist <- function(p1=.2,p2=.2,n1=100,n2=100,qq=NA)
{
	m	<- 10^4
	x1	<- rbinom(m,n1,p1)
	x2	<- rbinom(m,n2,p2)

	p1.hat<- x1/n1
	p2.hat<- x2/n2

	stat.diff		<- (p2.hat  - p1.hat) / sqrt( p2.hat*(1-p2.hat)/n2 + p1.hat*(1-p1.hat)/n1)
	mn.stat.diff	<- (p2  - p1) / sqrt( p2*(1-p2)/n2 + p1*(1-p1)/n1)

	theta.hat		<- (x2/(n2-x2)) / (x1/(n1-x1))
	var.theta.hat	<- (1/x1 + 1/(n1-x1) + 1/x2 + 1/(n2-x2))

	stat.fish		<- log(theta.hat) / sqrt(var.theta.hat)
	mn.stat.fish	<- log( (p2/(1-p2)) / (p1/(1-p1))) / sqrt( 1/(p1*n1) + 1/((1-p1)*n1) + 1/(p2*n2) + 1/((1-p2)*n2))

	pow.diff		<- mean(qq <= stat.diff)
	pval.fish		<- (phyper(x2,n2,n1,x1+x2,lower.tail = F) + phyper(x2-1,n2,n1,x1+x2,lower.tail = F))/2
	pow.fish		<- mean(qq <= qnorm(1 - pval.fish))

	
	if(is.na(qq))
		plot(density(stat.diff),col=3, ylim=c(0,.5))
	else
	{
		plot(density(stat.diff),col=3, ylim=c(0,.5),
			main = paste("Power: Diff.z = ",round(pow.diff,3),", OR.exact = ",round(pow.fish,3)))
		abline(v = qq,col=2)
	}

	browser()
	lines(density(stat.fish),col=4)
	curve(dnorm(x, mean=mn.stat.diff,sd=1), col = 3, lty = 2, lwd = 2, add = TRUE)
	curve(dnorm(x, mean=mn.stat.fish,sd=1), col = 4, lty = 2, lwd = 2, add = TRUE)
}


study.dist(p1=.3,p2=.3,n1=100,n2=100)
study.fisher.dist(p1=.3,p2=.3,n1=100,n2=100)


study.dist(p1=.15,p2=.4,n1=100,n2=100,qq = qnorm(.98))


study.dist(p1=.3,p2=.3,n1=100,n2=100)
study.dist(p1=.35,p2=.35,n1=100,n2=100)
study.dist(p1=.2,p2=.3,n1=100,n2=100)

study.dist(p1=.15,p2=.4,n1=100,n2=100)
study.dist(p1=.15,p2=.4,n1=100,n2=100,qq = qnorm(.98))
study.dist(p1=.15,p2=.4,n1=50,n2=50,qq = qnorm(.98))

study.dist(p1=.15,p2=.4,n1=60,n2=40,qq = qnorm(.98))
study.dist(p1=.15,p2=.4,n1=66,n2=44,qq = qnorm(.98))

#  Sample size of 140 patients given drugs 1 & 2 needed to ensure 80% power for both drugs!

study.dist(p1=.15,p2=.4,n1=84,n2=56,qq = qnorm(.98))