Comparison between normal and logistic priors for proportion Bayes factors

do_plots.R

devtools::source_gist("e49c345fcd6cfe8f32535eda4bd8c29b", filename='utility_functions.R')

# Prior Setup
p0 = .5
rscale = .5
interval = c(.5,1)
# Do not change shift unless you want different prior shift
# will cause deviation from BayesFactor package
# (not that there's anything wrong with that, but changes demo)
shift = qlogis(p0)

# data (a range of y values for comparison)
y = 75:125
N = 200

bf.norm = 1/BF01_norm(y,N,shift,rscale,interval,p0)

# Check against BayesFactor package
bf.logis = vecBF(y,N,p0,rscale,interval)


# make some plots
# red is normal prior
# black is logistic

par(mfrow=c(1,2))

plot.interval = c(-4,4)
## Show priors
xx = seq(plot.interval[1],plot.interval[2],len=100)
dens.ind = xx>qlogis(interval[1]) & xx<qlogis(interval[2])
# Logistic prior
plot(xx, dens.ind*dlogis(xx,shift,rscale)/normalize_logis(shift, rscale, qlogis(interval)),ty='l',xlab="Log Odds",ylab="Prior Density")
# Normal prior
lines(xx, dens.ind*dnorm(xx,shift,rscale*pi/sqrt(3))/normalize_norm(shift, rscale, qlogis(interval)),col="red")


plot(y,bf.logis,ty='l',log="y",ylab="Bayes factor",xlab=paste0("y (out of ",N,")"),las=1)
lines(y,bf.norm,col="red")
abline(h=1,col="gray",lty=2)

utility_functions.R

### Comparison between normal and logistic prior for Bayes factors for proportions
## Support functions

fullAlt_norm = Vectorize(function(theta, y, N, shift, scale){
  p = plogis(theta)
  exp(dbinom(y, N, p, log = TRUE) + dnorm(theta, shift, scale * pi/sqrt(3), log = TRUE))
},"theta")

normalize_norm = function(shift, scale, interval){
  diff(pnorm(interval, shift, scale * pi/sqrt(3)))
}

restrictedAlt_norm = function(theta,y,N,shift,scale,interval){
  fullAlt_norm(theta,y,N,shift,scale) / normalize_norm(shift, scale, interval) * (theta>interval[1] & theta<interval[2])
}

margLike_norm = function(y, N, shift, scale, interval){
  theta_interval = qlogis(sort(interval))
  integrate(restrictedAlt_norm, theta_interval[1], theta_interval[2], 
            y = y, N = N, shift = shift, scale = scale, interval = theta_interval)[[1]]
}

BF01_norm = Vectorize(function(y, N, shift, scale, interval, null.p){
  dbinom(y,N,null.p) / margLike_norm(y, N, shift, scale, interval)
},"y")


# For plotting the logistic prior
normalize_logis = function(shift, scale, interval){
  diff(plogis(interval, shift, scale))
}

# For comparing to the BayesFactor results
vecBF = Vectorize(function(y,...){as.vector(BayesFactor::proportionBF(y,...)[1])}, "y")