Study6.Results.txt: Study of the conservative bootstrap confidence interval for z-curve under full heterogeneity ---------------------------------------------------------------------------------- > ################################################################################## > # AnalyzeStudy6.txt > ################################################################################## > # Performance of the zcurve bootstrap confidence interval: Full heterogeneity > > rm(list=ls()) > options(scipen=999) # To avoid scientific notation! > simdata = scan("Study6.Data.txt") Read 1932147 items > # Calculate the number of simulations. It might not be exactly 10,000. There are > # 9 values of k = Number of studies > # 3 values of mean power after selection > # 9*3 = 27 lines of data for each simulation. > nsim = length(simdata)/(27*7) > cat("\nAnalyzing data from",nsim,"simulations \n\n") Analyzing data from 10223 simulations > nrows = length(simdata)/7; ncols = 7 > dim(simdata) = c(ncols,nrows) # Reversed: R reads by columns, rows > simdata = t(simdata) > colnames(simdata) = c("Simulation", "Nstudies", "PopMeanPower", "SampleMeanPower", + "Zcurve", "LowerCL", "UpperCL") > head(simdata) ; tail(simdata) Simulation Nstudies PopMeanPower SampleMeanPower Zcurve LowerCL UpperCL [1,] 1 25 0.25 0.1976075 0.192 0.010 0.487 [2,] 1 25 0.50 0.4797951 0.501 0.137 0.751 [3,] 1 25 0.75 0.7206065 0.831 0.684 0.922 [4,] 1 50 0.25 0.2501753 0.292 0.074 0.517 [5,] 1 50 0.50 0.4535663 0.650 0.443 0.768 [6,] 1 50 0.75 0.6930733 0.659 0.458 0.812 Simulation Nstudies PopMeanPower SampleMeanPower Zcurve LowerCL UpperCL [276016,] 36 3000 0.25 0.2484679 0.266 0.183 0.300 [276017,] 36 3000 0.50 0.4998756 0.476 0.430 0.532 [276018,] 36 3000 0.75 0.7498421 0.746 0.704 0.794 [276019,] 36 5000 0.25 0.2507172 0.226 0.186 0.291 [276020,] 36 5000 0.50 0.4932987 0.478 0.439 0.530 [276021,] 36 5000 0.75 0.7486494 0.734 0.694 0.777 > attach(data.frame(simdata)) # Now columns are available as variables > > # Experiment: A conservative correction (now part of the zboot function) > # delta = 0.02 > # UpperCL = UpperCL+delta; LowerCL = LowerCL-delta > Check1 = UpperCL > PopMeanPower; Check2 = PopMeanPower > LowerCL > InCI = Check1 * Check2 # Indicator for parameter in confidence interval > # cbind(Check1,Check2,InCI)[1:100,] # Check. > CIwidth = UpperCL-LowerCL > > ######################### Sample sizes ############################# > kounts = table(PopMeanPower, Nstudies) # Rows, cols ,panels > cat("\n Sample sizes \n\n"); print(kounts) Sample sizes Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 10223 10223 10223 10223 10223 10223 10223 10223 10223 0.5 10223 10223 10223 10223 10223 10223 10223 10223 10223 0.75 10223 10223 10223 10223 10223 10223 10223 10223 10223 > > ###################### Average Performance ###################### > estframe = aggregate(Zcurve,by=list(PopMeanPower, Nstudies),FUN=mean) > meanest = estframe$x; dim(meanest) = c(3,9) # Numbers of factor levels, in order > MeanEstimatedPower = kounts # To use the nice labels > for(j in 1:3) MeanEstimatedPower[j,] = meanest[j,] > cat("\n Mean Estimated Power \n\n") Mean Estimated Power > print(round(MeanEstimatedPower,2)) Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 0.28 0.26 0.25 0.24 0.23 0.23 0.23 0.23 0.23 0.5 0.54 0.53 0.52 0.50 0.50 0.49 0.49 0.49 0.48 0.75 0.77 0.77 0.76 0.76 0.75 0.74 0.74 0.74 0.73 > > ######################### Coverage ############################# > > ciframe = aggregate(InCI,by=list(PopMeanPower, Nstudies),FUN=mean) > meanci = round(100*ciframe$x,2); dim(meanci) = c(3,9) > ConfidenceIntervalCoverage = kounts > for(j in 1:3) ConfidenceIntervalCoverage[j,] = meanci[j,] > cat("\n CI Coverage \n\n") CI Coverage > ConfidenceIntervalCoverage Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 96.22 96.82 98.31 98.58 98.89 98.62 98.02 97.94 98.07 0.5 94.81 95.71 96.56 98.41 98.91 99.37 99.39 99.25 99.33 0.75 93.56 94.59 96.89 98.58 99.59 99.60 99.69 99.70 99.42 > > > ######################### CI Width ############################# > widframe = aggregate(CIwidth,by=list(PopMeanPower, Nstudies),FUN=mean) > meanwid = widframe$x; dim(meanwid) = c(3,9) # Numbers of factor levels, in order > Mean_CI_Width = kounts # To use the nice labels > for(j in 1:3) Mean_CI_Width[j,] = meanwid[j,] > cat("\n Mean CI Width \n\n") Mean CI Width > print(round(Mean_CI_Width,3)) Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 0.480 0.377 0.289 0.206 0.164 0.136 0.116 0.108 0.099 0.5 0.497 0.392 0.306 0.221 0.175 0.141 0.117 0.107 0.097 0.75 0.335 0.258 0.204 0.155 0.129 0.109 0.093 0.086 0.078 > > ######################### Margin of Error ############################# > MeanMarginOfError = Mean_CI_Width/2 > cat("\n Mean Margin of Error \n\n") Mean Margin of Error > print(round(MeanMarginOfError,3)) Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 0.240 0.188 0.145 0.103 0.082 0.068 0.058 0.054 0.049 0.5 0.249 0.196 0.153 0.111 0.087 0.071 0.059 0.054 0.048 0.75 0.167 0.129 0.102 0.077 0.064 0.054 0.047 0.043 0.039 > > ######################### Mean Lower Limit ############################# > > lowframe = aggregate(LowerCL,by=list(PopMeanPower, Nstudies),FUN=mean) > meanlow = lowframe$x; dim(meanlow) = c(3,9) # Numbers of factor levels, in order > MeanLowerCL = kounts # To use the nice labels > for(j in 1:3) MeanLowerCL[j,] = meanlow[j,] > cat("\n Mean Lower Limit \n\n") Mean Lower Limit > print(round(MeanLowerCL,2)) Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 0.06 0.09 0.12 0.14 0.16 0.17 0.17 0.18 0.18 0.5 0.26 0.32 0.36 0.39 0.41 0.42 0.43 0.43 0.44 0.75 0.56 0.61 0.65 0.67 0.68 0.69 0.69 0.70 0.70 > > ######################### Mean Upper Limit ############################# > > highframe = aggregate(UpperCL,by=list(PopMeanPower, Nstudies),FUN=mean) > meanhigh = highframe$x; dim(meanhigh) = c(3,9) # Numbers of factor levels, in order > MeanUpperCL = kounts # To use the nice labels > for(j in 1:3) MeanUpperCL[j,] = meanhigh[j,] > cat("\n Mean Upper Limit \n\n") Mean Upper Limit > print(round(MeanUpperCL,2)) Nstudies PopMeanPower 25 50 100 250 500 1000 2000 3000 5000 0.25 0.54 0.47 0.40 0.35 0.32 0.30 0.29 0.29 0.28 0.5 0.76 0.71 0.67 0.62 0.59 0.56 0.55 0.54 0.53 0.75 0.89 0.87 0.85 0.83 0.81 0.80 0.79 0.78 0.77