bayes_factors.R
## ----setup, include=FALSE------------------------------------------------
library(knitr)
library(rstanarm)
library(bayestestR)
library(ggplot2)
library(see)
library(emmeans)
library(lme4)
library(BayesFactor)
options(knitr.kable.NA = '')
opts_chunk$set(echo = TRUE)
opts_chunk$set(comment = ">")
knitr::opts_chunk$set(dpi=300)
theme_set(see::theme_modern())
options(digits = 2)
set.seed(5)
## ----deathsticks_fig, echo=FALSE, fig.cap="Bayesian analysis of the Students' (1908) Sleep data set.", fig.align='center', out.width="80%"----
knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/deathsticks.jpg")
## ----sleep_boxplot, echo=FALSE, message=FALSE, warning=FALSE-------------
ggplot(sleep, aes(x = group, y = extra, fill= group)) +
geom_boxplot() +
theme_classic()
## ---- echo=FALSE---------------------------------------------------------
null <- c(-.5,.5)
xrange <- c(-12,12)
ggplot() + aes(x = 0, y = 0) +
stat_function(aes(fill = "Null"),
fun = dnorm, args = list(sd = 2.5),
xlim = null, geom = "area") +
stat_function(aes(fill = "Alternative"),
fun = dnorm, args = list(sd = 2.5),
xlim = c(xrange[1],null[1]), geom = "area") +
stat_function(aes(fill = "Alternative"),
fun = dnorm, args = list(sd = 2.5),
xlim = c(null[2],xrange[2]), geom = "area") +
stat_function(fun = dnorm, args = list(sd = 2.5),
xlim = xrange, size = 1) +
scale_fill_flat(name = "") +
labs(x = "Drug effect", y = "Density") +
theme_modern() +
theme(legend.position = c(0.2, 0.8))
pnull <- diff(pnorm(null, sd = 2.5))
prior_odds <- (1 - pnull) / pnull
## ----rstanarm_disp, eval=FALSE, message=FALSE, warning=FALSE-------------
# library(rstanarm)
# model <- stan_glm(extra ~ group, data = sleep)
## ----rstanarm_fit, echo=FALSE, message=FALSE, warning=FALSE--------------
junk <- capture.output(model <- stan_glm(extra ~ group, data = sleep))
model_prior <- bayestestR:::.update_to_priors.stanreg(model)
posterior <- insight::get_parameters(model)$group2
prior <- insight::get_parameters(model_prior)$group2
f_post <- logspline::logspline(posterior)
dpost <- function(q){
logspline::dlogspline(q,f_post)
}
xrange <- c(-12,12)
ggplot() + aes(x = 0, y = 0) +
stat_function(aes(fill = "Null"),
fun = dpost,
xlim = null, geom = "area") +
stat_function(aes(fill = "Alternative"),
fun = dpost,
xlim = c(xrange[1],null[1]), geom = "area") +
stat_function(aes(fill = "Alternative"),
fun = dpost,
xlim = c(null[2],xrange[2]), geom = "area") +
stat_function(fun = dpost,
xlim = xrange, size = 1) +
scale_fill_flat(name = "") +
geom_vline(xintercept = point_estimate(posterior)$Median, size = 1, linetype = "dashed") +
labs(x = "Drug effect", y = "Density") +
theme_modern() +
theme(legend.position = c(0.2, 0.8))
My_first_BF <- bayesfactor_parameters(model, null = c(-.5,.5))
BF <- My_first_BF$BF[2]
post_odds <- prior_odds * BF
## ---- echo=TRUE, eval=FALSE----------------------------------------------
# My_first_BF <- bayesfactor_parameters(model, null = c(-.5,.5))
# My_first_BF
## ---- echo=FALSE, message=FALSE, warning=FALSE---------------------------
print(My_first_BF)
## ---- echo=TRUE, eval=FALSE----------------------------------------------
# plot(My_first_BF)
## ---- echo=FALSE, message=FALSE, warning=FALSE---------------------------
plot(bayesfactor_parameters(
data.frame(group2 = posterior),
data.frame(group2 = prior),
null = c(-.5,.5))) +
theme(legend.position = c(0.2, 0.8))
## ---- message=FALSE, warning=FALSE---------------------------------------
My_second_BF <- bayesfactor_parameters(model, null = 0)
My_second_BF
## ---- echo=TRUE, eval=FALSE----------------------------------------------
# plot(My_second_BF)
## ---- echo=FALSE, message=FALSE, warning=FALSE---------------------------
plot(bayesfactor_parameters(
data.frame(group2 = posterior),
data.frame(group2 = prior),
null = 0)) +
theme(legend.position = c(0.2, 0.8))
## ----prior_n_post_plot_one_sided, echo=FALSE, message=FALSE, warning=FALSE----
# Using "see"
bfsd <- bayesfactor_parameters(
data.frame(group2 = posterior),
data.frame(group2 = prior),
direction = ">"
)
plot(bfsd) +
theme(legend.position = c(0.8,0.8))
## ----savagedickey_one_sided, message=FALSE, warning=FALSE----------------
test_group2_right <- bayesfactor_parameters(model, direction = ">")
test_group2_right
## ------------------------------------------------------------------------
library(emmeans)
group_diff <- pairs(emmeans(model, ~ group))
group_diff
# pass the original model via prior
bayesfactor_parameters(group_diff, prior = model)
## ----brms_disp, eval=FALSE, message=FALSE, warning=FALSE-----------------
# library(brms)
#
# m0 <- brm(Sepal.Length ~ 1, # intercept only model
# data = iris, save_all_pars = TRUE)
# m1 <- brm(Sepal.Length ~ Petal.Length,
# data = iris, save_all_pars = TRUE)
# m2 <- brm(Sepal.Length ~ Species,
# data = iris, save_all_pars = TRUE)
# m3 <- brm(Sepal.Length ~ Species + Petal.Length,
# data = iris, save_all_pars = TRUE)
# m4 <- brm(Sepal.Length ~ Species * Petal.Length,
# data = iris, save_all_pars = TRUE)
## ----brms_models_disp, eval=FALSE----------------------------------------
# library(bayestestR)
# comparison <- bayesfactor_models(m1, m2, m3, m4, denominator = m0)
# comparison
## ----brms_models_print, echo=FALSE, message=FALSE, warning=FALSE---------
# dput(comparison)
comparison <- structure(
list(
Model = c(
"Petal.Length",
"Species",
"Species + Petal.Length",
"Species * Petal.Length",
"1"
),
BF = c(3.44736e+44, 5.628679e+29, 7.121386e+55, 9.149948e+55, 1)
),
class = c("bayesfactor_models", "see_bayesfactor_models", "data.frame"),
row.names = c(NA, -5L),
denominator = 5L,
BF_method = "marginal likelihoods (bridgesampling)"
)
comparison
## ----update_models1, message=FALSE, warning=FALSE------------------------
update(comparison, reference = 3)
## ----update_models2, message=FALSE, warning=FALSE------------------------
update(comparison, reference = 2)
## ----lme4_models, message=FALSE, warning=FALSE---------------------------
library(lme4)
m0 <- lmer(Sepal.Length ~ (1 | Species), data = iris)
m1 <- lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
m2 <- lmer(Sepal.Length ~ Petal.Length + (Petal.Length | Species), data = iris)
m3 <- lmer(Sepal.Length ~ Petal.Length + Petal.Width + (Petal.Length | Species), data = iris)
m4 <- lmer(Sepal.Length ~ Petal.Length * Petal.Width + (Petal.Length | Species), data = iris)
bayesfactor_models(m1, m2, m3, m4, denominator = m0)
## ----inclusion_brms, message=FALSE, warning=FALSE, eval=TRUE-------------
bayesfactor_inclusion(comparison)
## ----inclusion_brms2, message=FALSE, warning=FALSE, eval=TRUE------------
bayesfactor_inclusion(comparison, match_models = TRUE)
## ----JASP_all, message=FALSE, warning=FALSE, eval=TRUE-------------------
library(BayesFactor)
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
BF_ToothGrowth <- anovaBF(len ~ dose*supp, ToothGrowth)
bayesfactor_inclusion(BF_ToothGrowth)
## ----JASP_all_fig, echo=FALSE, message=FALSE, warning=FALSE--------------
knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP1.PNG")
## ----JASP_matched, message=FALSE, warning=FALSE, eval=TRUE---------------
bayesfactor_inclusion(BF_ToothGrowth, match_models = TRUE)
## ----JASP_matched_fig, echo=FALSE, message=FALSE, warning=FALSE----------
knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP2.PNG")
## ----JASP_Nuisance, message=FALSE, warning=FALSE, eval=TRUE--------------
BF_ToothGrowth_against_dose <- BF_ToothGrowth[3:4]/BF_ToothGrowth[2] # OR:
# update(bayesfactor_models(BF_ToothGrowth),
# subset = c(4, 5),
# reference = 3)
BF_ToothGrowth_against_dose
bayesfactor_inclusion(BF_ToothGrowth_against_dose)
## ----JASP_Nuisance_fig, echo=FALSE, message=FALSE, warning=FALSE---------
knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP3.PNG")
## ----plot_iris, echo=FALSE, message=FALSE, warning=FALSE-----------------
ggplot(iris, aes(Petal.Length, Sepal.Length, color = Species)) +
geom_point() +
scale_color_flat() +
theme(legend.position = c(0.2, 0.8))
## ---- message=FALSE, warning=FALSE, eval = FALSE-------------------------
# iris_model <- stan_glm(Sepal.Length ~ Species + Petal.Length,
# data = iris)
## ---- echo=FALSE, message=FALSE, warning=FALSE---------------------------
junk <- capture.output(iris_model <- stan_glm(Sepal.Length ~ Species + Petal.Length,
data = iris, refresh = 0))
model_prior <- bayestestR:::.update_to_priors.stanreg(iris_model)
priors <- insight::get_parameters(model_prior)
priors$`(Intercept)` <- NULL
ggplot(stack(priors), aes(values, fill = ind)) +
geom_density(color = NA) +
geom_vline(xintercept = 0, size = 1, linetype = "dashed") +
facet_grid(~ind) +
scale_fill_flat() +
theme(legend.position = "none")
# describe_posterior(iris_model)
## ------------------------------------------------------------------------
botanist_hypotheses <- c(
"Petal.Length > 0",
"(Speciesversicolor > 0) & (Speciesvirginica > 0)"
)
## ------------------------------------------------------------------------
botanist_BFs <- bayesfactor_restricted(iris_model, hypothesis = botanist_hypotheses)
botanist_BFs
## ---- eval=FALSE---------------------------------------------------------
# library(bayestestR)
# library(rstanarm)
# library(emmeans)
# library(ggplot2)
#
# contrasts(iris$Species) <- contr.sum
#
# fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris,
# prior = normal(0,0.1), # just to drive the point home
# family = gaussian())
# c_sum <- pairs(emmeans(fit_sum, ~ Species))
# c_sum
## ---- warning=FALSE, echo=FALSE------------------------------------------
contrasts(iris$Species) <- contr.sum
set.seed(5)
junk <- capture.output(
fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris,
# just to drive the point home, we'll use ultra-narrow priors
# (probably should not be used)
prior = normal(0,0.1),
family = gaussian())
)
c_sum <- pairs(emmeans(fit_sum, ~ Species))
c_sum
## ---- message=FALSE, echo=FALSE------------------------------------------
plot(bayesfactor_parameters(c_sum, fit_sum)) +
coord_cartesian(xlim = c(-2,1))
## ------------------------------------------------------------------------
contrasts(iris$Species) <- contr.bayes
## ---- eval=FALSE---------------------------------------------------------
# options(contrasts = c('contr.bayes', 'contr.poly'))
## ---- eval=FALSE---------------------------------------------------------
# fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris,
# prior = normal(0,0.1),
# family = gaussian())
# c_bayes <- pairs(emmeans(fit_bayes, ~ Species))
# c_bayes
## ---- warning=FALSE, echo=FALSE------------------------------------------
set.seed(5)
junk <- capture.output(
fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris,
prior = normal(0,0.1), # just to drive the point home
family = gaussian())
)
c_bayes <- pairs(emmeans(fit_bayes, ~ Species))
c_bayes
## ---- message=FALSE, echo=FALSE------------------------------------------
plot(bayesfactor_parameters(c_bayes, fit_bayes)) +
coord_cartesian(xlim = c(-2,1))
## ------------------------------------------------------------------------
hyp <- c(
# comparing 2 levels
"setosa > versicolor",
"setosa > virginica",
"versicolor > virginica",
# comparing 3 (or more) levels
"setosa < virginica & virginica < versicolor",
"virginica < setosa & setosa < versicolor",
"setosa < versicolor & versicolor < virginica"
)
## ------------------------------------------------------------------------
em_sum <- emmeans(fit_sum, ~Species)
em_sum # the posterior marginal means
bayesfactor_restricted(em_sum, fit_sum, hypothesis = hyp)
## ------------------------------------------------------------------------
em_bayes <- emmeans(fit_bayes, ~Species)
em_bayes # the posterior marginal means
bayesfactor_restricted(em_bayes, fit_bayes, hypothesis = hyp)