#		Simulate observational study 3*2 SNP association tables

Gen.tbl	<- function(n=1000,p.A = .6,mu.AA = -2,d.Aa.AA = .5,d.aa.Aa = 0)
{

	SNP.seq	<- sample(c("AA","Aa","aa"),n,replace=T,prob = c(p.A,2*p.A*(1-p.A),(1-p.A)^2))
	D.logit.vec	<- c(mu.AA,mu.AA + d.Aa.AA,mu.AA + d.Aa.AA + d.aa.Aa)
	D.prb.vec	<- exp(D.logit.vec) / (1 + exp(D.logit.vec))
	D.prb.seq	<- D.prb.vec[match(SNP.seq,c("AA","Aa","aa"))]
	D.seq		<- rbinom(n,1,D.prb.seq)
	return(table(SNP.seq,D.seq))
}

(a <- Gen.tbl(n=1000,p.A=.6,mu.AA=0,d.Aa.AA=.5,d.aa.Aa=0))
chisq.test(a)
chisq.test(a)$res


(a <- Gen.tbl(n=1000,p.A=.6,mu.AA=.5,d.Aa.AA=-1,d.aa.Aa=0))
chisq.test(a)
chisq.test(a)$res
	

(a <- Gen.tbl(n=1000,p.A=.6,mu.AA=-1,d.Aa.AA=-.5,d.aa.Aa=-.5))
chisq.test(a)
chisq.test(a)$res
	
(a <- Gen.tbl(n=1000,p.A=.6,mu.AA=-1,d.Aa.AA=0,d.aa.Aa=-.5))
chisq.test(a)
chisq.test(a)$res
	

#		Simulate 200K SNP study

#M			<-  2*10^5
#M.rec		<- 
#M.dom		<- 
#M.inc		<- 
#M.dec		<- 

#n.M			<- rep(1000,M)
#p.A.M		<- runif(M,min=0.5,max=0.8)
#mu.AA.M		<- rnorm(M,mean=-1,sd=1)
#d.Aa.AA.M	<- c(rep(0,M.rec),rnorm(M.dom,mean=0,sd=.5),runif(M.inc,min=0.1,max=0.4),runif(M.dec,min=-0.5,max=-0.1),rep(0,M-M.rec-M.dom-M.inc-M.dec))
#d.aa.Aa.M	<- c(rnorm(M.rec,mean=0,sd=1),rep(0,M.dom),runif(M.inc,min=0.1,max=0.4),runif(M.dec,min=-0.5,max=-0.1),rep(0,M-M.rec-M.dom-M.inc-M.dec))


#M.SNP.tbl			<- array(dim=c(3,2,M))
#dimnames(M.SNP.tbl)	<- list(c("AA","Aa","aa"),c("0","1"),paste("SNP",1:M))

for(i in 1:M)
		M.SNP.tbl[,,i]	<- Gen.tbl(n=n.M[i],p.A=p.A.M[i],mu.AA=mu.AA.M[i],d.Aa.AA=d.Aa.AA.M[i],d.aa.Aa=d.aa.Aa.M[i])



dput(M.SNP.tbl,"SNP data")