# Usage: 
# evidences(55,83,50,24,42,39) # The first 3 numbers refere to the genotypic frequencies
# of the case group; the last 3 refere to the genotypic frequencies of the control group

q1<-1/3
q2<-1/4
p2<-2/5
p1<-q1+1/2*q2-1/2*p2

# MC=1: use Monte Carlo method for pvalues
# a1, a2, a3, b1, b2, b3: parameters
# from the prior distribution (dirichlet); default: uniform  
# n1, n2, n3, m1, m2, m3: genotypic distribution
# B: number of iteractions     
evidences <-function(n1,n2,n3,m1,m2,m3,MC=0,
                     a1=1,a2=1,a3=1,b1=1,b2=1,b3=1,B=20000) {

  matriz<-matrix(c(n1,n2,n3,m1,m2,m3),2,3,byrow=T)
  n<-n1+n2+n3
  m<-m1+m2+m3
  a<-a1+a2+a3
  b<-b1+b2+b3

  # Allelic homogeneity: chi2
  r<-nlm(p=c(n2/n,m1/m,m2/m),f=MliklogH,n1=n1,n2=n2,n3=n3,m1=m1,m2=m2,m3=m3)$estimate #ML estimators
  p2c<-r[1]
  q1c<-r[2]
  q2c<-r[3]
  q3c<-1-q2c-q1c
  p1c<-q1c+1/2*q2c-1/2*p2c
  p3c<-1-p1c-p2c
  chiA<-(n1-p1c*n)^2/(p1c*n) + (n2-p2c*n)^2/(p2c*n) +
        (n3-p3c*n)^2/(p3c*n) +(m1-q1c*m)^2/(q1c*m) +
        (m2-q2c*m)^2/(q2c*m) +(m3-q3c*m)^2/(q3c*m)

  matrixAlleles<-matrix(c(2*n1+n2,2*n3+n2,2*m1+m2,2*m3+m2),2,2) # Allelic frequencies table

  # Samples from posterior distribution
  g11<-rgamma(B,shape=a1+n1, scale=1)
  g12<-rgamma(B,shape=a2+n2, scale=1)
  g13<-rgamma(B,shape=a3+n3, scale=1)
  g21<-rgamma(B,shape=b1+m1, scale=1)
  g22<-rgamma(B,shape=b2+m2, scale=1)
  g23<-rgamma(B,shape=b3+m3, scale=1)
  d11<-g11/(g11+g12+g13)
  d12<-g12/(g11+g12+g13)
  d13<-g13/(g11+g12+g13)
  d21<-g21/(g21+g22+g23)
  d22<-g22/(g21+g22+g23)
  d23<-g23/(g21+g22+g23)

  posteriors<-NULL;
  posteriors1<-NULL;
  posteriors2<-NULL;

  
  for (i in 1:B) {
    posteriors[i] <- Mliklog(c(d11[i],d12[i],d13[i],d21[i],d22[i],d23[i]),
                             n1+a1-1,n2+a2-1,n3+a3-1,m1+b1-1,m2+b2-1,m3+b3-1)
    posteriors1[i] <- Mliklog1(c(d11[i],d12[i],d13[i]),n1+a1-1,n2+a2-1,n3+a3-1)
    posteriors2[i] <- Mliklog1(c(d21[i],d22[i],d23[i]),m1+b1-1,m2+b2-1,m3+b3-1)
  }

  # Genotypic homogeneity: e-value 
  evH_G <- sum(posteriors>Mliklog(c((n1+a1-1+m1+b1-1)/(n+a+m+b-6),
                (n2+a2-1+m2+b2-1)/(n+a+m+b-6),
                (n3+a3-1+m3+b3-1)/(n+a+m+b-6),
                (n1+a1-1+m1+b1-1)/(n+a+m+b-6),
                (n2+a2-1+m2+b2-1)/(n+a+m+b-6),
                (n3+a3-1+m3+b3-1)/(n+a+m+b-6)),
                n1=n1+a1-1,n2=n2+a2-1,n3=n3+a3-1,m1=m1+b1-1,m2=m2+b2-1,m3=m3+b3-1))/B

  # Allelic homogeneity: e-value 
  evH_A <- sum(posteriors>nlm(p=c(n2/n,m1/m,m2/m),f=MliklogH,
               n1=n1+a1-1,n2=n2+a2-1,n3=n3+a3-1,m1=m1+b1-1,m2=m2+b2-1,m3=m3+b3-1)$minimum)/B

  # HW equilibrium: case group

  t<-(2*(n1+a1-1)+n2+a2-1)/(2*(n+a-3))

  evalueHWCase<-sum(posteriors1>Mliklog1(c(t^2,2*t*(1-t),(1-t)^2),
                    n1=n1+a1-1,n2=n2+a2-1,n3=n3+a3-1))/B

  pHW<-(2*n1+n2)/(2*n)
  evalueHWCase<-sum(posteriors1>Mliklog1(c(t^2,2*t*(1-t),(1-t)^2),n1=n1+a1-1,n2=n2+a2-1,n3=n3+a3-1))/B
  chiHWCase<-(n1-pHW^2*n)^2/(pHW^2*n) + (n2-2*pHW*(1-pHW)*n)^2/(2*pHW*(1-pHW)*n) + 
             (n3-(1-pHW)^2*n)^2/((1-pHW)^2*n)

  # HW equilibrium: control group
  t<-(2*(m1+b1-1)+m2+b2-1)/(2*(m+b-3))
  evalueHWControl<-sum(posteriors2>Mliklog1(c(t^2,2*t*(1-t),(1-t)^2),
                       n1=m1+b1-1,n2=m2+b2-1,n3=m3+b3-1))/B
  pHW<-(2*m1+m2)/(2*m)
  chiHWControl<-(m1-pHW^2*m)^2/(pHW^2*m) + (m2-2*pHW*(1-pHW)*m)^2/(2*pHW*(1-pHW)*m) + 
                (m3-(1-pHW)^2*m)^2/((1-pHW)^2*m)

  if(MC==0) {
     output<-list("pvalueH_G"=chisq.test(matriz,correct=F)$p.value,
                  "pvalueH_A"=1-pchisq(chiA,1),
                  "pvalueH_AUsual"=chisq.test(matrixAlleles,correct=F)$p.value,
                  "evalueH_G"=evH_G,"evalueH_A"=evH_A,
                  "pvalueHWCase"=1-pchisq(chiHWCase,1),"evalueHWCase"=evalueHWCase,
                  "pvalueHWControl"=1-pchisq(chiHWControl,1),"evalueHWControl"=evalueHWControl)
  }

  if(MC==1) {
     chi2<-sampleChiH_02(n1+n2+n3,m1+m2+m3,B=B);
     output<-list("pvalueH_G"=chisq.test(matriz,simulate.p.value=T,B=B,correct=F)$p.value,
                  "pvalueH_A"=sum(chi2>chiA)/length(chi2),
                  "pvalueH_AUsual"=chisq.test(matrixAlleles, simulate.p.value=T,B=B,correct=F)$p.value,
                  "evalueH_G"=evH_G,"evalueH_A"=evH_A,"pvalueHWCase"=1-pchisq(chiHWCase,1),
                  "evalueHWCase"=evalueHWCase,"pvalueHWControl"=1-pchisq(chiHWControl,1),
                  "evalueHWControl"=evalueHWControl)
  }

   return(output)
}


