library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(gridExtra)
library(metafor)  

#####################################################
## FUNCTIONS
#####################################################
#####################Functions for J used for bias correction in Mean and Variance
#FROM Hedges 1981 exactJ.r(df) = Je = c(m) = G(m/2)/[sqrt(m/2)*G((m-1)/2)], is exact function.
#FROM Hedges 1981 an approximation to Je: approxJ.r(df) = Ja = 1 - 3/(4m-1), where 
#G(x) = Gamma Function of x; m = df = nE + nC - 2 in our studies and nE, nC are sample sizes for endpoints.
#R can not handle Gamma functions with df > 300, so use approximate calculation.  Exceedingly negligible differences at this point.
approxJ.r=function(df){1-3/(4*df-1)}
exactJ.r=function(df){ifelse(df<300,gamma(df/2)/(sqrt(df/2)*gamma((df-1)/2)),approxJ.r(df))} #warning reflects that calculates gamma even when uses approxJ.r 

#nf.r(n1,n2) = ntilde = nE*nC/(nE + nC);
nf.r=function(n1,n2){n1*n2/(n1+n2)}

#Functions for calculating Var(d) using various estimation formulas.
#All are for var(du)  to get Var(dh) multiply by Je = exactJ.r(df).

#FROM Hedges 1981 Unadjusted variance:  exactVar.r(n1,n2,effSize) = Var(g) = m/[(m-2)*ntilde]*[1+ntilde*delta^2]-delta^2/c(m)^2
exactVarDu.r=function(n1,n2,effSize){
  nf=nf.r(n1,n2)
  df=n1+n2-2
  Je=exactJ.r(df)
  #  df/((df-2)*nf)*(1+nf*effSize^2)-(effSize/Je)^2
  df/((df-2)*nf)+effSize^2*(df/(df-2)-1/Je^2)
}

#Var_EX(Dh) = exactJ.r(df)^2*exactVarDu.r(n1,n2,effectSize)
#Approximate with mean delta is simply exact variance with mean delta
#exVarDu.mnDelta=exactVarDu.r(n1,n2,mean(effSize))

#FROM Hedges 1982 Large Sample Approximation to dH = gU.
#approxVar.r(n1,n2,effSize) = Var(gU) = (nE + nC)/(nE*nC) + delta^2/[2*(nE+nC)]
#This is for ADJUSTED d, dH = gU.  Var(dU) = Var(g) = Var(gU)/Je^2
approxVarDu.r=function(n1,n2,effSize){
  nf=nf.r(n1,n2)
  df=n1+n2-2
  Je=exactJ.r(df)
  (1/nf+effSize^2/(2*(n1+n2)))/Je^2
}

#Approximate with mean delta is simply large sample theory with mean delta
#approxVarDu.mnDelta=approxVarDu.r(n1,n2,mean(effSize))

#Hedges 1985 Calculate weights by dropping component with effect size; so only depends on N1 and N2.
noEffSizeVarDu.r=function(n1,n2){
  df=n1+n2-2
  df/(nf.r(n1,n2)*(df-2))
}

#White propose estimates based on variance of estimated variance.
#For dH
whiteVarDh.r=function(n1,n2,dh){
  nf=nf.r(n1,n2)
  df=n1+n2-2
  Je=exactJ.r(df)
  1/nf+dh^2*(1-(df-2)/(df*Je^2))
}

#For dU
whiteVarDu.r=function(n1,n2,du){
  nf=nf.r(n1,n2)
  df=n1+n2-2
  Je=exactJ.r(df)
  1/(Je^2*nf)+du^2*(1-((df-2)/(df*Je^2)))
}

#Calculate weighted average
calcWtMn.r=function(wt,x){sum(wt*x)/sum(wt)}

#Calculate variance using weights, 
#which assumes weights are reasonable estimates of underlying variance, so any old weights will not do...
calcSumWtsVar.r=function(x,nStudy){1/(sum(x))}

#Calculates weighted SSE for delta in simple weighted average when tau2 > 0.
#Generated from tau2 and weights.
calcWtdSS.r=function(obs,mn,tau2,wts){
  tWt=1/(tau2+1/wts)
  sum(tWt*(obs-mn)^2)
}

#Generates sample of Ns for the  M Studies of meta-analysis.
#Creates equal distribution over each value of N included in possible sample sizes.  
#(First Ns in vector get extra N when M is not a multiple of number unique N in meta-analysis)
genFixedSampleN.r=function(N,M){  #N should be a vector, M is scalar.
  Nvec=rep(floor(M/length(N)),length(N))
  if(M!=sum(Nvec)) {
    shortN=M-sum(Nvec)
    Nvec[1:shortN]=Nvec[1:shortN]+1
  }
  rep(N,Nvec)
}

