#################################
#load libraries, functions, and data needed for following code.
#################################
source("librariesAndFunctionsForMEE.R")

modFnFCd=read.csv("modFnFCd.csv")
#modFnFCd.wtsKnown=read.csv("modFnFCd_wtsKnown.csv")
modFnRSd=read.csv("modFnRSd.csv")
#modFnRSd2=read.csv("modFnRSd2.csv")
#modFnRSd2.wtsKnown=read.csv("modFnRSd2_wtsKnown.csv")
#modFnRSd2.wtsUnKnown=read.csv("modFnRSd2_wtsUnKnown.csv")  #Not needed for Figures.
modnbRSd=read.csv("modnbRSd.csv")
repData=read.csv("repData.csv")
#sigmaRatioData=read.csv("sigmaRatioData.csv")

#########################
#Data Generation
#########################
###THEORETICAL RESULTS for estimating effect size and Var(d).
#Data provided in "SigmaRatioData.txt".
sigmaRatioData=data.frame(tau2=rep(c(0.0,0.05,0.1,0.25,0.5,0.75,1),rep(600,7)),
                          effectSize=rep(rep((0:5)/2.5,rep(100,6)),7),
                          n=rep(3:102,42))%>%mutate(varEX=exactVarDu.r(n,n,effectSize),
                                                    varLS=approxVarDu.r(n,n,effectSize),
                                                    varNE=noEffSizeVarDu.r(n,n),
                                                    Je=exactJ.r(2*(n-1)),
                                                    dlb=effectSize-2*sqrt(tau2+Je^2*exactVarDu.r(n,n,effectSize)),
                                                    dub=effectSize+2*sqrt(tau2+Je^2*exactVarDu.r(n,n,effectSize)),
                                                    varEXlb=exactVarDu.r(n,n,dlb),
                                                    varEXub=exactVarDu.r(n,n,dub),
                                                    wlb=1/(tau2+Je^2*varEXlb),
                                                    wub=1/(tau2+Je^2*varEXub),
                                                    aN=Je^2*varNE,
                                                    bNd2=(2*(n-1)/(2*(n-1) - 2) - 1)*effectSize^2,
                                                    w=1/(tau2+Je^2*varEX),
                                                    ICC=tau2*w,
                                                    ICClb=tau2*wlb,
                                                    ICCub=tau2*wub,
                                                    invJe=(1-Je),
                                                    JeRatio=(Je/(1-Je)),
                                                    wtJratio=varNE*JeRatio,
                                                    Je2Ratio=Je^2/(1-Je),
                                                    wtJ2ratio=varNE*Je2Ratio,
                                                    effectSize=factor(effectSize))

theoryDesignData=expand.grid(tau2=c(0,0.05,0.1,0.25,0.5,0.75,1),
                             effectSize=(0:5)/2.5,
                             p=c(0,0.05,0.1,0.25,0.5,0.75,.95,1),
                             n=3:50,
                             kN=1:6)


ciTolRange=c(lb=qbinom(0.025,10000,0.95),ub=qbinom(.975,10000,0.95))

#######################################
#Negative Binomial Method Distribution of weights
######################################################
nbSim=10000
nbSimParamList=list(popnICCVec=c(0,1/10,1/3,0.5,2/3,0.9),  #within study Var relative to popnVar=1.
                    #nbSimParamList=list(tau2=c(0.0,0.05,0.1,0.25,0.5,0.75,1),
                    effectSizeVec=c(0,0.2,0.5,0.8,1.2))
nbNvec=c(4,9,16,25)  

nbParamDesign=expand.grid(popnICC=nbSimParamList$popnICCVec,effectSize=nbSimParamList$effectSizeVec)
#nbParamDesign=expand.grid(tau2C=nbSimParamList$tau2,effectSize=nbSimParamList$effectSizeVec)

nbParamDesign=nbParamDesign%>%mutate(tau2=popnICC/(1-popnICC))

tmpN1=data.frame()
gdN1=data.frame()

for(i in 1:4) {
  tmpN1=rbind(tmpN1,data.frame(mnN=nbNvec[i],N=rnbinom(nbSim,size=1.55,mu=nbNvec[i])))  #mu= 8 to 9 is considered average.  See Hillebrand & Gurevitch (2014).
  gdN1=rbind(gdN1,data.frame(mnN=nbNvec[i],N=nbNvec[i]*(1:8)))
  
}
tmpN1=tmpN1%>%mutate(simN1 = ifelse(N<3,3,N),simN2 = simN1)
gdN1=gdN1%>%mutate(simN1 = ifelse(N<3,3,N),simN2 = simN1)


