#####################################################################################################################
# This file uses the functions defined in "Functions_Rcpp.cpp" or "Functions_Rbase.r" to perform analyses presented in the manuscript.
# "Functions_Rcpp.cpp" are written in C++ language, which can be integrated with R by the 'Rcpp' package (our code also requires the "RcppArmadillo" package).
# We recommend to source "Functions_Rcpp.cpp", wihch is at least 10 times faster than sourcing "Functions_Rbase.r".
# The first part implements several simulation studies and saves the results for further use.
# The second part reproduces figures and tables in the manuscript and the Supplemental Material.


##########################################################################################
### Part I: Conduct simulation studies and save the results for further use
##########################################################################################

library(Rcpp)
sourceCpp("Functions_Rcpp.cpp")
# source("Functions_Rbase.r") 

#############################
### An auxiliary function ###
#############################

## Introduction of the function ##
# The 'fn_Simlt()' function simulates datasets and applies the general method, the GEI method and the RGEI method to the 
# generated datasets, for different combinations of parameter setting and sample size. 

## Input Variables ##
# Setting: A J*7 matrix that specifies J parameter settings, with each row being (\beta_{0},\beta_{e},\beta_{g},\beta_{eg},\pi_{e},\pi_{g},\rho).
# Size: A K*2 matrix that specifies K sample sizes, with the first column being the number of controls and the second being the number of cases.
# NREP: The number of datasets to be generated for each combination of parameter setting and sample size.
# Hypers: The prior hyper-parameters, (\sigma_{e},\sigma_{g},\sigma_{eg},\sigma_{\rho}), for the GEI/RGEI method. Note the GEI method fixes \sigma_{\rho}=0.

## Output ##
# A list 'ResList', with ResList[[j]][[k]] being the simulation result for the j-th parameter setting and k-th sample size
# ResList[[j]][[k]] is a list itself, containing the following elements:
#   GEN: a 3-column matrix for the general estimates. The first column is the point estimate, and the last two columns are the 95% CI
#   GEI: a 3-column matrix for the GEI estimates. The first column is the point estimate, and the last two columns are the 95% CI
#   RGEI: a 3-column matrix for the RGEI estimates. The first column is the point estimate, and the last two columns are the 95% CI
#   NumTheta: the number (0,1,or 2) of \theta that 'fit' the generated dataset under the GEI assumption
#   PropNonGEI: the proportion of mass that the general posterior distribution puts outside the GEI space

