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

modFnFCd=read.csv("modFnFCd.csv")
modFnFCd.wtsKnown=read.csv("modFnFCd_wtsKnown.csv")
modFnRSd=read.csv("modFnRSd.csv")
modFnRSd2.wtsKnown=read.csv("modFnRSd2_wtsKnown.csv")
modFnRSd2.wtsUnKnown=read.csv("modFnRSd2_wtsUnKnown.csv")  #Not needed for Figures.
modnbRSd=read.csv("modnbRSd.csv")  #Not needed for Figures.?
repData=read.csv("repData.csv")
sigmaRatioData=read.csv("sigmaRatioData.csv")  #This data is also provided in source code below.

source("TheoryAndSummaryDataForMEE.R")

########################################
#FUNCTIONS
########################################
addGGplot=function(ppReport,currentFig,title,height=6,width=9,left=0.5,top=1){
  ppReport= add_slide(ppReport,layout = "Title and Content", master = "Office Theme")
  
  ppReport= ph_with(x = ppReport, 
                    value = currentFig,
                    location=ph_location(height = height , width = width, left =left, top = top))
  
  ppReport= ph_with(x =ppReport, value = title, location = ph_location_type(type = "title"))
  
  ppReport
}

################################################
#GLOBAL VARIABLES
################################################
methodList=c("du.eq","du.xmDu","du.lmDu",
             "dh.eq","dh.n","dh.neDh","dh.xmDh", "dh.lmDh",
             "dh.exDh","dh.wxDh","dh.lsDh",
             "du.exDu","du.wxDu",
             "du.lsDu")


methodListN=length(methodList)

theFinalList=methodList

methodCalcList=c("du.eq","dh.eq","dh.n","dh.neDh",
                 "du.xmDu", "dh.xmDh", "du.lmDu", "dh.lmDh",
                 "du.exDu", "dh.exDh", "du.lsDu", "dh.lsDh","du.wxDu","dh.wxDh")

methodLabelList=c(expression(paste(d[U],".w(EQ)",sep="")),
                  expression(paste(d[U],".w(EX, ",bar(d)[U],")",sep="")),
                  expression(paste(d[U],".w(LS, ",bar(d)[U],")",sep="")),
                  expression(paste(d[H],".w(EQ)",sep="")),
                  expression(paste(d[H],".w(N)",sep="")),
                  expression(paste(d[H],".w(NE)",sep="")),
                  expression(paste(d[H],".w(EX, ",bar(d)[H],")",sep="")),
                  expression(paste(d[H],".w(LS, ",bar(d)[H],")",sep="")),
                  expression(paste(d[H],".w(EX, ",d[H],")",sep="")),
                  expression(paste(d[H],".w(WX, ",d[H],")",sep="")),
                  expression(paste(d[H],".w(LS, ",d[H],")",sep="")),
                  expression(paste(d[U],".w(EX, ",d[U],")",sep="")),
                  expression(paste(d[U],".w(WX, ",d[U],")",sep="")),
                  expression(paste(d[U],".w(LS, ",d[U],")",sep="")))

I2LabList=c(expression(paste(I[1]^2," = 0")),
           expression(paste(I[1]^2," = 1/10")),
           expression(paste(I[1]^2," = 1/3")),
           expression(paste(I[1]^2," = 1/2")),
           expression(paste(I[1]^2," = 2/3")),
           expression(paste(I[1]^2," = 9/10")))

NlabList=c(expression(paste(N[S]," = 4")),
           expression(paste(N[S]," = 9")),
           expression(paste(N[S]," = 16")),
           expression(paste(N[S]," = 25")))

NIv2LabList=c(expression(paste(n," = 3")),
              expression(paste(n," = 4")),
              expression(paste(n," = 6")),
              expression(paste(n," = 9")),
              expression(paste(n," = 16")),
              expression(paste(n," = 25")))

NaLabList=c(expression(paste(N[a]," = 4")),
           expression(paste(N[a]," = 9")),
           expression(paste(N[a]," = 16")),
           expression(paste(N[a]," = 25")))

NstarLabList=c(expression(paste(N["*"]," = 4")),
            expression(paste(N["*"]," = 9")),
            expression(paste(N["*"]," = 16")),
            expression(paste(N["*"]," = 25")))

SnLabList=c(expression(paste(S[N]," = 2")),
            expression(paste(S[N]," = 3")),
            expression(paste(S[N]," = 4")),
            expression(paste(S[N]," = 6")),
            expression(paste(S[N]," = 8")))

T2FracLabList=c("0","1/11","1/2","1","2","9")

T2simLabList=c(expression(paste(tau^2," = 0")),
            expression(paste(tau^2," = 1/11")),
            expression(paste(tau^2," = 1/2")),
            expression(paste(tau^2," = 1")),
            expression(paste(tau^2," = 2")),
            expression(paste(tau^2," = 9")))

T2labList=c(expression(paste(tau^2," = 0.05")),
            expression(paste(tau^2," = 0.10")),
            expression(paste(tau^2," = 0.25")),
            expression(paste(tau^2," = 0.50")),
            expression(paste(tau^2," = 0.75")),
            expression(paste(tau^2," = 1.00")))

tauLabList=c(expression(paste(tau^2," = 0.0")),
             expression(paste(tau^2," = 0.5")),
             expression(paste(tau^2," = 1.0")),
             expression(paste(tau^2," = 2.0")))

deltaLabList=c(expression(paste(delta," = 0.0")),
               expression(paste(delta," = 0.5")),
               expression(paste(delta," = 1.0")),
               expression(paste(delta," = 2.0")))

effSizeLabList=c(expression(paste(delta," = 0.0")),
                 expression(paste(delta," = 0.4")),
                 expression(paste(delta," = 0.8")),
                 expression(paste(delta," = 1.2")),
                 expression(paste(delta," = 1.6")),
                 expression(paste(delta," = 2.0")))

biasMethodList=c("dh.eq","dh.exDh","dh.lsDh","dh.wxDh","du.eq","du.exDu","du.lsDu","du.wxDu")
biasMethodLabelList=methodLabelList[c(4,9,11,10,1,12,14,13)]

nColorList=c("Red","Orange","Blue","Black")
mColorList=c("Black","Blue","DarkGreen","brown","Orange","Red")
kNcolorList=c("Black","Blue","DarkGreen","Tan","Orange","Red")
kNcolorShortList=c("Black","Orange","Red")
deltaColorList=c("Black","Blue","Green","brown","Red","Orange")
effectSizeColorList=c("Black","Blue","Green","Red","Orange")
popnIccColorList=c("Brown","Orange","Red","Blue","DarkGreen","Black")
methodColorList=c(rep("red",3),rep("Black",5),rep("purple",3),rep("darkOrange",2),"green")
biasedMethodColorList=c("Red","Black","darkorange4","blue","purple","green")

pColorList=c("Black","DarkGreen","Purple","Orange","Red","darkorange","Green","Blue","Brown","DarkGreen","Pink")

figEffectSize=2
facetFontSize=12
labFontSize=12

####################################
#FIGURES FOR MAIN TEXT
####################################

ppFileName=paste(Sys.Date(),"Figures for Hedges d Estimator MEE Paper.pptx",sep="")
ppReport=read_pptx("testFormat.pptx")

#Figure 1 is a graphic and produced elsewhere.
#Figure 2.  Theoretical weights adjusted to assess the range of potential weight values for different $\tau^2$.
nFig=2
currentFig=full_join(sigmaRatioData,sigmaRatioData%>%filter(n==3&effectSize==0),by=c("tau2"))%>%
  mutate(effectSize=factor(effectSize.x),
         ICC=ICC.x/ICC.y,
         ICC0=round(tau2/(tau2+1),2))%>%filter(tau2!=0&n.x<=50)%>%
  mutate(tau2lab=factor(tau2,labels=T2labList))%>%
  ggplot()+
  scale_y_continuous(trans = 'log2') +
  geom_line(aes(x=n.x,y=ICC,color=effectSize),size=0.75)+#xlim(0,60)+
  scale_color_manual(values=c("Black","Red","Blue","darkorange","DarkGreen","darkgoldenrod"))+
  facet_wrap(~tau2lab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste("Sample size (",n[T]," = ", n[C],")")),
       y=expression(omega[i]^0),
       color=expression(delta))

figLegText="Theoretical Weights: Rescaled"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 3.  Estimates of $\delta_P$ from 10,000 simulations for the listed estimators.  
#Each point represents the estimated result for a given scenario, 
#so points correspond to one of the parameter combinations of $\tau^2$, $M$, $k_N$, and $S_N$ with the $N_S$ and $\delta_P$ denoted in the plot.  Observations are jittered horizontally across a small range of $x$-axis to better observe results.

currentFig=simResults$modFnRSd%>%select(popnICC,effectSize,method,N,mn)%>%filter(popnICC<0.6)%>%
  mutate(ICC0=factor(popnICC,labels=c("0","1/10","1/3","1/2")))%>%
  ggplot(aes(x=N,y=mn,color=factor(method,levels=methodList,labels=methodLabelList)))+
  geom_jitter(width=0.5)+
  geom_hline(yintercept=c(0,0.2, 0.5, 0.8,1.2,2),color="darkorange")+
  scale_color_manual(values=methodColorList,labels = parse_format())+
  theme(legend.text = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y=expression(paste("Effect Size ", delta[P])),
       color="Estimator")

nFig=nFig+1
figLegText="E(d) for Simulated meta-analyses"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 4. Histograms of the percent bias found in the 2592 different simulation parameter scenarios for six representative estimators. 
#Estimators using individual d_(*i) in weight functions are biased, particularly when N_S=4 and to some extent when N_S=9.
#The results are restricted to case where $k_N = 4$. Results for $N = 25$ are similar to those seen for $N=16$.

