Cherry blossom

rm(list = ls())
library(LaplacesDemon)
library(tidyverse)
library(brms)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(coda)
library(mvtnorm)
library(loo)
library(dagitty)
library(MASS)
library(splines)
library(rethinking)
library(BayesFactor)
library(gRain)
library(bnlearn)
library(readxl)
library(mltools)
library(data.table)
library(lattice)
library(skimr)
library(lme4)
library(bayesrules)
library(tidyverse)
library(bayesplot)
library(rstanarm)
library(tidybayes)
library(broom.mixed)
library(janitor)
library(grid)
library(gridExtra)
library(ggstatsplot)
library(GGally)

### Example 15
data(cherry_blossom_sample)
skimr::skim(cherry_blossom_sample)
running <- cherry_blossom_sample
nrow(running)
ggplot(running, aes(x = age, y = net)) + 
  geom_point() + 
  facet_wrap(~ runner)+
  theme_bw()
ggplot(running, aes(x = runner, y = net)) + 
  geom_boxplot() + theme_bw()
ggplot(running, aes(x = age, y = net, group = runner)) + 
  geom_smooth(method = "lm", se = FALSE, color = "gray", size = 0.5) + 
  geom_abline(aes(intercept = 75.2, slope = 0.268), color = "blue") + theme_bw()
head(running, 2)
ggplot(running, aes(y = net, x = age)) + 
  geom_point() + theme_bw()

ggplot(running, aes(y = net, x = age, size=gun, color=as.factor(year))) + 
  geom_point() + theme_bw()

ggcorrmat(
  data     = running,
  title    = "Correlalogram for Cherry blossom dataset")

running$year<-as.factor(running$year)
running%>%
  dplyr::select(age, net, gun, year)%>%
  ggpairs(.,
          aes(color = year),
          upper = list(continuous = wrap('cor', size = 3)),
          lower = list(combo = wrap("facethist", bins = 30)),
          diag = list(continuous = wrap("densityDiag", alpha = 0.5)),
          title = "Scatterplot matrix of `Cherry Blossom` Grouped by Year")+theme_bw()

get_prior(net ~ age, 
          data = running, family = gaussian)
