Code to generate Bayes factors and plots for the probit models described in "What Are the Odds? Modern Relevance and Bayes Factor Solutions for MacAlister's Problem from the 1881 Educational Times."

create_probit_plots.R

#####################
# Code for computing probit model Bayes factors for
# 2x2 contingency table data
# Data and plots
# Code by Richard D. Morey and EJ Wagenmakers, July, 2015
# For:
# "What Are the Odds? Modern Relevance and Bayes 
# Factor Solutions for MacAlister's Problem from the 1881 Educational Times."
# Authors: Tahira Jamil, Maarten Marsman, Alexander Ly, 
# Richard D. Morey, Eric-Jan Wagenmakers
#
# See also the shiny applet at:
# https://richarddmorey.shinyapps.io/probitProportions
#####################

# These two lines of code read in the utility functions 
# that do all the work
require("devtools")
source_gist('https://gist.github.com/richarddmorey/d1ba3d4ff63a968944c3')


# MacAlister (1881) data
y = c(9,7)
N = c(14,10)

# limits on delta: one sided test
d.lim = c(0, Inf)

# Prior parameters
mu_mu = 0
sd_mu = sqrt(2)/2
mu_delta = 0
sd_delta = sqrt(2)

# Some helpful quantities
opt.par = c(mean(qnorm(y/N)), diff(qnorm(y/N)))
post.sd.delta = sqrt(sum((y/N*(1-y/N)/N)/dnorm(qnorm(y/N))^2))
post.lims = opt.par[2] + 2*c(-1,1)*sd_delta
dd = seq(max(c(post.lims[1],d.lim[1])),
         min(c(post.lims[2],d.lim[2])), len=100)
logConst = normConst(y, N, delta.lim = d.lim, mu_mu = mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

# This computes the Bayes factor in favor of the null
bf = bfProbit(y, N, log = FALSE, delta.lim = d.lim, opt.pars = opt.par, mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

# Create plot
yy = margDelta(dd-opt.par[2], y, N, logConst = logConst, delta.lim = d.lim, opt.pars = opt.par,mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    font.lab = 2, cex.axis = 1.3, bty = "n", las=1)

plot(dd,yy, ty='l',xlab = "",ylab="", lwd=2, axes=FALSE, xlim=c(0,2))
lines(dd, dnorm(dd, mu_delta, sd_delta)/diff(pnorm(d.lim, mu_delta, sd_delta)), lwd=2, lty=2)
abline(h=0, col="gray", lwd=2)
abline(v = 0, col="gray", lwd=2)
mtext(expression(paste("Effect (probit difference, ",delta,")",sep="")), side=1, line = 2.8, cex=1.5)
par(las=0)
mtext("Density", side=2, line = 0.8, cex=1.8)
axis(1)
text(0.85,1.2, "Posterior", cex=1.5)
text(1.6,0.45, "Prior", cex=1.5)


# Tuckman and Kennedy (2011) data
y = c(300,328)
N = c(351,351)
d.lim = c(0, Inf)
mu_mu = 0
sd_mu = .707
mu_delta = 0
sd_delta = sqrt(2)

opt.par = c(mean(qnorm(y/N)), diff(qnorm(y/N)))
post.sd.delta = sqrt(sum((y/N*(1-y/N)/N)/dnorm(qnorm(y/N))^2))
post.lims = opt.par[2] + 2*c(-1,1)*sd_delta
dd = seq(max(c(post.lims[1],d.lim[1])),
         min(c(post.lims[2],d.lim[2])), len=200)
logConst = normConst(y, N, delta.lim = d.lim, mu_mu = mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

# This computes the Bayes factor in favor of the null
bf = bfProbit(y, N, log = FALSE, delta.lim = d.lim, opt.pars = opt.par, mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

# Create plot
yy = margDelta(dd-opt.par[2], y, N, logConst = logConst, delta.lim = d.lim, opt.pars = opt.par,mu_mu, sd_mu = sd_mu, mu_delta = mu_delta, sd_delta = sd_delta)

par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5,
    font.lab = 2, cex.axis = 1.3, bty = "n", las=1)

plot(dd,yy, ty='l',xlab = "",ylab="", lwd=2, axes=FALSE, xlim=c(0,2))
lines(dd, dnorm(dd, mu_delta, sd_delta)/diff(pnorm(d.lim, mu_delta, sd_delta)), lwd=2, lty=2)
abline(h=0, col="gray", lwd=2)
abline(v = 0, col="gray", lwd=2)
mtext(expression(paste("Effect (probit difference, ",delta,")",sep="")), side=1, line = 2.8, cex=1.5)
par(las=0)
mtext("Density", side=2, line = 0.8, cex=1.8)
axis(1)
text(1,2.5, "Posterior", cex=1.5)
text(1.6,0.6, "Prior", cex=1.5)