currentFig=simResults$modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList))%>%
  filter(effectSize>0&(method=="dh.eq"|method=="dh.lsDh"|method=="dh.exDh"|method=="du.eq"|method=="du.lsDu"|method=="du.exDu"))%>%
  filter((100*bias/effectSize)>-20&(100*bias/effectSize)<15)%>%
  ggplot(aes(x=100*bias/effectSize,fill=factor(N)))+
  facet_wrap(~as.character(methodLab),nrow=2, labeller = label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_fill_manual(values=nColorList)+
  geom_histogram(bins=40)+xlim(-20,15)+
  labs(x="% Bias",y="Count",
       fill=expression(N[S]))

nFig=nFig+1
figLegText="Frequency of % Bias values"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 5.  Percent bias as a function of I_v^2, δ_P, and N_S. Observations are jittered horizontally across a small range of the x-axis to better observe results.  
#The results are restricted to the case where k_N=4. 
#Results for N_S=25 are similar to results for N_S=16.

###########
#Calculation of Iv^2
testGridData=expand.grid(totalVarVec=1,
                         popnICC=c(0,1/10,1/3,1/2,2/3,9/10),  #within study Var relative to popnVar=1.
                         effectSize=c(0,0.2,0.5,0.8,1.2,2), 
                         S=2,
                         M=c(5,10,15,25,50,100),  #number of endpoints
                         N=c(4,9,16,25),  #sample size for endpoint
                         kN= 2:4)%>%  #creates two tiers of sample size across endpoints in meta-analysis
  mutate(tau2=popnICC/(1-popnICC))%>%
  select(popnICC,effectSize,M,N,kN,tau2)  

gridS2=NULL

for(i in 1:dim(testGridData)[1]){
  gridResult=gridTypicalS2.r(testGridData$M[i],testGridData$N[i],testGridData$kN[i],testGridData$effectSize[i],testGridData$tau2[i])
  gridS2=rbind(gridS2,data.frame(typS2=gridResult[1],avgS2=gridResult[2]))
}

testGridData=data.frame(testGridData,gridS2)

#Plot of ICCv against bias
currentFig=right_join(simResults$modFnRSd%>%select(popnICC,effectSize,method,N,kN,M,bias),
                      testGridData,
                      by=c("popnICC","effectSize","M","N","kN"))%>%#filter(popnICC<0.6)%>%
  mutate(ICCv=tau2/(tau2+avgS2))%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NlabList),
         tau2=factor(tau2,labels=c("0","1/11","1/2","1","2","9")))%>%
  filter(effectSize>0&N < 20&(method=="dh.eq"|method=="dh.lsDh"|method=="du.lsDu"))%>%
  ggplot(aes(x=ICCv,y=100*bias/effectSize))+
  geom_point(aes(color=factor(N)))+
  geom_hline(yintercept=0)+
  facet_wrap( ~ as.character(methodLab),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=nColorList)+
  labs(x=expression(paste(I[v]^2,"(",tau^2,")",sep="")),
       y="Percent Bias",
       color=expression(N[S]))

nFig=nFig+1
figLegText="ICCv vs % Bias"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 6.  RMSE for a selected set of estimators relative to RMSE for d_H.w(EQ) as a function of I_v^2 and k_N. 
#Differences between N_S=16 and N_S=25 were negligible.
currentFig=full_join(relRmseData,testGridData,
                     by=c("popnICC","effectSize","M","N","kN"))%>%
  mutate(ICCv=tau2/(tau2+avgS2))%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NlabList))%>%
  filter((method=="dh.n"|method=="dh.neDh"|method=="du.lsDu"))%>%
  ggplot(aes(x=ICCv,y=relRMSE,color=factor(kN)))+
  facet_wrap(~methodLab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=kNcolorShortList)+
  geom_point()+
  labs(x=expression(paste(I[v]^2,"(",tau^2,")",sep="")),
       y="Relative RMSE",
       color=expression(k[N]))

nFig=nFig+1
figLegText="RMSE for Methods relative to RMSE of dh.w(EQ) Across ICCv"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 7.  RMSE of  d_U.w(LS,d_U) relative to RMSE of d_H.w(EQ) for a fixed effects model (i.e.  where τ^2is known to equal 0) as a function of S_N and k_N. Observations are jittered horizontally across a small range of the x-axis to better observe results.  Larger values than S_N=2 are expected to be observed in meta-analyses. Results for N_S>9 were similar to results for N_S=9 and are not shown. 
#The estimator d_U.w(LS,d_U) was selected for comparison as it had a smaller RMSE than any of the unbiased estimators.

currentFig=relRmseFeData%>%mutate(Nlab=factor(N,labels=NlabList),
                       SnLab=factor(Sn,labels=SnLabList))%>%
  filter(method=="du.lsDu"&Sn<8&N<10)%>%
  ggplot(aes(x=kN,y=relRMSE,color=factor(M)))+
  geom_jitter(width=0.1)+
  geom_hline(yintercept = 0.9,color="black")+
  facet_grid(SnLab~Nlab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=mColorList)+
  labs(x=expression(paste("Relative range of sample size ",k[N])),
       y="Relative RMSE",
       color="M")

nFig=nFig+1
figLegText="kN vs relRMSE"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#Figure 8.  Estimation of tau2.
#Bias in tau2 does not impact estimation of deltaP
mnPctBiasTau2=simResults$modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  group_by(N,method)%>%
  summarize(tauMN=100*mean(tau2.mn/trueTau2-1))%>%filter(N>9)

mnPctBiasTau2=simResults$modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  group_by(method)%>%
  summarize(tauPctBiasMn=100*mean(tau2.mn/trueTau2-1),
            tauAbsPctBiasMn=100*mean(abs(tau2.mn/trueTau2-1)),
            tauAbsPctBiasMax=100*max(abs(tau2.mn/trueTau2-1)))

#Figure 8a.
currentFig=simResults$modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NlabList))%>%
  ggplot(aes(x=trueTau2,y=tau2.mn,color=factor(N)))+
  scale_x_log10()+scale_y_log10()+
  geom_point()+
  facet_wrap(~methodLab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=nColorList)+
  geom_jitter(width= 0.05)+
  #  geom_hline(yintercept=0,color ="blue")+
  geom_abline(intercept=0,slope=1,color ="blue")+
  labs(x=expression(paste("True ", tau^2)),
       y=expression(paste( tau^2, "estimated ", " from simulations")), 
       color = expression(N[S]))

nFig=nFig+1
figLegText="True vs Estimated tau2"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. ",nFig,": ", figLegText,sep=""))

