Kuhberger Power Estimation

gistfile1.r

############DATA############
#All significant Z values: 1.9609,1.9611,1.9727,1.9859,2.0063,2.017,2.0188,2.0231,2.0316,2.0377,2.04,2.0444,2.0494,2.0579,2.062,2.0624,2.0662,2.0737,2.0802,2.0803,2.0821,2.0863,2.0895,2.095,2.0964,2.0983,2.1,2.1,2.1113,2.1219,2.124,2.1299,2.1314,2.1318,2.1444,2.1448,2.1481,2.1491,2.1551,2.1584,2.1599,2.1615,2.1622,2.163,2.1656,2.1711,2.1807,2.1836,2.1865,2.1879,2.1909,2.2068,2.2151,2.2205,2.2212,2.2237,2.2316,2.2369,2.2457,2.2572,2.2668,2.2716,2.2777,2.2892,2.2993,2.3001,2.3034,2.31,2.3171,2.3183,2.3232,2.3282,2.337,2.343,2.3455,2.3459,2.3473,2.3517,2.354,2.3543,2.36,2.3623,2.3697,2.3707,2.3767,2.3781,2.3833,2.3871,2.396,2.4221,2.43,2.4408,2.4463,2.4584,2.4714,2.4751,2.485,2.488,2.4889,2.5121,2.5141,2.5306,2.5387,2.5531,2.5698,2.5711,2.5726,2.5793,2.5803,2.5988,2.6322,2.6335,2.6396,2.6397,2.6468,2.6475,2.6664,2.6687,2.6788,2.702,2.7032,2.7212,2.7245,2.7289,2.729,2.7324,2.736,2.7371,2.7402,2.744,2.7537,2.7578,2.761,2.7775,2.7846,2.7929,2.7931,2.8054,2.8131,2.8474,2.8605,2.8657,2.8735,2.8782,2.8967,2.9005,2.9009,2.9228,2.9447,2.9468,2.9509,2.9531,2.9637,2.9776,2.9848,2.993,2.9977,3.0034,3.0133,3.0226,3.0304,3.0314,3.0359,3.0382,3.047,3.051,3.06,3.0681,3.1013,3.116,3.1176,3.1185,3.1248,3.1366,3.1585,3.18,3.2044,3.2151,3.2169,3.2265,3.2309,3.2388,3.2771,3.2857,3.3115,3.3149,3.3366,3.3478,3.3487,3.3555,3.3556,3.3566,3.4205,3.4205,3.4323,3.4349,3.4446,3.4472,3.4475,3.4627,3.4669,3.4673,3.4686,3.4944,3.4971,3.5251,3.5486,3.5487,3.5575,3.5713,3.5879,3.62,3.6453,3.6825,3.6851,3.6893,3.6949,3.7174,3.732,3.7378,3.7504,3.7613,3.7688,3.7836,3.822,3.8268,3.8297,3.8435,3.8658,3.8963,3.898,3.909,3.9165,3.9307,3.9321,3.9668,3.9687,3.9797,4.0175,4.0336,4.0451,4.0453,4.0788,4.1012,4.1099,4.1415,4.1483,4.1653,4.169,4.2156,4.2702,4.3392,4.3676,4.3715,4.3962,4.4001,4.454,4.4683,4.472,4.4976,4.5094,4.5243,4.5336,4.5352,4.5457,4.57,4.5881,4.6274,4.6429,4.7121,4.7251,4.7768,4.7975,4.87,4.8916,4.8916,4.979,5.0441,5.1213,5.1277,5.197,5.2666,5.276,5.298,5.3267,5.3267,5.3267,5.3267,5.3783,5.3896,5.4995,5.5384,5.7307,5.7307,5.7307,5.8245,6.1,6.1294,6.1648,6.1694,6.245,6.2517,6.2778,6.35,6.3613,6.439,6.4635,6.6893,6.7058,6.7446,6.7483,6.7619,6.9085,6.9676,6.9761,7.3755,7.4084,7.4861,7.5288,7.6026,7.9309,7.9887,8.064,8.7197,8.755,8.7832,9.0734,9.8566,9.8709,9.9599,9.9836,10.4862,10.6835,10.7068,10.9522,12.0523,13.6545,13.9459,14.1006,15.4901,16.7672,24.2452,25.4592,30.14
#All Z values: 0.0013,0.0114,0.0526,0.0682,0.0892,0.1395,0.1732,0.1969,0.2701,0.3753,0.4115,0.4339,0.4437,0.5194,0.63,0.6369,0.6492,0.6735,0.7858,0.8213,0.835,0.8363,0.8861,0.8894,0.8983,0.9016,0.9451,0.9565,0.9741,1.0246,1.0399,1.1222,1.1391,1.2222,1.2356,1.2987,1.3259,1.3914,1.5043,1.5295,1.5844,1.6145,1.6598,1.6751,1.684,1.6848,1.7013,1.7114,1.7278,1.7517,1.7669,1.7679,1.7723,1.808,1.8119,1.8292,1.8458,1.847,1.8493,1.906,1.9164,1.924,1.9435,1.9436,1.944,1.9582,1.9609,1.9611,1.9727,1.9859,2.0063,2.017,2.0188,2.0231,2.0316,2.0377,2.04,2.0444,2.0494,2.0579,2.062,2.0624,2.0662,2.0737,2.0802,2.0803,2.0821,2.0863,2.0895,2.095,2.0964,2.0983,2.1,2.1,2.1113,2.1219,2.124,2.1299,2.1314,2.1318,2.1444,2.1448,2.1481,2.1491,2.1551,2.1584,2.1599,2.1615,2.1622,2.163,2.1656,2.1711,2.1807,2.1836,2.1865,2.1879,2.1909,2.2068,2.2151,2.2205,2.2212,2.2237,2.2316,2.2369,2.2457,2.2572,2.2668,2.2716,2.2777,2.2892,2.2993,2.3001,2.3034,2.31,2.3171,2.3183,2.3232,2.3282,2.337,2.343,2.3455,2.3459,2.3473,2.3517,2.354,2.3543,2.36,2.3623,2.3697,2.3707,2.3767,2.3781,2.3833,2.3871,2.396,2.4221,2.43,2.4408,2.4463,2.4584,2.4714,2.4751,2.485,2.488,2.4889,2.5121,2.5141,2.5306,2.5387,2.5531,2.5698,2.5711,2.5726,2.5793,2.5803,2.5988,2.6322,2.6335,2.6396,2.6397,2.6468,2.6475,2.6664,2.6687,2.6788,2.702,2.7032,2.7212,2.7245,2.7289,2.729,2.7324,2.736,2.7371,2.7402,2.744,2.7537,2.7578,2.761,2.7775,2.7846,2.7929,2.7931,2.8054,2.8131,2.8474,2.8605,2.8657,2.8735,2.8782,2.8967,2.9005,2.9009,2.9228,2.9447,2.9468,2.9509,2.9531,2.9637,2.9776,2.9848,2.993,2.9977,3.0034,3.0133,3.0226,3.0304,3.0314,3.0359,3.0382,3.047,3.051,3.06,3.0681,3.1013,3.116,3.1176,3.1185,3.1248,3.1366,3.1585,3.18,3.2044,3.2151,3.2169,3.2265,3.2309,3.2388,3.2771,3.2857,3.3115,3.3149,3.3366,3.3478,3.3487,3.3555,3.3556,3.3566,3.4205,3.4205,3.4323,3.4349,3.4446,3.4472,3.4475,3.4627,3.4669,3.4673,3.4686,3.4944,3.4971,3.5251,3.5486,3.5487,3.5575,3.5713,3.5879,3.62,3.6453,3.6825,3.6851,3.6893,3.6949,3.7174,3.732,3.7378,3.7504,3.7613,3.7688,3.7836,3.822,3.8268,3.8297,3.8435,3.8658,3.8963,3.898,3.909,3.9165,3.9307,3.9321,3.9668,3.9687,3.9797,4.0175,4.0336,4.0451,4.0453,4.0788,4.1012,4.1099,4.1415,4.1483,4.1653,4.169,4.2156,4.2702,4.3392,4.3676,4.3715,4.3962,4.4001,4.454,4.4683,4.472,4.4976,4.5094,4.5243,4.5336,4.5352,4.5457,4.57,4.5881,4.6274,4.6429,4.7121,4.7251,4.7768,4.7975,4.87,4.8916,4.8916,4.979,5.0441,5.1213,5.1277,5.197,5.2666,5.276,5.298,5.3267,5.3267,5.3267,5.3267,5.3783,5.3896,5.4995,5.5384,5.7307,5.7307,5.7307,5.8245,6.1,6.1294,6.1648,6.1694,6.245,6.2517,6.2778,6.35,6.3613,6.439,6.4635,6.6893,6.7058,6.7446,6.7483,6.7619,6.9085,6.9676,6.9761,7.3755,7.4084,7.4861,7.5288,7.6026,7.9309,7.9887,8.064,8.7197,8.755,8.7832,9.0734,9.8566,9.8709,9.9599,9.9836,10.4862,10.6835,10.7068,10.9522,12.0523,13.6545,13.9459,14.1006,15.4901,16.7672,24.2452,25.4592,30.14


#R Code to compute average power 
#Simonsohn, Nelson and Simmons,  "P-Curve and Effect Size: Correcting for Publication Bias Using Only Significant Results"
#
#Prepared by Uri Simonsohn, uws@wharton.upenn.edu

#Last updated: 2014 05 13


#OVERVIEW
#
# (STEP 1) CREATE FUNCTIONS THAT FIND THE NONCENTRALITY PARAMETER FOR THE t,F,N,chisq, associated with a given power, for the observed d.f.
# (STEP 2) CREATE PP-VALUES FOR EACH OF THE FOUR DISTRIBUTIONS FOR HOW WELL A GIVEN POWER_EST FITS 
# (STEP 3) STACK-UP ALL THE PP-VALUES INTO A VECTOR
# (STEP 4) FIND AVERAGE POWER THAT BEST FITS

####################################################################################################################################
# (STEP 1) CREATE FUNCTIONS THAT FIND THE NONCENTRALITY PARAMETER FOR THE t,F,N,chisq, associated with a given power, for the observed d.f.
#
#
#
#SET OF FUNCTIONS 1. 
#COMPUTE GAP BETWEEN POWER AND DESIRED POWER FOR A GIVEN NCP 
# (minimize these in the next step to solve for the ncp that gives the desired power)
ncp_error.t = function(delta, power, x, df)      pt(x, df = df, ncp = delta) - (1-power)   #if this equals 0, we found the ncp.
ncp_error.f = function(delta, power, x, df1,df2) pf(x, df1 = df1, df2=df2, ncp = delta) - (1-power)   
ncp_error.c = function(delta, power, x, df)      pchisq(x, df = df, ncp = delta) - (1-power)   
ncp_error.z = function(delta, power, x)          pnorm(x, mean = delta,sd=1) - (1-power)   

