################################################################
#This code generates simulation parameter sets and runs the different parameter settings in blocks 
#as time to run is very long (4 - 10 days).

source("librariesAndFunctionsForMEE.R")

################################################################
#GLOBAL VARIABLES
################################################
#List of indicators for the different estimation methods.  See Tables 1 - 3 in paper for details.
methodList=c("du.eq","dh.eq","dh.n","dh.neDh",
             "dh.xmDh", "dh.lmDh","du.xmDu","du.lmDu",
             "dh.exDh","dh.lsDh","dh.wxDh","du.exDu","du.lsDu","du.wxDu")

methodListN=length(methodList)

theFinalList=methodList

nSim=10000 

#############################################################################
#Simulation Parameters and Variables
feSimParamList=list(totalVarVec=1,
                    popnICCVec=0,  #Population ICC = proportion of PopnVar to Total = Vp/(Vp + Ve).
                    effectSizeVec=c(0,0.2,0.5,0.8,1.2,2), 
                    Svec=c(2:4,6,8),
                    Mvec=c(5,10,15,25,50,100),  #number of endpoints
                    Nvec=c(4,9,16,25),  #sample size for endpoint
                    kNvec= 1:5)  #creates two tiers of sample size across endpoints in meta-analysis

reSimParamList=list(totalVarVec=1,
                    popnICCVec=c(1/10,1/3,1/2,2/3,9/10),  #within study Var relative to popnVar=1.
                    effectSizeVec=c(0,0.2,0.5,0.8,1.2), 
                    Svec=2,
                    Mvec=c(5,10,15,25,50),  #number of endpoints
                    Nvec=c(4,9,16,25),  #sample size for endpoint
                    kNvec= 2:4)  #creates two tiers of sample size across endpoints in meta-analysis

############################################################
#DATA DEVELOPMENT GRP1
############################################################
reSimParamList1=reSimParamList
reSimParamList1$popnICCVec=reSimParamList$popnICCVec[1]  

reSimParamList2=reSimParamList
reSimParamList2$popnICCVec=reSimParamList$popnICCVec[2]  

reSimParamList2.wtsKnwn=reSimParamList
reSimParamList2.wtsKnwn$popnICCVec=reSimParamList$popnICCVec[2]  
reSimParamList2.wtsKnwn$Mvec=reSimParamList$Mvec[3]  
reSimParamList2.wtsKnwn$effectSizeVec=reSimParamList$effectSizeVec[4]  

reSimParamList8=reSimParamList
reSimParamList8$effectSizeVec=c(reSimParamList$effectSizeVec,2)
reSimParamList8$Mvec=100  

reSimParamList8b=reSimParamList
reSimParamList8b$popnICCVec=reSimParamList$popnICCVec[c(3,5)]  
reSimParamList8b$effectSizeVec=c(reSimParamList$effectSizeVec,2)
reSimParamList8b$Mvec=100  

##################################################
#SIMULATIONS FOR FIXED AND RANDOM EFFECTS
#These calculations are time consuming, approx. 4 - 10 days for each on laptop.
##################################
#simGenMA_Methods.r is shell for running simulation in blocks based on subset descrbed in *ParamList
Sys.time()
modFnFCd.wtsKnwn=simGenMA_Methods.r(nSim,feSimParamList,genNtype="Fn",genDtype="FC",wtsKnown=TRUE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnFCd.wtsKnwn,file="modFnFCdwtsKnwn.Rdata")

Sys.time()
modFnRSd2.wtsKnwn=simGenMA_Methods.r(nSim,reSimParamList2.wtsKnwn,genNtype="Fn",genDtype="RS",wtsKnown=TRUE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd2.wtsKnwn,file="modFnRSd2wtsKnwn.Rdata")

Sys.time()
modFnFCd=simGenMA_Methods.r(nSim,feSimParamList,genNtype="Fn",genDtype="FC",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnFCd,file="modFnFCd.Rdata")

Sys.time()
modRnFCd=simGenMA_Methods.r(nSim,feSimParamList,genNtype="Rn",genDtype="FC",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modRnFCd,file="modRnFCd.Rdata")

Sys.time()
modFnRSd1=simGenMA_Methods.r(nSim,reSimParamList1,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd1,file="modFnRSd1.Rdata")

Sys.time()
modFnRSd2=simGenMA_Methods.r(nSim,reSimParamList2,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd2,file="modFnRSd2.Rdata")

Sys.time()
modFnRSd8=simGenMA_Methods.r(nSim,reSimParamList8,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd8,file="modFnRSd8.Rdata")

#Missed or incorrect popnICC
Sys.time()
modFnRSd8b=simGenMA_Methods.r(nSim,reSimParamList8b,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd8b,file="modFnRSd8b.Rdata")