#These alternatives are not used in the paper.
#Figure 8b.
simResults$modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  ggplot(aes(x=trueTau2,y=tau2.mn/trueTau2-1,color=factor(N)))+
  geom_point()+
  facet_wrap(~method)+
  scale_color_manual(values=nColorList)+
  geom_jitter(width= 0.05)+
  geom_abline(intercept=0,slope=1,color ="blue")+
  labs(x=expression(paste("True ", tau^2)),
       y=expression(paste("Relative error of ", tau^2, "and estimated ", tau^2, " from simulations")), 
       color = expression(N[S]))

#Figure 8c.
simResults$modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  ggplot(aes(x=trueTau2,y=tau2.mn/trueTau2-1,color=factor(N)))+
  geom_point()+
  facet_wrap(~method)+
  scale_color_manual(values=nColorList)+
  geom_jitter(width= 0.05)+
  geom_abline(intercept=0,slope=1,color ="blue")+
  labs(x=expression(paste("True ", tau^2)),
       y=expression(paste("Difference between ", tau^2, "and estimated ", tau^2, " from simulations")), 
       color = expression(N[S]))

####################################################################################
####################################
#TABLES FOR SUPPLEMENTARY MATERIALS
####################################
#Tables 1 -3 formed outside R.
#Table S1.  Range of weights from negative binomial method and grid method to generate study sample size distributions.
nbWtsRange=nbWtsData%>%filter(tau2>0&effectSize==0.5)%>%group_by(icc,effectSize,mnN)%>%summarise(minWt=min(nbWts.exDh),maxWt=max(nbWts.exDh))
gdWtsRange=gdWtsData%>%filter(tau2>0&effectSize==0.5)%>%group_by(icc,effectSize,mnN)%>%summarise(minWt=min(gdWts.exDh),maxWt=max(gdWts.exDh))

nbWtsRange=nbWtsData%>%filter(effectSize==0.5)%>%
  group_by(icc,effectSize,mnN)%>%
  summarise(minWt=min(nbWts.exDh),maxWt=max(nbWts.exDh))
gdWtsRange=gdWtsData%>%filter(effectSize==0.5)%>%
  group_by(icc,effectSize,mnN)%>%
  summarise(minWt=min(gdWts.exDh),maxWt=max(gdWts.exDh))


print(full_join(nbWtsRange,gdWtsRange%>%filter(mnN==4),by=c("icc","effectSize"))%>%
        mutate(pctCover=100*(maxWt.y-minWt.y)/(maxWt.x-minWt.x),
               n4Cover=100*(minWt.y-minWt.x)/(maxWt.x-minWt.x),
               n3Cover=pctCover+n4Cover))

write.table(full_join(nbWtsRange,gdWtsRange%>%filter(mnN==4),by=c("icc","effectSize"))%>%
              mutate(pctCover=100*(maxWt.y-minWt.y)/(maxWt.x-minWt.x),
                     n4Cover=100*(minWt.y-minWt.x)/(maxWt.x-minWt.x),
                     n3Cover=pctCover+n4Cover),paste(Sys.Date()," nbVsGdWtCover.csv",sep=","))

####################################
#FIGURES FOR SUPPLEMENTARY MATERIALS
####################################
##Figure SX.  Comparison of I_1^2(tau2) versus I_v^2(tau2).

iccCompData=expand.grid(N=c(3,4,6,9,16,25),  #ICC0 follows N = 4 and delta = 1.5 - 2.0 (approx).
                        effectSize=c(0,0.5,1,1.5,2,4),
                        ICC0=unique(simResults$modFnRSd$popnICC))%>%filter(ICC0>0)%>%
  mutate(v2=exactJ.r(N+N-2)^2*exactVarDu.r(N,N,effectSize),
         tau2=ICC0/(1-ICC0),
         ICC2=ICC0/(ICC0+2*(1-ICC0)),
         ICCv=ICC0/(ICC0+v2*(1-ICC0)),
         ICCva=tau2/(tau2+v2))

currentFig=iccCompData%>%
  mutate(NIv2Lab=factor(N,levels=c(3,4,6,9,16,25),labels=NIv2LabList))%>%
  ggplot()+
  geom_point(aes(x=ICC0,y=ICCv,color=factor(effectSize)))+
  geom_abline(intercept=0,slope=1,color="blue")+
  facet_wrap(~NIv2Lab,labeller=label_parsed)+
  scale_color_manual(values=c(effectSizeColorList,"hotpink"))+
  labs(x=expression(I[1]^2),y=expression(I[v]^2),color=expression(delta))

nFig=0
figLegText="Comparison of ICC1 and ICCv for same tau2"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S1.  Compare I^2 with Iv^2
testS2Data=unique(data.frame(N=simResults$modFnRSd$N,
                             M=simResults$modFnRSd$M,
                             effectSize=simResults$modFnRSd$effectSize,
                             ICC0=simResults$modFnRSd$popnICC))%>%
  filter(ICC0>0)%>%
  mutate(tau2=ICC0/(1-ICC0))%>%
  select(M,N,effectSize,tau2,ICC0)