#SET OF FUNCTIONS 2: MINIMIZE FUNCTIONS ABOVE
#t-test
getncp.t =function(df, power)   {      
  xc=qt(p=.975, df=df) # critical t-value
  return(uniroot(ncp_error.t, c(0, 37.62), x = xc, df =df, power=power)$root)   }  

#F-test
getncp.f =function(df1,df2, power)   {      
  xc=qf(p=.95, df1=df1,df2=df2) # critical F-value
  return(uniroot(ncp_error.f, c(0, 37.62), x = xc, df1 = df1,df2=df2, power=power)$root)  }


#chisq-test
getncp.c =function(df, power)   {      
  xc=qchisq(p=.95, df=df) # critical c-value
  return(uniroot(ncp_error.c, c(0, 37.62), x = xc, df = df, power=power)$root)   }

#Normal
getncp.z =function(power)   {      
  xc=qnorm(p=.975) # critical Z-value with df=1
  return(uniroot(ncp_error.z, c(0, 37.62), x = xc, power=power)$root)   }             


####################################################################################################################################
# (STEP 2) CREATE PP-VALUES FOR EACH OF THE FOUR DISTRIBUTIONS FOR HOW WELL A GIVEN POWER_EST FITS 
powerfit.t=function(t_obs, df_obs, power_est)    {
  ncp_est=mapply(getncp.t,df=df_obs,power=power_est)  #find ncp for each  that gives each test power.k
  p_larger=pt(t_obs,df=df_obs,ncp=ncp_est)            #prob t>tobs given ncp_est
  ppr=(p_larger-(1-power_est))/power_est              #condition on p<.05
  return(ppr)   }

