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:

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 http://gifmaker.me/).

	### plots of prior robustness check
	# gif created from PNG plots using http://gifmaker.me/
	# clear R's memory
	rm(list = ls())
	# load the Bayes factor package
	library(BayesFactor)
	# 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"), pt.bg = c("red", "cyan"), bty = "n")
	# save the current plot
	dev.off()
	}
	#------------------------------------------------------------------------------

view raw priorRobustnessCheck.R hosted with ❤ by GitHub

Recruitment Order & Sequential Bayes Factors

One advantage of using Bayes factors (or any other form of Bayesian analysis) is the ability to engage in optional stopping. That is, one can collect some data, perform the critical Bayesian test, and stop data collection once a pre-defined criterion has been obtained (e.g., until “strong” evidence has been found in support of one hypothesis over another). If the criterion has not been met, one can resume data collection until it has.

Optional stopping is problematic when using p-values, but is “No problem for Bayesians!” (Rouder, 2014; see also this blog post from Kruschke). A recent paper in press at Psychological Methods (Shoenbrodt et al., in press) shows you how to use optional stopping with Bayes factors—so-called “sequential Bayes factors” (SBF).

The figure below a recent set of data from my lab where we used SBF for optional stopping. Following the advice of the SBF paper, we set our stopping criteria as N > 20 and BF10 > 6 (evidence in favour of the alternative) or BF10 < 1/6 (evidence in favour of the null); these criteria are shown as horizontal dashed lines. The circles represent the progression of the Bayes factor (shown in log scale) as the sample size increased. We stopped at 76 participants, with “strong” evidence supporting the alternative hypothesis.

One thing troubled me when I was writing these results up: I could well imagine a skeptical reviewer noticing that the Bayes factor seemed to be in support of the null between sample sizes of 30 and 40, although the criterion for stopping was not reached. A reviewer might wonder, therefore, whether had some of the later participants who showed no effect (see e.g., participants 54–59) been recruited earlier, would we have stopped data collection and found evidence in favour of the null? Put another way, how robust was our final result to the order in which our sample were recruited to the experiment?

Assessing the Effect of Recruitment Order

I was not aware of any formal way to assess this question, so for the purposes of the paper I was writing I knocked up a quick simulation. I randomised the recruitment order of our sample, and plotted the sequential Bayes factor for this “new” recruitment order. I performed this 50 times, plotting the trajectory of the SBF. I reasoned that if our conclusion was robust against recruitment order, the stopping rule in favour of the null should be rarely (if ever) met. The results of this simulation are below. As can be seen, the stopping criterion for the null was never met, suggesting our results were robust to recruitment order.

 

Recruitment Order Simulation

Then, I became curious: How robust are SBFs to recruitment order when we know what the true effect is? I simulated a set of 100 participants for an effect known to be small (d = 0.3), and plotted the SBF as “participants” were recruited. I tried different random seeds until I found a pattern I was interested in: I wanted a set of data where the SBF comes very close to the stopping rule in favour of the null, even though we know the alternative is true. It took me a few tries, but eventually I got the plot below.

Note the dip toward the null criterion is even more extreme than in my experiment. Thus, this set of data seemed perfect to test my method for assessing the robustness of the final conclusion (“strong evidence for alternative”) against recruitment order. I followed the same procedure as before, simulating 50 random recruitment orders, and plotted the results below.

The null criterion is only met in one simulated “experiment”. This seems to be very pleasing, at least for the paper I am due to submit soon: the SBF (at least in these two small examples) seem robust against recruitment order.

R Code

If you would like to see the R code for the above plots, then you can download the html file for this post which includes R code.