set.seed(13254) # Import libraries library(tidyverse) library(magrittr) library(rstan) library(tidybayes) library(here) library(foreach) library(doParallel) library(loo) registerDoParallel() # library(future) # plan(multiprocess) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())
N = 50 G = 10 counts_source = "Acute_Myeloid_Leukemia_Primary_Blood_Derived_Cancer_-_Peripheral_Blood" counts_tidy <- load_test_data(N, G, counts_source) data_for_stan = data_for_stan_from_tidy(N, G, counts_tidy)
model_names = c( #"multinomial", "dirichlet_multinomial", #"negBinomial" "gamma_multinomial" )
models = foreach(model_name = model_names) %do% { library(rstan) #On Windows, %dopar%-ed code does not share the main session library(here) #On Windows, %dopar%-ed code does not share the main session stan_model(here("stan",sprintf("%s.stan", model_name))) } %>% setNames(model_names)
data_for_stan <- generate_data_gamma_multinomial(10, rep(1000, 30))$observed data_for_stan$generate_quantities <- 1 data_for_stan$my_prior <- c(0,5) data_for_stan$is_prior_asymetric <- 0 data_for_stan$exposure <- rowSums(data_for_stan$counts)
model_defs = tibble(model_name = model_names) models_list = models[model_defs$model_name] #Copy the data for all runs data_list = list() for(i in 1:length(models_list)) { data_list[[i]] = data_for_stan } fits = sampling_multi(models, data_list) %>% setNames(model_names)
loo_res <- fits %>% map(function(x) { log_lik <- extract_log_lik(x, merge_chains = FALSE) r_eff <- relative_eff(exp(log_lik)) loo(log_lik, r_eff = r_eff) }) loo_res %>% map(plot) compare(x = loo_res)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.