complete_pooled_model <- stan_glm(
  net ~ age, 
  data = running, family = gaussian, 
  prior_intercept = normal(0, 2.5, autoscale = TRUE),
  prior = normal(0, 2.5, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4, iter = 5000*2, seed = 84735)
tidy(complete_pooled_model, conf.int = TRUE, conf.level = 0.80)

examples <- running %>% 
  filter(runner %in% c("1", "20", "22"))
g1<-ggplot(examples, aes(x = age, y = net)) + 
  geom_point() + 
  facet_wrap(~ runner) + 
  geom_abline(aes(intercept = 75.2242, slope = 0.2678), 
              color = "blue") + theme_bw()
g2<-ggplot(examples, aes(x = age, y = net)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE, fullrange = TRUE) + 
  facet_wrap(~ runner) + 
  xlim(52, 62) + theme_bw()
grid.arrange(g1,g2)


model_summary<-tidy(complete_pooled_model, 
                    conf.int = TRUE, 
                    conf.level = 0.80)
model_summary
B0 <- model_summary$estimate[1];B0
B1 <- model_summary$estimate[2];B1
ggplot(running, aes(x = age, y = net)) + 
  geom_point() + 
  geom_abline(aes(intercept = B0, slope = B1))+
  theme_bw()
running_model_1_prior <- stan_glmer(
  net ~ age + (1 | runner), 
  data = running, family = gaussian,
  prior_intercept = normal(100, 10),
  prior = normal(2.5, 1), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735, 
  prior_PD = TRUE)
summary(running_model_1_prior)
set.seed(84735)
running %>% na.omit() %>%
  add_fitted_draws(running_model_1_prior, n = 9) %>%
  ggplot(aes(x = age, y = net)) +
  geom_line(aes(y = .value, group = paste(runner, .draw))) + 
  facet_wrap(~ .draw)+theme_bw()
running %>% na.omit() %>%
  add_predicted_draws(running_model_1_prior, n = 100) %>%
  ggplot(aes(x = net)) +
  geom_density(aes(x = .prediction, group = .draw)) +
  xlim(-100,300)+theme_bw()

running_model_1 <- update(running_model_1_prior, prior_PD = FALSE)
prior_summary(running_model_1)
summary(running_model_1)
mcmc_trace(running_model_1)
mcmc_dens_overlay(running_model_1)
mcmc_acf(running_model_1)
plot(neff_ratio(running_model_1))
plot(rhat(running_model_1))

tidy_summary_1 <- tidy(running_model_1, effects = "fixed",
                       conf.int = TRUE, conf.level = 0.80)
B0 <- tidy_summary_1$estimate[1]
B1 <- tidy_summary_1$estimate[2]
tidy_summary_1
running %>% na.omit() %>%
  add_fitted_draws(running_model_1, n = 200, re_formula = NA) %>%
  ggplot(aes(x = age, y = net)) +
  geom_line(aes(y = .value, group = .draw), alpha = 0.1) +
  geom_abline(intercept = B0, slope = B1, color = "blue") +
  lims(y = c(75, 110))+theme_bw()

runner_summaries_1 <- running_model_1 %>%
  spread_draws(`(Intercept)`, b[,runner]) %>% 
  mutate(runner_intercept = `(Intercept)` + b) %>% 
  dplyr::select(-`(Intercept)`, -b) %>% 
  median_qi(.width = 0.80) %>% 
  dplyr::select(runner, runner_intercept, .lower, .upper)
runner_summaries_1 %>% 
  filter(runner %in% c("runner:4", "runner:5"))

running %>% na.omit() %>%
  filter(runner %in% c("4", "5")) %>% 
  add_fitted_draws(running_model_1, n = 100) %>%
  ggplot(aes(x = age, y = net)) +
  geom_line(
    aes(y = .value, group = paste(runner, .draw), color = runner),
    alpha = 0.1) +
  geom_point(aes(color = runner))+
  theme_bw()

running %>% na.omit() %>%
  add_fitted_draws(running_model_1, n = 100) %>%
  ggplot(aes(x = age, y = net)) +
  geom_line(
    aes(y = .value, group = paste(runner, .draw), color = runner),
    alpha = 0.1) +
  geom_point(aes(color = runner))+
  facet_wrap(~runner, ncol=5, scales="free")+
  theme_bw()


# Plot runner-specific models with the global model
ggplot(running, aes(y = net, x = age, group = runner)) + 
  geom_abline(data = runner_summaries_1, color = "gray",
              aes(intercept = runner_intercept, slope = B1)) + 
  geom_abline(intercept = B0, slope = B1, color = "blue") + 
  lims(x = c(50, 61), y = c(50, 135))+
  theme_bw()

running %>% 
  filter(runner %in% c("4", "5", "20", "29")) %>% 
  ggplot(., aes(x = age, y = net)) + 
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE) + 
  facet_grid(~ runner)+
  theme_bw()
ggplot(running, aes(x = age, y = net, group = runner)) + 
  geom_smooth(method = "lm", se = FALSE, size = 0.5)+
  theme_bw()

options(mc.cores=parallel::detectCores())
running_model_2 <- stan_glmer(
  net ~ age + (age | runner),
  data = running, family = gaussian,
  prior_intercept = normal(100, 10),
  prior = normal(2.5, 1), 
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735, adapt_delta = 0.99999)

tidy(running_model_2, 
     effects = "fixed", 
     conf.int = TRUE, 
     conf.level = 0.80)
tidy(running_model_2, 
     effects = "ran_vals", 
     conf.int = TRUE, 
     conf.level = 0.80)
tidy(running_model_2, 
     effects = "ran_pars", 
     conf.int = TRUE, 
     conf.level = 0.80)