#genFixedSampleN.r(c(4,6,12),25)

############################
#SE Estimators
############################
#Don't think this is used in MEE
#trueWtMn.r=function(p,wt,mn){
#    sum(p*wt*mn)/sum(p*wt)  #generalize for pi i = 1, ... S, S > 2, sum pi = 1.
#}

#hauptmann-knapp estimate of variance(d) 
hkVarEst.r=function(d,w){
  wtAvg=sum(w*d)/sum(w)
  sum(w*(d-wtAvg)^2)/(sum(w)*(length(w)-1))  
}

#Confidence interval for observed var(d) given N and delta.  EX and LS are Exact and Large Sample approximation
varCI.r=function(n1,n2,delta,method="EX"){
  df=n1+n2-2
  Je=exactJ.r(df)
  if(method=="EX") result=Je*sqrt(exactVarDu.r(n1,n2,delta+qnorm(c(.025,0.5,.975))*Je*sqrt(exactVarDu.r(n1,n2,delta))))
  if(method=="LS") result=Je*sqrt(approxVarDu.r(n1,n2,(delta+qnorm(c(.025,0.5,.975))*Je*sqrt(exactVarDu.r(n1,n2,delta)))))
  if(method=="WH") result=sqrt(whiteVarDh.r(n1,n2,(delta+qnorm(c(.025,0.5,.975))*Je*sqrt(exactVarDu.r(n1,n2,delta)))))
  if(method=="WU") result=sqrt(whiteVarDu.r(n1,n2,(delta/Je+qnorm(c(.025,0.5,.975))*sqrt(exactVarDu.r(n1,n2,delta)))))
  
  result
}

#########################################
#THEORETICAL CALCULATIONS
#########################################

#################################################################################################
#CALCULATING MEAN FUNCTIONS
#####################################
#Code to estimate random effects model results.  Modified from Hamman 2018
#####################################
simRMA.r=function(x){ 
  n=length(x)/2
  mn=x[1:n]
  wt=x[(n+1):(2*n)]
  du_iv=rma(mn,wt,knha= TRUE,control=list(stepadj=0.5,maxiter = 10000)) # perform MA
  
  c(du_iv$b, # estimated effect size
    du_iv$se, # standard error of estimate
    du_iv$tau2,du_iv$ci.lb,du_iv$ci.ub)  #estimate of Var(p)
}

#simRMA.r(c(obsSimData$MnDu[,i],obsSimData$Wts.lmDh[,i]))


########################################################################################
#OUTER SHELL FOR SIMULATION RUN:
#########################################################################################
#Generalized version, so S > 2. pi = 1/S
#Function that drives entire simulation.
#This runs several different effects model, driven by 
#genNtype, which determines whether sample sizes are random "Rn" or fixed "Fn", and
#genDtype, which determines whether delta are random "RS" or fixed "FC".
# The function takes list of parameters (simParamList), creates full factorial of parameter combinations and
#runs simulation one parameter combination at a time 
simGenMA_Methods.r=function(nSim,simParamList,genNtype,genDtype,wtsKnown=FALSE){
  
  #  methodListN=length(methodList)  
  #range of parameter values examined
  
  totalVarVec=simParamList$totalVarVec
  #totalVarVec:  Everything in this system derives from of totalVarVec, which is assumed = 1 in simulations. 
  
  popnICCVec=simParamList$popnICCVec  #Population ICC = proportion of PopnVar to Total = Vpopn/(Vpopn + Verr).
  effectSizeVec=simParamList$effectSizeVec 
  
  Mvec=simParamList$Mvec  #number of endpoints
  
  Svec=simParamList$Svec  #number of unique sample sizes.
  Nvec=simParamList$Nvec  #sample size for endpoint
  kNvec= simParamList$kNvec  #creates S tiers of sample size between N and kN*N inclusive across endpoints in meta-analysis
  
  simDesign=expand.grid(totalVar=totalVarVec,popnICC=popnICCVec,effectSize=effectSizeVec,
                        M=Mvec,Sn=Svec,N=Nvec,kN=kNvec)
  
  resultsList=vector("list",dim(simDesign)[1])
  
  #Simulation parameters selected for each simulation run
  for(ii in 1:dim(simDesign)[1]){
    totalVar=simDesign$totalVar[ii]
    
    popnICC=simDesign$popnICC[ii]
    
    popnSD = sqrt(popnICC*totalVar)
    studySD=sqrt(totalVar-popnSD^2)
    
    effectSize=simDesign$effectSize[ii]
    
    popnMn = effectSize*studySD
    
    N1 = N2 = floor((0:(simDesign$Sn[ii]-1))*simDesign$N[ii]*(simDesign$kN[ii]-1)/(simDesign$Sn[ii]-1) + simDesign$N[ii])
    
    nStudy=simDesign$M[ii]
    
    #################################
    #Generate simulated study results
    obsSimData=simRun.r(nSim,nStudy,popnMn,popnSD,N1,studySD,genNtype,genDtype,wtsKnown)
    
    resultsList[[ii]]=simSummary.r(obsSimData,effectSize,popnICC,totalVar,popnMn,popnSD,studySD,simDesign$N[ii],simDesign$kN[ii],simDesign$Sn[ii],nStudy,genDtype)
    
  }  
  simResults = do.call(rbind, resultsList)
  
  simResults
  
}

