bayesian_penguins.R

library(palmerpenguins)
library(rstanarm)
library(tidybayes)
library(tidyverse)
library(brms)

minsmaxs = penguins %>%
           drop_na() %>% 
           group_by(species) %>% 
           summarise(low = min(flipper_length_mm),
                     high = max(flipper_length_mm), 
                     m = mean(flipper_length_mm),
                     s = sd(flipper_length_mm))


model_data = penguins %>% 
              drop_na() %>% 
              left_join(minsmaxs) %>% 
              mutate(x = (flipper_length_mm - m)/s)

model_1 = brm( bill_length_mm~x*species, data = model_data)
model_2 = brm( bill_length_mm~x+species, data = model_data)


bayes_lines = penguins %>%
  drop_na() %>% 
  group_by(species) %>% 
  modelr::data_grid(flipper_length_mm) %>% 
  left_join(minsmaxs) %>% 
  mutate(x = (flipper_length_mm - m)/s) %>% 
  add_fitted_draws(model, n = 100) %>% 
  ggplot()+
  geom_point(data = model_data, aes(flipper_length_mm, bill_length_mm, color = species, shape = species ))+
  geom_line(aes(flipper_length_mm, .value, group = interaction(.draw, species), color = species),alpha = 0.1)

  
    

bayes_lines +
  scale_color_brewer(palette = 'Set1')+
  labs(x = 'Flipper Length (mm)',
       y = 'Bill Length (mm)',
       title = 'Penguins, But Make It Bayesian!',
       subtitle = 'Posterior estimates of bill length as a function of flipper length')+
  theme_minimal()+
  theme(aspect.ratio = 0.62,
        legend.position = c((220-170)/(230-170), (35-30)/(60-30)))