## Introduction of the function ##
# The 'fn_GetBeta()' function computes the values of parameters that are identifiable in a case-control study without any additional assumption.

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.

## Output ##
# a vector of (\beta_{e}, \beta_{g}, \beta_{eg})

fn_GetBeta <- function(rctr, rcas)
{
  rctr <- rctr/sum(rctr)
  rcas <- rcas/sum(rcas)
  
  lrctr <- log(rctr)
  lrcas <- log(rcas)
  
  beta0 <- lrcas[1]-lrctr[1]
  beta1 <- lrcas[2]-lrctr[2]-beta0
  beta2 <- lrcas[3]-lrctr[3]-beta0
  beta3 <- lrcas[4]-lrctr[4]-beta0-beta1-beta2
  
  return( c(beta1 = beta1, beta2 = beta2, beta3 = beta3) )
}


## Introduction of the function ##
# The 'fn_LambdaToTheta()' function computes the value(s) of \theta given the case-control cell probabilities and the G-E odds ratio in the source population. 

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# a list containing zero, one, or two values of \theta. 

fn_LambdaToTheta <- function(rctr, rcas, rho=1)
{
  
  rctr <- rctr/sum(rctr)
  rcas <- rcas/sum(rcas)
  rdif <- rcas-rctr
  
  a <- rdif[1]*rdif[4]-rho*rdif[2]*rdif[3] 
  b <- rdif[1]*rctr[4]-rho*rdif[2]*rctr[3]+rdif[4]*rctr[1]-rho*rdif[3]*rctr[2]
  c <- rctr[1]*rctr[4]-rho*rctr[2]*rctr[3]
  d <- b*b-4*a*c
  
  ThetaList <- list()
  
  if(a==0) 
  { 
    theta <- -c/b; 
    if((0<theta)&&(theta<1)) { ThetaList <- c(ThetaList, list(theta)) }  
  }  
  else
  {
    if(d==0) 
    { 
      theta <- -b/(2*a) 
      if((0<theta)&&(theta<1)) { ThetaList <- c(ThetaList, list(theta)) } 
    }  
    if(d>0)  
    { 
      theta <- (-b+sqrt(d))/(2*a)
      if((0<theta)&&(theta<1)) { ThetaList <- c(ThetaList, list(theta)) } 
      
      theta <- (-b-sqrt(d))/(2*a)
      if((0<theta)&&(theta<1)) { ThetaList <- c(ThetaList, list(theta)) } 
    }
  }

  return(ThetaList)
}

## Introduction of the function ##
# The 'fn_LambdaToPhi()' determines the value(s) of \phi given the case-control cell probabilities and the G-E odds ratio in the source population

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). They can also be control/case cell counts if necessary.
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# a list containing zero, one, or two values of \phi

