MTCARS evaluation

rm(list = ls())
library(magrittr)
library(dplyr)
library(forcats)
library(modelr)
library(ggdist)
library(tidybayes)
library(ggplot2)
library(cowplot)
library(emmeans)
library(broom)
library(rstan)
library(rstanarm)
library(brms)
library(bayesplot)
library(MCMCglmm)
library(RColorBrewer)
library(ggpubr)
library(rstatix)
library(grid)
library(gridExtra)
require(bayesrules)
library(gganimate)
theme_set(theme_tidybayes() + panel_border())

mtcars$cyl_fc<-as.factor(mtcars$cyl)
ggplot(mtcars, aes(x=hp, 
                   y=mpg, 
                   col=cyl_fc))+
  geom_point(size=4)+
  labs(col="cyl")+
  theme_bw()
ggplot(mtcars, aes(x=hp, 
                   y=qsec, 
                   col=cyl_fc))+
  geom_point(size=4)+
  labs(col="cyl")+
  facet_wrap(~gear)+
  theme_bw()

get_prior(qsec ~ cyl + disp + drat + wt + hp + gear + carb + vs + am, 
          data = mtcars,
          family = gaussian())
prior1<-c(prior_string("student_t(3, 17.7, 2.5)", class = "Intercept"),
          prior_string("student_t(3, 0, 2.5)", class = "sigma"))
summary(mtcars$qsec)
fit <- brm(qsec ~ cyl + disp + drat + wt + hp + gear + carb + vs + am, 
           data = mtcars,
           prior = c(prior(normal(5, 2),      class = Intercept),
                     prior(normal(3, 0.5),     class = b, coef=cyl),
                     prior(normal(2, 0.2),     class = b, coef=hp),
                     prior(normal(1.2, 0.3),   class = b, coef=am),
                     prior(normal(0.8, 0.2),   class = b, coef=carb),
                     prior(normal(7.2, 5.6),   class = b, coef=disp),
                     prior(normal(3.5, 2.8),   class = b, coef=drat),
                     prior(normal(1.5, 0.3),   class = b, coef=gear),
                     prior(normal(0.6, 0.1),   class = b, coef=vs),
                     prior(normal(2.8, 0.4),   class = b, coef=wt),
                     prior(gamma(2, 1),        class = sigma)),
           family = gaussian())
fit$prior
summary(fit)
plot(fit)

pp_check(fit, ndraws=100)
pp_check(fit, type = "error_hist", ndraws = 11)
pp_check(fit, type = "scatter_avg", ndraws = 100)
pp_check(fit, type = "stat_2d")
pp_check(fit, type = "loo_pit")

yrep<-posterior_predict(fit)
y<-mtcars$qsec
group<-mtcars$cyl_fc
x<-mtcars$hp
ppc_dens_overlay(y, yrep[1:25,])
ppc_ecdf_overlay(y, yrep[sample(nrow(yrep), 25), ])
ppc_hist(y, yrep[1:8, ])
ppc_boxplot(y, yrep[1:8, ])
ppc_dens(y, yrep[200:202, ])
ppc_freqpoly(y, yrep[1:3,], alpha = 0.1, size = 1, binwidth = 5)
ppc_freqpoly_grouped(y, yrep[1:3,], group) + yaxis_text()
ppc_dens_overlay_grouped(y, yrep[1:25, ], group = group)
ppc_ecdf_overlay_grouped(y, yrep[1:25, ], group = group)
ppc_violin_grouped(y, yrep, group, size = 1.5)
ppc_violin_grouped(y, yrep, group, alpha = 0, y_draw = "points", y_size = 1.5)
ppc_violin_grouped(y, yrep, group, alpha = 0, y_draw = "both",
                   y_size = 1.5, y_alpha = 0.5, y_jitter = 0.33)

ppc_stat(y, yrep)
q25 <- function(y) quantile(y, 0.25)
ppc_stat(y, yrep, stat = "q25")
ppc_stat(y, yrep, stat = function(y) quantile(y, 0.25))
ppc_stat_grouped(y, yrep, group)
ppc_stat_grouped(y, yrep, group) + yaxis_text()

bayesplot_theme_set(ggplot2::theme_linedraw())
color_scheme_set("viridisE")
ppc_stat_2d(y, yrep, stat = c("mean", "sd"))

bayesplot_theme_set(ggplot2::theme_grey())
color_scheme_set("brewer-Paired")
ppc_stat_2d(y, yrep, stat = c("median", "mad"))

theme_set(theme_bw())
color_scheme_set("brightblue")
ppc_intervals(y, yrep)
ppc_intervals(y, yrep, size = 1.5, fatten = 0)
ppc_ribbon(y, yrep)
ppc_ribbon(y, yrep, y_draw = "points")
ppc_ribbon(y, yrep, y_draw = "both")

ppc_ribbon_grouped(y, yrep, x = x, group, y_draw = "both") +
  ggplot2::scale_x_continuous(breaks = pretty)

color_scheme_set("purple")
ppc_intervals(y = y, yrep = yrep, x = x, prob = 0.8) +
  panel_bg(fill="gray90", color = NA) +
  grid_lines(color = "white")

color_scheme_set("gray")
ppc_intervals(y, yrep, prob = 0.5) +
  ggplot2::scale_x_continuous(
    labels = rownames(mtcars),
    breaks = 1:nrow(mtcars)) +
  xaxis_text(angle = -70, vjust = 1, hjust = 0) +
  xaxis_title(FALSE)

ppc_error_hist(y, yrep[1:3, ])
ppc_error_hist_grouped(y, yrep[1:3, ], group)
ppc_error_scatter(y, yrep[10:14, ])
ppc_error_scatter_avg(y, yrep)
ppc_scatter_avg_grouped(y, yrep, group, facet_args = list(scales = "free_x"))

loo1 <- loo(fit, save_psis = TRUE, cores = 4)
psis1 <- loo1$psis_object
lw <- weights(psis1)
ppc_loo_pit_overlay(y, yrep, lw = lw)
ppc_loo_pit_qq(y, yrep, lw = lw)
ppc_loo_pit_qq(y, yrep, lw = lw, compare = "normal")
keep_obs <- 1:30
ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs)
ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs,
                  order = "median")