negBinS2=NULL
for(i in 1:dim(testS2Data)[1]){
  negBinResult=negBinTypicalS2.r(testS2Data$M[i],testS2Data$N[i],testS2Data$effectSize[i],testS2Data$tau2[i])
  negBinS2=rbind(negBinS2,data.frame(typS2=negBinResult[1],arthS2=negBinResult[2]))
}
testS2Data=data.frame(testS2Data,negBinS2)
names(testS2Data)
mean(testS2Data$tau2)
currentFig=testS2Data%>%mutate(Nlab=factor(N,labels=NaLabList))%>%
  mutate(tau2panelLab=factor(tau2,labels=T2FracLabList[-1]))%>%#filter(N < 20)%>%
  ggplot()+
  geom_point(aes(x=tau2/(tau2+typS2),y=tau2/(tau2+arthS2),color=tau2panelLab))+
  geom_abline(intercept=0,slope=1,color="blue")+
  geom_abline(intercept=0.33,slope=0,color="red")+
  geom_vline(xintercept=0.33,color="red")+
  facet_wrap(~Nlab,labeller=label_parsed)+
  scale_color_manual(values=popnIccColorList[-1])+
  labs(x=expression(I^2),
       y=expression(I[v]^2),color=expression(tau^2))

nFig=nFig+1
figLegText="Compare Iv^2 with I^2"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S2.  Comparison of variance functions $v(LS,d_*)$ and $v(EX.d_*)$ and $N_S$.  The purple line is the relationship between $J_e$ and $N_S$.

currentFig=data.frame(expand.grid(n=3:40,delta=(0:5)/2.5))%>%
  mutate(nf=n/2,df=2*n-2,
         Je=exactJ.r(df),    
         dh.lsDhVar=Je^2*approxVarDu.r(n,n,delta/Je),
         dh.exDhVar=Je^2*exactVarDu.r(n,n,delta))%>%
  mutate(delta=factor(delta))%>%rename(obsEffectSize=delta)%>%
  ggplot()+
  geom_line(aes(x=n,y=dh.lsDhVar/dh.exDhVar,colour=obsEffectSize))+
  geom_line(aes(x=n,y=Je),colour="purple",size=2)+
  geom_point(aes(x=n,y=dh.lsDhVar/dh.exDhVar,colour=obsEffectSize))+
  scale_color_manual(values=deltaColorList)+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste("Sample size (",n[T]," = ", n[C],")")),
       y=expression(paste("v(LS, ",d[Ui],")/v(EX, ",d[Hi],")")),
       color=expression(delta[i]))

nFig=nFig+1
figLegText="Function Approximation: Relative error of Var[LS](dH) to Var[EX](dH)"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S3.  Proportion of total variance attributed to $\tau^2$, $a(n_i)$, and $b(n)\delta_i^2$.  Observe that $\tau^2$ quickly becomes dominate source of variance.

termColorList=c("red","darkorange","black")
currentFig=expand.grid(n=3:40,effectSize=c(0,0.5,1,2),tau2=c(0,0.5,1,2))%>%
  mutate(nf=nf.r(n,n),
         df=n+n-2,
         Je=exactJ.r(df),
         deltaLab=factor(effectSize,levels=c(0,0.5,1,2),labels=deltaLabList),
         tau2a=factor(tau2,levels=c(0,0.5,1,2),labels=tauLabList),
         aN=Je^2*df/((df-2)*nf),bN=Je^2*effectSize^2*(df/(df-2)-1/Je^2))%>%
  filter(n<=20)%>%select(n,deltaLab,tau2a,aN,bN,tau2)%>%
  gather(key="term",value="value",aN,bN,tau2)%>%
  mutate(fillLab=factor(term,levels=c("bN","aN","tau2"),labels=c("b(n)","a(n)",expression(paste(tau^2)))))%>%
  ggplot(aes(fill=fillLab,y=value,x=n))+
  geom_bar(position="fill",stat="identity")+
  facet_grid(tau2a~deltaLab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_fill_manual(values=termColorList,labels = parse_format())+
  labs(x=expression(n[i]),y=expression(paste("Proportion of Total variance (", tau^2, " + ",upsilon[Hi]^2,")")),fill="Term")

nFig=nFig+1
figLegText="n vs Propn of total variance"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S4. Theoretical weights when $\tau^2$ = 0 for different $N_S$ and $\delta_i$..  

currentFig=sigmaRatioData%>%mutate(effectSize=factor(effectSize))%>%filter(tau2==0&n<=40)%>%
  ggplot()+
  geom_line(aes(x=n,y=1/varEX,color=effectSize))+
  scale_color_manual(values=c("Black","Red","Blue","darkorange","DarkGreen","darkgoldenrod"))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste("Sample size (",n[T]," = ", n[C],")")),
       y=expression(paste("w = 1/(",upsilon[i]^2,")")),
       color=expression(delta))

nFig=nFig+1
figLegText="Theoretical Weights: Fixed Effects"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S5. Weight function $\omega_i^T = \tau^2/(\tau^2 + \upsilon_i^2)$ against $N_S$ for different $\delta_i$ and $\tau^2$.

currentFig=sigmaRatioData%>%mutate(effectSize=factor(effectSize))%>%filter(tau2!=0&n<=40)%>%
  mutate(tau2lab=factor(tau2,labels=T2labList))%>%
  ggplot()+
  geom_line(aes(x=n,y=ICC,color=effectSize))+
  scale_color_manual(values=c("Black","Red","Blue","darkorange","DarkGreen","darkgoldenrod"))+
  facet_wrap(~tau2lab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste("Sample size (",n[T]," = ", n[C],")")),
       y=expression(paste(w[i], " = ",tau^2, "/(",tau^2, " + ",upsilon[i]^2,")")),
       color=expression(delta[i]))

nFig=nFig+1
figLegText="Theoretical Weights: A comparison"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S6.  Distribution of negative binomial method with grid method results for comparison.

#Scale weights for negBin (nb) and grid (gd) methods to case with n = 3 and delta = 0.5 
nbWtsScData=full_join(nbWtsData%>%filter(effectSize==0.5),
                      unique(nbWtsData%>%filter(N==3,effectSize==0.5)%>%select(-N,-effectSize)),
                      by=c("icc","mnN"))