fn_LambdaToPhi <- function(rctr, rcas, rho=1)
{

  rctr <- rctr/sum(rctr)
  rcas <- rcas/sum(rcas)
  rdif <- rcas-rctr
  
  a <- rdif[1]*rdif[4]-rho*rdif[2]*rdif[3]
  b <- rdif[1]*rctr[4]-rho*rdif[2]*rctr[3]+rdif[4]*rctr[1]-rho*rdif[3]*rctr[2]
  c <- rctr[1]*rctr[4]-rho*rctr[2]*rctr[3]
  d <- b*b-4*a*c
  
  PhiList <- list()
  
  if(a==0) 
  { 
    theta <- -c/b
    if((0<theta)&&(theta<1)) 
    { 
      pctr  <- rctr*(1-theta)
      pcas  <- rcas*theta
      lpctr <- log(pctr)
      lpcas <- log(pcas)
      
      prvle <- pctr[2]+pctr[4]+pcas[2]+pcas[4]
      prvlg <- pctr[3]+pctr[4]+pcas[3]+pcas[4]
      beta0 <- lpcas[1]-lpctr[1]
      beta1 <- lpcas[2]-lpctr[2]-beta0
      beta2 <- lpcas[3]-lpctr[3]-beta0
      beta3 <- lpcas[4]-lpctr[4]-beta0-beta1-beta2

      Phi <- c(beta0 = beta0, 
               beta1 = beta1, 
               beta2 = beta2,
               beta3 = beta3, 
               prvle = prvle, 
               prvlg = prvlg, 
               rho   = rho )

      PhiList <- c(PhiList, list(Phi))
    }  
  }  
  else
  {
    if(d==0) 
    { 
      theta <- -b/(2*a)
      if((0<theta)&&(theta<1)) 
      { 
        pctr  <- rctr*(1-theta)
        pcas  <- rcas*theta
        lpctr <- log(pctr)
        lpcas <- log(pcas)
        
        prvle <- pctr[2]+pctr[4]+pcas[2]+pcas[4]
        prvlg <- pctr[3]+pctr[4]+pcas[3]+pcas[4]
        beta0 <- lpcas[1]-lpctr[1]
        beta1 <- lpcas[2]-lpctr[2]-beta0
        beta2 <- lpcas[3]-lpctr[3]-beta0
        beta3 <- lpcas[4]-lpctr[4]-beta0-beta1-beta2
        
        Phi <- c(beta0 = beta0, 
                 beta1 = beta1, 
                 beta2 = beta2,
                 beta3 = beta3, 
                 prvle = prvle, 
                 prvlg = prvlg, 
                 rho   = rho )

        PhiList <- c(PhiList, list(Phi)) 
      }   
    }  
    if(d>0)  
    { 
      theta <- (-b+sqrt(d))/(2*a)
      if((0<theta)&&(theta<1)) 
      { 
        pctr  <- rctr*(1-theta)
        pcas  <- rcas*theta
        lpctr <- log(pctr)
        lpcas <- log(pcas)
        
        prvle <- pctr[2]+pctr[4]+pcas[2]+pcas[4]
        prvlg <- pctr[3]+pctr[4]+pcas[3]+pcas[4]
        beta0 <- lpcas[1]-lpctr[1]
        beta1 <- lpcas[2]-lpctr[2]-beta0
        beta2 <- lpcas[3]-lpctr[3]-beta0
        beta3 <- lpcas[4]-lpctr[4]-beta0-beta1-beta2
        
        Phi <- c(beta0 = beta0, 
                 beta1 = beta1, 
                 beta2 = beta2,
                 beta3 = beta3, 
                 prvle = prvle, 
                 prvlg = prvlg, 
                 rho   = rho )

        PhiList <- c(PhiList, list(Phi)) 
      }
      
      theta <- (-b-sqrt(d))/(2*a)
      if((0<theta)&&(theta<1)) 
      { 
        pctr  <- rctr*(1-theta)
        pcas  <- rcas*theta
        lpctr <- log(pctr)
        lpcas <- log(pcas)
        
        prvle <- pctr[2]+pctr[4]+pcas[2]+pcas[4]
        prvlg <- pctr[3]+pctr[4]+pcas[3]+pcas[4]
        beta0 <- lpcas[1]-lpctr[1]
        beta1 <- lpcas[2]-lpctr[2]-beta0
        beta2 <- lpcas[3]-lpctr[3]-beta0
        beta3 <- lpcas[4]-lpctr[4]-beta0-beta1-beta2
        
        Phi <- c(beta0 = beta0, 
                 beta1 = beta1, 
                 beta2 = beta2,
                 beta3 = beta3, 
                 prvle = prvle, 
                 prvlg = prvlg, 
                 rho   = rho )

        PhiList <- c(PhiList, list(Phi)) 
      }
    }
  }

  return(PhiList)
}

## Introduction of the function ##
# The 'fn_PhiToCellprb()' function performs the transformation from \phi to cell probabilities (Pr[D,G,E]).

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# rho:   the GE odds ratio in the source population, with the default of 1 (GEI assumption)

## Output ##
# a list with two elements:
#     pctr: the vector of the control cell probabilities, Pr[G,E,D=0].
#     pacs: the vector of the case cell probabilities, Pr[G,E,D=1].

