R/bf_mws.R

In tomfaulkenberry/anovaBFcalc: ANOVA Bayes Factor Calculator

#' Compute MWS Bayes factors from ANOVA summary statistics.
#' 
#' This function employs the MWS method, based on the Pearson
#' Type VI priors recommended by Maruyama (2009) and Wang and 
#' Sun (2014).
#' 
#' Minimally, the user must provide three inputs:
#' F = the observed F statistic
#' df1 = the "between-treatments" degrees of freedom
#' df2 = the residual degrees of freedom
#' 
#' The function outputs the value of BF_01 by default,
#' though this can be changed to BF_10 by the user (see
#' below).

#' 
#' Additional options may be specified -- see below
#' @param F the observed F statistic
#' @param df1 the between-treatments degrees of freedom
#' @param df2 the residual degrees of freedom
#' @param report.as a string indicating whether to report Bayes factor
#'     as support for null ("BF01") or alternative ("BF10"). Defaults
#'     to "BF01"
#' @param alpha hyperparameter for scale of Pearson Type VI prior. Wang and Sun
#'     recommend choosing alpha between -0.5 and 0. Defaults to -0.5, which 
#'     provides asymptotic approximation to multivariate Cauchy prior.  
#' @export
#' @author Thomas J. Faulkenberry <faulkenberry@tarleton.edu>
#' @references Maruyama, Y. (2009). A Bayes factor with reasonable model selection consistency for 
#'     ANOVA model. arXiv preprint arXiv:0906.4329.
#'     
#'     Wang, M., & Sun, X. (2014). Bayes Factor Consistency for One-way Random Effects Model. 
#'     Communications in Statistics - Theory and Methods, 43(23), 5072–5090. 
#'     doi:10.1080/03610926.2012.739252
#' 
#' @examples
#' ## using MWS method with alpha = -0.5 (default)
#' ## observed F(2,15) = 7.16
#' bf_mws(F=7.16, df1=2, df2=15)
#' 
#' ## using alpha=0 and reporting as BF10
#' bf_mws(F=7.16, df1=2, df2=15, report.as="BF10", alpha=0)

bf_mws <- function(F, df1, df2, report.as="BF01", alpha=-0.5) {
  bf = gamma(df1/2+alpha+1)*gamma(df2/2)/gamma(df1/2+df2/2)/gamma(alpha+1)*(df2/(df2+df1*F))^(alpha-df2/2+1)
  
  if (report.as == "BF10") {
     return(c(BF10 = bf))
  }
  else {
    return(c(BF01 = 1/bf))
  }
}