library(rstan) library(rstanarm) library(loo) library(tidyverse) library(bayesplot) library(shinystan) library(here) library(lme4) library(gridExtra) library(readr) rstan_options(auto_write = TRUE) #cache model options(mc.cores = parallel::detectCores()) #multiple cores for multiple chains ###set wd setwd("D:/Google Drive/R/RStan/sales_marketing_pipeline_revisited") ###import data pipelinedata <- read_csv("RStan/sales_marketing_pipeline_revisited/data/pipelinedata.csv") ###select key salespeople to allow convergence d <- pipelinedata %>% dplyr::filter(dsp_id %in% c(2, 4, 6, 11, 18, 23)) d$dsp_id <- as.numeric(recode(d$dsp_id, '2' = '1', '4' = '2', '6' = '3', '11' = '4', '18' = '5', '23' = '6')) ###lme4 m_lme <- lme4::glmer(formula = close ~ (1|year_id) + (1|dsp_id) - 1, data = d, family = binomial) # (1+year_id/dsp_id) coef(m_lme)$year_id #coef(m_lme)$`dsp_id:year_id`[[1]][1:5] ###RStanARM m_rstanarm <- stan_glmer(formula = close ~ (1|year_id) + (1|dsp_id) - 1, data = d, family = binomial, prior = normal(location = 0, scale = 1, autoscale = FALSE), prior_intercept = normal(location = 0, scale = 1, autoscale = FALSE), prior_aux = exponential(rate = 1, autoscale = FALSE), cores = 4, chains = 4, iter = 250, warmup = 50, seed = 1, refresh = 1) coef(m_rstanarm)$year_id #coef(m_rstanarm)$`dsp_id:year_id`[[1]][1:5] ###Stan ####prep data N_obs = nrow(d) #level 1 n close = d$close year_id = d$year_id N_year = length(unique(year_id)) dsp_id = d$dsp_id N_dsp = length(unique(dsp_id)) stan_rdump(c("N_obs", "close", "year_id", "dsp_id", "N_year", "N_dsp"), file="data/stan_data_dump.R") data <- read_rdump("data/stan_data_dump.R") #import data ####prep model m1 <- stanc(file = "model/stan_model.stan") # Check Stan file m1 <- stan_model(stanc_ret = m1) # Compile Stan code ####fit model fit1 <- sampling(m1, #use sampling function, rather than stan() warmup=50, # iter=250, #total iterations including warmup seed=1, #ensures reproducibility data=data, chains = 1, #each chain has different initial condition which helps to find pathological neighborhood of posterior cores = getOption("mc.cores", 1L), verbose = FALSE, refresh = 100) ####coefficients list_of_draws <- rstan::extract(fit1) year_df <- data.frame(list_of_draws$year) year_df <- year_df - median(list_of_draws$year_hyp) apply(year_df, 2, median)