In stemangiola/RNAseq-noise-model:

devtools::load_all()
library(here)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(tidyverse)
library(cowplot)

model_multinomial = stan_model(here("stan","multinomial.stan"))
model_lognormal_multinomial = stan_model(here("stan","logNormal_multinomial.stan"))
model_nb_multinomial_approx = stan_model(here("stan","nb_multinomial_approx.stan"))

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
N_steps = 1000

generator <- function() { generate_data_multinomial(G, sums) }
sbc_res_multinomial <- sbc(model_multinomial, generator, N_steps = N_steps, control = list(adapt_delta = 0.9))
saveRDS(sbc_res_multinomial, here("local_temp_data",paste0("sbc_multinomial_",N_steps,".rds")))
plot_sbc_params(sbc_res_multinomial$params)
summarise_sbc_diagnostics(sbc_res_multinomial)

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
N_steps = 50

generator <- function() { generate_data_lognormal_multinomial(G, sums) }
sbc_res_lognormal_multinomial <- sbc(model_lognormal_multinomial, generator, N_steps = N_steps, control = list(adapt_delta = 0.9, max_treedepth = 12), iter = 500)
saveRDS(sbc_res_lognormal_multinomial, here("local_temp_data",paste0("sbc_lognormal_multinomial_",N_steps,".rds")))
plot_sbc_params(sbc_res_lognormal_multinomial$params %>% filter(!grepl("theta_z",param_name)))
plot_sbc_params(sbc_res_lognormal_multinomial$params %>% filter(grepl("theta_z\\[2,",param_name)))
summarise_sbc_diagnostics(sbc_res_lognormal_multinomial)

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 

generator <- function() { generate_data_nb_multinomial(G, sums, lambda_raw_sigma = 2)}
sbc_res <- sbc(model_nb_multinomial_approx, generator, 1000, control = list(adapt_delta = 0.9))
plot_sbc_params(sbc_res$params)
summarise_sbc_diagnostics(sbc_res)

model_gamma_multinomial = stan_model(here("stan","gamma_multinomial.stan"))

G = 4
N = 3
sums = G * 10 + rnbinom(N, mu = G * 500, size = 1) 
N_steps = 100

generator <- function() { generate_data_gamma_multinomial(G, sums)}
sbc_res <- sbc(model_gamma_multinomial, generator, N_steps, control = list(adapt_delta = 0.95))
plot_sbc_params(sbc_res$params, n_bins = 10)
summarise_sbc_diagnostics(sbc_res)

true_mu <- sbc_res$true_values %>% map_dbl(function(x) { x$mu}) 
true_phi <- sbc_res$true_values %>% map_dbl(function(x) { x$phi}) 
n_zeroes <- sbc_res$data %>% map_dbl(function(x) {sum(x$counts == 0)})
plot(true_mu, sbc_res$diagnostics$n_divergent)

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