[R] Bayes Factor from 2 x 2 contingency table

bayes_factor.cpp

// [[Rcpp::plugins("cpp11")]]
#include <Rcpp.h>
using namespace Rcpp;

double bf(NumericVector v, int n) {
  int h = sum(rbeta(n, 1 + v[0], 1 + v[1]) > rbeta(n, 1 + v[2], 1 + v[3]));

  if (h == n) {
    return R_PosInf;
  } else {
    return (double)h / (n - h);
  }
}

// [[Rcpp::export]]
NumericVector bayes_factor(NumericVector v, int n = 1000) {
  return NumericVector::create(
    _["a"] = v[0],
    _["b"] = v[1],
    _["c"] = v[2],
    _["d"] = v[3],
    _["bf"] = bf(v, n)
  );
}

bayes_factor_2x2.R

#!/usr/bin/env Rscript
#
#               outcome
#              +       -
#          +-------+-------+
#        + |   a   |   b   | a + b
#  group   +-------+-------+
#        - |   c   |   d   | c + d
#          +-------+-------+
#            a + c   b + d
#

sapply(c('dplyr', 'data.table', 'Rcpp', 'rbenchmark'),
       function(p) require(p, character.only = TRUE))

# C++
sourceCpp('bayes_factor.cpp') # bayes_factor()

# R
bayes_factor_r <- function(v_abcd, iter = 1000) {
  sim <- function(v, n) return(rbeta(n, 1 + v['a'], 1 + v['b']) > rbeta(n, 1 + v['c'], 1 + v['d']))
  odds <- function(h, n) return(h / (n - h))
  return(c(v_abcd, bf = odds(sum(sim(v_abcd, iter)), iter)))
}

# Fisher's exact test
fisher_test <- function(v_abcd, alt = 'two.sided', cnf = 0.95) {
  f <- fisher.test(matrix(v_abcd, nrow = 2), alternative = alt, conf.level = cnf)
  return(c(v_abcd,
           p_val = f$p.value,
           or = f$estimate[[1]],
           or_cil = f$conf.int[1],
           or_ciu = f$conf.int[2]))
}

# generate sample data for test
dt_test <- setnames(data.table(matrix(abs(ceiling(rnorm(400) * 100)), ncol = 4)),
                    c('a', 'b', 'c', 'd'))

# run
print(benchmark(dt_cpp <- as.data.table(t(apply(dt_test, 1, bayes_factor))),
                dt_r <- as.data.table(t(apply(dt_test, 1, bayes_factor_r))),
                dt_fex <- as.data.table(t(apply(dt_test, 1, fisher_test)))))
print(dt_rcpp <- dt_cpp %>%
                   inner_join(dt_r, by = c('a', 'b', 'c', 'd')) %>%
                   inner_join(dt_fex, by = c('a', 'b', 'c', 'd')))