#R Code behind Figure 2 in
#Simonsohn, Simmons and Nelson, "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
##########################################################
#####################################################################################################
#SIMPLIFIED LOSS FUNCTION, USED THROUGHOUT BELOW
#see appendix in paper for explanations and more robust version (robust to user input with t<0 and/or p>.05)
loss=function(t_obs,df_obs,d_est) {
ncp_est=sqrt((df_obs+2)/4)*d_est
tc=qt(.975,df_obs)
power_est=1-pt(tc,df_obs,ncp_est)
p_larger=pt(t_obs,df=df_obs,ncp=ncp_est)
ppr=(p_larger-(1-power_est))/power_est
KSD=ks.test(ppr,punif)$statistic
return(KSD) }
######################################################################################
library(metafor) #this is the meta-analysis package that contains the trim-and-fill routine
f2=function(n0,n1,ki,d_mean,d_sd) {
#SYNTAX
# n0 is the smallest sample size considered
# n1 is the largest sample size in the set
# ki is the number of p<.05 studies of a given sample size, e.g., if ki=100 and n0=20 and n1=25, there will 100 n=20 studies, 100 n=21, 100 n=22, etc.
# d_mean is the average of the true effect size being simulated, cohen-d
# d_sd is the stadnard deviation of the true effect size being simulated.
# to simulate studies from the exact same effect size we set d_sd=0
#EXPLANATION:
# This program generates Figure 2 which shows naive, trim-and-fill corrected, and p-curve based effect size estimates
# It generates significant studies with samples sizes between n=n0 and n=n1, the same number of each, where the true effect size is d~N(d_mean,d_sd)
#
#STEPS
# 1) Generate the sample sizes to be used
# 2) Draw true effect sizes from a population distribution N(d_mean,d_sd), one effect size is drawn per simulated study
# 3) Generate only significant t-tests, under the noncentrality parameter (ncp) implied by the sample sizes and true ds
# 4) Compute p-curve's estimate by minimizing loss for those t-tests
# 5) Convert ts into ds
# 6) Use metaanalysis to compute the naive d and the trim and fill correction
#SET SEED TO ALWAYS GET SAME RESULT
set.seed(150*d_mean) #set seed as a function of d_mean so that the results are always the same, but across different settings the seed varies
#1) Generate the sample sizes to be used
ni=rep(n0:n1,each=ki) #vector ni has one scalar per study, the first ki elements are n0, the next ki are n0+1....
df=2*ni-2 #degrees of freedom, is a vector
#2) true d
di=rnorm(n=ki*(n1-n0+1),mean=d_mean,sd=d_sd)
#3) Generate significant t-tests
tc=qt(.975,df) #3.1 critical t-test that is significant with the sample size
ncpi=sqrt(ni/2)*di #3.2 noncentarlity parameter used to draw from noncentral t, vector as different ncp for each n
poweri=1-pt(tc,df=df,ncp=ncpi) #3.3 power is computed to use in the next line to generate only p<.05 t-values
rp=runif(n=(n1-n0+1)*ki,min=1-poweri,max=1) #3.4 generates the randomly draw "p"-values from a uniform (p-values but computed with the noncentral
# This latter step is a bit tricky if you are not familiar with noncentral distributions
# Once we know power, we know that every t-value above the power_th percentile of them, under the ncp,
# will give a p<.05 if we evaluate it with the central t so we draw only from that range.
# Example. Say that for a given n=20 & d=5-->power =33%. To draw significant t-values under the null of d=.5, then,
# we compute the inverse of randomly drawn uniform values between .67 and 1 (the biggest third) form the qt(rt(min=.67,max=1),ncp=sqrt(20/2)*.5, df=38). .
t_publish=qt(p=rp,df=df,ncp=sqrt(ni/2)*di) #t-values associated with those uniform cdf values
#4) Find d-hat for p-curve
dp=optimize(loss,c(-.3,2),df_obs=df,t_obs=t_publish)$minimum #optimize misbehaves very far from the truth, so constrained to -.3,2
#5) #Convert t-values into d-values using formula: d=2*t/sqrt(N) -shown in main paper, see footnote
d_publish=(2*t_publish)/sqrt(df+2)
#6) Meta-analysis estimates
#6.1) Compute the variance of each d-value in the set.
#the variance of a standardized coefficient depends only on sample size and the coefficient.
#I use the formula provided in the Meta-analysis book by Hedges and Olkin (1985), "Statistical Methods for Meta-Analysis", page 80, Equation 8
#var(g=1/n~ + (d**2)/(2*(N-3.94))
# Where:
# g is what I refer to as d here, the standardized difference of means
# n~ is n1*n2/(n1+n2) and n1,n2 refer to sample sizes of each samples, if n1=n2 then n~=n**2/2*n=n/2 [sorry, ~ is not approximate, but tilde]
# N is the total n, if n1=n2 then N=2n
vard_publish= 2/ni +(d_publish**2)/(2*(ni-3.94)) #this is var(d)
#6.2) Run the standard meta-analysis routine
naive=rma(yi=d_publish,vi=vard_publish,method="FE") #once again, you need to load metafor() to run RMA
#this will run a fixed-effect meta-analysis (hence method="FE"
#using as the d.v. the d_publish vector and as the variance, its variance vard_publish)
#the result, naive, contains many values, not just the point estimate, we need that full set
#for trim-and-fill correction
#6.3) Apply the trim and fill correction to the estimate
tf=trimfill(naive)
#Report results as they are being computed
cat("# tests ",length(t_publish))
cat("\n naive estimate: ",naive$b) #the $b indicates that from the vector naive we just want the point estimate, $b
cat("\n trim and fill: ",tf$b)
cat("\n p-curve: ",dp)
cat("\n mean(di)=",mean(di)," sd(di)=",sd(di) ) #this verifies the simulation obtained the desired distribution of true-effect sizes
cat("\n =========================================== \n")
return(c(naive$b,tf$b,dp)) #Return results as vector of point estimate by the navie method of average p<.05, the trim-and-fill correction, and p-curves
}
#mat() is a function that runs the tf() function for d=[0,.8] and saves results into a matrix
mat=function(ki,n0,n1,d_sd) {
#Run the tf() for d=(0,.2,.4,.6,.8)
r1=f2(ki=ki,n0=n0,n1=n1,d_mean=.0,d_sd=d_sd)
r2=f2(ki=ki,n0=n0,n1=n1,d_mean=.2,d_sd=d_sd)
r3=f2(ki=ki,n0=n0,n1=n1,d_mean=.4,d_sd=d_sd)
r4=f2(ki=ki,n0=n0,n1=n1,d_mean=.6,d_sd=d_sd)
r5=f2(ki=ki,n0=n0,n1=n1,d_mean=.8,d_sd=d_sd)
r=c(r1,r2,r3,r4,r5)
return(matrix(r,5,3,byrow=TRUE))
}
#Fig 2A: n=20, sd(d)=0
m20=mat(ki=5000,n0=20,n1=20,d_sd=0)
#Fig 2B: n=U[5-35], sd(d)=0
m5_35=mat(ki=161,n0=5,n1=35,d_sd=0) #NOte, with 161 studies per sample-size, the total number of studies is ~5000
#Fig 2C: n=U[5-35], sd(d)=.1
m5_35sd2=mat(ki=161,n0=5,n1=35,d_sd=.2)
#I copy pasted the results into excel to plot Figure 2
#Note: R gives a lot of warnings when using the non-central student distribution. These warnings point to trivial discrepancies
# in precision way outside the range we may care about. Ignoring them is, apparently, the best thing to do.