Iris

rm(list = ls())
library(rethinking)
library(dplyr)
library(ggplot2)
library(grid)
library(gridExtra)
library(GGally)
library(MVN)
library(data.table)
library(mltools)
library(brms)
library(tidybayes)
library(bayesplot)

dev.off()
skimr::skim(iris)
ggplot(iris, aes(y= Sepal.Length, x = Petal.Length, 
                 colour=Species, 
                 size=Sepal.Width)) + 
  geom_point(alpha=0.8)+
  theme_bw()
model1<-lm(Sepal.Length~Petal.Length + Species, data=iris)
summary(model1)
par(mfrow = c(2, 2))
plot(model1)

model2<- rstanarm::stan_lm(
  Sepal.Length~Petal.Length + Species, data=iris,
  prior = rstanarm::R2(0.75))
rstanarm::R2(0.75)
summary(model2)
plot(model2)
bayesplot::color_scheme_set("red")
bayesplot::mcmc_areas(model2, prob = 0.8,area_method = "equal area") + 
  geom_vline(xintercept = 0)
neff_ratio(model2)
rhat(model2)
bayesplot::mcmc_trace(model2, size = 0.1)
bayesplot::mcmc_dens_overlay(model2)
bayesplot::mcmc_acf(model2)
bayesplot::pp_check(model2)

iris %>%
  tidybayes::add_fitted_draws(model2, n = 100) %>%
  ggplot(aes(x = Petal.Length, 
             y = Sepal.Width)) +
  geom_line(
    aes(y = .value, 
        group = paste(Species, .draw), 
        color = Species),
    alpha = 0.1) +
  geom_point(aes(color = Species)) + theme_bw()
predictions <- posterior_predict(model2, newdata = iris)
bayesplot::ppc_intervals(iris$Sepal.Width, 
                         yrep = predictions, 
                         x = iris$Petal.Length, 
                         prob = 0.5, prob_outer = 0.95)
bayesplot::ppc_intervals_grouped(iris$Sepal.Width, 
                                 yrep = predictions, 
                                 x = iris$Petal.Length,
                                 group=iris$Species,
                                 prob = 0.5, prob_outer = 0.95)
bayesplot::ppc_ribbon_grouped(iris$Sepal.Width, 
                              yrep = predictions, 
                              x = iris$Petal.Length,
                              group=iris$Species,
                              prob = 0.5, prob_outer = 0.95)
bayesrules::prediction_summary(model2, data = iris)
cv_procedure <- bayesrules::prediction_summary_cv(model = model2,
                                                  group = "Species",
                                                  data  = iris, 
                                                  k = 10)
cv_procedure$cv
model_elpd <- loo(model2, save_psis = TRUE)
print(model_elpd)
model_elpd$estimates
plot(model_elpd)
yrep <- posterior_predict(model2)
bayesplot::ppc_loo_pit_overlay(
  y = iris$Sepal.Width,
  yrep = yrep,
  lw = weights(model_elpd$psis_object)
)
bayesplot::ppc_loo_pit_qq(
  y = iris$Sepal.Width,
  yrep = yrep,
  lw = weights(model_elpd$psis_object)
)

posterior_interval(model2, prob = 0.80)
iris %>% 
  group_by(Species)%>%
  tidybayes::add_epred_draws(model2, ndraws = 20) %>%
  ggplot(aes(x = Petal.Length, y = Sepal.Width)) +
  geom_point(aes(color=Species), alpha=0.5)+
  geom_line(aes(y = .epred, group = paste(Species, .draw), color=Species),lwd=1) + 
  facet_wrap(~ .draw) +theme_bw()
iris %>%
  tidybayes::add_predicted_draws(model2, ndraws = 100) %>%
  ggplot(aes(x = Petal.Length)) +
  geom_density(aes(x = .prediction, group = .draw), alpha=0.5) +
  facet_wrap(~Species)+theme_bw()

predict_flower <- posterior_predict(
  model2, 
  newdata = data.frame(Species = c("setosa", "versicolor","virginica"),
                       Petal.Length = c(1,1,1)))