tidy(running_model_2, 
     effects = "auxiliary", 
     conf.int = TRUE, 
     conf.level = 0.80)

runner_chains_2 <- running_model_2 %>%
  spread_draws(`(Intercept)`, b[term, runner], `age`) %>% 
  pivot_wider(names_from = term, names_glue = "b_{term}",
              values_from = b) %>% 
  mutate(runner_intercept = `(Intercept)` + `b_(Intercept)`,
         runner_age = age + b_age)
runner_summaries_2 <- runner_chains_2 %>% 
  group_by(runner) %>% 
  summarize(runner_intercept = median(runner_intercept),
            runner_age = median(runner_age))
head(runner_summaries_2, 3)
ggplot(running, aes(y = net, x = age, group = runner)) + 
  geom_point(alpha=0.2)+
  geom_abline(data = runner_summaries_2, color = "gray",
              aes(intercept = runner_intercept, slope = runner_age)) + 
  lims(x = c(50, 61), y = c(50, 135))+theme_bw()

runner_summaries_2 %>% 
  filter(runner %in% c("runner:1", "runner:10"))
tidy(running_model_2, effects = "ran_pars")
p1<-pp_check(complete_pooled_model) + 
  labs(x = "net", title = "complete pooled model")
p2<-pp_check(running_model_1) + 
  labs(x = "net", title = "running model 1")
p3<-pp_check(running_model_2) + 
  labs(x = "net", title = "running model 2")
grid.arrange(p1,p2,p3,ncol=3)

running %>% 
  na.omit() %>% 
  prediction_summary(model = complete_pooled_model, data = .)
running %>% 
  na.omit() %>% 
  prediction_summary(model = running_model_1, data = .)
running %>% 
  na.omit() %>% 
  prediction_summary(model = running_model_2, data = .)
running %>% 
  na.omit() %>% 
prediction_summary_cv(model = running_model_2, data =.,
                      k = 10, group = "runner")
elpd_hierarchical_1 <- loo(running_model_1)
elpd_hierarchical_2 <- loo(running_model_2)
loo_compare(elpd_hierarchical_1, elpd_hierarchical_2)

running %>% 
  filter(runner %in% c("1", "10")) %>% 
  ggplot(., aes(x = age, y = net)) + 
  geom_point() + 
  facet_grid(~ runner) + 
  lims(x = c(54, 61))

predict_next_race_1 <- posterior_predict(
  running_model_1, 
  newdata = data.frame(runner = c("1", "Miles", "10"),
                       age = c(61, 61, 61)))
class(predict_next_race_1)
predict_next_race_2 <- posterior_predict(
  running_model_2, 
  newdata = data.frame(runner = c("1", "Miles", "10"),
                       age = c(61, 61, 61)))


p1<-mcmc_areas(predict_next_race_1, prob = 0.8) +
  ggplot2::scale_y_discrete(labels = c("runner 1", "Miles", "runner 10"))
p2<-mcmc_areas(predict_next_race_2, prob = 0.8) +
  ggplot2::scale_y_discrete(labels = c("runner 1", "Miles", "runner 10"))
grid.arrange(p1,p2,ncol=2)

predict_next_race_1<-as.data.frame(predict_next_race_1);predict_next_race_1
predict_next_race_1$model="Model 1"
predict_next_race_2<-as.data.frame(predict_next_race_2);predict_next_race_2
predict_next_race_2$model="Model 2"

combined<-rbind(predict_next_race_1, predict_next_race_2)
combined_melt<-reshape2::melt(combined)
variable.labs<-c("runner 1", "Miles", "runner 10")
names(variable.labs) <- c("1", "2", "3")
ggplot(combined_melt,
       aes(x=value, fill=as.factor(model)))+
  geom_density(alpha=0.5)+
  theme_bw()+
  facet_wrap(vars(variable),
             labeller = labeller(variable = variable.labs),
             scales="free")+
  labs(x="Estimate P | age = 61", 
       y="Density", 
       fill="Model", 
       title="Comparing Bayesian Model estimates")