gdWtsScData=full_join(gdWtsData%>%filter(effectSize==0.5),
                      unique(nbWtsData%>%filter(N==3,effectSize==0.5)%>%select(-N,-effectSize)),
                      by=c("tau2","icc","mnN"))

#The official entry
names(nbWtsScData)
table(nbWtsScData$icc)
currentFig=nbWtsScData%>%
  mutate(tau2=factor(icc,labels=T2simLabList))%>%
  filter(effectSize==0.5&icc<0.8)%>%
  ggplot()+
  geom_histogram(aes(x=nbWts.exDh.x/nbWts.exDh.y))+
  facet_grid(factor(mnN,labels=NstarLabList)~tau2,scales="free",labeller=label_parsed)+
  geom_vline(data=gdWtsScData%>%mutate(kN=N/mnN)%>%
               mutate(tau2=factor(icc,labels=T2simLabList))%>%
               filter((effectSize==0.5&icc<0.8)&(kN==1|kN==4|kN==8)),
             aes(xintercept = gdWts.exDh/nbWts.exDh,color=factor(kN)))+
  scale_color_manual(values=c("Blue","orange","red"))+
  labs(x="Weights",y="Count",color=expression(k[N]))

nFig=nFig+1
figLegText="Study Sample Size Distn from NegBin method"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S7.  Relative efficiency using negative binomial method for study sample size generation. 
#Results for relative RMSE using negative binomial distribution to generate study sample sizes for meta-analysis simulations.  
#RMSE for a selected set of estimators relative to RMSE for d_H.w(EQ) as a function of I_0^2 and N_S. 
#Differences between N_S=16 and N_S=25 were negligible.

currentFig=relRmseNbData%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NaLabList[1:3]),
         tau2=factor(popnICC/(1-popnICC),labels=c("0","1/11","1/2","1","2","9")))%>%
  filter(N<20&(method=="dh.n"|method=="dh.neDh"|method=="du.lsDu"))%>%
  ggplot(aes(x=tau2,y=relRMSE,color=factor(M)))+
  facet_grid(Nlab~methodLab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=mColorList)+
  geom_jitter(width= 0.05)+
  labs(x=expression(tau^2),y=expression(paste("RMSE(.)/RMSE(",d[H],".w(EQ)")),
       color="M")

nFig=nFig+1
figLegText="Relative RMSE for weighted versus unweighted averages"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S8.  Ratio of Relative Efficiencies of the two study sample size distribution methods.
currentFig=right_join(relRmseNbData%>%select(method,M,popnICC,effectSize,N,ICC0,relRMSE),
           relRmseData%>%filter(N==4&kN==4)%>%select(method,M,popnICC,effectSize,ICC0,relRMSE),
           by=c("method","M","popnICC","effectSize"))%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NaLabList[1:3]),
         tau2=factor(popnICC/(1-popnICC),labels=c("0","1/11","1/2","1","2","9")))%>%
  filter(N<20&(method=="dh.n"|method=="dh.neDh"|method=="du.lsDu"))%>%
  #  mutate(ICC0=factor(popnICC,labels=c("0","1/10","1/3")))%>%
  ggplot(aes(x=tau2,y=relRMSE.y/relRMSE.x,color=factor(M)))+
  facet_grid(Nlab~methodLab,labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=mColorList)+
  geom_jitter(width= 0.05)+
  labs(x=expression(tau^2),y=expression(paste(e[g],"/",e[n]),sep=""),
       color="M") 

nFig=nFig+1
figLegText="Ratio of relative efficiencies for negBin and n = (4, 16) grid"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S9.  Percent bias for fixed effects, i.e., $I_0^2 = 0$ is known and weights are also known. 

currentFig=modFnFCd.wtsKnwn%>%
  filter(effectSize>0&(method%in%biasMethodList))%>%
  mutate(methodLab=factor(as.character(method),levels=biasMethodList,labels=biasMethodLabelList),
         method=factor(as.character(method),levels=biasMethodList))%>%
  ggplot(aes(x=N,y=100*bias/effectSize,color=factor(kN)))+
  geom_jitter(width=0.5)+
  facet_wrap(~methodLab,nrow=2, labeller = label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=kNcolorList)+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y="Percent Bias",color=expression(k[N]))

nFig=nFig+1
figLegText="N vs % Bias when tau2 = 0, Weights Known"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S10.  Percent bias for fixed effects, i.e., $I_0^2 = 0$ is known, but weights are estimated.

currentFig=modFnFCd%>%
  filter(effectSize>0&(method%in%biasMethodList))%>%
  mutate(methodLab=factor(as.character(method),levels=biasMethodList,labels=biasMethodLabelList),
         method=factor(as.character(method),levels=biasMethodList))%>%
  ggplot(aes(x=N,y=100*bias/effectSize,color=factor(kN)))+
  geom_jitter(width=0.5)+
  facet_wrap(~methodLab,nrow=2, labeller = label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=kNcolorList)+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y="Percent Bias",color=expression(k[N]))

nFig=nFig+1
figLegText="N vs % Bias when tau2 = 0, Weights estimated"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S11.  Percent bias for random effects, $I_0^2 = 1/3$ and weights are treated as known in the simulations.

