Bayes_reproducibility.R

## from the reproducibility project code here https://osf.io/vdnrb/
#make sure this file is in your working directory
info <- GET('https://osf.io/fgjvw/?action=download', write_disk('rpp_data.csv', overwrite = TRUE)) #downloads data file from the OSF
MASTER <- read.csv("rpp_data.csv")[1:167, ]
colnames(MASTER)[1] <- "ID" # Change first column name to ID to be able to load .csv file

studies<-MASTER$ID[!is.na(MASTER$T_r..O.) & !is.na(MASTER$T_r..R.)] ##to keep track of which studies are which
studies<-studies[-31]##remove the problem studies (46 and 139)
studies<-studies[-80]

orig<-MASTER$T_r..O.[studies] ##read in the original rs that have matching rep rs
rep<-MASTER$T_r..R.[studies]  ##read in the rep rs that have matching original rs

N.R<-MASTER$T_N_R_for_tables[studies] ##n of replications for analysis
N.O<-MASTER$T_N_O_for_tables[studies] ##n of original studies for analysis

p<-MASTER$T_pval_USE..R.[studies] #extract p-values for the studies

bfRep<- numeric(length=95) #prepare for running replications against original study posterior

#download the code for replication functions from here https://osf.io/v7nux/ and load functions into globale environment
for(i in 1:95){
        bfRep[i]<- repBfR0(nOri=N.O[i],rOri=orig[i],nRep=N.R[i],rRep=rep[i])
}

#bfRep Remove hash to disply BFs

#create bin numbers
bfstrength<-seq(1,8,1)
#create labels for the bins
barlabels<-c("BF<1/100","1/100<BF<1/10","1/10<BF<1/3","1/3<BF<1", "1<BF<3","3<BF<10","10<BF<100","BF>100")
#create category labels to add to bins
g<-c("Very strong","Strong","Moderate","Insensitive","Insensitive", "Moderate","Strong","Very strong")
barlabels<-as.vector(barlabels)


#create new dummy variables for different BFs in the categories
bf<-numeric(length=95)
for(i in 1:95){
        if(bfRep[i]<.001){ #BF<1/100
                bf[i]<-bfstrength[1]
        }
        if(bfRep[i]<.1 & bfRep[i]>=.001){ #1/100<BF<1/10
                bf[i]<-bfstrength[2]
        } 
        if(bfRep[i]<.3333 & bfRep[i]>=.1){ #1/10<BF<1/3
                bf[i]<-bfstrength[3]
        }
        if(bfRep[i]<1 & bfRep[i]>=.3333){ #1/3<BF<1
                bf[i]<-bfstrength[4]
        }
        if(bfRep[i]>1 & bfRep[i]<=3){ #1<BF<3
                bf[i]<-bfstrength[5]
        }
        if(bfRep[i]>3 & bfRep[i]<=10){ #3<BF<10
                bf[i]<-bfstrength[6]
        }
        if(bfRep[i]>10 & bfRep[i]<=100){ #10<BF<100
                bf[i]<-bfstrength[7]
        }
        if(bfRep[i]>100){ #BF>100
                bf[i]<-bfstrength[8]
        }
}

table(bf) #shows counts for each bin

#plot the bins
barplot(table(bf),names.arg=barlabels,border="gray80",col="paleturquoise1",
        xlab="Replication Bayes Factor Categories",las=1,
        main="Bayesian Replication Outcomes from the Reproducibility Project: Psychology",
        ylim=c(0,25), sub="BFs > 1 are evidence in favor of the original effect")
arrows(x0=4.8,x1=-.1,y0=21,y1=21,col="red",lwd=5) #replication strength arrow
arrows(x0=5,x1=9.8,y0=21,y1=21,col="darkorchid1",lwd=5) #ditto
text("STRONGER REPLICATION 'FAILURE'",y=22,x=2.5,cex=1.5,col="red") #add text above arrows
text("STRONGER REPLICATION 'SUCCESS'",y=22,x=7.5,cex=1.5,col="darkorchid1") #ditto
text(g,x=c(.7,1.9,3.1,4.3,5.5,6.7,7.9,9.1),y=rep(.5,8)) #add verbal labels to the bottom of bins (e.g., very strong)

k=c(2/95,16/95,20/95,17/95,9/95,8/95,8/95,16/95) #percent of reps in each bin

tab2<-cbind(studies,N.O,N.R,orig,rep,bfRep,p,bf,g) #create matrix of values