###############################################################################
### Compute Bayes factors after a replication study (t-test) 
### to find out whether the original findings are replicated
### Last revised: 19-7-2013
### Author: Josine Verhagen, A.J.Verhagen@uva.nl
### www.josineverhagen.com
### INPUT:  t original, n original, t replication, n replication
### OUTPUT: Replication Bayes factor, JZS Bayes factors one and two sided,
### Bayes factor Bayarri and Mayoral,  Meta-analysis Bayes factor, 
### prior and posterior distributions for effect size
###############################################################################

# clears workspace:  
rm(list=ls(all=TRUE)) 

# Set work directory  (replace with own path)
#dir <- "C:/Users/Josine/Documents/My Dropbox/UVA/JosineEJ/Replicationprior/Code/"
#WB <- "C:/Users/Josine/WinBUGS14"

dir <- "C:/Users/averhag2/Dropbox/UVA/JosineEJ/ReplicationPrior/Code/"

setwd(dir)

# Load and if necessary install the required libraries
library(R2WinBUGS)  # Only necessary for posterior plots for JZS prior
library(MCMCpack)

# Load the functions necessary to compute the replication prior (file saved in work directory)
source(paste(dir,'Repfunctionspack.R', sep = "")) 

####################################################################
#From raw data, first compute trep, the t value from the replication
####################################################################

x <- rnorm(100,.5,1)
y <- rnorm(100,0,1)

# Make sure to input the vectors in the order of the original study. 


#########################
#For a one sample t-test: 
#########################

# For a one sample test: 
trep <- t.test(x,y, paired = TRUE)$statistic
n2 <- length(x) 

#Example: get data from original study, t and n  
tobs <- 2  # t observed in the original study  
n1 <- 50   # number of subjects in the original study 

### Function to compute all relevant Bayes factors for replication: ###
# Independent JZS Bayes factors in Original and Replication studies, one and two sided. BF10, H1 (effect) vs H0 (no effect) 
# Replication Bayes factors: Is the effect the same as in the original study or zero? BF10, H1 (effect size = posterior original study) vs H0 (effect szie = 0)   
# Bayarri and Mayoral Bayes Factor: is the replicated effect size equal to the original or different? BF01, H0 (same effect size) vs H1 (different effect size) 
# Meta-analysis Bayes factor: combining all effect sizes. BF10, H1 (there is an overall effect) vs H0 (there is no overall effect).    

BFSALL(tobs,trep,n1,n2, Type = "ALL")

### Plots for one and two-sided independent JZS Bayes factors: normalized ###
### Savage Dickey Bayes factor output ###
# To make these plots, you need WinBUGS. You can download it: 
# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml
# Check where the program folder is and insert in as WB= 
# The default is: WB ="c:/Program Files/WinBUGS14"

WB ="c:/Program Files/WinBUGS14"

# Original 
PlotJZS(tobs,n1, WB=WB)
# Replication two-sided 
PlotJZS(trep,n2, WB=WB)
# Replication one-sided Upper 
PlotJZS(trep,n2, side=2, WB=WB)
# Replication one-sided Lower 
#PlotJZS(trep,n2, side=1, WB=WB)

# Replication Bayes Factor + a plot of prior and posterior: 
out <- ReplicationBayesfactorNCT(tobs, trep, n1, n2, plot=1, post =1)

# save plot (you can do this for every plot) 
title <-  paste ("Replication Bayes factor")
dev.copy(device=png, file= paste(dir,title,".png", sep=""))
dev.off()


############################
#for a two sample t test: 
############################

# for a two sample test: 
trep <- t.test(x,y, paired = FALSE)$statistic
n2 <- length(x) 
m2 <- length(y) 

#Example data: 
tobs <- 2  # t observed original study  
n1 <- 50   # number of subjects in group 1 in the original study 
m1 <- 50   # number of subjects in group 2 in the original study 
sample  <- 2

# Independent JZS Bayes factors in Original and Replication studies, one and two sided: 
# Replication Bayes factors: Is the effect the same as in the original study or zero? 
# Bayarri and Mayoral Bayes Factor: is the replicated effect size equal to the original or different
# Meta-analysis Bayes factor: combining all effect sizes  

BF <- BFSALL(tobs, trep, n1, n2, m1=m1,m2=m2, sample=sample, Type = "ALL")
BF

### Plots for one and two-sided independent JZS Bayes factors: normalized ###
### Savage Dickey Bayes factor output ###
# To make these plots, you need WinBUGS. You can download it: 
# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml
# Check where the program folder is and insert in as WB= 
# The default is: WB ="c:/Program Files/WinBUGS14"

WB ="c:/Program Files/WinBUGS14"

# Original two-sided
PlotJZS(tobs,n1,m1=m1, sample=sample, WB=WB)
#Replication two-sided 
PlotJZS(trep, n2, m1 = m2, sample=sample, WB=WB)
#Replication one-sided Upper
PlotJZS(trep, n2, m1 = m2, sample=sample, side=2, WB=WB)
#Replication one-sided 
#PlotJZS(trep, n2, m2, sample=sample, side=2, WB=WB)