currentFig=modFnRSd2.wtsKnwn%>%
  filter(effectSize>0&(method%in%biasMethodList))%>%
  mutate(methodLab=factor(as.character(method),levels=biasMethodList,labels=biasMethodLabelList),
         method=factor(as.character(method),levels=biasMethodList))%>%
  ggplot(aes(x=N,y=100*bias/effectSize,color=factor(kN)))+
  geom_point()+
  facet_wrap(~methodLab,nrow=2, labeller = label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=kNcolorShortList)+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y="Percent Bias",color=expression(k[N]))

nFig=nFig+1
figLegText="N vs % Bias when ICC0 = 1/3, Weights Known"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S12.  Percent bias for random effects, $I_0^2 = 1/3$ and weights are estimated in the simulations.

currentFig=modFnRSd2%>%
  filter(popnICC==(1/3))%>%
  filter(effectSize>0&(method%in%biasMethodList))%>%
  mutate(methodLab=factor(as.character(method),levels=biasMethodList,labels=biasMethodLabelList),
         method=factor(as.character(method),levels=biasMethodList))%>%
  ggplot(aes(x=N,y=100*bias/effectSize,color=factor(kN)))+
  geom_jitter(width=0.5)+
  facet_wrap(~methodLab,nrow=2, labeller = label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=kNcolorShortList)+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y="Percent Bias",color=expression(k[N]))

nFig=nFig+1
figLegText="N vs % Bias when ICC0 = 1/3, Weights Estimated"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S13.  Coverage of distribution for estimated Var(d) by plotting $Var(d_q)$, where $d_q$ is the $2.5^{th}$ and $97.5^{th}$ percentile of distribution of $d$.
#Data not provided.  See TheoryAndSummaryDataForMEE.txt for code to calculate this data.
currentFig=full_join(varCI.exData,varCI.lsData,by=c("n1","n2","effectSize"),suffix=c(".ex",".ls"))%>%
  mutate(varLBs.ex=min(varLB.ex/varMed.ex,varMed.ex/varMed.ex),
         varMeds.ex=varMed.ex/varMed.ex,
         varUBs.ex=varUB.ex/varMed.ex,
         varLBs.ls=min(varLB.ls/varMed.ex,varMed.ls/varMed.ex),
         varUBs.ls=varUB.ls/varMed.ex,
         effectSizeLab=factor(effectSize,levels=c(0,0.4,0.8,1.2,1.6,2),labels=effSizeLabList))%>%
  ggplot()+
  geom_line(aes(x=n1,y=varLBs.ls),color="red")+
  geom_line(aes(x=n1,y=varMeds.ex),color="blue")+
  geom_line(aes(x=n1,y=varUBs.ls),color="red")+
  geom_line(aes(x=n1,y=varLBs.ex),color="black")+
  geom_line(aes(x=n1,y=varUBs.ex),color="black")+
  facet_wrap(~effectSizeLab,labeller=label_parsed,scales="free")+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste("Sample size (",n[T]," = ", n[C],")")),y=expression(paste("Standardized Confidence bounds for (v(f,d)/v(EX,",delta[H],")")))

nFig=nFig+1
figLegText="95% CI for observed SE(dH) with Estimators EX and LS"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S14.  Proportion of simulation runs with $\delta_P$ found to be in 95% confidence intervals against $M$.  This covers the set of unbiased estimators and the nearly unbiased estimator $d_U.w(LS,d_U)$.

currentFig=simResults$modFnRSd%>%select(method,N,M,popnICC,pctSig)%>%
  filter(method=="dh.eq"|method=="dh.n"|method=="dh.neDh"|method=="dh.xmDh"|method=="dh.lmDh"|method=="du.lsDu")%>%
  mutate(tau2=factor(popnICC/(1-popnICC),labels=T2FracLabList))%>%
  ggplot(aes(x=M,y=pctSig,color=tau2))+
  geom_jitter(width=1)+
  scale_color_manual(values=c(popnIccColorList,"purple"))+
   facet_wrap(~factor(method,levels=methodList,label=methodLabelList),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=paste("M = number of endpoints"),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(tau^2))

nFig=nFig+1
figLegText="M vs % Significant tests in Simulation"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S15.  Proportion of simulation runs with $\delta_P$ found to be in 95% confidence intervals against $I_0^2$.  This covers the set of unbiased estimators and the nearly unbiased estimator $d_U.w(LS,d_U)$.

currentFig=full_join(simResults$modFnRSd%>%select(method,N,kN,M,effectSize,popnICC,pctSig),testGridData,
                     by=c("popnICC","effectSize","M","N","kN"))%>%
  mutate(ICCv=tau2/(tau2+avgS2))%>%
  filter(method=="dh.eq"|method=="dh.n"|method=="dh.neDh"|method=="dh.xmDh"|method=="dh.lmDh"|method=="du.lsDu")%>%
  mutate(tau2=factor(popnICC/(1-popnICC),labels=T2FracLabList))%>%
  ggplot(aes(x=ICCv,y=pctSig,color=factor(N)))+
  geom_point()+
  geom_abline(intercept=0.95,slope=0,color="blue")+
  scale_color_manual(values=nColorList)+
  facet_grid(scales = "free")+
  facet_wrap(~factor(method,levels=methodList,label=methodLabelList),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste(I[v]^2)),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(N[S]))

nFig=nFig+1
figLegText="Proportion where true effect Size found in 95% CI"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S16.  Proportion of simulation runs with $\delta_P$ found to be in 95% confidence intervals against $N_S$.  This covers the set of biased estimators.  

