Code for animating robustness check of Bayesian prior.

priorRobustnessCheck.R

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