fn_PhiToCellprb <- function(Betas, prvle, prvlg, rho=1)
{
  intTerm  <- sqrt((1+(prvle+prvlg)*(rho-1))^2-4*rho*(rho-1)*prvle*prvlg)
  addTerm  <- (4*rho*prvle*prvlg-(rho-1)*(prvle+prvlg)^2-2*(prvle+prvlg))/(2+2*intTerm);
  logitPr0 <- Betas[1]
  logitPr1 <- Betas[1]+Betas[2]
  logitPr2 <- Betas[1]+Betas[3]
  logitPr3 <- Betas[1]+Betas[2]+Betas[3]+Betas[4]
  
  pctr <- c( (1+addTerm-(prvle+prvlg)/2)*(1/(1+exp(logitPr0))),
             ( -addTerm+(prvle-prvlg)/2)*(1/(1+exp(logitPr1))),
             ( -addTerm-(prvle-prvlg)/2)*(1/(1+exp(logitPr2))),
             (  addTerm+(prvle+prvlg)/2)*(1/(1+exp(logitPr3))) )
  
  pcas <- c( (1+addTerm-(prvle+prvlg)/2)*(exp(logitPr0)/(1+exp(logitPr0))),
             ( -addTerm+(prvle-prvlg)/2)*(exp(logitPr1)/(1+exp(logitPr1))),
             ( -addTerm-(prvle-prvlg)/2)*(exp(logitPr2)/(1+exp(logitPr2))),
             (  addTerm+(prvle+prvlg)/2)*(exp(logitPr3)/(1+exp(logitPr3))) )

  return( list(pctr=pctr, pcas=pcas))
}


## Introduction of the function ##
# The 'fn_LogCondtPri()' function computes the logarithm of the prior density for \phi|\rho, which is assumed to be independent of \rho in the manuscript.

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# hypr1: the hyper-parameter, \sigma_{e}, for the prior of \beta_{e}
# hypr2: the hyper-parameter, \sigma_{g}, for the prior of \beta_{g}
# hypr3: the hyper-parameter, \sigma_{eg}, for the prior of \beta_{eg}

## Output ##
# log of the prior density for \phi|\rho

fn_LogCondtPri <- function(Betas, prvle, prvlg, hypr1, hypr2, hypr3)
{
  dns <- ((exp(Betas[1]))/((1+exp(Betas[1]))^2))*
    (exp(-((Betas[2]/hypr1)^2)/2)/(hypr1*sqrt(2*pi)))*
    (exp(-((Betas[3]/hypr2)^2)/2)/(hypr2*sqrt(2*pi)))*
    (exp(-((Betas[4]/hypr3)^2)/2)/(hypr3*sqrt(2*pi)))
  
  return( log(dns) )
}

## Introduction of the function ##
# The 'fn_LogDetJacobian()' function computes log of the absolute value of the Jacobian determinant as given in Appendix A of the manuscript.

## Input Variables ##
# Betas: the vector of (\beta_{0}, \beta_{e}, \beta_{g}, \beta_{eg})
# prvle: the prevalence of environmental exposure Pr(E=1)
# prvlg: the prevalence of genetic variant Pr(G=1)
# rho:   the GE odds ratio in the source population, with the default of 1 (GEI assumption)

## Output ##
# log of the absolute value of the Jacobian determinant

