## ----setup, include=FALSE----------------------------------------------------- if (!(requireNamespace("rstanarm", quietly = TRUE) && requireNamespace("BayesFactor", quietly = TRUE) && requireNamespace("emmeans", quietly = TRUE) && requireNamespace("logspline", quietly = TRUE) && requireNamespace("lme4", quietly = TRUE) && requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("see", quietly = TRUE) )) { knitr::opts_chunk$set(eval = FALSE) } else { library(knitr) library(insight) library(bayestestR) library(rstanarm) library(BayesFactor) library(emmeans) library(ggplot2) library(see) options(knitr.kable.NA = '', digits = 2) opts_chunk$set(echo = TRUE, comment = ">", message = FALSE, warning = FALSE, dpi = 150) theme_set(theme_modern()) set.seed(4) } ## ----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------------------------------------------------ ggplot(sleep, aes(x = group, y = extra, fill= group)) + geom_boxplot() + theme_classic() ## ----rstanarm_model, eval = FALSE--------------------------------------------- # library(rstanarm) # # model <- stan_glm(extra ~ group, data = sleep, # prior = normal(0, 3, autoscale = FALSE)) ## ---- echo=FALSE-------------------------------------------------------------- model <- stan_glm(extra ~ group, data = sleep, prior = normal(0, 3, autoscale = FALSE), refresh = 0) ## ---- echo=FALSE-------------------------------------------------------------- null <- c(-1,1) xrange <- c(-10,10) x_vals <- seq(xrange[1], xrange[2], length.out = 400) d_vals <- dnorm(x_vals, sd = 3) in_null <- null[1] < x_vals & x_vals < null[2] range_groups <- rep(0, length(x_vals)) range_groups[!in_null & x_vals < 0] <- -1 range_groups[!in_null & x_vals > 0] <- 1 ggplot(mapping = aes(x_vals, d_vals, fill = in_null, group = range_groups)) + geom_area(color = "black", size = 1) + scale_fill_flat(name = "", labels = c("Alternative", "Null")) + 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_fit, echo=FALSE------------------------------------------------- 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) d_vals_post <- logspline::dlogspline(x_vals,f_post) ggplot(mapping = aes(x_vals, d_vals_post, fill = in_null, group = range_groups)) + geom_area(color = "black", size = 1) + scale_fill_flat(name = "", labels = c("Alternative", "Null")) + labs(x = "Drug effect", y = "Density") + theme_modern() + theme(legend.position = c(0.2, 0.8)) My_first_BF <- bayesfactor_parameters(model, null = c(-1,1), prior = model_prior) BF <- My_first_BF$BF[2] post_odds <- prior_odds * BF med_post <- point_estimate(posterior)$Median ## ---- eval=FALSE-------------------------------------------------------------- # My_first_BF <- bayesfactor_parameters(model, null = c(-1, 1)) # My_first_BF ## ---- echo=FALSE-------------------------------------------------------------- print(My_first_BF) ## ----------------------------------------------------------------------------- library(see) plot(My_first_BF) ## ---- eval=FALSE-------------------------------------------------------------- # My_second_BF <- bayesfactor_parameters(model, null = 0) # My_second_BF ## ---- echo=FALSE-------------------------------------------------------------- My_second_BF <- bayesfactor_parameters( data.frame(group2 = posterior), data.frame(group2 = prior), null = 0) print(My_second_BF) ## ----------------------------------------------------------------------------- plot(My_second_BF) ## ----savagedickey_one_sided, eval=FALSE--------------------------------------- # test_group2_right <- bayesfactor_parameters(model, direction = ">") # test_group2_right ## ----prior_n_post_plot_one_sided, echo=FALSE---------------------------------- test_group2_right <- bayesfactor_parameters( data.frame(group2 = posterior), data.frame(group2 = prior), null = 0, direction = ">" ) BF <- test_group2_right$BF print(test_group2_right) ## ----------------------------------------------------------------------------- plot(test_group2_right) ## ---- eval=FALSE-------------------------------------------------------------- # my_first_si <- si(model, BF = 1) # my_first_si ## ---- echo=FALSE-------------------------------------------------------------- my_first_si <- si(data.frame(group2 = posterior), data.frame(group2 = prior), BF = 1) print(my_first_si) ## ----------------------------------------------------------------------------- plot(my_first_si) ## ----brms_disp, eval=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-------------------------------------------- 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)", unsupported_models = FALSE ) print(comparison) ## ----update_models1----------------------------------------------------------- update(comparison, reference = 3) ## ----update_models2----------------------------------------------------------- update(comparison, reference = 2) ## ----lme4_models, eval=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) ## ---- echo=FALSE-------------------------------------------------------------- structure(list(Model = c( "Petal.Length + (1 | Species)", "Petal.Length + (Petal.Length | Species)", "Petal.Length + Petal.Width + (Petal.Length | Species)", "Petal.Length * Petal.Width + (Petal.Length | Species)", "1 + (1 | Species)"), BF = c(8.24027869011648e+24, 4.7677519818206e+23, 1.52492156042604e+22, 5.93045520305254e+20, 1)), class = c("bayesfactor_models", "see_bayesfactor_models", "data.frame"), row.names = c(NA, -5L), denominator = 5L, BF_method = "BIC approximation", unsupported_models = FALSE) ## ---- eval=FALSE-------------------------------------------------------------- # iris_model <- stan_glm(Sepal.Length ~ Species + Petal.Length, # data = iris, # prior = normal(0, c(2, 1.2, 1.2), autoscale = FALSE)) ## ---- echo=FALSE-------------------------------------------------------------- iris_model <- stan_glm(Sepal.Length ~ Species + Petal.Length, data = iris, prior = normal(0, c(2, 1.2, 1.2), autoscale = FALSE), refresh = 0) ## ----------------------------------------------------------------------------- botanist_hypotheses <- c( "Petal.Length > 0", "(Speciesversicolor > 0) & (Speciesvirginica > 0)" ) ## ---- eval=FALSE-------------------------------------------------------------- # botanist_BFs <- bayesfactor_restricted(iris_model, hypothesis = botanist_hypotheses) # botanist_BFs ## ---- echo=FALSE-------------------------------------------------------------- model_prior <- bayestestR:::.update_to_priors.stanreg(iris_model) botanist_BFs <- bayesfactor_restricted(iris_model, prior = model_prior, hypothesis = botanist_hypotheses) print(botanist_BFs) ## ----plot_iris, echo=FALSE---------------------------------------------------- ggplot(iris, aes(Petal.Length, Sepal.Length, color = Species)) + geom_point() + scale_color_flat() + theme(legend.position = c(0.2, 0.8)) ## ----inclusion_brms----------------------------------------------------------- bayesfactor_inclusion(comparison) ## ----inclusion_brms2---------------------------------------------------------- bayesfactor_inclusion(comparison, match_models = TRUE) ## ----JASP_all----------------------------------------------------------------- library(BayesFactor) data(ToothGrowth) ToothGrowth$dose <- as.factor(ToothGrowth$dose) BF_ToothGrowth <- anovaBF(len ~ dose*supp, ToothGrowth, progress = FALSE) bayesfactor_inclusion(BF_ToothGrowth) ## ----JASP_all_fig, echo=FALSE------------------------------------------------- knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP1.PNG") ## ----JASP_matched------------------------------------------------------------- bayesfactor_inclusion(BF_ToothGrowth, match_models = TRUE) ## ----JASP_matched_fig, echo=FALSE--------------------------------------------- knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP2.PNG") ## ----JASP_Nuisance------------------------------------------------------------ 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-------------------------------------------- knitr::include_graphics("https://github.com/easystats/easystats/raw/master/man/figures/bayestestR/JASP3.PNG") ## ---- eval=FALSE-------------------------------------------------------------- # mod <- stan_glm(mpg ~ wt + am, # data = mtcars, # prior = normal(0, c(10,10), autoscale = FALSE), # diagnostic_file = file.path(tempdir(), "df1.csv")) # # mod_carb <- stan_glm(mpg ~ wt + am + carb, # data = mtcars, # prior = normal(0, c(10,10,20), autoscale = FALSE), # diagnostic_file = file.path(tempdir(), "df0.csv")) # # bayesfactor_models(mod_carb, denominator = mod) ## ---- echo=FALSE-------------------------------------------------------------- mod <- stan_glm(mpg ~ wt + am, data = mtcars, prior = normal(0, c(10,10), autoscale = FALSE), diagnostic_file = file.path(tempdir(), "df1.csv"), refresh = 0) mod_carb <- stan_glm(mpg ~ wt + am + carb, data = mtcars, prior = normal(0, c(10,10,20), autoscale = FALSE), diagnostic_file = file.path(tempdir(), "df0.csv"), refresh = 0) BF_carb <- bayesfactor_models(mod_carb, denominator = mod, verbose = FALSE) BF <- BF_carb$BF[1] print(BF_carb) ## ----------------------------------------------------------------------------- hdi(mod_carb, ci = .95) ## ----------------------------------------------------------------------------- BMA_draws <- weighted_posteriors(mod, mod_carb) BMA_hdi <- hdi(BMA_draws, ci = .95) BMA_hdi plot(BMA_hdi) ## ---- echo=FALSE-------------------------------------------------------------- set.seed(1) ## ----------------------------------------------------------------------------- library(emmeans) groups <- emmeans(model, ~ group) group_diff <- pairs(groups) (groups_all <- rbind(groups, group_diff)) # pass the original model via prior bayesfactor_parameters(groups_all, prior = model) ## ---- echo=FALSE-------------------------------------------------------------- set.seed(1) ## ----------------------------------------------------------------------------- library(modelbased) estimate_contrasts(model, test = "bf") ## ---- eval=FALSE-------------------------------------------------------------- # contrasts(iris$Species) <- contr.sum # # fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris, # prior = normal(0, c(1, 1), autoscale = FALSE), # prior_PD = TRUE, # sample priors # family = gaussian()) # # pairs_sum <- pairs(emmeans(fit_sum, ~ Species)) # pairs_sum ## ---- echo=FALSE-------------------------------------------------------------- contrasts(iris$Species) <- contr.sum fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris, prior = normal(0, c(1, 1), autoscale = FALSE), prior_PD = TRUE, # sample priors family = gaussian(), refresh = 0) pairs_sum <- pairs(emmeans(fit_sum, ~ Species)) em_pairs_samples <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(pairs_sum, names = FALSE))) print(pairs_sum) ggplot(stack(em_pairs_samples), aes(x = values, fill = ind)) + geom_density(size = 1) + facet_grid(ind ~ .) + labs(x = "prior difference values") + theme(legend.position = "none") ## ----------------------------------------------------------------------------- contrasts(iris$Species) <- contr.bayes ## ---- eval=FALSE-------------------------------------------------------------- # options(contrasts = c('contr.bayes', 'contr.poly')) ## ---- eval=FALSE-------------------------------------------------------------- # contrasts(iris$Species) <- contr.bayes # # fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris, # prior = normal(0, c(1, 1), autoscale = FALSE), # prior_PD = TRUE, # sample priors # family = gaussian()) # # pairs_bayes <- pairs(emmeans(fit_bayes, ~ Species)) # pairs_bayes ## ---- echo=FALSE-------------------------------------------------------------- contrasts(iris$Species) <- contr.bayes fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris, prior = normal(0, c(1, 1), autoscale = FALSE), prior_PD = TRUE, # sample priors family = gaussian(), refresh = 0) pairs_bayes <- pairs(emmeans(fit_bayes, ~ Species)) em_pairs_samples <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(pairs_bayes, names = FALSE))) print(pairs_bayes) ggplot(stack(em_pairs_samples), aes(x = values, fill = ind)) + geom_density(size = 1) + facet_grid(ind ~ .) + labs(x = "prior difference values") + theme(legend.position = "none") ## ----------------------------------------------------------------------------- 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" ) ## ---- eval=FALSE-------------------------------------------------------------- # contrasts(iris$Species) <- contr.sum # # fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris, # prior = normal(0, c(1, 1), autoscale = FALSE), # family = gaussian()) # # em_sum <- emmeans(fit_sum, ~ Species) # the posterior marginal means # # bayesfactor_restricted(em_sum, fit_sum, hypothesis = hyp) ## ---- echo=FALSE-------------------------------------------------------------- contrasts(iris$Species) <- contr.sum fit_sum <- stan_glm(Sepal.Length ~ Species, data = iris, prior = normal(0, c(1, 1), autoscale = FALSE), family = gaussian(), refresh = 0) em_sum <- emmeans(fit_sum, ~ Species) # the posterior marginal means bayesfactor_restricted(em_sum, fit_sum, hypothesis = hyp) ## ---- eval=FALSE-------------------------------------------------------------- # contrasts(iris$Species) <- contr.bayes # # fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris, # prior = normal(0, c(1, 1), autoscale = FALSE), # family = gaussian()) # # em_bayes <- emmeans(fit_sum, ~ Species) # the posterior marginal means # # bayesfactor_restricted(em_bayes, fit_sum, hypothesis = hyp) ## ---- echo=FALSE-------------------------------------------------------------- contrasts(iris$Species) <- contr.bayes fit_bayes <- stan_glm(Sepal.Length ~ Species, data = iris, prior = normal(0, c(1, 1), autoscale = FALSE), family = gaussian(), refresh = 0) em_bayes <- emmeans(fit_bayes, ~ Species) # the posterior marginal means bayesfactor_restricted(em_bayes, fit_bayes, hypothesis = hyp)