# Replication Bayes Factor + a plot of prior and posterior 
out <- ReplicationBayesfactorNCT(tobs, trep, n1, n2, m1=m1,m2=m2, sample=sample, plot=1, post=1)

# save plot (you can do this for every plot) 
title <-  paste ("Replication Bayes factor")
dev.copy(device=png, file= paste(dir,title,".png", sep=""))
dev.off()

###########################
#Examples
###########################


### Elliot et al. (2010):  Red, rank and romance ###
# Original study: t value: 2.18, n1: 10, m1: 11 
# 3 Replication studies: t values: 3.06,.25,2.44; n2: 27,27,16; m1: 27,27,17 

# Independent JZS Bayes factors in Original and Replication studies, one and two sided: 
# Replication Bayes factors: Is the effect the same as in the original study or zero? 
# Bayarri and Mayoral Bayes Factor: is the replicated effect size equal to the original or different
# Meta-analysis Bayes factor: combining all effect sizes  

BFelliot <- BFSALL(2.18,c(3.06,.25,2.44),10, c(27,27,16), 11, c(27,27,17), sample=2, Type = "ALL")
BFelliot

# Without the male study: 
BFelliot2 <- BFSALL(2.18,c(3.06,2.44),10, c(27,16), 11, c(27,17), sample=2, Type = "ALL")
BFelliot2

# with the later replication in Elliot and Maier (2013)
#The full study yielded the following results: Women in the red
#condition (M 
 6.29, SD 
 1.35) perceived the target man to be
#more attractive than women in the gray condition (M 
 5.93,
#SD 
 1.49), but this difference fell short of statistical significance
#(t 
 1.51, p 
 .13, d 
 0.25). 

effectsizes(2.18,c(3.06,.25,2.44,1.51),10, c(27,27,16,75), 11, c(27,27,17,69), sample=2)

### Plots for one and two-sided independent JZS Bayes factors: normalized ###
### Savage Dickey Bayes factor output ###
# To make these plots, you need WinBUGS. You can download it: 
# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/contents.shtml
# Check where the program folder is and insert in as WB= 
# The default is: WB ="c:/Program Files/WinBUGS14"

WB ="c:/Program Files/WinBUGS14"

# Original JZS plot
PlotJZS(2.18,10,11, sample = 2, WB=WB)

# save plot (you can do this for every plot) 
title <-  paste ("JZS original Elliot")
dev.copy(device=png, file= paste(dir,title,".png", sep=""))
dev.off()

#Replication JZS plots two-sided 
PlotJZS(3.06,27,27, sample = 2, WB=WB)
PlotJZS(.25,27,27, sample = 2, WB=WB)
PlotJZS(2.44,16,17, sample = 2, WB=WB)

# Replication JZS plots one-sided (Upper) (yhigh = 2 makes sure all plots have equal height)
PlotJZS(3.06,27,27, sample = 2, side=2, WB=WB, yhigh = 2)
PlotJZS(.25,27,27, sample = 2, side=2, WB=WB, yhigh = 2)
PlotJZS(2.44,16,17, sample = 2, side=2, WB=WB, yhigh = 2)

# Plots Replication Bayes Factor (yhigh = 2 makes sure all plots have equal height)
ReplicationBayesfactorNCT(2.18, 3.06, 10, 27, 11, 27, sample =2, plot=1 , post = 1, yhigh = 2)
ReplicationBayesfactorNCT(2.18, .25, 10, 27, 11, 27, sample =2, plot=1 , post = 1, yhigh = 2)
ReplicationBayesfactorNCT(2.18, 2.44, 10,16 ,11, 17, sample =2, plot=1 , post = 1, yhigh = 2)

#####################################################################################
### Shanks et al. (2013) replicating Professor priming (LeBoeuff and Estes, 2004) ###
# original study: t values  tPS2, n1: 22 , m1: 22 
#tPS <- sqrt(7.12) # Compute t value
tPS2 <- (56.2- 45.2) / (sqrt((20*11.1^2 + 22* 13^2)/ 42)  *  sqrt(1/21 + 1/23) )  # t-value computed form the reults table


# first replication study: -.25 - 25, 24
# the second  study should have effect size: t =  -1.25 p = .11 16 , 16 

#If all are counted as replications of the same effec: 
BFshanksall <-  BFSALL(tPS2,c(-.25, -1.25),22, c(25,16), 22,c(24,16),sample=2, Type="ALL")
BFshanksall

effectsizes(tPS2,c(-.25, -1.25),22, c(25,16), 22,c(24,16),sample=2)

# Original
PlotJZS(tPS2,22,22, sample = 2, WB=WB)

#Replication plots two-sided 
PlotJZS(-.25,25,24, sample = 2, WB=WB)
PlotJZS(-1.25,16,16, sample = 2, WB=WB)


# Replication plots one-sided (Upper) (yhigh = 2 makes sure all plots have equal height)
PlotJZS(-.25,25,24, sample = 2, side = 2, WB=WB)
PlotJZS(-1.25,16,16, sample = 2, side=2, WB=WB)

