MTCARS

#### MTCARS ####
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)
library(ggstatsplot)
theme_set(theme_tidybayes() + panel_border())

data(mtcars)
dim
skimr::skim(mtcars)
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()
ggscatterstats(
  data  = mtcars,
  x     = hp,
  y     = mpg,
  xlab  = "Horsepower",
  ylab  = "Miles per gallon (MPG)",
  title = "Relationship between miles per gallon and the amount of horsepower of a car")

mtcars$cyl_fac<-as.factor(mtcars$cyl)
grouped_ggscatterstats(
  data  = mtcars,
  x     = hp,
  y     = mpg,
  grouping.var     = cyl_fac,
  xlab  = "Horsepower",
  ylab  = "Miles per gallon (MPG)",
  ggtheme          = ggplot2::theme_grey(),
  plotgrid.args    = list(nrow = 1),
  annotation.args  = list(title = "Relationship between miles per gallon and the amount of horsepower of a car across number of cylinders"))

ggcorrmat(
  data     = mtcars,
  title    = "Correlalogram for mtcars dataset")
grouped_ggcorrmat(
  data         = mtcars,
  type         = "robust",
  grouping.var = cyl_fac,
  matrix.type  = "lower")

get_prior(mpg~hp*cyl, data=mtcars)
m_brm = brm(
  mpg ~ hp * cyl, 
  data = mtcars)
summary(m_brm)
plot(m_brm)
get_variables(m_brm)
ndraws = 4000
m_brm %>%
  recover_types(mtcars)%>%
  spread_draws(b_hp, b_cyl)%>%
  ggplot(., aes(x=b_hp, fill=as.factor(.chain)))+
  geom_density(alpha=0.5)+
  labs(fill="Chain")+
  theme_bw()
m_brm %>%
  recover_types(mtcars)%>%
  spread_draws(b_hp, b_cyl)%>%
  ggplot(., aes(x=b_hp, 
                y=b_cyl, 
                color=as.factor(.chain)))+
  geom_point(alpha=0.7)+
  labs(color="Chain")+
  theme_bw()
m_brm %>%
  recover_types(mtcars)%>%
  spread_draws(b_hp, b_Intercept)%>%
  ggplot(., aes(x=b_hp, 
                y=b_Intercept, 
                color=as.factor(.chain)))+
  geom_point(alpha=0.7)+
  labs(color="Chain")+
  theme_bw()
m_brm %>%
  recover_types(mtcars)%>%
  spread_draws(b_hp, sigma)%>%
  ggplot(., aes(x=b_hp, 
                y=sigma, 
                color=as.factor(.chain)))+
  geom_point(alpha=0.7)+
  labs(color="Chain")+
  theme_bw()

summary(mtcars$mpg)
get_prior(mpg~hp*cyl, data=mtcars)
m_brm = brm(
  mpg ~ hp * cyl, 
  data = mtcars, 
  prior = c(prior(normal(15, 4),      class = Intercept),
            prior(normal(3, 0.5),     class = b, coef=cyl),
            prior(normal(2, 0.2),     class = b, coef=hp),
            prior(normal(1.7, 0.6),   class = b, coef=hp:cyl),
            prior(gamma(4, 1),        class = sigma)),
  chains=4, 
  cores=6)
m_brm$prior
summary(m_brm)
plot(m_brm)
pp_check(m_brm, ndraws=100)


ndraws = 50
p = mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 101)) %>%
  add_epred_draws(m_brm, ndraws = ndraws) %>%
  ggplot(aes(x = hp, y = mpg, color = ordered(cyl))) +
  geom_line(aes(y = .epred, group = paste(cyl, .draw))) +
  geom_point(data = mtcars) +
  scale_color_brewer(palette = "Dark2") +
  transition_states(.draw, 0, 1) +
  shadow_mark(future = TRUE, color = "gray50", alpha = 1/20)
animate(p, nframes = ndraws, fps = 2.5, 
        width = 800, height = 500, res = 150, dev = "png", type = "cairo")
getwd()
anim_save("mtcars.gif")

mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 101)) %>%
  add_predicted_draws(m_brm) %>%
  ggplot(aes(x = hp, y = mpg, color = ordered(cyl), fill = ordered(cyl))) +
  stat_lineribbon(aes(y = .prediction), .width = c(.95, .80, .50), alpha = 1/4) +
  geom_point(data = mtcars) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Dark2")
mtcars %>%
  group_by(cyl) %>%
  data_grid(hp = seq_range(hp, n = 101)) %>%
  add_predicted_draws(m_brm) %>%
  ggplot(aes(x = hp, y = mpg)) +
  stat_lineribbon(aes(y = .prediction), .width = c(.99, .95, .8, .5), color = brewer.pal(5, "Blues")[[5]]) +
  geom_point(data = mtcars) +
  scale_fill_brewer() +
  facet_grid(. ~ cyl, space = "free_x", scales = "free_x")


mtcars_clean = mtcars %>%
  mutate(cyl = ordered(cyl))
head(mtcars_clean)
get_prior(cyl ~ mpg, 
          data = mtcars_clean, 
          family = cumulative)
summary(mtcars$cyl)
m_cyl = brm(
  cyl ~ mpg, 
  data = mtcars_clean, 
  family = cumulative,
  seed = 58393,
  prior=c(prior(normal(1, 0.5),      class = Intercept, coef=1),
          prior(normal(-1,0.5),      class = Intercept, coef=2),
          prior(normal(1.7, 0.6),   class = b, coef=mpg)),
  chains=4, 
  cores=6)
plot(m_cyl)
summary(m_cyl)
pp_check(m_cyl, ndraws=100)

tibble(mpg = 21) %>%
  add_epred_draws(m_cyl) %>%
  median_qi(.epred)
data_plot = mtcars_clean %>%
  ggplot(aes(x = mpg, y = cyl, color = cyl)) +
  geom_point() +
  scale_color_brewer(palette = "Dark2", name = "cyl")
fit_plot = mtcars_clean %>%
  data_grid(mpg = seq_range(mpg, n = 101)) %>%
  add_epred_draws(m_cyl, value = "P(cyl | mpg)", category = "cyl") %>%
  ggplot(aes(x = mpg, y = `P(cyl | mpg)`, color = cyl)) +
  stat_lineribbon(aes(fill = cyl), alpha = 1/5) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2")
plot_grid(ncol = 1, align = "v",
          data_plot,
          fit_plot)

ndraws = 100
p = mtcars_clean %>%
  data_grid(mpg = seq_range(mpg, n = 101)) %>%
  add_epred_draws(m_cyl, value = "P(cyl | mpg)", category = "cyl") %>%
  ggplot(aes(x = mpg, y = `P(cyl | mpg)`, color = cyl)) +
  # we remove the `.draw` column from the data for stat_lineribbon so that the same ribbons
  # are drawn on every frame (since we use .draw to determine the transitions below)
  stat_lineribbon(aes(fill = cyl), alpha = 1/5, color = NA, data = . %>% select(-.draw)) +
  # we use sample_draws to subsample at the level of geom_line (rather than for the full dataset
  # as in previous HOPs examples) because we need the full set of draws for stat_lineribbon above
  geom_line(aes(group = paste(.draw, cyl)), size = 1, data = . %>% sample_draws(ndraws)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  transition_manual(.draw)
animate(p, nframes = ndraws, fps = 2.5, 
        width = 800, height = 500, res = 150, dev = "png", type = "cairo")
getwd()
anim_save("mtcars2.gif")

tibble(mpg = 20) %>%
  add_epred_draws(m_cyl, value = "P(cyl | mpg = 20)", category = "cyl") %>%
  ungroup() %>%
  select(.draw, cyl, `P(cyl | mpg = 20)`) %>%
  gather_pairs(cyl, `P(cyl | mpg = 20)`, triangle = "both") %>%
  filter(.row != .col) %>%
  ggplot(aes(.x, .y)) +
  geom_point(alpha = 1/50) +
  facet_grid(.row ~ .col) +
  ylab("P(cyl = row | mpg = 20)") +
  xlab("P(cyl = col | mpg = 20)")