Skip to content

Instantly share code, notes, and snippets.

@davharris
Last active September 2, 2015 06:10
Show Gist options
  • Save davharris/9c49c875898a38ecb0a5 to your computer and use it in GitHub Desktop.
Save davharris/9c49c875898a38ecb0a5 to your computer and use it in GitHub Desktop.
Code for a kernel density plot showing the Bayes Factors associated with the "Bayesian Reproducibility Project" blog post by Alex Etz.
# 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)
@davharris
Copy link
Author

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment