In stemangiola/RNAseq-noise-model:

This notebook tests that the implemented models can be run on simulated data and roughly recover parameters, the tests are meant to be quick and not really rigorous

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

model_dirichlet_multinomial = stan_model(here("stan","dirichlet_multinomial_base.stan"))
model_generalized_dirichlet_multinomial = stan_model(here("stan","generalized_dirichlet_multinomial.stan"))
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"))
model_nb_deseq2 = stan_model(here("stan","negBinomial_deseq2.stan"))

Multinomial

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
data = generate_data_multinomial(G, sums)
fit = sampling(model_multinomial, data = data$observed)
evaluation_summary(fit, data$true)

Neg. binomial multinomial

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
data = generate_data_nb_multinomial(G, sums, lambda_raw_sigma = 2)
fit = sampling(model_nb_multinomial_approx, data = data$observed)
evaluation_summary(fit, data$true)

Dirichlet multinomial - Martin's variant

G = 5
sums = c(100,200,300)
lambda_raw_prior = 2 
data = generate_data_dirichlet_multinomial_base(G, sums, lambda_raw_prior)
fit = sampling(model_dirichlet_multinomial, data = data$observed)
evaluation_summary(fit, data$true)

LogNormal Multinomial

G = 4
N = 6
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
data = generate_data_lognormal_multinomial(G, sums)
fit = sampling(model_lognormal_multinomial, data = data$observed, control = list(max_treedepth = 12, adapt_delta = 0.9), iter = 500)
evaluation_summary(fit, data$true)

G = 5
N = 10
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
data = generate_data_lognormal_multinomial(G, sums, is_prior_assymetric = TRUE)
fit = sampling(model_lognormal_multinomial, data = data$observed, control = list(max_treedepth = 12, adapt_delta = 0.9), iter = 500)
evaluation_summary(fit, data$true)

Gamma multinomial

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

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

data = generate_data_gamma_multinomial(G, sums)
fit = sampling(model_gamma_multinomial, data = data$observed)
evaluation_summary(fit, data$true)

Gamma multinomial 2

model_gamma_multinomial_2 = stan_model(here("stan","gamma_multinomial_2.stan"))

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

data = generate_data_gamma_multinomial_2(G, sums, asymptDisp = 0.5, extraPois = 1.3)
data$observed$N_samples_log_lik <- 1
fit = sampling(model_gamma_multinomial_2, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit, data$true)

Gamma multinomial - NB approx

model_gamma_multinomial_2_nb = stan_model(here("stan","gamma_multinomial_2_nb.stan"))

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

data = generate_data_gamma_multinomial_2(G, sums, asymptDisp = 0.5, extraPois = 1.3)
fit = sampling(model_gamma_multinomial_2_nb, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit, data$true)

data$observed$N_samples_log_lik <- 1
data$observed$generate_log_lik = 0
fit_g = sampling(model_gamma_multinomial, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit_g, data$true)

model_dirichlet_multinomial = stan_model(here("stan","dirichlet_multinomial_base.stan"))

G = 10
N = 6
sums = G * 10 + rnbinom(N, mu = G * 50, size = 1) 
data <- generate_data_dirichlet_multinomial_base(G, sums, lambda_raw_prior = 2)

data$observed$variant = 0
fit_d <- sampling(model_dirichlet_multinomial, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit_d, data$true)

data$observed$variant = 1
fit_n <- sampling(model_dirichlet_multinomial, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit_n, data$true)

data$observed$variant = 2
fit_g <- sampling(model_dirichlet_multinomial, data = data$observed, control = list(adapt_delta = 0.99))
evaluation_summary(fit_g, data$true)

Negative binomial, DeSEQ2 style

#TODO This crashes (caused by G and N values)
G = 10
N = 15

data = generate_data_negbinomial_deseq2(G, N, asymptDisp = 0.5, extraPois = 1.3)
fit = sampling(model_nb_deseq2, data = data$observed)
evaluation_summary(fit, data$true)

Generalized dirichlet multinom

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

#Right now, I don't have direct simulator
data = generate_data_gamma_multinomial(G, sums)
fit = sampling(model_generalized_dirichlet_multinomial, data = data$observed)
fit