#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.