fn_LogDetJacobian <- function(Betas, prvle, prvlg, rho=1)
{
  cellprb <- fn_PhiToCellprb(Betas, prvle, prvlg, rho)
  
  pctr  <- cellprb$pctr
  pcas  <- cellprb$pcas

  theta <- sum(pcas)
  peg   <- pctr+pcas
  
  d1 <- ((prvle+prvlg)*(rho-1)+1)^2-4*rho*(rho-1)*prvle*prvlg
  d2 <- (rho-1)*(prvle+prvlg)^2+2*(prvle+prvlg)-4*rho*prvle*prvlg
  ke <- ((rho+1)*prvlg-(rho-1)*prvle-1)/(2*sqrt(d1))
  kg <- ((rho+1)*prvle-(rho-1)*prvlg-1)/(2*sqrt(d1))
  kr <- d2*(d2+2*prvle*prvlg-prvle-prvlg)/(2*sqrt(d1)*(1+sqrt(d1))^2)-
        (prvle-prvlg)^2/(2+2*sqrt(d1))
  
  u <- c(-1,1,-1,1)
  v <- c(-1,-1,1,1)
  w <- c(1,-1,-1,1)
  
  DevPctr <- matrix(0, 4, 7)
  DevPctr[ ,1] <- (pctr/peg)*(w*ke+u/2)
  DevPctr[ ,2] <- (pctr/peg)*(w*kg+v/2)
  DevPctr[ ,3] <- (pctr/peg)*(w*kr)
  DevPctr[ ,4] <- -pcas*pctr/peg
  DevPctr[2,5] <- DevPctr[2,4]
  DevPctr[4,5] <- DevPctr[4,4]
  DevPctr[3,6] <- DevPctr[3,4]
  DevPctr[4,6] <- DevPctr[4,4]
  DevPctr[4,7] <- DevPctr[4,4]
  
  DevPcas <- matrix(0, 4, 7)
  DevPcas[ ,1] <- (pcas/peg)*(w*ke+u/2)
  DevPcas[ ,2] <- (pcas/peg)*(w*kg+v/2)
  DevPcas[ ,3] <- (pcas/peg)*(w*kr)
  DevPcas[ ,4] <- pcas*pctr/peg
  DevPcas[2,5] <- DevPcas[2,4]
  DevPcas[4,5] <- DevPcas[4,4]
  DevPcas[3,6] <- DevPcas[3,4]
  DevPcas[4,6] <- DevPcas[4,4]
  DevPcas[4,7] <- DevPcas[4,4]
  
  DevTheta <- apply(DevPcas, 2, sum)
    
  Jcob <- matrix( nrow = 7, ncol = 7, byrow = TRUE,
            data = c( (DevPctr[1, ]*(1-theta)+DevTheta*pctr[1])/((1-theta)^2),
                      (DevPctr[2, ]*(1-theta)+DevTheta*pctr[2])/((1-theta)^2),
                      (DevPctr[3, ]*(1-theta)+DevTheta*pctr[3])/((1-theta)^2),
                      (DevPcas[1, ]*(  theta)+DevTheta*pcas[1])/((  theta)^2),
                      (DevPcas[2, ]*(  theta)+DevTheta*pcas[2])/((  theta)^2),
                      (DevPcas[3, ]*(  theta)+DevTheta*pcas[3])/((  theta)^2),
                       DevTheta) )
  
  return( log(abs(det(Jcob))) ) 
}


## Introduction of the function ##
# The 'fn_LogInverseProp()' function computes the factor that accounts for the transformation from (\lambda, \rho) to \xi = (\lambda, \theta), on a logarithmic scale

## Input Variables ##
# rctr/rcas: the control/case cell probabilities (Order {GE: 00, 01, 10, 11}). 
# rho: the GE odds ratio in the source population, with the default being 1 (the GEI assumption)

## Output ##
# the factor that accounts for the transformation from (\lambda, \rho) to \xi = (\lambda, \theta), on a logarithmic scale

fn_LogInverseProp <- function(rctr, rcas, rho=1)
{ 
  rctr <- rctr/sum(rctr)
  rcas <- rcas/sum(rcas)
  rdif <- rcas-rctr
  
  InvProp   <- 0
  ThetaList <- fn_LambdaToTheta(rctr, rcas, rho)
  nsol      <- length(ThetaList)
  
  for(i in 1:nsol)
  {
    theta <- ThetaList[[i]]
    
    tm1 <- (2*rdif[1]*rdif[4]*theta+rdif[1]*rctr[4]+rctr[1]*rdif[4])*(rdif[2]*theta+rctr[2])*(rdif[3]*theta+rctr[3])
    tm2 <- (2*rdif[2]*rdif[3]*theta+rdif[2]*rctr[3]+rctr[2]*rdif[3])*(rdif[1]*theta+rctr[1])*(rdif[4]*theta+rctr[4])
    div <- ((rdif[2]*theta+rctr[2])*(rdif[3]*theta+rctr[3]))^2
    
    InvProp = InvProp+abs(div/(tm1-tm2))
  }
  
  return( log(InvProp) )
}



#################################################
###                MAIN FUNCTIONS             ###
#################################################

## Introduction of the function ##
# The 'fn_MainGeneral()' function samples from the general posterior distribution of \beta_{eg}, which assumes a non-informative conjugate prior.

## Input Variables ##
# dctr/dcas: the control/case cell counts. (Order {GE: 00, 01, 10, 11})
# n: the size of the sample from the general posterior distributon