powerfit.f=function(f_obs, df1_obs, df2_obs, power_est)    {
  ncp_est=mapply(getncp.f,df1=df1_obs, df2=df2_obs,power=power_est)  #find ncp for each  that gives each test power.k
  p_larger=pf(f_obs,df1=df1_obs,df2=df2_obs, ncp=ncp_est)        #prob t>tobs given ncp_est
  ppr=(p_larger-(1-power_est))/power_est          #condition on p<.05
  return(ppr)   }

powerfit.z=function(z_obs, power_est)    {
  ncp_est=mapply(getncp.z,power=power_est)  
  p_larger=pnorm(z_obs,mean=ncp_est)        
  ppr=(p_larger-(1-power_est))/power_est          
  return(ppr)     }


powerfit.c=function(c_obs, df_obs, power_est)    {
  ncp_est=mapply(getncp.c,df=df_obs,power=power_est)  
  p_larger=pchisq(c_obs,df=df_obs,ncp=ncp_est)        
  ppr=(p_larger-(1-power_est))/power_est          
  return(ppr)   }

####################################################################################################################################
# (STEP 3) STACK-UP ALL THE PP-VALUES INTO A VECTOR

powerfit.all=function(power_est)
{
  ppr.all=c()
  #for each kind of test, check if there are any significant values, if there are, add ppr to overall ppr
  if (length(t.value.sig)>0) ppr.all=c(ppr.all, powerfit.t(t_obs=t.value.sig, df_obs=t.df.sig, power_est=power_est))
  if (length(f.value.sig)>0) ppr.all=c(ppr.all, powerfit.f(f_obs=f.value.sig, df1_obs=f.df1.sig, df2_obs=f.df2.sig, power_est=power_est))
  if (length(z.value.sig)>0) ppr.all=c(ppr.all, powerfit.z(z_obs=z.value.sig, power_est=power_est))
  if (length(c.value.sig)>0) ppr.all=c(ppr.all, powerfit.c(c_obs=c.value.sig, df_obs=c.df.sig, power_est=power_est))
  
  #outlier robust pprs  - WINSORIZE AT 1%
  #ppr.all=pmax(.01,ppr.all)
  #ppr.all=pmin(.99,ppr.all)
  KSD=ks.test(ppr.all,punif)$statistic                #KS test on the resulting pprs
  return(KSD)
}

####################################################################################################################################
# (STEP 4) PLOT FIT OF EACH POSSIBLE POWER VALUE 6%-99%

#FITTING POWER WITH P-CURVE
plotfit=function()
{
  # Fit will be evaluated at every possible value of power between 6% and 99% in steps of 1%, stored in fit()
  fit=c()
  # Evalaute fit for every value
  for (i in 6:99)   fit=c(fit,powerfit.all(i/100))
  # Find the minimum
  mini=match(min(fit),fit)       #which ith power level considered leads to best estimate
  hat=(5+mini)/100               #convert that into the power level, the ith value considered is (6+ith)/1000
  #Plot results
  #create the x-axis
  x.power=seq(from=6,to=99)/100 
  #Draw teh line
  plot(x.power,fit,xlab="Underlying Power", ylab="Loss (D stat in KS test)",ylim=c(min(fit)-.05,max(fit)+.05), 
       main="How well does each power level fit? \n(lower is better)")  
  #Make red dot at the estimate
  points(hat,min(fit),pch=19,col="red",cex=2)    
  #Put a label with the estimate value
  text(max(.15,hat),min(fit)-.04,paste0("P-curve's Estimate: Average Power=",hat*100,"%"),col="red")
}


