
#R Code behind Supplement 3 
#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 09 04
##########################################################
n.small=10
n.big=30
d.small=0
d.big=0
tf_share=0
tf_pmax=1
tf_pmin=0
k.small=10000
k.big=10000

#This program simulates is a modification of the R Code for supplement 2.
#It allows simulating studies with two underlying effect sizes and sample size.
#It also allows adding to the meta-analysis a pre-specified subset of n.s. findings


library(metafor)  #this is the meta-analysis package that contains the trim-and-fill routine.

  #Function tf_sees_ns allows simulations trim-and-fill observing non-significant findings.

    tf_sees_ns=function(n.small, n.big, k.small, k.big, d.small, d.big, tf_share, tf_pmax, tf_pmin) 
      #Syntax
      #n.small: sample size of set of studies with small sample size. In the paper we used n.small=10
      #n.big:   sample size of set of studies with big sample size. In the paper we used n.big=30
      #d.small: underlying true effect size of studies with a small sample size, in the paper we used d.small=0
      #d.small: underlying true effect size of studies with a big sample size, in the paper we used d.big=0
      #tf_share: number of n.s. studies that are observed, expressed as a % of the number of 
      #          significant studies that are observed. In the paper we varied between 0 and 100%
      #tf_pmax: highest n.s. p-value that's observable. we ended up not using it. It allows considering a case where only '
      #         p-values below this number are observable. For example, fi only p<.25 would be shared with meta-analyst
      #         we set it to 1 so that all p-values are observable
      #tf_pmin:  analgous with smallest p-value (imagine skeptical replicators only sharing if their p>.2). We also did not use
      #         and set it at .025 throughout (recall these are one sided p-values for a two-sided test, to p>.025 is n.s.)
    {
  
    set.seed(1975+d.small*1000+length(n.small))    #seed in a way that same paramters for simulation lead to same result but different
                                                   #parameters use other random numbers, to avoid having the same draws influence 
                                                   #all simulations.
  
  #1) Generate the sample sizes to be used
        ni.small=rep(n.small,each=(k.small/length(n.small)))    #sample size of the small studies
        ni.big=rep(n.big,each=(k.big/length(n.big)))            #sample size of the big studies
        di=rep(c(d.small,d.big),c(k.small,k.big))               #effect size of the underlying studies
        ni=c(ni.small,ni.big)                                   #merge sample sizes
        df_all=2*ni-2                                              #express as d.f. in a t-test
              
  #2) Generate  t-tests
     t_all=rt(n=(k.big+k.small),df=df_all,ncp=sqrt(ni/2)*di)       #draw t-values under null
     pi=1-pt(t_all,df=df_all)                                      #p-value of this simulated study (one sided)
     see_ns=ifelse(pi>.025 & pi>=tf_pmin & pi<=tf_pmax,1,0)     #binary 1/0 variable for whether a study is observable 
    
  
  #P<.05
    t_sig=subset(t_all,pi<.025)[1:500]                          #selects 500 p<.05 studies
    df_sig=subset(df_all,pi<.025)[1:500]                        #and their d.f.
    
  #n.s.
    t_ns0=subset(t_all,pi>=tf_pmin & pi<=tf_pmax)               #selects the corresponding share of n.s. studies
    df_ns0=subset(df_all,pi>=tf_pmin & pi<=tf_pmax)
  
    #how many nns
      nns=round(tf_share*length(t_sig)/2)                      #scalar use to select the subset of n.s. studies
    #keep only those  
        t_ns= c(t_ns0[1:nns],  t_ns0[(length(t_ns0)-nns+1):(length(t_ns0))])   #keep nns/2 of the small studies, keep nns/2 of the large studies
        df_ns=c(df_ns0[1:nns],df_ns0[(length(t_ns0)-nns+1):(length(t_ns0))])   #and their d.f.
  
    
        if (tf_share==0)                       #of 0 n.s. studies are observable, create an empty vector
            { 
            t_ns=c()
            df_ns=c()
        }
    
  
  
  #3) Combined set of p<.05 and n.s. studies
     t_obs=c(t_sig,t_ns)
     df_obs=c(df_sig,df_ns)
     n_obs=c((df_sig+2)/2,(df_ns+2)/2)
     d_obs=(2*t_obs)/sqrt(df_obs+2)
     vard_obs= 2/n_obs +(d_obs**2)/(2*(n_obs-3.94))  
  
  #4) Metaanalysis estimates
  
  #4.1) Run the standard meta-analysis routine                                    
      naive=rma(yi=d_obs,vi=vard_obs,method="FE")
      
  #4.2) Apply the trim and fill correction to the estimate
      tf=trimfill(naive)
      return(c(naive$b,tf$b,length(t_sig),length(t_ns)))    #output the two point estimates and the underlying # of studies
  #Report results as they are being computed
      cat("\n# p<.05 studies       ",length(t_sig))
      cat("\n# observable studies  ",length(t_ns))
      cat("\n naive estimate:   ",naive$b)                  
      cat("\n trim and fill:    ",tf$b)
  }


  #We run 10000 studies to get about 250 p<.05 with n=10, and about 250 with n=30. Start with 0 n.s. studies, then add 50 at a time
      r1=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=0, tf_pmax=1, tf_pmin=.025) 
      r2=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.1, tf_pmax=1, tf_pmin=.025) 
      r3=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.2, tf_pmax=1, tf_pmin=.025) 
      r4=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.3, tf_pmax=1, tf_pmin=.025) 
      r5=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.4, tf_pmax=1, tf_pmin=.025) 
      r6=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.5, tf_pmax=1, tf_pmin=.025) 
      r7=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.6, tf_pmax=1, tf_pmin=.025) 
      r8=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.7, tf_pmax=1, tf_pmin=.025) 
      r9=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.8, tf_pmax=1, tf_pmin=.025) 
      r10=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=.9, tf_pmax=1, tf_pmin=.025) 
      r11=tf_sees_ns(n.small=10, n.big=30, k.small=10000, k.big=10000, d.small=0, d.big=0, tf_share=1, tf_pmax=1, tf_pmin=.025) 

  rs=cbind(r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11)
  

 #the matrix rs is then copy-pasted onto Excel to plot the chart.
 



  