## Output ##
# a list of two elements:
#     beta3GEN: a sample of size n from the general posterior distribution of \beta_{eg}
#     propNonGEI: the proportion of mass that the general posterior puts outside the GEI space

fn_MainGeneral <- function(dctr, dcas, n)
{
  beta3GEN = nonGEI <- rep(NA, n)

  for(mnlp in 1:n)
  {
    mdctr <- rgamma(4, 1+dctr); mdctr <- mdctr/sum(mdctr)
    mdcas <- rgamma(4, 1+dcas); mdcas <- mdcas/sum(mdcas)

    beta3GEN[mnlp] <- fn_GetBeta(mdctr, mdcas)[3]

    ThetaList = fn_LambdaToTheta(mdctr, mdcas, 1)
    if (length(ThetaList)==0) { nonGEI[mnlp] = 1 } else { nonGEI[mnlp] = 0; } 
  }

  return( list(beta3GEN=beta3GEN, propNonGEI=mean(nonGEI)) )
}

## Introduction of the function ##
# The 'fn_MainRGEI()' function obtains an importance weighted sample to numerically represent the RGEI posterior distribution of \phi.
# Note that the GEI posterior distribution is just a special case where \sigma_{\rho} = 0.

## Input Variables ##
# dctr/dcas: control/case cell counts
# ESS: the desired effective sample size of the importance sample, with the default being 10000
# hypr1: the hyper-parameter, \sigma_{e}, for the prior of \beta_{e}
# hypr2: the hyper-parameter, \sigma_{g}, for the prior of \beta_{g}
# hypr3: the hyper-parameter, \sigma_{eg}, for the prior of \beta_{eg}
# hyprr: the hyper-parameter, \sigma_{\rho}, for the prior of \rho, with the default being 0 (GEI)

## Output ##
# a list of two elements:
#     Phi: A 7-column matrix, with each row being a sample point of \phi
#     Weight: the importance weight associated with the corresponding sample point

fn_MainRGEI <- function(dctr, dcas, ESS=10000, hypr1=2, hypr2=2, hypr3=2, hyprr=0)
{  
  nlp <- 0
  ess <- 0    
  LogWeight <- NULL
  MatrixPhi <- NULL
  
  # dynamically increase the sample size until reach the desired effective sample size
  while( ess<ESS ) 
  { 
    LogWeight <- c(LogWeight, rep(0, 5000))
    MatrixPhi <- rbind(MatrixPhi, matrix(0, 5000, 7))
    
    for (i in 1:5000) 
    {
      nlp <- nlp + 1

      rho  <- exp(rnorm(1, 0, hyprr))       
      rctr <- rgamma(4, dctr+1, 1)
      rctr <- rctr/sum(rctr)     
      rcas <- rgamma(4, dcas+1, 1)
      rcas <- rcas/sum(rcas)
      
      PhiList <- fn_LambdaToPhi(rctr, rcas, rho)
      
      while(length(PhiList)==0)
      {
        rho  <- exp(rnorm(1, 0, hyprr))            
        rctr <- rgamma(4, dctr+1, 1)
        rctr <- rctr/sum(rctr)
        rcas <- rgamma(4, dcas+1, 1)
        rcas <- rcas/sum(rcas)
        
        PhiList <- fn_LambdaToPhi(rctr, rcas, rho) 
      }
      
      nsol <- length(PhiList)
      Phi  <- PhiList[[sample.int(nsol,1)]]      

      MatrixPhi[nlp,] <- Phi
      
      Betas <- Phi[1:4]
      prvle <- Phi[5]
      prvlg <- Phi[6]

      LogWeight[nlp ] <- fn_LogCondtPri(Betas, prvle, prvlg, hypr1, hypr2, hypr3)-
                         fn_LogDetJacobian(Betas, prvle, prvlg, rho)+
                         fn_LogInverseProp(rctr, rcas, rho)
    }

    rWeight <- exp(LogWeight-max(LogWeight))
    Weight  <- rWeight/sum(rWeight)
    ess     <- (sum(Weight))^2/sum(Weight^2)
  }
  
  return(   list( Phi=MatrixPhi, Weight=Weight) )
}