fn_Simlt <- function(Setting, Size, NREP, Hypers)
{
  J <- nrow(Setting)
  K <- nrow(Size) 
  ResList <- NULL

  for(j in 1:J)
  {
    Betas <- Setting[j,1:4]
    prvle <- Setting[j,5]
    prvlg <- Setting[j,6]
    rho   <- Setting[j,7]

    Cellprb <- fn_PhiToCellprb(Betas, prvle, prvlg, rho)
    rctr <- Cellprb$pctr/sum(Cellprb$pctr)
    rcas <- Cellprb$pcas/sum(Cellprb$pcas)
    
    tmpList <- NULL
    for(k in 1:K)
    {
      nctr <- Size[k,1]
      ncas <- Size[k,2]
      
      GEN  <- matrix(NA, NREP, 3)
      GEI  <- matrix(NA, NREP, 3)
      RGEI <- matrix(NA, NREP, 3)
      NumTheta   <- rep(NA, NREP)
      PropNonGEI <- rep(NA, NREP)
 
      for(nlp in 1:NREP)
      {
        print(c(j,k,nlp))
        
        dctr <- c(rmultinom(1, nctr, rctr))
        dcas <- c(rmultinom(1, ncas, rcas))
        NumTheta[nlp] <- length(fn_LambdaToTheta(dctr, dcas, rho))
        
        # General method
        resGEN <- fn_MainGeneral(dctr, dcas, 10000)
        PropNonGEI[nlp] <- resGEN$propNonGEI
        
        beta3GEN <- resGEN$beta3GEN
        GEN[nlp, ] <- c(mean(beta3GEN), quantile(beta3GEN, probs=c(0.025, 0.975)))
        
        # GEI method
        resGEI <- fn_MainRGEI(dctr, dcas, 10000, Hypers[1], Hypers[2], Hypers[3], 0)
        
        beta3GEI <- with(resGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
        GEI[nlp, ] <- c(mean(beta3GEI), quantile(beta3GEI, probs=c(0.025, 0.975)))
        
        # RGEI method    
        resRGEI <- fn_MainRGEI(dctr, dcas, 10000, Hypers[1], Hypers[2], Hypers[3], Hypers[4])
        
        beta3RGEI <- with(resRGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
        RGEI[nlp, ] <- c(mean(beta3RGEI), quantile(beta3RGEI, probs=c(0.025, 0.975)))

        tmpList[[k]] <- list(GEN = GEN, GEI = GEI, RGEI = RGEI, NumTheta = NumTheta, PropNonGEI = PropNonGEI)
      }

      ResList[[j]] <- tmpList
    } 
  }
  
  return(ResList) 
}


#############################
### Simulation Studies    ###
#############################

# This simulation study concerns seven parameter settings and there sample sizes. 
# Its results are presented in the Supplemental Material
Setting <- rbind( c(log(0.05/0.95), log(1.2), 0       , log(1.5), 0.1, 0.1, 1),
                  c(log(0.05/0.95), log(1.2), 0       , log(1.5), 0.1, 0.3, 1),
                  c(log(0.05/0.95), log(1.2), 0       , log(1.5), 0.3, 0.1, 1),
                  c(log(0.05/0.95), log(1.2), 0       , log(1.5), 0.3, 0.3, 1),
                  c(log(0.05/0.95), log(1.2), log(1.2), log(1.5), 0.1, 0.1, 1),
                  c(log(0.05/0.95), log(1.2), log(1.2), log(1.5), 0.1, 0.3, 1),
                  c(log(0.05/0.95), log(1.2), log(1.2), log(1.5), 0.3, 0.3, 1))
Size <- matrix(rep(c(400,2000,4000), 2), 3, 2)
ResList <- fn_Simlt(Setting, Size, NREP=2000, Hypers=c(2.5, 2.5, 3.5, 0.1))
save(ResList, file ="Simlt_AllSettings.RData")

# This simulation study focuses on the fourth parameter setting of the previous simulation study, which is the one used in the manuscript.
# It also considers three more larger sample sizes (20000/200000/2000000) to adress the asymptotic behavior of the GEI method.
# For consistency, the result for the three smaller sample sizes (800/4000/8000) are directly extracted from the previous simulation.
load("Simlt_AllSettings.RData")
tmp1 <- ResList[[4]]
rm(ResList)
Setting <- rbind( c(log(0.05/(1-0.05)), log(1.2), 0, log(1.5), 0.3, 0.3, 1) )
Size <- matrix(rep(c(10000,100000,1000000), 2), 3, 2)
tmp2 <- fn_Simlt(Setting, Size, NREP=2000, Hypers=c(2.5, 2.5, 3.5, 0.1))[[1]]
ResList <- c(tmp1, tmp2)
save(ResList, file="Simlt_MainSetting.RData")

# This simulation study investigates the impact of the size of \beta_{eg}
Setting <- rbind( c(log(0.05/0.95), log(1.2), 0,  log(1.0), 0.3, 0.3, 1),
                   c(log(0.05/0.95), log(1.2), 0,  log(1.5), 0.3, 0.3, 1),                  
                   c(log(0.05/0.95), log(1.2), 0,  log(3.0), 0.3, 0.3, 1) )
Size <- matrix(rep(c(400,2000,4000), 2), 3, 2)
ResList <- fn_Simlt(Setting, Size, NREP=2000, Hypers=c(2.5, 2.5, 3.5, 0.1))
save(ResList, file="Simlt_beta3Strength.RData")

# This simulation study compares the general method, the GEI method and the RGEI method when the GEI assumption is violated.
Setting <- rbind( c(log(0.05/0.95), log(1.2), 0,  log(1.5), 0.3, 0.3, 1.1) )
Size <- matrix(rep(4000, 2), 1, 2)
ResList <- fn_Simlt(Setting, Size, NREP=2000, Hypers=c(2.5, 2.5, 3.5, 0.1))
save(ResList, file="Simlt_nonGEI.RData")




##########################################################################################
### Part II: Reproduce figures and tables in the manuscript and the Supplemental Material
##########################################################################################

#######################################
###  Tables in the main manuscript  ###
#######################################

################
### Table 2  ###
################

# setwd(Path4RData)
load("Simlt_MainSetting.RData") 
MSE = CPB <- NULL

for(k in 1:3)
{ 
  res <- ResList[[k]]    
  MSE <- with(res, rbind(MSE, c(mean((GEN[,1]-log(1.5))^2),
                                mean((GEI[,1]-log(1.5))^2),
                                mean((RGEI[,1]-log(1.5))^2) ) ) )
  CPB <- with(res, rbind(CPB, c(mean(GEN[,2]<log(1.5)&log(1.5)<GEN[,3]),
                                mean(GEI[,2]<log(1.5)&log(1.5)<GEI[,3]),
                                mean(RGEI[,2]<log(1.5)&log(1.5)<RGEI[,3]) ) ) ) 
}
print(MSE) # Mean square error of the point estimate
print(CPB) # Coverage probability of the interval estimate

################
### Table 3  ###
################

# setwd(Path4RData)
load("Simlt_nonGEI.RData") 
res <- ResList[[1]][[1]]
MSE <- with(res, c(mean((GEN[,1]-log(1.5))^2), 
                   mean((GEI[,1]-log(1.5))^2),
                   mean((RGEI[,1]-log(1.5))^2) ) ) 
CPB <- with(res, c(mean(GEN[,2]<log(1.5)&log(1.5)<GEN[,3]),
                   mean(GEI[,2]<log(1.5)&log(1.5)<GEI[,3]),
                   mean(RGEI[,2]<log(1.5)&log(1.5)<RGEI[,3]) ) ) 
print(MSE) # Mean square error of the point estimate
print(CPB) # Coverage probability of the interval estimate

################
### Table 5  ###
################

# If package 'CGEN' has not been installed, run the following two lines:
#   source("http://bioconductor.org/biocLite.R")
#   biocLite("CGEN")
library(CGEN) 
library(Rcpp)
sourceCpp("Functions_Rcpp.cpp")
# source("Functions_Rbase.r") 

#############
### Dataset 1 
datCtr <- c(167, 34, 69, 11)
datCas <- c(60, 12, 32, 9)

# General
resGEN <- fn_MainGeneral(datCtr, datCas, 10000)
beta3GEN <- resGEN$beta3GEN
GEN <- c(mean(beta3GEN), quantile(beta3GEN, probs=c(0.025, 0.975)))
# GEI 
resGEI <- fn_MainRGEI(datCtr, datCas, 10000, 2.5, 2.5, 3.5, 0)
beta3GEI <- with(resGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
GEI <- c(mean(beta3GEI), quantile(beta3GEI, probs=c(0.025, 0.975)))
# RGEI (\sigma_{\rho} = 0.1)
resRGEI <- fn_MainRGEI(datCtr, datCas, 10000, 2.5, 2.5, 3.5, 0.01)
beta3RGEI <- with(resRGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
RGEI <- c(mean(beta3RGEI), quantile(beta3RGEI, probs=c(0.025, 0.975)))
# Empirical Bayes
data <- data.frame(D = c(rep(c(0,0,0,0), datCtr), rep(c(1,1,1,1), datCas)),
                   G = c(rep(c(0,0,1,1), datCtr), rep(c(0,0,1,1), datCas)),
                   E =c(rep(c(0,1,0,1), datCtr), rep(c(0,1,0,1), datCas)) )
resCGEN <- getSummary(snp.logistic(data, "D", "G", main.vars=c("E"), int.vars=c("E")))
EB <- resCGEN$EB[4,1]+c(0,-1.96,1.96)*resCGEN$EB[4,2]
# Case-Only
CO <- log(datCas[1]*datCas[4]/datCas[2]/datCas[3])+c(0,-1.96,1.96)*sqrt(sum(1/datCas))
  
print(rbind(GEN, GEI, RGEI, EB, CO))

#############
### Dataset 2
datCtr <- c(172, 230, 58, 52)
datCas <- c(106, 156, 83, 157)
# General
resGEN <- fn_MainGeneral(datCtr, datCas, 10000)
beta3GEN <- resGEN$beta3GEN
GEN <- c(mean(beta3GEN), quantile(beta3GEN, probs=c(0.025, 0.975)))
# GEI 
resGEI <- fn_MainRGEI(datCtr, datCas, 10000, 2.5, 2.5, 3.5, 0)
beta3GEI <- with(resGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
GEI <- c(mean(beta3GEI), quantile(beta3GEI, probs=c(0.025, 0.975)))
# RGEI (\sigma_{\rho} = 0.1)
resRGEI <- fn_MainRGEI(datCtr, datCas, 10000, 2.5, 2.5, 3.5, 0.01)
beta3RGEI <- with(resRGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
RGEI <- c(mean(beta3RGEI), quantile(beta3RGEI, probs=c(0.025, 0.975)))
# Empirical Bayes
data <- data.frame(D = c(rep(c(0,0,0,0), datCtr), rep(c(1,1,1,1), datCas)),
                   G = c(rep(c(0,0,1,1), datCtr), rep(c(0,0,1,1), datCas)),
                   E =c(rep(c(0,1,0,1), datCtr), rep(c(0,1,0,1), datCas)) )
resCGEN <- getSummary(snp.logistic(data, "D", "G", main.vars=c("E"), int.vars=c("E")))
EB <- resCGEN$EB[4,1]+c(0,-1.96,1.96)*resCGEN$EB[4,2]
# Case-Only
CO <- log(datCas[1]*datCas[4]/datCas[2]/datCas[3])+c(0,-1.96,1.96)*sqrt(sum(1/datCas))
  
print(rbind(GEN, GEI, RGEI, EB, CO))


#######################################
###  Figures in the main manuscript ###
#######################################

################
### Figure 1 ###
################

# setwd(Path4RData)
options(scipen=100)
load("Simlt_MainSetting.RData")
# setwd(Path4Figure)
postscript("Main-Figure_CompareLength.eps", height=5, width=7, horizontal=FALSE)
op <- par(mfrow = c(2,3), mar = c(4,4,4,2)+0.1)

for(k in 1:6)
{
  res <- ResList[[k]]
  GEN <- res$GEN[,3] - res$GEN[,2]
  GEI <- res$GEI[,3] - res$GEI[,2]
  
  col <- rep("black", length(GEN))
  col[res$NumTheta!=0] <- "gray70"
  xylim <- range(GEN, GEI)
  
  plot(GEN ~ GEI, col=col, pch=20, xlim=xylim, ylim=xylim, 
       main = paste("n =", c(800, 4000, 8000, 20000, 200000, 2000000)[k]) )
  abline(a=0,b=1)
}

dev.off()

################
### Figure 2 ###
################

# setwd(Path4RData)
options(scipen=100)
load("Simlt_MainSetting.RData")
# setwd(Path4Figure)
postscript("Main-Figure_LengthRatio.eps", height=3, width=7, horizontal=FALSE)
op <- par(mfrow = c(1,3), mar = c(4,4,6,2)+0.1)

for(k in 1:3)
{ 
  ### CI length ratio: GEI/GEN
  res <- ResList[[k]]
  LengthRatio <- with(res, (GEI[,3]-GEI[,2])/(GEN[,3]-GEN[,2]))
  PropNonGEI  <- 100*res$PropNonGEI    
  
  col <- rep("black", length(LengthRatio))
  col[res$NumTheta!=0] <- "gray70"
  
  plot(LengthRatio ~ PropNonGEI, col=col, pch=20, xlim = c(0,100), ylim = c(0.5,2), log="y",
       xlab = "% mass outside GEI-space", ylab = "Length Ratio", main = "", cex = 0.5) 
  mtext(paste("n = ", c(800, 4000, 8000)[k]), side=3, line=1)
  if(k == 2) { mtext("Length Ratio vs % mass outside GEI-space", side=3, line=4) }
}

dev.off()

################
### Figure 3 ###
################

# setwd(Path4RData)
options(scipen=100)
load("Simlt_beta3Strength.RData")
# setwd(Path4Figure)
postscript("Main-Figure_BetaStrength.eps", height=8, width=9, horizontal=FALSE)
op = par(mfrow = c(3,3), mar = c(4,4,4,2)+0.1)

for(j in 1:3)
{  
  for(k in 1:3)
  {   
    ### CI length ratio: GEI/GEN
    res <- ResList[[j]][[k]]    
    LengthRatio <- with(res, (GEI[,3] - GEI[,2])/(GEN[,3] - GEN[,2]) )
    PropNonGEI  <- 100*res$PropNonGEI   
    
    col <- rep("black", length(LengthRatio))
    col[res$NumTheta!=0] <- "gray70"
  
    plot(LengthRatio ~ PropNonGEI, pch=20, col=col, ylim = c(0.5,2), xlim = c(0,100), log="y",
         xlab = "% mass outside GEI-space", ylab = "Length Ratio", main = "", cex = 0.5)
    mtext(substitute(paste(beta[eg], "= log(", v1, "), n = ", v2, sep=""),
                     list(v1=c(1.0, 1.5, 3.0)[j], v2=c(800, 4000, 8000)[k])),
          side=3, line=1)
  }
}

dev.off()

################
### Figure 4 ###
################

library(Rcpp)
sourceCpp("Functions_Rcpp.cpp")
# source("Functions_Rbase.r") 

Cellprb <- fn_PhiToCellprb(c(log(0.05/0.95), log(1.2), 0, log(1.5)), 0.3, 0.3, 1)
rctr <- Cellprb$pctr/sum(Cellprb$pctr)
rcas <- Cellprb$pcas/sum(Cellprb$pcas)

set.seed(210)
dctr <- c(rmultinom(1, 400, rctr))
dcas <- c(rmultinom(1, 400, rcas))

# General
resGEN <- fn_MainGeneral(dctr, dcas, 10000)
beta3GEN <- resGEN$beta3GEN
# GEI
resGEI <- fn_MainRGEI(dctr, dcas, 10000, 2.5, 2.5, 3.5, 0)
beta3GEI <- with(resGEI, sample(Phi[,4], nrow(Phi), TRUE, Weight))
# RGEI-1 (\sigma_{\rho} = 0.01)
resRGEI_1 <- fn_MainRGEI(dctr, dcas, 10000, 2.5, 2.5, 3.5, 0.01)
beta3RGEI_1 <- with(resRGEI_1, sample(Phi[,4], nrow(Phi), TRUE, Weight))
# RGEI-2 (\sigma_{\rho} = 0.1)
resRGEI_2 <- fn_MainRGEI(dctr, dcas, 10000, 2.5, 2.5, 3.5, 0.1)
beta3RGEI_2 <- with(resRGEI_2, sample(Phi[,4], nrow(Phi), TRUE, Weight))
# RGEI-3 (\sigma_{\rho} = 4)
resRGEI_3 <- fn_MainRGEI(dctr, dcas, 10000, 2.5, 2.5, 3.5, 4)
beta3RGEI_3 <- with(resRGEI_3, sample(Phi[,4], nrow(Phi), TRUE, Weight))

# setwd(Path4Figure)
postscript("Main-Figure_FiveDensity.eps", height=4, width=7, horizontal=FALSE)
op <- par(mar=c(2.5,4,4,2)+0.1)
plot(density(beta3GEN, adjust=2), col="grey", lwd=2,  lty=1, xlab=NA,
     main=expression(paste("5 Posterior Density Estimates for ", beta[eg])))
lines(density(beta3GEI, adjust=2), col="grey", lwd=2, lty=3)
lines(density(beta3RGEI_1, adjust=2), col="black", lwd=2, lty=3)
lines(density(beta3RGEI_2, adjust=2), col="black", lwd=2, lty=2)
lines(density(beta3RGEI_3, adjust=2), col="black", lwd=2, lty=1)
legend("topright", legend=c("GEI","GEN","RGEI-1","RGEI-2", "RGEI-3"), lty=c(3,1,3,2,1),
       lwd=2, col=c("grey","grey", "black", "black", "black"), cex=0.7)
dev.off()


#########################################
###  Tables/Figures in the Supplement ###
#########################################

################
### Table 1  ###
################

# setwd(Path4RData)
load("Simlt_AllSettings.RData") 
MSE = CPB <- NULL
for(j in 1:7)
{ 
  for(k in 1:3)
  { 
    res <- ResList[[j]][[k]]    
    MSE <- with(res, rbind(MSE, c(mean((GEN[,1]-log(1.5))^2),
                                  mean((GEI[,1]-log(1.5))^2),
                                  mean((RGEI[,1]-log(1.5))^2) ) ) )
    CPB <- with(res, rbind(CPB, c(mean(GEN[,2]<log(1.5)&log(1.5)<GEN[,3]),
                                  mean(GEI[,2]<log(1.5)&log(1.5)<GEI[,3]),
                                  mean(RGEI[,2]<log(1.5)&log(1.5)<RGEI[,3]) ) ) ) 
  }
}
print(MSE) # mean square error of the point estimate
print(CPB) # Coverage Probabilities of the interval estimate


################
### Figure 1 ###
################

# setwd(Path4RData)
load("Simlt_AllSettings.RData") 
# setwd(Path4Figure)
png(filename = "Sup-Figure_CompareLength.png", width = 3150, height = 1500)
op = par(mfcol = c(3,7), mar = c(4,4,4,2)+0.1, cex=1.5)

for(j in 1:7)
{ 
  for(k in 1:3)
  { 
    res <- ResList[[j]][[k]]
    GEN <- res$GEN[,3] - res$GEN[,2]
    GEI <- res$GEI[,3] - res$GEI[,2]
    
    col <- rep("black", length(GEN))
    col[res$NumTheta!=0] <- "gray70"
    xylim <- range(GEN, GEI)
    
    plot(GEN ~ GEI, col=col, pch=20, xlim=xylim, ylim=xylim, 
         main = paste("n =", c(800, 4000, 8000)[k]) )
    abline(a=0,b=1)
  }  
}
dev.off()

################
### Figure 2 ###
################

# setwd(Path4RData)
load("Simlt_AllSettings.RData") 
# setwd(Path4Figure)
png(filename = "Sup-Figure_LengthRatio.png", width = 3150, height = 1500)
op = par(mfcol = c(3,7), mar = c(4,4,2,2)+0.1, cex=1.5)

for(j in 1:7)
{ 
  for(k in 1:3)
  { 
    res <- ResList[[j]][[k]]    
    LengthRatio <- with(res, (GEI[,3]-GEI[,2])/(GEN[,3]-GEN[,2]))
    PropNonGEI  <- 100*res$PropNonGEI    

    col <- rep("black", length(LengthRatio))
    col[res$NumTheta!=0] <- "gray70"
    
    plot(LengthRatio ~ PropNonGEI, col=col, pch=20, xlim = c(0,100), ylim = c(0.5,2), log="y",
         xlab = "% mass outside GEI-space", ylab = "Length Ratio", main = "")
  }
}

dev.off()