R/inferBF.R

In mattansb/BFEffect: Analyse models and effects estimated with BayesFactor

#' Create hypotheses for inferBF
#'
#' @param ... (Optionally named) hypotheses to be tested. The arguments in ... are automatically quoted and evaluated in the context of the prior/posterior samples.
#'
#' @return A list of named quoted hypotheses.
#'
#' @examples
#' \dontrun{
#'
#' hyp(ROPE = -.2 < b & b < .2,
#'     b - 0)
#' ## 2 Hypotheses:
#' ##        Hypothesis.specification
#' ## ROPE         -0.2 < b & b < 0.2
#' ## Test_2                    b - 0
#' }
#'
#'
#' @importFrom rlang enquos
hyp <- function(...) {
  # NOT EXPORTED
  X <- enquos(...)

  empty_names <- map_lgl(names(X),~identical(.x, ""))
  names(X)[empty_names] <- paste0('Test_',seq_along(empty_names)[empty_names])
  names(X) <- make.unique(names(X),sep = '_')
  attr(X,'class') <- c('hypBF','quosures')
  X
}

#' Compute BFs for hypotheses
#'
#' Compute BFs for hypotheses. A kin to `brms::hypothesis`.
#'
#' Please note that `inferBF` is in alpha stages of development, and though it is suited
#' for estimating effects (contrasts / slopes), it should not be used for estimating predicted
#' values (due to `BayesFactor` currently \href{https://github.com/richarddmorey/BayesFactor/issues/133}{not supporting "intercept" estimation}).
#'
#' @author Mattan S. Ben-Shachar
#' @param BF an object returned by any of the BayesFactor modeling functions
#' @param hypothesis specified hypotheses, returned from `hyp`.
#' @param hdi.level credible level
#' @param iterations number of posterior iterations
#' @param seed A single numeric value passed to set.seed to make results reproducible.
#' @param ... passed to posterior sampling function
#' @param nice.names whould the coeff names be a little nicer? (for \code{BFBayesFactor} objects)
#' @param index which BF mddel to test? (for \code{BFBayesFactor} objects)
#'
#' @return A list of data.frames, for each hypothesis, that prints nicely.
#'
#' @details
#' For directed hypotheses:
#' \enumerate{
#'   \item \code{BF.v.Prior} - The change in probability of hypothesis from the prior distribution
#'   to the posterior. Can also be thought as the BF of the restricted model vs. the
#'   unrestricted model.
#'   \item \code{BF.v.Null} - The evidence ratio between the denominator of the BF object
#'   and the restricted mode. Calculated as BF.v.Prior times the models BF.
#'   \item \code{BF.v.Opposite} - The BF of `a>b`/`b>a`
#'   \item \code{Post.Prob} - Posterior probability
#' }
#' For undirected (point) hypotheses:
#' \enumerate{
#'   \item \code{Estimate} - The median of the posterior distribution
#'   \item \code{Est.Error} - Standard error of the posterior distribution
#'   \item \code{HDI.Lower} - upper bound of HDI credible interval
#'   \item \code{HDI.Upper} - lower bound of HDI credible interval
#'   \item \code{BF} - The Savage-Dickey density ratio
#' }
#'
#' @seealso \code{\link[BFEffects]{as.data.frame.inferBF}}
#'
#'
#' @examples
#' \dontrun{
#' # =======
#' # anovaBF
#' # =======
#' library(BayesFactor)
#' library(BFEffects)
#' data(puzzles)
#' BF <- anovaBF(RT ~ shape*color + ID, data=puzzles, whichRandom = "ID", progress = FALSE)
#'
#' # If no hypotheses are passed, `inferBF` returns the paramater names:
#' inferBF(BF[4])
#' ## [1] "mu"                     "round"                  "square"                 "color"                  "monochromatic"
#' ## [6] "X1"                     "X2"                     "X3"                     "X4"                     "X5"
#' ## [11] "X6"                     "X7"                     "X8"                     "X9"                     "X10"
#' ## [16] "X11"                    "X12"                    "round...color"          "round...monochromatic"  "square...color"
#' ## [21] "square...monochromatic" "sig2"                   "g_shape"                "g_color"                "g_ID"
#' ## [26] "g_shape.color"
#' ## Test a model that has two simple effects of color in each level of shape.
#'
#' puzzle_hypotheses <- hyp(directed_h = (color + round...color < monochromatic + round...monochromatic) &
#'                            (color + square...color < monochromatic + square...monochromatic),
#'                          point_h = color - monochromatic)
#' inferBF(BF[4],puzzle_hypotheses)
#' ## Hypotheses tested based on the model:
#' ##         RT ~ shape + color + shape:color + ID
#' ##
#' ## Directional Tests:
#' ##            BF.v.Prior BF.v.Null BF.v.Opposite Post.Prob
#' ## directed_h      2.757    11.661        11.484      0.92
#' ##
#' ## Point Tests:
#' ##         Estimate Est.Error HDI.Lower HDI.Upper    BF
#' ## point_h    -0.86     0.383    -1.576    -0.079 3.128
#' ## ---
#' ## HDI level: 0.95
#' ## Point BF calculated using the Savage-Dickey method
#'
#'
#' # ============
#' # regressionBF
#' # ============
#' data(attitude)
#' BF <- regressionBF(rating ~ ., data = attitude, progress = FALSE)
#' # test for a positive slope for complaints and a positive slope for learning.
#' attitude_hypo <- hyp(complaints = complaints > 0,
#'                      attitude   = learning > 0,
#'                      both       = learning > 0 & complaints > 0)
#' inferBF(head(BF)[2], attitude_hypo)
#' ## Hypotheses tested based on the model:
#' ##   rating ~ complaints + learning
#' ##
#' ## Directional Tests:
#' ##            BF.v.Prior BF.v.Null BF.v.Opposite Post.Prob
#' ## complaints      2.008  416208.6           Inf     1.000
#' ## attitude        1.861  385670.5        12.699     0.927
#' ## both            3.730  773203.4        12.699     0.927
#'
#' # =================
#' # Advanced analyses
#' # =================
#' BF <- regressionBF(rating ~ ., data = attitude, progress = FALSE)
#'
#' advance_hypo <- hyp(
#'   # privileges is not between [-0.2,0.2]
#'   # Post.Prob is the % of posterior mass outside the ROPE
#'   ROPE_privileges = !(privileges > -0.2 & privileges < 0.2),
#'   # Estimate privileges
#'   estimate_privileges = privileges,
#'   # Estimate privileges based on a one sided prior
#'   estimate_privileges2 = ifelse(privileges > 0,privileges,NA)
#' )
#'
#' inferBF(BF[2], advance_hypo)
#' ## Hypotheses tested based on the model:
#' ##         rating ~ privileges
#' ##
#' ## Directional Tests:
#' ##                 BF.v.Prior BF.v.Null BF.v.Opposite Post.Prob
#' ## ROPE_privileges       1.23     3.907         4.118     0.805
#' ##
#' ## Point Tests:
#' ##                      Estimate Est.Error HDI.Lower HDI.Upper    BF
#' ## estimate_privileges     0.339     0.172     0.017     0.679 3.394
#' ## estimate_privileges2    0.342     0.162     0.037     0.655 5.940
#' ## ---
#' ## HDI level: 0.95
#' ## Point BF calculated using the Savage-Dickey method
#' }
#'
#' @references Wagenmakers, E. J., Lodewyckx, T., Kuriyal, H., & Grasman, R. (2010).
#' Bayesian hypothesis testing for psychologists: A tutorial on the Savage-Dickey
#' method. Cognitive psychology, 60(3), 158-189.
#'
#' Morey, R. D., & Wagenmakers, E. J. (2014). Simple relation between Bayesian
#' order-restricted and point-null hypothesis tests. Statistics & Probability
#' Letters, 92, 121-124.
#'
#' Morey, R. D. (2015). Multiple Comparisons with BayesFactor, Part 2 - order
#' restrictions [Blog post]. https://richarddmorey.org/2015/01/multiple-comparisons-with-bayesfactor-part-2-order-restrictions/
#'
#'
#' @importFrom HDInterval hdi
#' @importFrom rlang eval_tidy
inferBF <- function(BF, hypothesis = NULL, hdi.level = .95,
                    iterations = 10000,seed = NULL,
                    index = 1, nice.names = TRUE,...){
  # NOT EXPORTED

  if (length(hypothesis) == 0) iterations <- 3
  set.seed(seed)

  BF <- BF[index]

  # get posteriors
  post_samples <- as.data.frame(BayesFactor::posterior(BF, iterations = iterations, progress = FALSE,...))

  # get priors
  prior_samples <- sample_from_priors(BF, post_samples)

  # nice names
  if (nice.names) {
    post_samples <- nicify_names(post_samples)
    prior_samples <- nicify_names(prior_samples)
  }

  if (length(hypothesis) == 0) {
    return(names(post_samples))
  }

  fin_BFs <- vector('list',length = 0L)

  for (h in seq_along(hypothesis)) {
    resA <- eval_tidy(hypothesis[[h]],data = post_samples)
    res0 <- eval_tidy(hypothesis[[h]],data = prior_samples)

    if (is.logical(resA)) {
      # 2 directional BF
      Post.Prob <- mean(resA, na.rm = TRUE)
      Prior.Prob <- mean(res0, na.rm = TRUE)
      BF.v.Prior <- Post.Prob / Prior.Prob

      fin_BFs[[names(hypothesis)[h]]] <- data.frame(
        BF.v.Prior    = BF.v.Prior,
        BF.v.Null     = BF.v.Prior * unname(as.vector(BF)),
        BF.v.Opposite = Post.Prob/(1-Post.Prob),
        Post.Prob     = Post.Prob,
        row.names = names(hypothesis)[h]
      )
    } else {
      X0 <- na.omit(res0)
      XA <- na.omit(resA)

      HDI_int <- unname(hdi(XA,credMass = hdi.level))

      dp0 <- get_dfun(X0)(0)
      dpA <- get_dfun(XA)(0)
      dp0 <- ifelse(is.na(dp0),0,dp0)
      dpA <- ifelse(is.na(dpA),0,dpA)

      fin_BFs[[names(hypothesis)[h]]] <- data.frame(
        Estimate  = median(resA,na.rm = TRUE),
        Est.Error = sd(resA,na.rm = TRUE),
        HDI.Lower = HDI_int[1],
        HDI.Upper = HDI_int[2],
        BF        = dp0 / dpA,
        row.names = names(hypothesis)[h]
      )
    }
  }

  attr(fin_BFs,'class') <- c('BFEffect','inferBF','list')
  attr(fin_BFs,'model') <- BF@numerator[[1]]@longName
  attr(fin_BFs,'level') <- hdi.level
  attr(fin_BFs,'hyp') <- purrr::map_chr(hypothesis,quo_name)

  return(fin_BFs)
}