date()
m<-10
## m is the matching times
nc<-300
## nc is the number of case 
pl<-c(0.05,0.1,0.15)
## pl is the allele frequency for each subpopulation
g0l<-(1-pl)^2
g1l<-2*pl*(1-pl)
g2l<-pl^2
## gil=P(Gi/Cl)
f0l<-c(0.01,0.05,0.02)
gama<-2
f2l<-gama*f0l
#f1l<-f0l
#f1l<-(f0l+f2l)/2
f1l<-f2l
## fil=P(case/Gi,Cl)
sl<-c(0.25,0.5,0.25)
## sl is the sampled propotion of each subpopulation
p0l<-f0l*g0l/(f0l*g0l+f1l*g1l+f2l*g2l)
p1l<-f1l*g1l/(f0l*g0l+f1l*g1l+f2l*g2l)
p2l<-f2l*g2l/(f0l*g0l+f1l*g1l+f2l*g2l)
## pil=P(Gi/case,Cl)
q0l<-(1-f0l)*g0l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
q1l<-(1-f1l)*g1l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
q2l<-(1-f2l)*g2l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
## qil=P(Gi/control,cl)
rafa<-0.05
## rafa is the type I error

## Power of single MTT#########################################################################
smual<-function(x){
  a<-sum(sl*(p2l+x*p1l-q2l-x*q1l))
  return(a)
}
ssigmaal<-function(x){
  a<-m*sum(sl*(x^2*p1l+p2l-x^2*p1l^2-p2l^2-2*x*p2l*p1l))+sum(sl*(x^2*q1l+q2l-x^2*q1l^2-q2l^2-2*x*q1l*q2l))
  return(a)
}
ssigma0l<-function(x){
  a<-sum(sl*((x^2*p1l+p2l)+m*(x^2*q1l+q2l)-(m-1)*(x*q1l+q2l)^2-2*(x*p1l+p2l)*(x*q1l+q2l)))
  return(a)
}

power1<-1-pnorm((qnorm(1-rafa/2)*sqrt(ssigma0l(0))-sqrt(nc)*sqrt(m)*smual(0))/sqrt(ssigmaal(0)))+pnorm((-qnorm(1-rafa/2)*sqrt(ssigma0l(0))-sqrt(nc)*sqrt(m)*smual(0))/sqrt(ssigmaal(0)))
power2<-1-pnorm((qnorm(1-rafa/2)*sqrt(ssigma0l(0.5))-sqrt(nc)*sqrt(m)*smual(0.5))/sqrt(ssigmaal(0.5)))+pnorm((-qnorm(1-rafa/2)*sqrt(ssigma0l(0.5))-sqrt(nc)*sqrt(m)*smual(0.5))/sqrt(ssigmaal(0.5)))
power3<-1-pnorm((qnorm(1-rafa/2)*sqrt(ssigma0l(1))-sqrt(nc)*sqrt(m)*smual(1))/sqrt(ssigmaal(1)))+pnorm((-qnorm(1-rafa/2)*sqrt(ssigma0l(1))-sqrt(nc)*sqrt(m)*smual(1))/sqrt(ssigmaal(1)))
## Power of MAX##################################################################################
theta<-function(fol,f1l,f2l){
  p0l<-f0l*g0l/(f0l*g0l+f1l*g1l+f2l*g2l)
  p1l<-f1l*g1l/(f0l*g0l+f1l*g1l+f2l*g2l)
  p2l<-f2l*g2l/(f0l*g0l+f1l*g1l+f2l*g2l)
  q0l<-(1-f0l)*g0l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  q1l<-(1-f1l)*g1l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  q2l<-(1-f2l)*g2l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  smual<-NULL
  ssigma0l<-NULL
  x<-c(0,0.5,1)
  for(i in 1:3){
  smual[i]<-sum(sl*(p2l+x[i]*p1l-q2l-x[i]*q1l))
  ssigma0l[i]<-sum(sl*((x[i]^2*p1l+p2l)+m*(x[i]^2*q1l+q2l)-(m-1)*(x[i]*q1l+q2l)^2-2*(x[i]*p1l+p2l)*(x[i]*q1l+q2l)))
  }
  return(nc*m*smual/sqrt(nc*m*ssigma0l))
}
xigma<-function(f0l,f1l,f2l){
  p0l<-f0l*g0l/(f0l*g0l+f1l*g1l+f2l*g2l)
  p1l<-f1l*g1l/(f0l*g0l+f1l*g1l+f2l*g2l)
  p2l<-f2l*g2l/(f0l*g0l+f1l*g1l+f2l*g2l)
  q0l<-(1-f0l)*g0l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  q1l<-(1-f1l)*g1l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  q2l<-(1-f2l)*g2l/((1-f0l)*g0l+(1-f1l)*g1l+(1-f2l)*g2l)
  sigma0<-function(x){
    a<-nc*m*sum(sl*((x^2*p1l+p2l)+m*(x^2*q1l+q2l)-(m-1)*(x*q1l+q2l)^2-2*(x*p1l+p2l)*(x*q1l+q2l)))
    return(a)
  }
  cova<-function(x1,x2){
    b<-nc*sum(sl*(m^2*(x1*x2*p1l+p2l-x1*x2*p1l^2-(x1+x2)*p1l*p2l-p2l^2)+m*(x1*x2*q1l+q2l-x1*x2*q1l^2-(x1+x2)*q1l*q2l-q2l^2)))
    return(b) 
  } 
  c<-c(cova(0,0)/sigma0(0),cova(0.5,0)/sqrt(sigma0(0.5)*sigma0(0)),cova(1,0)/sqrt(sigma0(1)*sigma0(0)),cova(0,0.5)/sqrt(sigma0(0)*sigma0(0.5)),cova(0.5,0.5)/sigma0(0.5),cova(1,0.5)/sqrt(sigma0(1)*sigma0(0.5)),cova(0,1)/sqrt(sigma0(0)*sigma0(1)),cova(0.5,1)/sqrt(sigma0(0.5)*sigma0(1)),cova(1,1)/sigma0(1))
  return(matrix(c,nrow=3)) 
}
trinorm<-function(n,a,b){
c<-eigen(b)
d<-c$values
e<-c$vectors
f<-e%*%diag(sqrt(abs(d)))%*%t(e)
g<-NULL
for(i in 1:n){
h<-rnorm(3)
i<-f%*%as.matrix(h)
i<-as.vector(i)+a
g<-c(g,i)
}
j<-matrix(data=g,nrow=3,ncol=n)
return(t(j))
}

