################################# #load libraries, functions, and data needed for following code. ################################# source("librariesAndFunctionsForMEE.R") modFnFCd=read.csv("modFnFCd.csv") #modFnFCd.wtsKnown=read.csv("modFnFCd_wtsKnown.csv") modFnRSd=read.csv("modFnRSd.csv") #modFnRSd2=read.csv("modFnRSd2.csv") #modFnRSd2.wtsKnown=read.csv("modFnRSd2_wtsKnown.csv") #modFnRSd2.wtsUnKnown=read.csv("modFnRSd2_wtsUnKnown.csv") #Not needed for Figures. modnbRSd=read.csv("modnbRSd.csv") repData=read.csv("repData.csv") #sigmaRatioData=read.csv("sigmaRatioData.csv") ######################### #Data Generation ######################### ###THEORETICAL RESULTS for estimating effect size and Var(d). #Data provided in "SigmaRatioData.txt". sigmaRatioData=data.frame(tau2=rep(c(0.0,0.05,0.1,0.25,0.5,0.75,1),rep(600,7)), effectSize=rep(rep((0:5)/2.5,rep(100,6)),7), n=rep(3:102,42))%>%mutate(varEX=exactVarDu.r(n,n,effectSize), varLS=approxVarDu.r(n,n,effectSize), varNE=noEffSizeVarDu.r(n,n), Je=exactJ.r(2*(n-1)), dlb=effectSize-2*sqrt(tau2+Je^2*exactVarDu.r(n,n,effectSize)), dub=effectSize+2*sqrt(tau2+Je^2*exactVarDu.r(n,n,effectSize)), varEXlb=exactVarDu.r(n,n,dlb), varEXub=exactVarDu.r(n,n,dub), wlb=1/(tau2+Je^2*varEXlb), wub=1/(tau2+Je^2*varEXub), aN=Je^2*varNE, bNd2=(2*(n-1)/(2*(n-1) - 2) - 1)*effectSize^2, w=1/(tau2+Je^2*varEX), ICC=tau2*w, ICClb=tau2*wlb, ICCub=tau2*wub, invJe=(1-Je), JeRatio=(Je/(1-Je)), wtJratio=varNE*JeRatio, Je2Ratio=Je^2/(1-Je), wtJ2ratio=varNE*Je2Ratio, effectSize=factor(effectSize)) theoryDesignData=expand.grid(tau2=c(0,0.05,0.1,0.25,0.5,0.75,1), effectSize=(0:5)/2.5, p=c(0,0.05,0.1,0.25,0.5,0.75,.95,1), n=3:50, kN=1:6) ciTolRange=c(lb=qbinom(0.025,10000,0.95),ub=qbinom(.975,10000,0.95)) ####################################### #Negative Binomial Method Distribution of weights ###################################################### nbSim=10000 nbSimParamList=list(popnICCVec=c(0,1/10,1/3,0.5,2/3,0.9), #within study Var relative to popnVar=1. #nbSimParamList=list(tau2=c(0.0,0.05,0.1,0.25,0.5,0.75,1), effectSizeVec=c(0,0.2,0.5,0.8,1.2)) nbNvec=c(4,9,16,25) nbParamDesign=expand.grid(popnICC=nbSimParamList$popnICCVec,effectSize=nbSimParamList$effectSizeVec) #nbParamDesign=expand.grid(tau2C=nbSimParamList$tau2,effectSize=nbSimParamList$effectSizeVec) nbParamDesign=nbParamDesign%>%mutate(tau2=popnICC/(1-popnICC)) tmpN1=data.frame() gdN1=data.frame() for(i in 1:4) { tmpN1=rbind(tmpN1,data.frame(mnN=nbNvec[i],N=rnbinom(nbSim,size=1.55,mu=nbNvec[i]))) #mu= 8 to 9 is considered average. See Hillebrand & Gurevitch (2014). gdN1=rbind(gdN1,data.frame(mnN=nbNvec[i],N=nbNvec[i]*(1:8))) } tmpN1=tmpN1%>%mutate(simN1 = ifelse(N<3,3,N),simN2 = simN1) gdN1=gdN1%>%mutate(simN1 = ifelse(N<3,3,N),simN2 = simN1) nbWtsData=data.frame() gdWtsData=data.frame() for(i in 1:dim(nbParamDesign)[1]){ nbWts.exDu=1/(nbParamDesign$tau2[i]+exactVarDu.r(tmpN1$simN1,tmpN1$simN2,nbParamDesign$effectSize[i])) nbWts.exDh=1/(nbParamDesign$tau2[i]+exactJ.r(tmpN1$simN1+tmpN1$simN2-2)^2*exactVarDu.r(tmpN1$simN1,tmpN1$simN2,nbParamDesign$effectSize[i])) nbWtsData=rbind(nbWtsData,data.frame(tau2=nbParamDesign$tau2[i],icc=nbParamDesign$popnICC[i],effectSize=nbParamDesign$effectSize[i], mnN=tmpN1$mnN,N=tmpN1$simN1,nbWts.exDu=nbWts.exDu,nbWts.exDh=nbWts.exDh)) gdWts.exDu=1/(nbParamDesign$tau2[i]+exactVarDu.r(gdN1$simN1,gdN1$simN2,nbParamDesign$effectSize[i])) gdWts.exDh=1/(nbParamDesign$tau2[i]+exactJ.r(gdN1$simN1+gdN1$simN2-2)^2*exactVarDu.r(gdN1$simN1,gdN1$simN2,nbParamDesign$effectSize[i])) gdWtsData=rbind(gdWtsData,data.frame(tau2=nbParamDesign$tau2[i],icc=nbParamDesign$popnICC[i],effectSize=nbParamDesign$effectSize[i], mnN=gdN1$mnN,N=gdN1$simN1,gdWts.exDu=gdWts.exDu,gdWts.exDh=gdWts.exDh)) } ################################################## #RMSE Summary calculations ################################################## relRmseData=full_join(modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse), modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>% filter(method=="dh.eq"), by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>% mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2), relRMSE=rmse/rmse.dh.eq, ICC0=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10"))) rmseSumData= relRmseData%>%group_by(method,ICC0,N)%>% summarize(RMSE.mn=sqrt(mean(relRMSE^2)), RMSE.min=min(relRMSE), RMSE.max=max(relRMSE),.groups = 'drop') relRmseFeData=full_join(modFnFCd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse), modFnFCd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>% filter(method=="dh.eq"), by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>% mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2), relRMSE=rmse/rmse.dh.eq) #NegBinom Data relRmseNbData=full_join(modnbRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse), modnbRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>% filter(method=="dh.eq"), by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dh.eq"))%>% mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2), relRMSE=rmse/rmse.dh.eq, ICC0=factor(popnICC,labels=c("0","1/10","1/3","0.5","2/3","9/10"))) rmseSumNbData= relRmseNbData%>%group_by(method,ICC0,N)%>% summarize(RMSE.mn=sqrt(mean(relRMSE^2)), RMSE.min=min(relRMSE), RMSE.max=max(relRMSE),.groups = 'drop') ###################################### #% Bias Summary calculations ###################################### maxPctBias=unlist(modFnRSd%>%filter(effectSize>0)%>%summarize(max(abs(bias/effectSize)))) maxPctBias.duLS=unlist(modFnRSd%>%filter(effectSize>0&method=="du.lsDu")%>%summarize(max(abs(bias/effectSize)))) pctNoteBias=prop.table(with(modFnRSd%>%filter(effectSize>0),table(method,abs(bias/effectSize)> 0.05)),1) pctNoteBiasByN=prop.table(with(modFnRSd%>%filter(effectSize>0),table(method,N,abs(bias/effectSize)> 0.05)),c(1,3))[,,2] sumBiasByN=modFnRSd%>% filter(effectSize>0)%>%group_by(method,N)%>% summarise(mnBias=mean(bias/effectSize), medBias=median(bias/effectSize), sdBias=sqrt(var(bias/effectSize)), maxBias=max(abs(bias)/effectSize),.groups = 'drop') ###By M biasVsVarSumData=full_join(modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse), modFnRSd%>%select(method,N,kN,Sn,M,popnICC,effectSize,bias,rmse)%>% filter(method=="dh.eq")%>% select(N,kN,Sn,M,popnICC,effectSize,rmse), by=c("N","kN","Sn","M","popnICC","effectSize"),suffix=c("",".dHeq"))%>% mutate(avgN=N*(1+kN)/2,bN2=bias^2/(rmse^2-bias^2), relRMSE=rmse/rmse.dHeq)%>% group_by(method,M)%>%filter(effectSize>0)%>% summarize(RMSE.mn=sqrt(mean(relRMSE^2)), bias.mn=100*mean(abs(bias/effectSize)), bN.mn=sqrt(mean(bN2)), RMSE.max=sqrt(max(relRMSE^2)), #RMSE.max=ifelse(abs(min(log(relRMSE^2)))>abs(max(log(relRMSE^2))),sqrt(min(relRMSE^2)),sqrt(max(relRMSE^2))), bias.max=100*max(abs(bias/effectSize)), bN.max=sqrt(max(bN2)),.groups = 'drop') ################################# #VARIOUS SUMMARIES ################################# #Mean Bias mnBiasResult=sumBiasByN%>%dplyr::filter(method=="du.lsDu"&N==25)%>%select(mnBias) #Not in data list for MEE #Notable differences between observed q and 100(1-alpha) infPctNoteBias=prop.table(with(modFnRSd,table(method,abs(pctSig-0.95)> 0.025)),1) infPctNoteBiasByN=prop.table(with(modFnRSd,table(method,N,abs(pctSig-0.95)> 0.025)),c(1,3))[,,2] #Notable improvement in relative RMSE noteRelRmse=prop.table(with(relRmseData,table(method,abs(relRMSE-1)> 0.1)),1) noteRelRmseByN=prop.table(with(relRmseData,table(method,kN,abs(relRMSE-1)> 0.1)),c(1,3))[,,2] noteRelRmse5=prop.table(with(relRmseData,table(method,abs(relRMSE-1)> 0.05)),1) #Neg Binom results noteNbRelNbRmse=prop.table(with(relRmseNbData,table(method,abs(relRMSE-1)> 0.1)),1) noteNbRelNbRmseByN=prop.table(with(relRmseNbData,table(method,kN,abs(relRMSE-1)> 0.1)),c(1,3))[,,2] noteRelNbRmse5=prop.table(with(relRmseNbData,table(method,abs(relRMSE-1)> 0.05)),1) #Mean bias of effect size for random (> 0) and fixed (==0) effects mnBiasTbl=with(modFnRSd%>%filter(effectSize>0),tapply(bias/effectSize,list(method,N),mean)) mnBias0Tbl=with(modFnRSd%>%filter(effectSize==0),tapply(bias,list(method,N),mean)) #Not in data list for MEE #Notable bias mnInfBiasTbl=with(modFnRSd%>%filter(effectSize>0),tapply(abs(pctSig-0.95)> 0.025,list(method,N,M),mean)) #Not in data list for MEE #Mean RMSE mnRMSEtbl=with(modFnRSd,sqrt(tapply(rmse^2,list(method,N),mean))) #Not in data list for MEE ############################################ #Calculate CI for variance estimate based on uncertainty in delta. ############################################ varCI.exData=theoryDesignData%>%mutate(n1=n,n2=n)%>%select(n1,n2,effectSize)%>%distinct()%>%rowwise%>% mutate(varLB=varCI.r(n1,n2,effectSize,method="EX")[1], varMed=varCI.r(n1,n2,effectSize,method="EX")[2], varUB=varCI.r(n1,n2,effectSize,method="EX")[3]) varCI.lsData=theoryDesignData%>%mutate(n1=n,n2=n)%>%select(n1,n2,effectSize)%>%distinct()%>%rowwise%>% mutate(varLB=varCI.r(n1,n2,effectSize,method="LS")[1], varMed=varCI.r(n1,n2,effectSize,method="LS")[2], varUB=varCI.r(n1,n2,effectSize,method="LS")[3]) ###################################### #Replicate data for Bland-Altman analysis. #This data is provided in repData.csv. Code shows how it was created. #repData=full_join(gather(modFnRSd2%>%select(popnICC,M,N,kN,Sn,method,effectSize,bias,pctSig,simSD), # "metric","value",bias,pctSig,simSD), # gather(modFnRSd3%>%select(popnICC,M,N,kN,Sn,method,effectSize,bias,pctSig,simSD), # "metric","value",bias,pctSig,simSD), # by=c("popnICC","M","N","kN","Sn","effectSize","method","metric"),suffix=c("",".rep")) # repSD=repData%>%group_by(metric)%>%summarize(sqrt(var(value-value.rep)/2),.groups = 'drop')