#################################
#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.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")

################################################
#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[0]^2," = 0")),
            expression(paste(I[0]^2," = 1/10")),
            expression(paste(I[0]^2," = 1/3")),
            expression(paste(I[0]^2," = 1/2")),
            expression(paste(I[0]^2," = 2/3")),
            expression(paste(I[0]^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")))

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")))

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("Red","Orange","Yellow","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
####################################
#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$.

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),
       #       y=expression(paste(omega[i]^0," are weights standardized by ", omega,"(d = 0|n = 3)")),
       color=expression(delta))

#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.

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(axis.title = element_text(size = labFontSize))+
  labs(x=expression(paste(N[S]," = Smallest sample size ")),
       y=expression(paste("Effect Size ", delta[P])),
       color="Estimator")


#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$.

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",
       fill=expression(N[S]))

#Figure 5.  Percent bias as a function of I_0^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.

modFnRSd%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NlabList))%>%
  filter(effectSize>0&N < 20&(method=="dh.eq"|method=="dh.lsDh"|method=="du.lsDu")&kN==4)%>%
  ggplot(aes(x=popnICC,y=100*bias/effectSize))+
  geom_jitter(aes(color=factor(effectSize)),width=0.02)+
  geom_hline(yintercept=0)+
  facet_grid(Nlab ~ 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=effectSizeColorList)+
  labs(x=expression(paste(I[0]^2)),y="Percent Bias",color=expression(delta[P]))

#Figure 6.  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.
relRmseData%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NlabList))%>%
  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=ICC0,y=relRMSE,color=factor(kN)))+
  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=kNcolorShortList)+
  geom_jitter(width= 0.05)+
  labs(x=expression(I[0]^2),y=expression(paste("RMSE(.)/RMSE(",d[H],".w(EQ)")),
       color=expression(k[N]))

#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.

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=expression(paste("RMSE(",d[U],".w(LS, ",d[U],"))/RMSE(",d[H],".w(EQ))")),
       color="M")

#Figure 8.  Estimation of tau2.
#Bias in tau2 does not impact estimation of deltaP
#Figure 8a. true vs. estimated tau2
modFnRSd%>%mutate(trueTau2=popnICC/(1-popnICC))%>%
  filter(method!="du.eq"&method!="dh.eq")%>%filter(popnICC>0)%>%
  ggplot(aes(x=trueTau2,y=tau2.mn,color=factor(N)))+
  scale_x_log10()+scale_y_log10()+
  geom_point()+
  facet_wrap(~method)+
  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]))


####################################
#TABLES FOR SUPPLEMENTARY MATERIALS
####################################
#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))

####################################
#FIGURES FOR SUPPLEMENTARY MATERIALS
####################################
#Figure S1.  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$.

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]))

#Figure S2.  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")
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")

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

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))

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

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]))

#Figure S5.  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"))

nbWtsScData%>%filter(effectSize==0.5&icc<0.8)%>%
  ggplot()+
  geom_histogram(aes(x=nbWts.exDh.x/nbWts.exDh.y))+
  facet_grid(factor(mnN,labels=NstarLabList)~factor(icc,labels=I2LabList[1:5]),scales="free",labeller=label_parsed)+
  geom_vline(data=gdWtsScData%>%mutate(kN=N/mnN)%>%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]))

#Figure S6.  Relative efficiency using negative binomial method for study sample size generation. 

relRmseNbData%>%
  mutate(methodLab=factor(method,levels=methodList,labels=methodLabelList),
         method=factor(method,levels=methodList),
         Nlab=factor(N,labels=NaLabList[1:3]))%>%
  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=ICC0,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(I[0]^2),y=expression(paste("RMSE(.)/RMSE(",d[H],".w(EQ)")),
       color="M")

#Figure S7.  Ratio of Relative Efficiencies of the two study sample size distribution methods.

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]))%>%
  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=ICC0.x,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(I[0]^2),y=expression(paste(e[g],"/",e[n]),sep=""),
       color="M") 

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

modFnFCd.wtsKnown%>%
  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]))

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

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]))

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

modFnRSd2.wtsKnown%>%
  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]))

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

modFnRSd2.wtsUnKnown%>%
  filter(popnICC<0.4&popnICC>0.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]))

#Figure S12.  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.
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],")")))

#Figure S13.  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)$.

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")%>%
  ggplot(aes(x=M,y=pctSig,color=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10"))))+
  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(paste(I[0]^2)))

#Figure S14.  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)$.

modFnRSd%>%select(method,N,kN,M,popnICC,pctSig)%>%#filter(N<5)%>%
  filter(method=="dh.eq"|method=="dh.n"|method=="dh.neDh"|method=="dh.xmDh"|method=="dh.lmDh"|method=="du.lsDu")%>%
  ggplot(aes(x=popnICC,y=pctSig,color=factor(N)))+
  geom_jitter(width=0.02)+
  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[0]^2)),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(N[S]))

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

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")

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

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")%>%
  ggplot(aes(x=M,y=pctSig,color=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10"))))+
  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(paste(I[0]^2)))

#Figure S17.  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.

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=="du.lmDu"|method=="du.xmDu")%>%
  ggplot(aes(x=popnICC,y=pctSig,color=factor(N)))+
  geom_jitter(width=0.02)+
  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[0]^2),
       y="Proportion of Sim Distn in 95% Confidence Interval",
       color=expression(N[S]))


#Figure S18. 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.

rmseSumData%>%
  #  mutate(ICC0=factor(popnICC)),labels=c("0","1/10","1/3","2/3","9/10")))%>%
  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(ICC0~Summary,scales="free")+
  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")


#Figure S19.  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.  

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")