In OlivierGranacher/modelStan: Helper Functions for Stan Models

knitr::opts_chunk$set(echo = T, include = T, eval = T, message = F, warning = F, fig.asp = 0.4, fig.path = 'Figs/')
  library(rethinking)
  library(rstan)
  library(data.table)
  library(modelStan)
  library(ggplot2)
  library(magrittr)

  #devtools::install_github("rmcelreath/rethinking", ref = "Experimental")
  options(mc.cores = parallel::detectCores())
  rstan::rstan_options(auto_write = TRUE)
  pal_disc_qual <- 'harmonic' # discrete palette from colorspace qualitative
  scaleCustom <- colorspace::scale_fill_discrete_qualitative(palette = pal_disc_qual,
                                                             aesthetics = c('color', 'fill'),
                                                             name = "")
  theme_set(ggthemes::theme_tufte())
  theme_set(theme_minimal())

  set.seed(4)
  simDT <- modelStan::simPotData(N = 10,
                                 np = c(3, 5),
                                 Y = c(1, 0),
                                 sigma = c(2, 2),
                                 sigma_g = c(.2, .2),
                                 time_effect = T,
                                 p = 5,
                                 q = 5,
                                 e = 5
                                 )
  # plot of Y mean evolution with time
  simDT[, .(Y_mean = mean(Y)), by = id_dat] %>% 
  {ggplot(., (aes(x = id_dat))) +
    geom_line(aes(y = Y_mean), color = "grey50", size = .5, alpha = .5) +
    geom_point(aes(y = Y_mean), size = 1) +
    ggtitle("Valeur moyenne de Y dans le temps") +
    scale_x_continuous(minor_breaks = NULL, breaks = NULL, name = NULL) +
    scale_y_continuous(breaks = c(min(.$Y_mean), mean(.$Y_mean), max(.$Y_mean)),
                       name = NULL,
                       labels = scales::comma_format(.01)) +
        ggthemes::geom_rangeframe(inherit.aes = F, aes(y = Y_mean))
      }

  d <- tidybayes::compose_data(simDT[, j =  !c("dat")],
                          GC = tidybayes::x_at_y(G, C),
                          BC = tidybayes::x_at_y(B, C),
                          .n_name = tidybayes::n_prefix("N")
  )
  d$N_dat <- length(unique(simDT[, id_dat]))
  d$S <- 8 #  sd for mean effect 
  d$E <- 0  # mean expected effect

lm basic model

  mod1 <- lm(formula = Y ~ G, data = simDT)
  summary(mod1)

Modele Stan base non-centered without time effect

  mod.0nc <- rstan::stan(file = "simPotData_0nc.stan",
                       data = d,
                       chains = 4,
                       iter = 2000)
 precis(mod.0nc, 2)
 modelStan::plotStanParam(mod.0nc, "alpha_diff") +
   ggtitle("non-centered WITHOUT time effect")

  modelStan::addStanPred(mod.0nc, data = simDT, "Y_sim") %>% 
   {ggplot(., aes(x = C)) + 
   geom_point(aes(y = Y), color = "grey50") +
   geom_linerange(aes(y = pred, ymin = q10, ymax = q90),
                   color = "indianred",
                   shape = 21,
                   fill = NA,
                   fatten = .4) +
   stat_summary(aes(y = Y),fun.y = mean, geom = "point", color = "black", size = 3) +
    ggtitle("non-centered WITHOUT time effect") + ggthemes::geom_rangeframe(aes(y = Y))}

