###############################################################################
### 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)