In stemangiola/RNAseq-noise-model:

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)

stemangiola/RNAseq-noise-model documentation built on Oct. 18, 2019, 1:47 p.m.