############################################################
#DATA DEVELOPMENT GRP2
############################################################

reSimParamList3=reSimParamList
reSimParamList3$popnICCVec=reSimParamList$popnICCVec[3]  

reSimParamList4=reSimParamList
reSimParamList4$popnICCVec=reSimParamList$popnICCVec[4]  

reSimParamList5=reSimParamList
reSimParamList5$popnICCVec=reSimParamList$popnICCVec[5]  

reSimParamList0=reSimParamList
reSimParamList0$popnICCVec=reSimParamList$popnICCVec[1]  

reSimParamList8a=reSimParamList
reSimParamList8a$popnICCVec=reSimParamList$popnICCVec[1]  
reSimParamList8a$effectSizeVec=c(reSimParamList$effectSizeVec,2)
reSimParamList8a$Mvec=100 


##################################################
#SIMULATIONS FOR RANDOM EFFECTS (PopnVar > 0)
#These calculations are time consuming, approx. 3 days on laptop.
##################################
#TESt rUN USING Lin PARAMETER SETTINGS
Sys.time()
modFnRSd3=simGenMA_Methods.r(nSim,reSimParamList3,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd3,file="modFnRSd3.Rdata")

Sys.time()
modFnRSd4=simGenMA_Methods.r(nSim,reSimParamList4,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd4,file="modFnRSd4.Rdata")

Sys.time()
modFnRSd5=simGenMA_Methods.r(nSim,reSimParamList5,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd5,file="modFnRSd5.Rdata")

Sys.time()
modFnRSd0=simGenMA_Methods.r(nSim,reSimParamList0,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd0,file="modFnRSd0.Rdata")

Sys.time()
modFnRSd8a=simGenMA_Methods.r(nSim,reSimParamList8a,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd8a,file="modFnRSd8a.Rdata")

############################################################
#DATA DEVELOPMENT GRP3
############################################################

reSimParamList6=reSimParamList
reSimParamList6$popnICCVec=reSimParamList$popnICCVec[6]  

reSimParamList7=reSimParamList
reSimParamList7$effectSizeVec=2  

##################################################
#SIMULATIONS FOR RANDOM EFFECTS (PopnVar > 0)
#These calculations are time consuming, approx. 4 - 10 days on laptop.
##################################

Sys.time()
modFnRSd6=simGenMA_Methods.r(nSim,reSimParamList6,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd6,file="modFnRSd6.Rdata")

Sys.time()
modFnRSd7=simGenMA_Methods.r(nSim,reSimParamList7,genNtype="Fn",genDtype="RS",wtsKnown=FALSE)%>%
  mutate(method=factor(method,levels=theFinalList))
Sys.time()
save(modFnRSd7,file="modFnRSd7.Rdata")


############################################################
#These load in simulation results.  
#*FCd is for fixed effects results and can be run relatively quickly.
#*RSd# is for random effects and these took 4 - 10 days each on a laptop.
load("finalData/modFnFCd.Rdata")
load("finalData/modFnRSd0.Rdata")  #This is popnICC = 0
load("finalData/modFnRSd1.Rdata")  #This is popnICC = 0.1
load("finalData/modFnRSd2.Rdata")  #This is a rep of modFnRSd3
load("finalData/modFnRSd3.Rdata")   #This is popnICC = 1/3
load("finalData/modFnRSd4.Rdata")   #This is popnICC = 0.5
load("finalData/modFnRSd5.Rdata")  #This is popnICC = 2/3
load("finalData/modFnRSd6.Rdata")  #This is popnICC = 0.9
load("finalData/modFnRSd7.Rdata")  #This adds effect size = 2
load("finalData/modFnRSd8.Rdata")  #This adds M = 100
load("finalData/modFnRSd8a.Rdata")  #This adds M = 100 in part (missed this in 8)
load("finalData/modFnRSd8b.Rdata")  #This adds M = 100 in part (missed this in 8)
load("finalData/modnbRSd.Rdata")  #This is simulation data when negative binomial distribution is used to derive sample size n.

simResults=vector("list",2)
names(simResults)=c("modFnFCd","modFnRSd")

simResults$modFnFCd=modFnFCd
simResults$modFnRSd=rbind(modFnRSd0,modFnRSd1,modFnRSd3,modFnRSd4,modFnRSd5,modFnRSd6,modFnRSd7,
                          modFnRSd8a,modFnRSd8,modFnRSd8b)%>%filter(popnICC<=0.9)

#Only modFnFCd, modFnRSd, modFnRSd2, modnbRSd, modFnFC.wtsKnown, modFnRSd2.wtsKnown, repData are given in data archive.