bayesplot::mcmc_areas(predict_flower,prob = 0.8)+
  ggplot2::scale_y_discrete(labels = c("setosa", "versicolor","virginica"))

get_prior(Sepal.Length~Petal.Length+Sepal.Width + Petal.Width + Species, data=iris)
fit<-brm(Sepal.Length~Petal.Length+Sepal.Width + Petal.Width + Species, 
         data=iris, 
         prior=c(
           prior(normal(5, 1),  class = Intercept),
           prior(normal(0, 1),    class = b, coef=Petal.Length),
           prior(normal(0, 1),    class = b, coef=Petal.Width),
           prior(normal(2, 1),  class = b, coef=Speciesversicolor),
           prior(normal(3, 1),  class = b, coef=Speciesvirginica),
           prior(gamma(1, 1),     class = sigma)),
         chains=4, cores=6)
fit$prior
summary(fit)
plot(fit)
pp_check(fit, ndraws = 100)
stancode(fit)
pp_check(fit, 
         type = "error_hist", 
         ndraws = 50)
pp_check(fit, 
         type = "scatter_avg", 
         ndraws = 100)
pp_check(fit, 
         type = "stat_2d")
pp_check(fit, type = "rootogram")
pp_check(fit, type = "loo_pit")

ce<-conditional_effects(fit, surface = TRUE)
pl_ce<-plot(ce, rug=TRUE, theme=theme_bw())
pl_ce




posterior_interval(fit, prob = 0.80)
iris %>% 
  group_by(Species)%>%
  tidybayes::add_epred_draws(fit, ndraws = 20) %>%
  ggplot(aes(x = Sepal.Width, y = Sepal.Length)) +
  geom_point(aes(color=Species), alpha=0.5)+
  geom_line(aes(y = .epred, group = paste(Species, .draw), color=Species),lwd=1) + 
  facet_wrap(~ .draw) +
  theme_bw()
iris %>%
  tidybayes::add_predicted_draws(fit, ndraws = 100) %>%
  ggplot(aes(x = Sepal.Length)) +
  geom_density(aes(x = .prediction, group = .draw), alpha=0.5) +
  facet_wrap(~Species)+theme_bw()

iris %>%
  add_epred_draws(fit,ndraws=150)%>%
  ggplot(aes(x = Sepal.Width, 
             y = Sepal.Length, 
             color = as.factor(Species))) +
  geom_line(aes(y = .epred, group = paste(Species, .draw)), alpha = 0.25) +
  geom_point(data = iris, color="black")+
  facet_wrap(~Species)+
  labs(color="Species")+
  theme_bw()+
  theme(legend.position = "none")


iris %>%
  add_epred_draws(fit,ndraws=150) %>%
  ggplot(aes(x = Sepal.Length,
             y = .epred,
             color = as.factor(Species), 
             group=.draw))+
  geom_point()+
  geom_abline(aes(intercept=0, slope=1), color="black", lty=2)+
  facet_wrap(~Species)+
  labs(color="Species")+
  theme_bw()+
  theme(legend.position = "none")

iris %>%
  add_epred_draws(fit,ndraws=150) %>%
  ggplot(aes(x = Sepal.Length,
             y = .epred,
             color = as.factor(Species),
             fill = as.factor(Species),
             group=.draw))+
  geom_smooth(alpha=0.1, se=FALSE)+
  geom_abline(aes(intercept=0, slope=1), color="black", lty=2)+
  facet_wrap(~Species, scales="free")+
  labs(color="Species")+
  theme_bw()+
  theme(legend.position = "none")

iris %>%
  add_epred_draws(fit,ndraws=150) %>%
  ggplot(aes(x = Sepal.Width, 
             y = Sepal.Length, 
             color = as.factor(Species))) +
  geom_density_2d(aes(y = .epred, group = paste(Species, .draw)), alpha = 0.25) +
  geom_point(data = iris, color="black")+
  facet_wrap(~Species)+
  labs(color="Species")+
  theme_bw()+
  theme(legend.position = "none")