cv<-function(n,rafa,a,b){
	t<-trinorm(n,a,b)
	c<-rep(0,n)
	for(i in 1:n){
		c[i]=max(abs(t[i,1]),abs(t[i,2]),abs(t[i,3]))
	}
	d<-sort(c)
	result<-d[(1-rafa)*n+1]
	return(result)
}
cvn<-cv(10000,rafa,c(0,0,0),xigma(f0l,f0l,f0l))
## cvn is the critical value under null hypothesis
powermax<-function(n,a,b,c){
  d<-trinorm(n,a,b)
  e<-NULL
  for(i in 1:n ){
    e[i]<-max(abs(d[i,1]),abs(d[i,2]),abs(d[i,3])) 
  }
  f<-mean(e>c)
  return(f)
}

cvn


## Power of MERT######################################################################
erou1<-sum(sl*(p2l*p0l+m*(p2l*q0l+q2l*p0l)+m^2*q2l*q0l))
erou2<-sum(sl*(p2l*(p0l+p1l)+m*(p2l*(q0l+q1l)+q2l*(p0l+p1l))+m^2*q2l*(q0l+q1l)))
erou3<-sum(sl*(p0l*(p1l+p2l)+m*(p0l*(q1l+q2l)+q0l*(p1l+p2l))+m^2*q0l*(q1l+q2l)))
erou<-erou1/sqrt(erou2*erou3)
thetamert<-(theta(f0l,f1l,f2l)[1]+theta(f0l,f1l,f2l)[3])/sqrt(2*(1+erou))
varmert<-(xigma(f0l,f1l,f2l)[1,1]+xigma(f0l,f1l,f2l)[3,3]+2*xigma(f0l,f1l,f2l)[1,3])/(2*(1+erou))
powermert<-1-pnorm((qnorm(1-rafa/2)-thetamert)/sqrt(varmert))+pnorm((-qnorm(1-rafa/2)-thetamert)/sqrt(varmert))

## Power of 2-degree####################################################################
power2d<-function(n,a,b){
  d<-trinorm(n,a,b)
  rou12<-nc*sum(sl*(m*p2l+m^2*q2l-m*(p2l*q1l+p1l*q2l+2*p2l*q2l)-(m^2-m)*(q2l*q1l+q2l^2)))
  rou11<-nc*sum(sl*(m*p2l+m^2*q2l-m*(m-1)*q2l^2-2*m*p2l*q2l))
  rou22<-nc*sum(sl*(m*(p1l+p2l)+m^2*(q1l+q2l)-m*(m-1)*(q1l+q2l)^2-2*m*(p1l+p2l)*(q1l+q2l)))
  rou2d<-rou12/sqrt(rou11*rou22)
  e<-NULL
  for(i in 1:n){
    e[i]<-(d[i,1]^2+d[i,3]^2-2*rou2d*d[i,1]*d[i,3])/(1-rou2d^2)
  }
  f<-mean(e>qchisq(1-rafa,df=2))
  return(f)
}


powermax(10000,theta(f0l,f1l,f2l),xigma(f0l,f1l,f2l),cvn)
#powermert
power2d(10000,theta(f0l,f1l,f2l),xigma(f0l,f1l,f2l))
power1
power2
power3
date()