#Plots Replication Bayes Factor (yhigh = 2.2 makes sure all plots have equal height: search for a good number if you want this)
ReplicationBayesfactorNCT(tPS2,-.25,22,25,22,24,sample=2, plot=1 , post = 1, yhigh = 2.2)
ReplicationBayesfactorNCT(tPS2, -1.25,22, 16, 21,16,sample=2, post = 1, plot=1, yhigh=2.2)


####################################################################################
### Neill and Kahan (1999)  replicating Negative priming (Milliken et al., 1998) ###
# Original study: t value: 3.29, n1: 20  
# 2 Replication studies: t values: 2.06,-2.4 n2: 30,43 

BFNeill <-  BFSALL(3.29,c(2.06,-2.4), 20, c(30,43), Type = "ALL")
BFNeill

effectsizes(3.29,c(2.06,-2.4), 20, c(30,43))

# Original
PlotJZS(3.29,20,sample = 1, WB=WB)

#Replication plots two-sided 
PlotJZS(2.06,30, sample = 1, WB=WB)
PlotJZS(-2.4,43, sample = 1, WB=WB)

# Replication plots one-sided (Upper) (yhigh = 2 makes sure all plots have equal height)
PlotJZS(2.06,30, sample = 1, side = 2, WB=WB)
PlotJZS(-2.4,43, sample = 1, side = 2, WB=WB)

# Replication plots one-sided (lower) (yhigh = 2 makes sure all plots have equal height)
PlotJZS(2.06,30, sample = 1, side = 1, WB=WB)
PlotJZS(-2.4,43, sample = 1, side = 1, WB=WB)

#Plots Replication Bayes Factor (yhigh = 3.5 makes sure all plots have equal height: search for a good number if you want this)
ReplicationBayesfactorNCT(3.29, 2.06,20, 30,sample=1, post = 1, plot=1, yhigh=3.5)
ReplicationBayesfactorNCT(3.29,-2.4,20,43,sample=1, plot=1 , post = 1, yhigh = 3.5)

##########################################################################################


###########################
#Normal Approximation 
##########################

PD <- matrix(0,6,4)
nt <- c(10,20,50,100,200,1000)

for (m in 1:6) {
  PD[m,1] <- PosteriorDeltaLC(2,nt[m])[[1]]
  PD[m,2] <- PosteriorDelta(2,nt[m])[[1]]
  PD[m,3] <- PosteriorDeltaLC(2,nt[m])[[2]]
  PD[m,4] <- PosteriorDelta(2,nt[m])[[2]]
}

ns <- c(10,20,50,100,200,1000)
vxlim <- c(2,2,1,1,.5,.5)
vylim <- c(2,3,5,6,8,9)

rownames(PD) <- ns
colnames(PD) <- c("mean Lecoutre", "mean Cumming", "sd Lecoutre", "sd Cumming") 

library(xtable)

print(xtable(PD), file =  paste(dir,"Normapproxcomp.tex", sep=""))
      

m<- 3

par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,font.lab = 2, cex.axis = 1.3, bty = "n", las=1)

for (m in 1:5){
plot(density(rnorm(1000000,PD[m,1],PD[m,3])), ylim = c(0,vylim[m]), xlim = c(-.1,vxlim[m]), lty = 2, xlab= "",  main = "", col = 'gray', lwd=2)
par(new=T)
plot(density(rnorm(1000000,PD[m,2],PD[m,4])), ylim = c(0,vylim[m]), xlim = c(-.1,vxlim[m]), lty = 1, xlab= "",  main = "", col = 'gray', lwd=2 )
par(new=T)

seqD <-  seq(-1,7,.01)
densD <- dt(2,(ns[m]-1),seqD)
#plot(seqD/sqrt(ns[m]), densD*sqrt(ns[m]), ylim = c(0,vylim[m]), xlim = c(-.1,vxlim[m]), lty = 1, xlab= "", ylab = "",  main = paste("Posterior t=",2,"n=", ns[m]) , lwd=2, type= 'l', col = 'black')
plot(seqD/sqrt(ns[m]), densD*sqrt(ns[m]), ylim = c(0,vylim[m]), xlim = c(-.1,vxlim[m]), lty = 1, xlab= "", ylab = "",  main = "", lwd=2, type= 'l', col = 'black')
par(new=F)

axis(1)
axis(2)
mtext(expression(Effect ~ size ~ delta), side=1, line = 2.8, cex=2)

legend(-.1,vylim[m], c("Posterior distribution","Lambda prime approximation","Iterative approximation"), lwd= c(1,1,1),  lty = c(1,1,2), col = c(1,'gray','gray',1), bty = 'n' , cex= 1.2)

legend((vxlim[m] - ((vxlim[m]- (-.1))/1.9)),vylim[m]/2 , paste("n = ", ns[m]),  bty = 'n' , cex= 2)

title <-  paste ("Normapproxt2n", ns[m], sep = "") 

dev.copy(device=png, file= paste(dir,title,".png", sep=""))
dev.off()
}

mtext(expression(t[obs] ~ = ~ 2, n = 10), side=1, line = 2.8, cex=2)