Modele Stan base non-centered with time effect

  mod.0nc_time <- rstan::stan(file = "simPotData_0nc_time.stan",
                       data = d,
                       chains = 4,
                       iter = 2000)
 precis(mod.0nc_time, 2)
 modelStan::plotStanParam(mod.0nc_time, c("alpha_diff")) +
    ggtitle("non-centered with time effect")

  modelStan::addStanPred(mod.0nc_time, data = simDT, "Y_sim") %>% 
   ggplot(aes(x = C)) + 
   geom_point(aes(y = Y)) +
   geom_pointrange(aes(y = pred, ymin = q10, ymax = q90), color ="indianred", shape = 21, fill = NA) +
   stat_summary(aes(y = Y),fun.y = mean, geom = "point", color = "black", size = 3) +
    ggtitle("non-centered with time effect")

  # Graphique de l'evolution temporelle de Y
 simDT %>% 
    dplyr::group_by(id_dat) %>% 
    dplyr::summarise(Y_mean = mean(Y),
                     Y_q10  = quantile(Y, .1),
                     Y_q90  = quantile(Y, .9)
                     ) %>% 
    modelStan::addStanPred(mod.0nc_time, data = ., "Y_dat") %>% 
    {ggplot(., (aes(x = id_dat))) +
    geom_pointrange(aes(y = Y_mean, ymin = Y_q10, ymax = Y_q90), color = "grey80", fatten  = .3) +
    geom_point(aes(y = pred), color = "indianred4", size = 2, shape = 21) +
    ggtitle("Valeur Moyenne Y dans le temps") +
    scale_x_continuous(minor_breaks = NULL, breaks = NULL, name = NULL) +
    scale_y_continuous(breaks = c(min(.$Y_mean), mean(.$Y_mean), max(.$Y_mean)),
                       name = NULL,
                       labels = scales::comma_format(.01)) +
        ggthemes::geom_rangeframe(inherit.aes = F, aes(y = Y_mean))
      }

Modele Stan base non-centered with hierarchical time effect

  mod.0nc_Htime <- rstan::stan(file = "simPotData_0nc_Htime.stan",
                       data = d,
                       chains = 4,
                       iter = 2000)
 precis(mod.0nc_Htime, 2)
 modelStan::plotStanParam(mod.0nc_time, c("alpha_diff")) +
    ggtitle("non-centered with H time effect")

  modelStan::addStanPred(mod.0nc_Htime, data = simDT, "Y_sim") %>% 
   ggplot(aes(x = C)) + 
   geom_point(aes(y = Y)) +
   geom_pointrange(aes(y = pred, ymin = q10, ymax = q90), color ="indianred", shape = 21, fill = NA) +
   stat_summary(aes(y = Y),fun.y = mean, geom = "point", color = "black", size = 3) +
    ggtitle("non-centered with H time effect")

  # Graphique de l'evolution temporelle de Y
 simDT %>% 
    dplyr::group_by(id_dat) %>% 
    dplyr::summarise(Y_mean = mean(Y),
                     Y_q10  = quantile(Y, .1),
                     Y_q90  = quantile(Y, .9)
                     ) %>% 
    modelStan::addStanPred(mod.0nc_Htime, data = ., "Y_dat") %>% 
    {ggplot(., (aes(x = id_dat))) +
    geom_pointrange(aes(y = Y_mean, ymin = Y_q10, ymax = Y_q90), color = "grey80", fatten  = .3) +
    geom_point(aes(y = pred), color = "indianred4", size = 2, shape = 21) +
    ggtitle("Valeur Moyenne Y dans le temps") +
    scale_x_continuous(minor_breaks = NULL, breaks = NULL, name = NULL) +
    scale_y_continuous(breaks = c(min(.$Y_mean), mean(.$Y_mean), max(.$Y_mean)),
                       name = NULL,
                       labels = scales::comma_format(.01)) +
        ggthemes::geom_rangeframe(inherit.aes = F, aes(y = Y_mean))
      }

Getting parameter values with tidybayes package

# table of results for 2 modeles 
dplyr::bind_rows(mod.0nc %>% 
  tidybayes::recover_types(simDT) %>% 
  tidybayes::spread_draws(alpha_diff[G]) %>% 
  tidybayes::mean_qi(), 
  mod.0nc_time %>% 
  tidybayes::recover_types(simDT) %>% 
  tidybayes::spread_draws(alpha_diff[G]) %>% 
  tidybayes::mean_qi(),
  mod.0nc_Htime %>% 
  tidybayes::recover_types(simDT) %>% 
  tidybayes::spread_draws(alpha_diff[G]) %>% 
  tidybayes::mean_qi()
) %>% dplyr::bind_cols(model = c("mod.0nc", "mod.0nc_time", "mod.0nc_Htime")) %>% 
  ggplot(aes( x = model)) + geom_pointrange(aes(y = alpha_diff, ymin = .lower, ymax= .upper)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = 2)
# Putting names back in data  
# mod.0nc %<>% 
#   tidybayes::recover_types(simDT)

# plot of results for alpha top
 # mod.0nc %>% 
 #   tidybayes::spread_draws(alpha_top[G]) %>% 
 #   ggplot(aes(x = alpha_top, y = G)) +
 #   tidybayes::stat_pointintervalh()