#ENTERING THE
#t-tests
t.value.sig=c()
t.df.sig=c()
#Ftest
f.value.sig=c()
f.df1.sig=c()
f.df2.sig=c()
#Normal
z.value.sig=c(1.9609,1.9611,1.9727,1.9859,2.0063,2.017,2.0188,2.0231,2.0316,2.0377,2.04,2.0444,2.0494,2.0579,2.062,2.0624,2.0662,2.0737,2.0802,2.0803,2.0821,2.0863,2.0895,2.095,2.0964,2.0983,2.1,2.1,2.1113,2.1219,2.124,2.1299,2.1314,2.1318,2.1444,2.1448,2.1481,2.1491,2.1551,2.1584,2.1599,2.1615,2.1622,2.163,2.1656,2.1711,2.1807,2.1836,2.1865,2.1879,2.1909,2.2068,2.2151,2.2205,2.2212,2.2237,2.2316,2.2369,2.2457,2.2572,2.2668,2.2716,2.2777,2.2892,2.2993,2.3001,2.3034,2.31,2.3171,2.3183,2.3232,2.3282,2.337,2.343,2.3455,2.3459,2.3473,2.3517,2.354,2.3543,2.36,2.3623,2.3697,2.3707,2.3767,2.3781,2.3833,2.3871,2.396,2.4221,2.43,2.4408,2.4463,2.4584,2.4714,2.4751,2.485,2.488,2.4889,2.5121,2.5141,2.5306,2.5387,2.5531,2.5698,2.5711,2.5726,2.5793,2.5803,2.5988,2.6322,2.6335,2.6396,2.6397,2.6468,2.6475,2.6664,2.6687,2.6788,2.702,2.7032,2.7212,2.7245,2.7289,2.729,2.7324,2.736,2.7371,2.7402,2.744,2.7537,2.7578,2.761,2.7775,2.7846,2.7929,2.7931,2.8054,2.8131,2.8474,2.8605,2.8657,2.8735,2.8782,2.8967,2.9005,2.9009,2.9228,2.9447,2.9468,2.9509,2.9531,2.9637,2.9776,2.9848,2.993,2.9977,3.0034,3.0133,3.0226,3.0304,3.0314,3.0359,3.0382,3.047,3.051,3.06,3.0681,3.1013,3.116,3.1176,3.1185,3.1248,3.1366,3.1585,3.18,3.2044,3.2151,3.2169,3.2265,3.2309,3.2388,3.2771,3.2857,3.3115,3.3149,3.3366,3.3478,3.3487,3.3555,3.3556,3.3566,3.4205,3.4205,3.4323,3.4349,3.4446,3.4472,3.4475,3.4627,3.4669,3.4673,3.4686,3.4944,3.4971,3.5251,3.5486,3.5487,3.5575,3.5713,3.5879,3.62,3.6453,3.6825,3.6851,3.6893,3.6949,3.7174,3.732,3.7378,3.7504,3.7613,3.7688,3.7836,3.822,3.8268,3.8297,3.8435,3.8658,3.8963,3.898,3.909,3.9165,3.9307,3.9321,3.9668,3.9687,3.9797,4.0175,4.0336,4.0451,4.0453,4.0788,4.1012,4.1099,4.1415,4.1483,4.1653,4.169,4.2156,4.2702,4.3392,4.3676,4.3715,4.3962,4.4001,4.454,4.4683,4.472,4.4976,4.5094,4.5243,4.5336,4.5352,4.5457,4.57,4.5881,4.6274,4.6429,4.7121,4.7251,4.7768,4.7975,4.87,4.8916,4.8916,4.979,5.0441,5.1213,5.1277,5.197,5.2666,5.276,5.298,5.3267,5.3267,5.3267,5.3267,5.3783,5.3896,5.4995,5.5384,5.7307,5.7307,5.7307,5.8245,6.1,6.1294,6.1648,6.1694,6.245,6.2517,6.2778,6.35,6.3613,6.439,6.4635,6.6893,6.7058,6.7446,6.7483,6.7619,6.9085,6.9676,6.9761,7.3755,7.4084,7.4861,7.5288,7.6026,7.9309,7.9887,8.064,8.7197,8.755,8.7832,9.0734,9.8566,9.8709,9.9599,9.9836,10.4862,10.6835,10.7068,10.9522,12.0523,13.6545,13.9459,14.1006,15.4901,16.7672,24.2452,25.4592,30.14)
#Chisq
c.value.sig=c()
c.df.sig=c()


#FIND BETS-FITTING POWER VIA P-CURVE
plotfit()                                         #plot how well each level of power fits (in 1% increments)
optimize(powerfit.all,c(.06,.999))$minimum        #estimate best fitting value