#Simulations
# clears workspace:  
rm(list=ls(all=TRUE)) 

# Set work directory  (replace with own path)

#dir <- "C:/Users/Josine/Documents/My Dropbox/UVA/JosineEJ/Replicationprior/Code/"
dir <- "c:/Dropbox/UVA/JosineEJ/Replicationprior/Code/"
setwd(dir)


# Load and if necessary install the required libraries
library(R2WinBUGS)  # Only necessary for posterior plots for JZS prior
library(MCMCpack)

# Load the functions necessary to compute the replication prior (file saved in work directory)

source(paste(dir,'Repfunctions.R', sep = "")) 


#Prepare simulations: 

torig <- c(1,2,3,4)
trep5050 <- c(-1,0,1,2,3,4)
trep10050 <- c(-.1,0,.1,.2,.3,.4) * sqrt(50)
trep50100 <- c(-.14,0,.14,.28,.42,.57) * sqrt(100)
trep100100 <- c(-1,0,1,2,3,4)

drep5050 <- c(-.14,0,.14,.28,.42,.57)
drep10050 <- c(-.1,0,.1,.2,.3,.4) 
drep50100 <- c(-.14,0,.14,.28,.42,.57) 
drep100100 <- c(-1,0,1,2,3,4)

drep <- c(drep5050,drep10050,drep50100,drep100100)


trep <- c(trep5050,trep10050,trep50100,trep100100)
norig  <- c(rep(50,6),rep(100,6),rep(50,6),rep(100,6))
nrep  <-  c(rep(50,6),rep(50,6),rep(100,6),rep(100,6))

##################################################
##New replication Bayes factor 
#################################################


SIM <- list()

k<-0
for (or in 1:4) { 
for (rep in 1:24){
k <- k+1
SIM[[k]] <- ReplicationBayesfactorNCT(torig[or], trep[rep], norig[rep], nrep[rep], plot=0 , post = 0)
}}

tableSIMnew <- matrix(0,8,12)

for ( k in 1:12){
tableSIMnew[1,k] <- SIM[[k]]$BF
tableSIMnew[5,k] <- SIM[[k+12]]$BF 
tableSIMnew[2,k] <- SIM[[k+24]]$BF
tableSIMnew[6,k] <- SIM[[k+36]]$BF
tableSIMnew[3,k] <- SIM[[k+48]]$BF
tableSIMnew[7,k] <- SIM[[k+60]]$BF
tableSIMnew[4,k] <- SIM[[k+72]]$BF
tableSIMnew[8,k] <- SIM[[k+84]]$BF
}

tableSIMnew2 <- tableSIMnew
tableSIMnew2[1:4,7:12] <-  tableSIMnew[5:8,1:6]
tableSIMnew2[5:8,1:6] <-  tableSIMnew[1:4,7:12]

round(tableSIMnew2,2)

##################################################
##Meta-analysis Bayes factor
#################################################


# and now for the other tests 
SIMmeta <- list()

k<-0
for (or in 1:4) { 
  for (rep in 1:24){
    k <- k+1
    SIMmeta[[k]] <- BFmeta(torig[or], trep[rep],1, norig[rep], nrep[rep])
  }}


tableSIMmeta <- matrix(0,8,12)

for ( k in 1:12){
  tableSIMmeta[1,k] <- SIMmeta[[k]]
  tableSIMmeta[5,k] <- SIMmeta[[k+12]]
  tableSIMmeta[2,k] <- SIMmeta[[k+24]]
  tableSIMmeta[6,k] <- SIMmeta[[k+36]]
  tableSIMmeta[3,k] <- SIMmeta[[k+48]]
  tableSIMmeta[7,k] <- SIMmeta[[k+60]]
  tableSIMmeta[4,k] <- SIMmeta[[k+72]]
  tableSIMmeta[8,k] <- SIMmeta[[k+84]]
}

tableSIMmeta2 <- tableSIMmeta
tableSIMmeta2[1:4,7:12] <-  tableSIMmeta[5:8,1:6]
tableSIMmeta2[5:8,1:6] <-  tableSIMmeta[1:4,7:12]

round(tableSIMmeta2,2)

# to print to latex: 
#library(xtable)
#print(xtable(round(tableSIMmeta2,2) ), only.contents=TRUE, ,include.colnames = F,
#      hline.after = NULL, file = "tableSIMmeta2.tex") 
#print(xtable(round(tableSIMmeta2E,2),display =c('s','E','E','E','E','E', 'E','E','E','E','E','E')  ), only.contents=TRUE, ,include.colnames = F,
#      hline.after = NULL, file = "tableSIMmeta2.tex") 


##################################################
##Equality of effect sizes Bayes factor
#################################################


SIMMB <- list()

k<-0
for (or in 1:4) { 
  for (rep in 1:24){
    k <- k+1
    SIMMB[[k]] <- BayMayBF(torig[or], trep[rep], norig[rep], nrep[rep])
  }}


tableSIMMB <- matrix(0,8,12)

for ( k in 1:12){
  tableSIMMB[1,k] <- SIMMB[[k]]
  tableSIMMB[5,k] <- SIMMB[[k+12]] 
  tableSIMMB[2,k] <- SIMMB[[k+24]]
  tableSIMMB[6,k] <- SIMMB[[k+36]]
  tableSIMMB[3,k] <- SIMMB[[k+48]]
  tableSIMMB[7,k] <- SIMMB[[k+60]]
  tableSIMMB[4,k] <- SIMMB[[k+72]]
  tableSIMMB[8,k] <- SIMMB[[k+84]]
}

tableSIMMB2 <- tableSIMMB
tableSIMMB2[1:4,7:12] <-  tableSIMMB[5:8,1:6]
tableSIMMB2[5:8,1:6] <-  tableSIMMB[1:4,7:12]

round(tableSIMMB2,2)

#library(xtable)
#print(xtable(round(tableSIMMB2,2)), only.contents=TRUE, ,include.colnames = F,
#hline.after = NULL, file = "tableSIMMB2.tex") 
#print(xtable(round(tableSIMMB2[1:4,1:6],2)), only.contents=TRUE, ,include.colnames = F,
#hline.after = NULL, file = "tableSIMMB25050.tex") 
#print(xtable(round(tableSIMMB2[5:8,1:6],2)), only.contents=TRUE, ,include.colnames = F,
#hline.after = NULL, file = "tableSIMMB210050.tex") 
#print(xtable(round(tableSIMMB2[1:4,7:12],2)), only.contents=TRUE, ,include.colnames = F,
#      hline.after = NULL, file = "tableSIMMB250100.tex") 
#print(xtable(round(tableSIMMB2[5:8,7:12],2)), only.contents=TRUE, ,include.colnames = F,
#      hline.after = NULL, file = "tableSIMMB2100100.tex") 

##################################################
##For the default JZS bayes factor (two-sided)
#################################################

torig2 <- rep(torig,6)
norig2 <- rep(norig,6)

SIMdef <- list()

k<-0
  for (rep in 1:24){
   k <- k+1
    SIMdef[[k]] <- JZSBF(torig2[rep], trep[rep], norig2[rep], nrep[rep])
  }

tableSIMdef <- matrix(0,1,24)

for ( k in 1:24){
  tableSIMdef[k] <- SIMdef[[k]][2]
}

simdf <- round(t(cbind(trep[1:12],drep[1:12],matrix(round(tableSIMdef,2),,2))),2)
simdf

#print(xtable(simdf), only.contents=TRUE, ,include.colnames = F,
#      hline.after = NULL, file = "tableSIMdef.tex")