# minus likelihood log.
Mliklog <- function(prob,n1,n2,n3,m1,m2,m3) {

  x<-c(n1,n2,n3,m1,m2,m3)
  if(any(prob==0) & any(x[prob==0]!=0)) {
    return(Inf);
  } else {
    return(-(n1*log(prob[1]) + n2*log(prob[2]) + n3*log(prob[3]) +
           m1*log(prob[4]) + m2*log(prob[5]) + m3*log(prob[6])));
  }
}

 

# minus likelihood log under H_A
MliklogH <- function(prob,n1,n2,n3,m1,m2,m3) {
  p2<-prob[1];
  q1<-prob[2];
  q2<-prob[3];

  if((n1==0 & (q1+1/2*q2-1/2*p2)==0) | (n2==0 & p2==0) | 
     (n3==0 & (1-q1-1/2*q2-1/2*p2)==0) | (m1==0 & q1==0) | 
     (m2==0 & q2==0) | (m3==0 & (1-q1-q2)==0)) {
    return(Inf);
  } else {
    return(-(n1*log(q1+1/2*q2-1/2*p2) + n2*log(p2) + n3*log(1-q1-1/2*q2-1/2*p2) + 
           m1*log(q1) + m2*log(q2) + m3*log(1-q1-q2)));
  }
}

sampleChiH_02 <-function(n,m,B=B) {

  chi<-NULL;
  for (i in 1:B) {
    x<-runif(n,0,1)
    y<-runif(m,0,1)
    n1<-sum(x<p1)
    n3<-sum(x>p1+p2)
    n2<-n-n1-n3
    m1<-sum(y<q1)
    m3<-sum(y>q1+q2)
    m2<-m-m1-m3

    r<-nlm(p=c(n2/n,m1/m,m2/m),f=MliklogH,n1=n1,n2=n2,n3=n3,m1=m1,m2=m2,m3=m3)$estimate
    p2c<-r[1]
    q1c<-r[2]
    q2c<-r[3]
    q3c<-1-q2c-q1c
    p1c<-q1c+1/2*q2c-1/2*p2c
    p3c<-1-p1c-p2c

    chi[i]<-(n1-p1c*n)^2/(p1c*n) + (n2-p2c*n)^2/(p2c*n) +(n3-p3c*n)^2/(p3c*n) +
            (m1-q1c*m)^2/(q1c*m) +(m2-q2c*m)^2/(q2c*m) +(m3-q3c*m)^2/(q3c*m)
   }
   return(chi);

}

 

 

Mliklog1 <- function(prob,n1,n2,n3) {

  x<-c(n1,n2,n3)
  if(any(prob==0) & any(x[prob==0]!=0)) {
    return(Inf);
  } else {
    return(-(n1*log(prob[1]) + n2*log(prob[2]) + n3*log(prob[3])));
  }

}