iris %>%
  add_epred_draws(fit, dpar = c("mu", "sigma")) %>%
  sample_draws(10) %>%
  ggplot(aes(y = Species)) +
  stat_dist_slab(aes(dist = "norm", arg1 = mu, arg2 = sigma), 
                 slab_color = "gray65", alpha = 1/10, fill = NA) +
  geom_point(aes(x = Sepal.Length), 
             data = iris, 
             shape = 21, 
             fill = "#9ECAE1", size = 2)+
  theme_bw()

iris %>%
  add_epred_draws(fit, dpar = c("mu", "sigma")) %>%
  ggplot(aes(x = Species)) +
  stat_dist_slab(aes(dist = "norm", arg1 = mu, arg2 = sigma), 
                 slab_color = "gray65", alpha = 1/10, fill = NA, 
                 data = . %>% sample_draws(10), scale = .5) +
  stat_halfeye(aes(y = .epred), side = "bottom", scale = .5) +
  geom_point(aes(y = Sepal.Length), data = iris, shape = 21, 
             fill = "#9ECAE1", size = 2, position = position_nudge(x = -.2))+
  theme_bw()

summary(iris$Sepal.Width)
summary(iris$Petal.Length)
summary(iris$Petal.Width)
nd <- tibble(Sepal.Width = seq(from = 2, to = 6, length.out = 20), 
             Petal.Length = seq(from = 1, to = 6, length.out = 20), 
             Petal.Width = seq(from = 0.1, to = 2, length.out = 20),
             Species = c("setosa", "versicolor", "virginica","setosa", "versicolor", "virginica",
                         "setosa", "versicolor", "virginica","setosa", "versicolor", "virginica",
                         "setosa", "versicolor", "virginica","setosa", "versicolor", "virginica",
                         "setosa", "versicolor"))
nd
p1<-predict(fit, newdata = nd) %>%
  data.frame() %>%
  bind_cols(nd) %>% 
  ggplot(aes(x = Petal.Width)) +
  geom_pointrange(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
                  color = "blue", shape = 20, lty=2) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "blue", alpha=0.1) +
  geom_line(aes(y = Estimate),
            color = "blue", lwd=1) +
  geom_point(data = iris, 
             aes(y = Sepal.Length), color="black") +
  labs(title="Lets see what the posterior looks like",
       subtitle = "Data with the percentile-based 95% intervals and the means of the posterior predictions",
       x = "Petal Width",
       y = "Sepal Length") +
  theme_bw()+
  facet_wrap(~Species, scales="free")
p2<-predict(fit, newdata = nd) %>%
  data.frame() %>%
  bind_cols(nd) %>% 
  ggplot(aes(x = Petal.Length)) +
  geom_pointrange(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
                  color = "blue", shape = 20, lty=2) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "blue", alpha=0.1) +
  geom_line(aes(y = Estimate),
            color = "blue", lwd=1) +
  geom_point(data = iris, 
             aes(y = Sepal.Length), color="black") +
  labs(x = "Petal Length",
       y = "Sepal Length") +
  theme_bw()+
  facet_wrap(~Species, scales="free")
p3<-predict(fit, newdata = nd) %>%
  data.frame() %>%
  bind_cols(nd) %>% 
  ggplot(aes(x = Sepal.Width)) +
  geom_pointrange(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
                  color = "blue", shape = 20, lty=2) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "blue", alpha=0.1) +
  geom_line(aes(y = Estimate),
            color = "blue", lwd=1) +
  geom_point(data = iris, 
             aes(y = Sepal.Length), color="black") +
  labs(x = "Sepal Width",
       y = "Sepal Length") +
  theme_bw()+
  facet_wrap(~Species, scales="free")
grid.arrange(p1,p2,p3)


fit%>%
  recover_types(iris) %>%
  emmeans::emmeans( ~ Species, data = iris) %>%
  emmeans::contrast(method = "pairwise") %>%
  gather_emmeans_draws() %>%
  ggplot(aes(x = .value, y = contrast)) +
  stat_halfeye()+
  theme_bw()