currentFig=simResults$modFnRSd%>%select(method,N,kN,M,popnICC,pctSig)%>%#filter(N<5)%>%
  filter(method=="dh.exDh"|method=="dh.lsDh"|method=="du.eq"|method=="du.exDu"|method=="dh.wxDh"|method=="du.wxDu")%>%
  ggplot(aes(x=N,y=pctSig,color=factor(M)))+
  geom_jitter(width=0.5)+
  scale_color_manual(values=mColorList)+
  facet_grid(scales = "free")+
  facet_wrap(~factor(method,levels=methodList,label=methodLabelList),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color="M")

nFig=nFig+1
figLegText="N vs % significant tests in simulation"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S17.  Proportion of simulation runs with $\delta_P$ found to be in 95% confidence intervals against $M$.  This covers the set of biased estimators.

currentFig=simResults$modFnRSd%>%select(method,N,M,popnICC,pctSig)%>%
  filter(method=="dh.exDh"|method=="dh.lsDh"|method=="du.eq"|method=="du.exDu"|method=="du.lmDu"|method=="du.xmDu")%>%
  mutate(tau2=factor(popnICC/(1-popnICC),labels=T2FracLabList))%>%
  ggplot(aes(x=M,y=pctSig,color=factor(tau2,labels=T2FracLabList)))+
  geom_jitter(width=1)+
  scale_color_manual(values=c(popnIccColorList,"purple"))+
  facet_grid(scales = "free")+
  facet_wrap(~factor(method,levels=methodList,label=methodLabelList),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=paste("M = number of endpoints"),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(tau^2))

nFig=nFig+1
figLegText="M vs % significant tests in simulation"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S18.  Proportion of simulation runs with $\delta_P$ found to be in 95% confidence intervals against $I_0^2$.  This covers the set of biased estimators.
currentFig=full_join(simResults$modFnRSd%>%select(method,N,kN,M,effectSize,popnICC,pctSig),testGridData,
          by=c("popnICC","effectSize","M","N","kN"))%>%
  mutate(ICCv=tau2/(tau2+avgS2))%>%
  filter(method=="dh.exDh"|method=="dh.lsDh"|method=="du.eq"|method=="du.exDu"|method=="du.lmDu"|method=="du.xmDu")%>%
  mutate(tau2=factor(popnICC/(1-popnICC),labels=T2FracLabList))%>%
  ggplot(aes(x=ICCv,y=pctSig,color=factor(N)))+
  geom_point()+
  scale_color_manual(values=nColorList)+
  facet_grid(scales = "free")+
  facet_wrap(~factor(method,levels=methodList,label=methodLabelList),labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x=expression(I[v]^2),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(N[S]))

nFig=nFig+1
figLegText="Proportion where true effect Size found in 95% CI"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S19. Minimum, mean, and maximum RMSE across different parameter scenarios for each of the biased estimators.  While under certain conditions other biased estimators have lower RMSE than $d_U.w(LS,d_U)$, the benefit is minimal and under other conditions the other biased estimators do considerably worse than $d_U.w(LS,d_U)$.
#Data not provided.  See TheoryAndSummaryDataForMEE.txt for code to calculate this data.
factor(rmseSumData$ICC0,labels=I2LabList)
levels(rmseSumData$ICC0)

names(rmseSumData)
levels(rmseSumData%>%
         mutate(tau2=factor(as.character(ICC0),levels=levels(ICC0),labels=T2simLabList))%>%select(tau2))

currentFig=rmseSumData%>%
  mutate(tau2=factor(as.character(ICC0),levels=levels(ICC0),labels=T2simLabList))%>%
  gather(key="metric",value="value",RMSE.min,RMSE.mn,RMSE.max)%>%
  separate(metric,c("Stat","Summary"))%>%
  mutate(Summary=factor(Summary,levels=c("min","mn","max"),labels=c("Min","Mean","Max")))%>%
  filter(ICC0 == "0" | ICC0 == "1/10"| ICC0 == "1/3"|ICC0 == "1/2")%>%
  filter(method=="dh.exDh"|method=="dh.lsDh"|method=="du.eq"|method=="du.exDu"|method=="du.lsDu"|method=="dh.wxDh")%>%
  ggplot(aes(x=N,y=value,color=factor(method,levels=methodList,label=methodLabelList)))+
  facet_grid(tau2~Summary,scales="free",labeller=label_parsed)+
  theme(strip.text.x = element_text(size = facetFontSize),
        strip.text.y = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  scale_color_manual(values=biasedMethodColorList,labels = parse_format())+
  geom_point()+
  geom_line()+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y=expression(paste("RMSE(.)/RMSE(",d[H],"w(EQ))")),
       color="Estimator")

nFig=nFig+1
figLegText="RMSE for Methods relative to RMSE(dh.w(EQ))"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

#Figure S20.  Bland-Altman plot of difference between replicated simulations against average of the two reps.  Each point denotes a replication pair for each of the 2592 simulation parameter settings.  

currentFig=repData%>%
  mutate(metricLab=factor(metric,levels=c("bias","pctSig","simSD"),labels=c("Bias","Percent Coverage","Simulation Distn SD")))%>%
  ggplot()+
  geom_point(aes(x=(value+value.rep)/2,y=value.rep-value))+
  facet_wrap(~metricLab,scales="free_x")+
  theme(strip.text.x = element_text(size = facetFontSize))+
  theme(axis.title = element_text(size = labFontSize))+
  labs(x="Mean of estimator trait over two simulation replications",
       y="Difference in estimator trait between two simulation replications")

nFig=nFig+1
figLegText="Bland-Altman Plot"
ppReport=addGGplot(ppReport,currentFig,paste("Fig. S",nFig,": ", figLegText,sep=""))

print(ppReport,target = ppFileName)