##############################################
#MAIN FUNCTIONS FOR SIMULATION
##############################################
#generates "raw" data, a vector of nSim*nStudy endpoints across nSim simulations of meta-analysis with nStudy endpoints. 
simRun.r=function(nSim,nStudy,popnMn,popnSD,N1,studySD,genNtype="Rn",genDtype="RS",wtsKnown=FALSE){
  
  simIter=rep(1:nSim,rep(nStudy,nSim))
  
  if(genNtype=="Rn"){
    simN1 = simN2 = sample(N1,nStudy*nSim,replace=TRUE,prob=rep(1/length(N1),length(N1)))} else 
      if(genNtype=="Fn"){
        simN1 = simN2 = rep(genFixedSampleN.r(N1,nStudy),nSim) } else {
          tmpN1=rnbinom(nStudy*nSim,size=1.55,mu=N1)
          simN1 = simN2 = ifelse(tmpN1<3,3,tmpN1)} # determine number of replicates.  Assuming N1 = N2...
  
  #Calculate Je    
  Je=exactJ.r(simN1+simN2-2)
  
  #Create true study means.
  if(genDtype=="RS") studyMn=rnorm(nStudy*nSim,popnMn,popnSD) else
    studyMn=rep(rnorm(nStudy,popnMn,popnSD),nSim) #Fixed effects model sets popnSD = 0. Changes for each trial in sim only for RS.
  
  #These are the observed values from each study in each simulation. 
  errMn=rnorm(nStudy*nSim,studyMn,studySD*sqrt(1/simN1+1/simN2))
  errSD=sqrt(studySD^2*rchisq(nStudy*nSim,simN1+simN2-2)/(simN1+simN2-2))
  
  #generate endpoint effect size values for dU and dH estimators
  MnDu=errMn/errSD
  MnDh=Je*MnDu
  
  obsMnDu=rep(tapply(MnDu,simIter,mean),rep(nStudy,nSim))
  obsMnDh=rep(tapply(MnDh,simIter,mean),rep(nStudy,nSim))
  
  obsStudyMn=rep(tapply(studyMn,simIter,mean),rep(nStudy,nSim))
  
  #Calculate Weights for different methods.
  #weights are for each endpoint based on N and estimated delta, MnDu or MnDh if weights unknown, studyMn/studySD if known.
  
  
  Wts.eq=rep(1/nStudy,nStudy*nSim)
  Wts.n=simN1
  Wts.neDh=1/(Je^2*noEffSizeVarDu.r(simN1,simN2))
  
  #Hedges & Olkin 1985 method: best method if assume fixed effect model.
  #Large sample approximation of variance, du and dH
  if(!wtsKnown) {
    Wts.lmDu=1/approxVarDu.r(simN1,simN2,obsMnDu)
    Wts.lmDh=1/(Je^2*approxVarDu.r(simN1,simN2,obsMnDh))
    
    Wts.xmDu=1/exactVarDu.r(simN1,simN2,obsMnDu)      
    Wts.xmDh=1/(Je^2*exactVarDu.r(simN1,simN2,obsMnDh))
    
    Wts.lsDu=1/approxVarDu.r(simN1,simN2,MnDu)
    Wts.lsDh=1/(Je^2*approxVarDu.r(simN1,simN2,MnDh))
    
    Wts.exDu=1/exactVarDu.r(simN1,simN2,MnDu)
    Wts.exDh=1/(Je^2*exactVarDu.r(simN1,simN2,MnDh))
    
    Wts.wxDu=1/whiteVarDu.r(simN1,simN2,MnDu)
    Wts.wxDh=1/whiteVarDh.r(simN1,simN2,MnDh)
  }
  
  if(wtsKnown) {
    Wts.lmDu=1/approxVarDu.r(simN1,simN2,obsStudyMn/studySD)
    Wts.lmDh=1/(Je^2*approxVarDu.r(simN1,simN2,obsStudyMn/studySD))
    
    Wts.xmDu=1/exactVarDu.r(simN1,simN2,obsStudyMn/studySD)
    Wts.xmDh=1/(Je^2*exactVarDu.r(simN1,simN2,obsStudyMn/studySD))
    
    Wts.lsDu=1/approxVarDu.r(simN1,simN2,studyMn/studySD)
    Wts.lsDh=1/(Je^2*approxVarDu.r(simN1,simN2,studyMn/studySD))
    
    Wts.exDu=1/exactVarDu.r(simN1,simN2,studyMn/studySD)
    Wts.exDh=1/(Je^2*exactVarDu.r(simN1,simN2,studyMn/studySD))
    
    Wts.wxDu=1/whiteVarDu.r(simN1,simN2,studyMn/studySD)
    Wts.wxDh=1/whiteVarDh.r(simN1,simN2,studyMn/studySD)
  }  
  
  #  }
  
  data.frame(simIter,simN1,simN2,studyMn,errMn,errSD,MnDu,MnDh,
             Wts.eq,Wts.n,Wts.neDh,
             Wts.xmDu,Wts.xmDh,Wts.lmDu,Wts.lmDh,
             Wts.exDu,Wts.exDh,Wts.lsDu,Wts.lsDh,
             Wts.wxDu,Wts.wxDh)
}


