set.seed(444)
library(bayestestR)
prior <- data.frame(
A = rnorm(500),
B = rnorm(500),
C = rnorm(500)
)
posterior <- data.frame(
A = rnorm(500, .4, 0.7),
B = rnorm(500, -.2, 0.4),
C = rnorm(500, 0, 0.5)
)
hyps <- c(
"A > B & B > C",
"A > B & A > C",
"C > A"
)
(b <- bayesfactor_restricted(posterior, hypothesis = hyps, prior = prior))
bool <- as.logical(b, which = "posterior")
head(bool)
if (FALSE) { # require("see") && require("patchwork")
see::plots(
plot(estimate_density(posterior)),
# distribution **conditional** on the restrictions
plot(estimate_density(posterior[bool[, hyps[1]], ])) + ggplot2::ggtitle(hyps[1]),
plot(estimate_density(posterior[bool[, hyps[2]], ])) + ggplot2::ggtitle(hyps[2]),
plot(estimate_density(posterior[bool[, hyps[3]], ])) + ggplot2::ggtitle(hyps[3]),
guides = "collect"
)
}
if (FALSE) { # require("rstanarm")
# \donttest{
# rstanarm models
# ---------------
data("mtcars")
fit_stan <- rstanarm::stan_glm(mpg ~ wt + cyl + am,
data = mtcars, refresh = 0
)
hyps <- c(
"am > 0 & cyl < 0",
"cyl < 0",
"wt - cyl > 0"
)
bayesfactor_restricted(fit_stan, hypothesis = hyps)
# }
}
if (FALSE) { # require("rstanarm") && require("emmeans")
# \donttest{
# emmGrid objects
# ---------------
# replicating http://bayesfactor.blogspot.com/2015/01/multiple-comparisons-with-bayesfactor-2.html
data("disgust")
contrasts(disgust$condition) <- contr.equalprior_pairs # see vignette
fit_model <- rstanarm::stan_glm(score ~ condition, data = disgust, family = gaussian())
em_condition <- emmeans::emmeans(fit_model, ~condition, data = disgust)
hyps <- c("lemon < control & control < sulfur")
bayesfactor_restricted(em_condition, prior = fit_model, hypothesis = hyps)
# > # Bayes Factor (Order-Restriction)
# >
# > Hypothesis P(Prior) P(Posterior) BF
# > lemon < control & control < sulfur 0.17 0.75 4.49
# > ---
# > Bayes factors for the restricted model vs. the un-restricted model.
# }
}
Run the code above in your browser using DataLab