Animating Robustness-Check of Bayes Factor Priors

Today I submitted a paper which contained some Bayesian analysis using Bayes factors (a default Bayesian t-test). As the test used a default prior on the effect size (a Cauchy distribution centred on zero with rate [r] = 0.707), I wanted to appease any reviewer concern that the final Bayes factor reported was some peculiarity of the prior rather than due to the data.

So, I performed a standard robustness check on the prior where I plot the resulting Bayes factor from a wide range of prior value (in this case, across a wide range of Cauchy rate parameters). The plot below is the result:paperRobust

The plot shows the Bayes factor with the default prior as the black dot; the line represents how the Bayes factor changes as a function of the Cauchy prior width. This analysis suggests the final conclusion is not a peculiarity of the prior, and as such the reported Bayes factor is quite robust.

What this plot doesn’t show, though, is what the prior distribution looks like as the rate (r) changes. Ideally, I wanted to show what the Cauchy distribution looks like across the entire range of r values explored in the plot above; this helps me (and I guess other people) understand how the prior is influencing the Bayes factor. Visualisation is very important sometimes.

So, I knocked up some R code which plots the robustness check above, but plots next to it the density of the Cauchy distribution with the rate parameter currently being explored. The animation cycles through increasing values for rate r, and shows the density for this rate on the left, and the Bayes factor for my data with this prior on the right. For me, this animation provides a nice insight to how the Bayes factor is being influenced by the prior.


R Code

The code for generating the individual plots is below. Note that all of these plots then need to be fed into your favourite gif-making device! (I used

### plots of prior robustness check
# gif created from PNG plots using
# clear R's memory
rm(list = ls())
# load the Bayes factor package
# set working directory (will be used to save files here, so make sure
# this is where you want to save your plots!)
setwd <- "D:/Work//Blog_YouTube code//Blog/Prior Robust Visualisation/plots"
### declare some variables for the analysis
# what is the t-value for the data?
tVal <- 3.098
# how many points in the prior should be explored?
nPoints <- 1000
# what Cauchy rates should be explored?
cauchyRates <- seq(from = 0.01, to = 1.5, length.out = nPoints)
# what effect sizes should be plotted?
effSize <- seq(from = -2, to = 2, length.out = nPoints)
# get the Bayes factor for each prior value
bayesFactors <- sapply(cauchyRates, function(x)
exp(ttest.tstat(t = tVal, n1 = 76, rscale = x)[['bf']]))
### do the plotting
# how many plots do we want to produce?
nPlots <- 50
plotWidth <- round(seq(from = 1, to = nPoints, length.out = nPlots), 0)
# loop over each plot
for(i in plotWidth){
# set up the file
currFile <- paste(getwd(), "/plot_", i, ".png", sep = "")
# initiate the png file
png(file = currFile, width = 1200, height = 1000, res = 200)
# change the plotting window so plots appear side-by-side
par(mfrow = c(1, 2))
# do the prior density plot
d <- dcauchy(effSize, scale = cauchyRates[i])
plot(effSize, d, type = "l", ylim = c(0, 5), xlim = c(-2, 2),
ylab = "Density", xlab = "Effect Size (d)", lwd = 2, col = "gray48",
main = paste("Rate (r) = ", round(cauchyRates[i], 3), sep = ""))
# do the Bayes factor plot
plot(cauchyRates, bayesFactors, type = "l", lwd = 2, col = "gray48",
ylim = c(0, max(bayesFactors)), xaxt = "n",
xlab = "Cauchy Prior Width (r)", ylab = "Bayes Factor (10)")
abline(h = 0, lwd = 1)
abline(h = 6, col = "black", lty = 2, lwd = 2)
axis(1, at = seq(0, 1.5, 0.25))
# add the BF at the default Cauchy point
points(0.707, 9.97, col = "black", cex = 1.5, pch = 21, bg = "red")
# add the BF for the Cauchy prior currently being plotted
points(cauchyRates[i], bayesFactors[i], col = "black", pch = 21, cex = 1.3,
bg = "cyan")
# add legend
legend(x = 0.25, y = 3, legend = c("r = 0.707", paste("r = ",
round(cauchyRates[i], 3),
sep = "")),
pch = c(21, 21), lty = c(NA, NA), lwd = c(NA, NA), pt.cex = c(1, 1),
col = c("black", "black"), = c("red", "cyan"), bty = "n")
# save the current plot

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s