#######################################################
#Summarizes simulation Data
simSummary.r=function(obsSimData,effectSize,popnICC,totalVar,popnMn,popnSD,studySD,N,kN,Sn,nStudy,genDtype){
  
  nSim=max(obsSimData$simIter)
  
  if(genDtype=="RS") {
    reResult=simMnCalc(obsSimData,genDtype) 
    obsSimMn=reResult$mn
    reSE.mn=sqrt(apply(reResult$se^2,2,mean))
    remCIpctSig=c(NA,NA,apply(effectSize<reResult$ci.ub&effectSize>reResult$ci.lb,2,mean))
    pctTau2eq0=c(NA,NA,apply(reResult$estTau2<0.001,2,mean))
    tau2.mn=c(NA,NA,apply(reResult$estTau2,2,mean))  #Result for du.eq and dh.eq not computed, so reported as NA.
    tau2.sd=c(NA,NA,sqrt(apply(reResult$estTau2,2,var)))      #Result for du.eq and dh.eq not computed, so reported as NA.
    
    #Expands mean values to have same length as simulation*nStudy vector for use in calculations below.  
    obsEstTau2=as.data.frame(apply(reResult$estTau2,2,function(x,nStudy,nSim){rep(x,rep(nStudy,nSim))},nStudy,nSim))
    names(obsEstTau2)=paste("tau2",names(obsEstTau2),sep=".")
  } else {
    obsSimMn=simMnCalc(obsSimData,genDtype)
    obsEstTau2=NULL
  }
  
  #Expands mean values to have same length as simulation*nStudy vector for use in calculations below.  
  obsSimMnRep=as.data.frame(apply(obsSimMn,2,function(x,nStudy,nSim){rep(x,rep(nStudy,nSim))},nStudy,nSim))
  
  #Estimates of variance based on inverse of sum of weights formula.    
  sumWtsVar=simSumWtsVarCalc(obsSimData,genDtype,obsEstTau2)  
  
  wtdSS=simWtdSSCalc(obsSimData,obsSimMnRep,genDtype,obsEstTau2) #sum weighted deviations for H-K method.
  
  hkSEobs=sqrt(wtdSS*sumWtsVar/(nStudy-1))  #H-K method:  method is equivalent to rma(method="FE",knha=TRUE)
  
  modWtdSS=apply(wtdSS,2,function(x){ifelse(x<1,1,x)})  #modified H-K
  modhkSEobs=sqrt(modWtdSS*sumWtsVar/(nStudy-1))
  
  #Generate HK-like estimate of variance, i.e., weighted variance calculation, using study variance only, no tau2.
  if(genDtype=="RS"){
    eWtdSS=simWtdSSCalc(obsSimData,obsSimMnRep,"FC") #sum weighted deviations for H-K method. "FC" used to force use of err var only.
    sumeWtsVar=simSumWtsVarCalc(obsSimData,"FC")  
    eWthkSEobs=sqrt(eWtdSS*sumeWtsVar/(nStudy-1))  #H-K method:  method is equivalent to rma(method="FE",knha=TRUE)
    eWthkSEobs.mn=sqrt(apply(eWthkSEobs^2,2,mean))
    eWthkSEobs.sd=sqrt(apply(eWthkSEobs^2,2,var))
  }
  
  ####Calculate unweighted and weighted average with effectSize known, but frequency of N varies around propn p.
  ### variance for simple average based on true parameter values ###
  simN1mat=matrix(obsSimData$simN1,nrow=nStudy,ncol=nSim,byrow=FALSE)
  exactMnDuVar=apply(simN1mat,2,
                     function(x,effectSize,M){sum(exactVarDu.r(x,x,effectSize))/M^2},effectSize,nStudy)
  
  exactMnDhVar=apply(simN1mat,2,
                     function(x,effectSize,M){sum(exactJ.r(2*x-2)^2*exactVarDu.r(x,x,effectSize))/M^2},effectSize,nStudy)
  
  ###variance for weighted average using inversve sum of weights, which are themsevles the inverse of true variance values ###
  exactWtMnDuVar=apply(simN1mat,2,
                       function(x,effectSize){1/sum(1/(exactVarDu.r(x,x,effectSize)))},effectSize)
  
  ###variance for weighted average using inverse sum of weights, which are themselves the inverse of true variance values ###
  exactWtMnDhVar=apply(simN1mat,2,
                       function(x,effectSize){1/sum(1/(exactJ.r(2*x-2)^2*exactVarDu.r(x,x,effectSize)))},effectSize)
  
  
  ############################################################################
  #Calculation of weighted and unweighted means for nSim simulations
  methodList=names(obsSimMn)
  methodListN=length(methodList)
  
  #Estimates of Population parameters based on Simulated Meta-analyses with nStudy endPoints
  simResult=data.frame(method=factor(methodList,levels=methodList),
                       effectSize=rep(effectSize,methodListN),
                       popnICC=rep(popnICC,methodListN),
                       totalVar=rep(totalVar,methodListN),
                       M=rep(nStudy,methodListN),
                       popnMn=rep(popnMn,methodListN),
                       popnSD=rep(popnSD,methodListN),
                       Sn=rep(Sn,methodListN),
                       N=rep(N,methodListN),
                       kN=rep(kN,methodListN),
                       studySD=rep(studySD,methodListN),
                       
                       #Simulated estimates of popn means from meta-analysis means across simulations
                       ###   ASSESSMENT OF MEAN ESTIMATOR    ########
                       mn=apply(obsSimMn,2,mean),  #Sim distn mean
                       bias=apply(obsSimMn,2,mean)-effectSize, #Bias from true mean
                       
                       #Simulated StdDev
                       simSD=sqrt(apply(obsSimMn,2,var)),   #sim distn SD
                       
                       
                       ####                    Assessment of Variance Estimators       ############################
                       #              Calculating mean and SD of inverse sum of weights variance calculation from simulated distn of sumWtsVar
                       sumWtsSE.mn=sqrt(apply(sumWtsVar,2,mean)),  
                       sumWtsSE.sd=sqrt(apply(sumWtsVar,2,var)),   
                       
                       hkSE.mn=sqrt(apply(hkSEobs^2,2,mean)),
                       hkSE.sd=sqrt(apply(hkSEobs^2,2,var)),
                       
                       modhkSE.mn=sqrt(apply(modhkSEobs^2,2,mean)),
                       
                       ######   Theoretical SE when delta is known, but N may vary (genNtype=="Rn")  Not sure of its purpose.                                
                       du.eqWtsSE.mn=sqrt(mean(exactMnDuVar)),  
                       du.exDhWtsSE.mn=sqrt(mean(exactWtMnDuVar)),  
                       dh.eqWtsSE.mn=sqrt(mean(exactMnDhVar)),  
                       dh.exDhWtsSE.mn=sqrt(mean(exactWtMnDhVar)),  
                       
                       #              Assessment of Mean and SD in tandem.
                       #Simulated RMSE: Summary of Bias and Var
                       rmse=sqrt(apply((obsSimMn-effectSize)^2,2,mean)), # RMSE or Sqrt(var+bias^2)
                       
                       #simulated distribution of effect size
                       simEffSize.mn=apply(obsSimMn/hkSEobs,2,mean),
                       simEffSize.sd=sqrt(apply(obsSimMn/hkSEobs,2,var)),
                       
                       #Check for bias in Confidence intervals (CI) Pct Significant
                       pctSig=apply(effectSize<(obsSimMn+qt(.975,nStudy-1)*hkSEobs)&effectSize>(obsSimMn-qt(.975,nStudy-1)*hkSEobs),2,mean),
                       modPctSig=apply(effectSize<(obsSimMn+qt(.975,nStudy-1)*modhkSEobs)&effectSize>(obsSimMn-qt(.975,nStudy-1)*modhkSEobs),2,mean))#,
  
  if(genDtype=="RS") {
    simResult=data.frame(simResult,
                         reSE.mn=reSE.mn,
                         remCIpctSig=remCIpctSig,
                         tau2.mn=tau2.mn,tau2.sd=tau2.sd,pctTau2eq0=pctTau2eq0,
                         eWthkSEobs.mn=eWthkSEobs.mn,eWthkSEobs.sd=eWthkSEobs.sd)
  }
  
  simResult
}

###############################################################
#CALCULATION OF BASIC STATISTICS ACROSS DIFFERENT WEIGHTING METHODS
###############################################################
#Calculate mean via random or fixed effects (determined by genDtype)
simMnCalc=function(simData,genDtype){
  if(genDtype=="RS"){
    
    nSim=max(simData$simIter)
    
    du.eq=apply(matrix(simData$MnDu,ncol=nSim,byrow=FALSE),2,mean)
    dh.eq=apply(matrix(simData$MnDh,ncol=nSim,byrow=FALSE),2,mean)
    
    dh.n=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                  vi=matrix(1/Wts.n,ncol=nSim,byrow=FALSE)),2,simRMA.r))    
    dh.neDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.neDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    du.xmDu=with(simData,apply(rbind(yi=matrix(MnDu,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.xmDu,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    dh.xmDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.xmDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    du.lmDu=with(simData,apply(rbind(yi=matrix(MnDu,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.lmDu,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    dh.lmDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.lmDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    du.exDu=with(simData,apply(rbind(yi=matrix(MnDu,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.exDu,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    dh.exDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.exDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    du.lsDu=with(simData,apply(rbind(yi=matrix(MnDu,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.lsDu,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    dh.lsDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.lsDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    du.wxDu=with(simData,apply(rbind(yi=matrix(MnDu,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.wxDu,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    dh.wxDh=with(simData,apply(rbind(yi=matrix(MnDh,ncol=nSim,byrow=FALSE),
                                     vi=matrix(1/Wts.wxDh,ncol=nSim,byrow=FALSE)),2,simRMA.r))
    
    list(mn=data.frame(du.eq=du.eq,dh.eq=dh.eq,dh.n=dh.n[1,],dh.neDh=dh.neDh[1,],
                       du.xmDu=du.xmDu[1,],dh.xmDh=dh.xmDh[1,],du.lmDu=du.lmDu[1,],dh.lmDh=dh.lmDh[1,],
                       du.exDu=du.exDu[1,],dh.exDh=dh.exDh[1,],du.lsDu=du.lsDu[1,],dh.lsDh=dh.lsDh[1,],
                       du.wxDu=du.wxDu[1,],dh.wxDh=dh.wxDh[1,]),
         se=data.frame(du.eq=sqrt(apply(matrix(simData$MnDu,ncol=nSim,byrow=FALSE),2,var)),
                       dh.eq=sqrt(apply(matrix(simData$MnDh,ncol=nSim,byrow=FALSE),2,var)),
                       dh.n=dh.n[2,],dh.neDh=dh.neDh[2,],
                       du.xmDu=du.xmDu[2,], dh.xmDh=dh.xmDh[2,], du.lmDu=du.lmDu[2,], dh.lmDh=dh.lmDh[2,],
                       du.exDu=du.exDu[2,], dh.exDh=dh.exDh[2,], du.lsDu=du.lsDu[2,], dh.lsDh=dh.lsDh[2,],
                       du.wxDu=du.wxDu[2,],dh.wxDh=dh.wxDh[2,]),
         estTau2=data.frame(dh.n=dh.n[3,],dh.neDh=dh.neDh[3,],
                            du.xmDu=du.xmDu[3,], dh.xmDh=dh.xmDh[3,], du.lmDu=du.lmDu[3,], dh.lmDh=dh.lmDh[3,],
                            du.exDu=du.exDu[3,], dh.exDh=dh.exDh[3,], du.lsDu=du.lsDu[3,], dh.lsDh=dh.lsDh[3,],
                            du.wxDu=du.wxDu[3,],dh.wxDh=dh.wxDh[3,]),
         ci.lb=data.frame(dh.n=dh.n[4,],dh.neDh=dh.neDh[4,],
                          du.xmDu=du.xmDu[4,], dh.xmDh=dh.xmDh[4,], du.lmDu=du.lmDu[4,], dh.lmDh=dh.lmDh[4,],
                          du.exDu=du.exDu[4,], dh.exDh=dh.exDh[4,], du.lsDu=du.lsDu[4,], dh.lsDh=dh.lsDh[4,],
                          du.wxDu=du.wxDu[4,],dh.wxDh=dh.wxDh[4,]),
         ci.ub=data.frame(dh.n=dh.n[5,],dh.neDh=dh.neDh[5,],
                          du.xmDu=du.xmDu[5,], dh.xmDh=dh.xmDh[5,], du.lmDu=du.lmDu[5,], dh.lmDh=dh.lmDh[5,],
                          du.exDu=du.exDu[5,], dh.exDh=dh.exDh[5,], du.lsDu=du.lsDu[5,], dh.lsDh=dh.lsDh[5,],
                          du.wxDu=du.wxDu[5,],dh.wxDh=dh.wxDh[5,]))
    
  } else {
    simData%>%group_by(simIter)%>%
      summarize(du.eq=calcWtMn.r(Wts.eq,MnDu),
                dh.eq=calcWtMn.r(Wts.eq,MnDh),
                dh.n=calcWtMn.r(Wts.n,MnDh),
                dh.neDh=calcWtMn.r(Wts.neDh,MnDh),
                
                du.xmDu=calcWtMn.r(Wts.xmDu,MnDu),
                dh.xmDh=calcWtMn.r(Wts.xmDh,MnDh),
                
                du.lmDu=calcWtMn.r(Wts.lmDu,MnDu),
                dh.lmDh=calcWtMn.r(Wts.lmDh,MnDh),
                
                #these assign endpoint delta estimates to weights; doesn't make sense for fixed effects.  mean delta better estimate.
                
                du.exDu=calcWtMn.r(Wts.exDu,MnDu),
                dh.exDh=calcWtMn.r(Wts.exDh,MnDh),
                
                du.lsDu=calcWtMn.r(Wts.lsDu,MnDu),
                dh.lsDh=calcWtMn.r(Wts.lsDh,MnDh),
                
                du.wxDu=calcWtMn.r(Wts.wxDu,MnDu),
                dh.wxDh=calcWtMn.r(Wts.wxDh,MnDh),.groups="drop_last")%>%select(-simIter)
  }
}
#simMnCalc(obsSimData,genDtype)
######################################
#Calculate estimate of Var(d) from simulated data
######################################
simSumWtsVarCalc=function(simData,genDtype,obsEstTau2=NULL){
  
  #Estimates of variance for each endpoint derived for each study.    
  if(genDtype=="RS"){
    cbind(simData,obsEstTau2)%>%group_by(simIter)%>%
      summarize(du.eq=calcSumWtsVar.r(Wts.eq),
                dh.eq=calcSumWtsVar.r(Wts.eq),
                dh.n=calcSumWtsVar.r(1/(tau2.dh.n+1/Wts.n)),
                dh.neDh=calcSumWtsVar.r(1/(tau2.dh.neDh+1/Wts.neDh)),
                
                du.xmDu=calcSumWtsVar.r(1/(tau2.du.xmDu+1/Wts.xmDu)),
                dh.xmDh=calcSumWtsVar.r(1/(tau2.dh.xmDh+1/Wts.xmDh)),
                
                du.lmDu=calcSumWtsVar.r(1/(tau2.du.lmDu+1/Wts.lmDu)),
                dh.lmDh=calcSumWtsVar.r(1/(tau2.dh.lmDh+1/Wts.lmDh)),
                
                #these assign endpoint delta estimates to weights; doesn't make sense for fixed effects.  mean delta better estimate.
                
                du.exDu=calcSumWtsVar.r(1/(tau2.du.exDu+1/Wts.exDu)),
                dh.exDh=calcSumWtsVar.r(1/(tau2.dh.exDh+1/Wts.exDh)),
                
                du.lsDu=calcSumWtsVar.r(1/(tau2.du.lsDu+1/Wts.lsDu)),
                dh.lsDh=calcSumWtsVar.r(1/(tau2.dh.lsDh+1/Wts.lsDh)),
                
                du.wxDu=calcSumWtsVar.r(1/(tau2.du.wxDu+1/Wts.wxDu)),
                dh.wxDh=calcSumWtsVar.r(1/(tau2.dh.wxDh+1/Wts.wxDh)),.groups="drop_last")%>%select(-simIter)  
  } else {
    simData%>%group_by(simIter)%>%
      summarize(du.eq=calcSumWtsVar.r(Wts.eq),
                dh.eq=calcSumWtsVar.r(Wts.eq),
                dh.n=calcSumWtsVar.r(Wts.n),
                dh.neDh=calcSumWtsVar.r(Wts.neDh),
                
                du.xmDu=calcSumWtsVar.r(Wts.xmDu),
                dh.xmDh=calcSumWtsVar.r(Wts.xmDh),
                
                du.lmDu=calcSumWtsVar.r(Wts.lmDu),
                dh.lmDh=calcSumWtsVar.r(Wts.lmDh),
                
                #these assign endpoint delta estimates to weights; doesn't make sense for fixed effects.  mean delta better estimate.
                
                du.exDu=calcSumWtsVar.r(Wts.exDu),
                dh.exDh=calcSumWtsVar.r(Wts.exDh),
                
                du.lsDu=calcSumWtsVar.r(Wts.lsDu),
                dh.lsDh=calcSumWtsVar.r(Wts.lsDh),
                
                du.wxDu=calcSumWtsVar.r(Wts.wxDu),
                dh.wxDh=calcSumWtsVar.r(Wts.wxDh),.groups="drop_last")%>%select(-simIter)  
    
  }
}

###############################################
simWtdSSCalc=function(simData,mnResults,genDtype,obsEstTau2=NULL){
  if(genDtype=="RS"){
    
    cbind(simData,mnResults,obsEstTau2)%>%group_by(simIter)%>%
      summarize(du.eq=sum(Wts.eq*(MnDu-du.eq)^2),
                dh.eq=sum(Wts.eq*(MnDh-dh.eq)^2),
                
                dh.n=calcWtdSS.r(MnDh,dh.n,tau2.dh.n,Wts.n),
                
                dh.neDh=calcWtdSS.r(MnDh,dh.neDh,tau2.dh.neDh,Wts.neDh),
                
                du.xmDu=calcWtdSS.r(MnDu,du.xmDu,tau2.du.xmDu,Wts.xmDu),
                dh.xmDh=calcWtdSS.r(MnDh,dh.xmDh,tau2.dh.xmDh,Wts.xmDh),
                
                du.lmDu=calcWtdSS.r(MnDu,du.lmDu,tau2.du.lmDu,Wts.lmDu),
                dh.lmDh=calcWtdSS.r(MnDh,dh.lmDh,tau2.dh.lmDh,Wts.lmDh),
                
                #These don't make sense for fixed effects as assume delta same for all endpoints,
                
                du.exDu=calcWtdSS.r(MnDu,du.exDu,tau2.du.exDu,Wts.exDu),
                dh.exDh=calcWtdSS.r(MnDh,dh.exDh,tau2.dh.exDh,Wts.exDh),
                
                du.lsDu=calcWtdSS.r(MnDu,du.lsDu,tau2.du.lsDu,Wts.lsDu),
                dh.lsDh=calcWtdSS.r(MnDh,dh.lsDh,tau2.dh.lsDh,Wts.lsDh),
                
                du.wxDu=calcWtdSS.r(MnDu,du.wxDu,tau2.du.wxDu,Wts.wxDu),
                dh.wxDh=calcWtdSS.r(MnDh,dh.wxDh,tau2.dh.wxDh,Wts.wxDh),.groups="drop_last")%>%select(-simIter)
  } else {
    cbind(simData,mnResults)%>%group_by(simIter)%>%
      summarize(du.eq=sum(Wts.eq*(MnDu-du.eq)^2),
                dh.eq=sum(Wts.eq*(MnDh-dh.eq)^2),
                dh.n=sum(Wts.n*(MnDh-dh.n)^2),
                
                dh.neDh=sum(Wts.neDh*(MnDh-dh.neDh)^2),
                
                du.xmDu=sum(Wts.xmDu*(MnDu-du.xmDu)^2),
                dh.xmDh=sum(Wts.xmDh*(MnDh-dh.xmDh)^2),
                
                du.lmDu=sum(Wts.lmDu*(MnDu-du.lmDu)^2),
                dh.lmDh=sum(Wts.lmDh*(MnDh-dh.lmDh)^2),
                
                #These don't make sense for fixed effects as assume delta same for all endpoints,
                du.exDu=sum(Wts.exDu*(MnDu-du.exDu)^2),
                dh.exDh=sum(Wts.exDh*(MnDh-dh.exDh)^2),
                
                du.lsDu=sum(Wts.lsDu*(MnDu-du.lsDu)^2),
                dh.lsDh=sum(Wts.lsDh*(MnDh-dh.lsDh)^2),
                
                du.wxDu=sum(Wts.wxDu*(MnDu-du.wxDu)^2),
                dh.wxDh=sum(Wts.wxDh*(MnDh-dh.wxDh)^2),.groups="drop_last")%>%select(-simIter)
    
  } 
}