Functions for computing the Bayes factor t test for an arbitrary prior on effect size

arbitrary_prior.R

## This function computes the (normalized) marginal likelihood.
## You will not use this function directly.
p.y.alt = function(t.stat, N, log.prior.dens,lo=-Inf,up=Inf,...){
  normalize = integrate(function(delta,...){
    exp(log.prior.dens(delta, ...))
  }, lower = lo, upper = up, ...)[[1]]
  py = integrate(function(delta,t.stat,N,...){
    exp(
      dt(t.stat, N - 1, ncp = delta*sqrt(N), log = TRUE) +
      log.prior.dens(delta, ...)  
    )
  },lower = lo, upper = up, t.stat = t.stat, N = N, ...)[[1]]
  py/normalize
}

## This function is what you will use.
## The ... argument is for the prior function arguments.
B01 = function(t.stat, N, log.prior.dens, lo = -Inf, up = Inf, ...){
  dt(t.stat, N - 1)/p.y.alt(t.stat,N,log.prior.dens,lo,up,...)
}

## This is a prior function. You can define whatever prior function you like.
## The function above will try to automatically normalize it, so
## it only has to integrate to a finite number (not necessarily to 1)
## But it must return the log density!
cauchy.prior = function(delta, rscale){
  dcauchy(delta, scale = rscale, log = TRUE)
}

# Normal prior function
normal.prior = function(delta, mu=0, sd=1){
  dnorm(delta, mu, sd, log = TRUE)
}

# Gamma prior function
gamma.prior = function(delta, scale = 1, shape = 1){
  dgamma(abs(delta), shape=shape, scale=scale, log = TRUE)
}

#######
# Demonstrate use of code
#######

t = 1
N = 10

# Full Cauchy 
B01(t,N,cauchy.prior,-Inf,Inf,rscale=sqrt(2)/2)
# Test against BayesFactor
1/BayesFactor::ttest.tstat(t, n1=N,simple=TRUE)

# Half Cauchy 
B01(t,N,cauchy.prior,0,Inf,rscale=sqrt(2)/2)
# Test against BayesFactor
1/BayesFactor::ttest.tstat(t, n1=N,nullInterval = c(0,Inf),simple=TRUE)

# Half Normal
B01(t,N,normal.prior,0,Inf,mu=0,sd=1)

# Gamma (chosen so median is sqrt(2)/2)
B01(t,N,gamma.prior,0,Inf,scale=0.2644319,shape=3)