nbWtsData=data.frame()
gdWtsData=data.frame()

for(i in 1:dim(nbParamDesign)[1]){
  nbWts.exDu=1/(nbParamDesign$tau2[i]+exactVarDu.r(tmpN1$simN1,tmpN1$simN2,nbParamDesign$effectSize[i]))
  nbWts.exDh=1/(nbParamDesign$tau2[i]+exactJ.r(tmpN1$simN1+tmpN1$simN2-2)^2*exactVarDu.r(tmpN1$simN1,tmpN1$simN2,nbParamDesign$effectSize[i]))
  
  nbWtsData=rbind(nbWtsData,data.frame(tau2=nbParamDesign$tau2[i],icc=nbParamDesign$popnICC[i],effectSize=nbParamDesign$effectSize[i],
                                       mnN=tmpN1$mnN,N=tmpN1$simN1,nbWts.exDu=nbWts.exDu,nbWts.exDh=nbWts.exDh))
  
  gdWts.exDu=1/(nbParamDesign$tau2[i]+exactVarDu.r(gdN1$simN1,gdN1$simN2,nbParamDesign$effectSize[i]))
  gdWts.exDh=1/(nbParamDesign$tau2[i]+exactJ.r(gdN1$simN1+gdN1$simN2-2)^2*exactVarDu.r(gdN1$simN1,gdN1$simN2,nbParamDesign$effectSize[i]))
  
  gdWtsData=rbind(gdWtsData,data.frame(tau2=nbParamDesign$tau2[i],icc=nbParamDesign$popnICC[i],effectSize=nbParamDesign$effectSize[i],
                                       mnN=gdN1$mnN,N=gdN1$simN1,gdWts.exDu=gdWts.exDu,gdWts.exDh=gdWts.exDh))
  
}


##################################################
#RMSE Summary calculations
##################################################
relRmseData=full_join(modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse),
                      modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>%
                        filter(method=="dh.eq"),
                      by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>%
  mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2),
         relRMSE=rmse/rmse.dh.eq,
         ICC0=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10")))

rmseSumData= relRmseData%>%group_by(method,ICC0,N)%>%
  summarize(RMSE.mn=sqrt(mean(relRMSE^2)),
            RMSE.min=min(relRMSE),
            RMSE.max=max(relRMSE),.groups = 'drop')

relRmseFeData=full_join(modFnFCd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse),
                        modFnFCd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>%
                          filter(method=="dh.eq"),
                        by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>%
  mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2),
         relRMSE=rmse/rmse.dh.eq)

#NegBinom Data
relRmseNbData=full_join(modnbRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse),
                        modnbRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>%
                          filter(method=="dh.eq"),
                        by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>%
  mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2),
         relRMSE=rmse/rmse.dh.eq,
         ICC0=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10")))

rmseSumNbData= relRmseNbData%>%group_by(method,ICC0,N)%>%
  summarize(RMSE.mn=sqrt(mean(relRMSE^2)),
            RMSE.min=min(relRMSE),
            RMSE.max=max(relRMSE),.groups = 'drop')

######################################
#% Bias Summary calculations
######################################
maxPctBias=unlist(modFnRSd%>%filter(effectSize>0)%>%summarize(max(abs(bias/effectSize))))
maxPctBias.duLS=unlist(modFnRSd%>%filter(effectSize>0&method=="du.lsDu")%>%summarize(max(abs(bias/effectSize))))

pctNoteBias=prop.table(with(modFnRSd%>%filter(effectSize>0),table(method,abs(bias/effectSize)> 0.05)),1)
pctNoteBiasByN=prop.table(with(modFnRSd%>%filter(effectSize>0),table(method,N,abs(bias/effectSize)> 0.05)),c(1,3))[,,2]

sumBiasByN=modFnRSd%>%
  filter(effectSize>0)%>%group_by(method,N)%>%
  summarise(mnBias=mean(bias/effectSize),
            medBias=median(bias/effectSize),
            sdBias=sqrt(var(bias/effectSize)),
            maxBias=max(abs(bias)/effectSize),.groups = 'drop')

