Code for a kernel density plot showing the Bayes Factors associated with the "Bayesian Reproducibility Project" blog post by Alex Etz.

bayes-reproducibility-density.R

# Start by running the code on Alex Etz's blog post
# http://alexanderetz.com/2015/08/30/the-bayesian-reproducibility-project/

# Cap the Bayes Factor at 10^12 to allow zooming in further
bfRepCapped = pmin(bfRep, 1E12)


# Kernel density estimate of the log-transformed Bayes Factors
D = density(log(bfRepCapped))

# Plot the kernel density with the original scale along the x axis
plot(D$x, D$y, yaxs = "i", ylim = c(0, 1.04 * max(D$y)), type = "l", axes = FALSE, ylab = "", xlab = "Bayes Factor (log scale)")
axis(1, log(10^seq(-10, 100, 2)), 10^seq(-10, 100, 2))
abline(v = log(1), lwd = 2) # Vertical line at BF=1
rug(log(bfRepCapped)) # Add a small line along the x axis for each data point


text(log(1E12), .02, "(10^12 or greater)", cex = .7)

An updated version with more labels:

# Cap the Bayes Factor at 10^12 to allow zooming in further
bfRepCapped = pmin(bfRep, 1E12)

# Kernel density estimate of the log-transformed Bayes Factors
D = density(log(bfRepCapped))

# Plot the kernel density with the original scale along the x axis
plot(D$x, D$y, yaxs = "i", ylim = c(0, .18), type = "l", axes = FALSE, ylab = "", xlab = "Bayes Factor (log scale)")
axis(1, log(10^seq(-10, 100, 2)), 10^seq(-10, 100, 2))
abline(v = log(1), lwd = 2) # Vertical line at BF=1
abline(v = log(10^c(-1, 1)), col = "darkgray")

rug(log(bfRepCapped)) # Add a small line along the x axis for each data point


text(log(1E12), .02, "(10^12 or greater)", cex = .7)
arrows(log(1E1), .15, log(1E7), .15)
arrows(log(1E-1), .15, log(1E-7), .15)
text(log(1E-1), .165, "Strong replication failure", cex = .8, pos = 2)
text(log(1E1), .165, "Strong replication success", cex = .8, pos = 4)

text(log(1), .2, "Weak/moderate\nreplication\nevidence", xpd = TRUE, cex = .8)