###By M 
biasVsVarSumData=full_join(modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse),
                           modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>%
                             filter(method=="dh.eq")%>%
                             select(N,kN,Sn,M,popnICC,effectSize,rmse),
                           by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dHeq"))%>%
  mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2),
         relRMSE=rmse/rmse.dHeq)%>%
  group_by(method,M)%>%filter(effectSize>0)%>%
  summarize(RMSE.mn=sqrt(mean(relRMSE^2)),
            bias.mn=100*mean(abs(bias/effectSize)),
            bN.mn=sqrt(mean(bN2)),
            RMSE.max=sqrt(max(relRMSE^2)),
            #RMSE.max=ifelse(abs(min(log(relRMSE^2)))>abs(max(log(relRMSE^2))),sqrt(min(relRMSE^2)),sqrt(max(relRMSE^2))),
            bias.max=100*max(abs(bias/effectSize)),
            bN.max=sqrt(max(bN2)),.groups = 'drop')

#################################
#VARIOUS SUMMARIES
#################################
#Mean Bias
mnBiasResult=sumBiasByN%>%dplyr::filter(method=="du.lsDu"&N==25)%>%select(mnBias) #Not in data list for MEE

#Notable differences between observed q and 100(1-alpha)
infPctNoteBias=prop.table(with(modFnRSd,table(method,abs(pctSig-0.95)> 0.025)),1)
infPctNoteBiasByN=prop.table(with(modFnRSd,table(method,N,abs(pctSig-0.95)> 0.025)),c(1,3))[,,2]

#Notable improvement in relative RMSE
noteRelRmse=prop.table(with(relRmseData,table(method,abs(relRMSE-1)> 0.1)),1)
noteRelRmseByN=prop.table(with(relRmseData,table(method,kN,abs(relRMSE-1)> 0.1)),c(1,3))[,,2]

noteRelRmse5=prop.table(with(relRmseData,table(method,abs(relRMSE-1)> 0.05)),1)

#Neg Binom results
noteNbRelNbRmse=prop.table(with(relRmseNbData,table(method,abs(relRMSE-1)> 0.1)),1)
noteNbRelNbRmseByN=prop.table(with(relRmseNbData,table(method,kN,abs(relRMSE-1)> 0.1)),c(1,3))[,,2]

noteRelNbRmse5=prop.table(with(relRmseNbData,table(method,abs(relRMSE-1)> 0.05)),1)

#Mean bias of effect size for random (> 0) and fixed (==0) effects
mnBiasTbl=with(modFnRSd%>%filter(effectSize>0),tapply(bias/effectSize,list(method,N),mean))
mnBias0Tbl=with(modFnRSd%>%filter(effectSize==0),tapply(bias,list(method,N),mean)) #Not in data list for MEE

#Notable bias
mnInfBiasTbl=with(modFnRSd%>%filter(effectSize>0),tapply(abs(pctSig-0.95)> 0.025,list(method,N,M),mean)) #Not in data list for MEE

#Mean RMSE
mnRMSEtbl=with(modFnRSd,sqrt(tapply(rmse^2,list(method,N),mean))) #Not in data list for MEE

############################################
#Calculate CI for variance estimate based on uncertainty in delta.  
############################################
varCI.exData=theoryDesignData%>%mutate(n1=n,n2=n)%>%select(n1,n2,effectSize)%>%distinct()%>%rowwise%>%
  mutate(varLB=varCI.r(n1,n2,effectSize,method="EX")[1],
         varMed=varCI.r(n1,n2,effectSize,method="EX")[2],
         varUB=varCI.r(n1,n2,effectSize,method="EX")[3])

varCI.lsData=theoryDesignData%>%mutate(n1=n,n2=n)%>%select(n1,n2,effectSize)%>%distinct()%>%rowwise%>%
  mutate(varLB=varCI.r(n1,n2,effectSize,method="LS")[1],
         varMed=varCI.r(n1,n2,effectSize,method="LS")[2],
         varUB=varCI.r(n1,n2,effectSize,method="LS")[3])

######################################
#Replicate data for Bland-Altman analysis.
#This data is provided in repData.csv.  Code shows how it was created.
#repData=full_join(gather(modFnRSd2%>%select(popnICC,M,N,kN,Sn,method,effectSize,bias,pctSig,simSD),
#                         "metric","value",bias,pctSig,simSD),
#                  gather(modFnRSd3%>%select(popnICC,M,N,kN,Sn,method,effectSize,bias,pctSig,simSD),
#                         "metric","value",bias,pctSig,simSD),
#                  by=c("popnICC","M","N","kN","Sn","effectSize","method","metric"),suffix=c("",".rep"))
#
repSD=repData%>%group_by(metric)%>%summarize(sqrt(var(value-value.rep)/2),.groups = 'drop')