FishersExactStanTest

library(tidyverse)
library(rstan)
library(bayesplot)

#code
model_code <- "
data {
  int<lower = 2> K; // # of categories
  vector[K] priorprobs; //prior
  int<lower=0> n_total_1; //total in group 1
  int<lower=0> n_total_2; //total in group 2
  int<lower=0> n_counts_1[K]; // counts 1
  int<lower=0> n_counts_2[K]; // counts 2
}
parameters { 
  simplex[K] theta_1;
  simplex[K] theta_2;
}
model {
  //priors
  theta_1 ~ dirichlet(priorprobs); // group 1 thetas
  theta_2 ~ dirichlet(priorprobs); // group 2 thetas
  
  
  //likelihood
  n_counts_1 ~ multinomial(theta_1);
  n_counts_2 ~ multinomial(theta_2);
}
generated quantities {
  vector[K] diff;
  
  diff = theta_1 - theta_2;
}


"
n <- 500
group1counts<- as.vector(rmultinom(1, n, c(40,5,6,7,8,9,8,7,6,5,80)))
group2counts<- as.vector(rmultinom(1, n, c(1,2,3,4,5,6,5,4,3,2,1)))

standata <- list(K = 11,
                 priorprobs = rep(1,11),
                 n_total_1 = 200,
                 n_total_2 = 200,
                 n_counts_1 = group1counts,
                 n_counts_2 = group2counts)

fishersExact <- stan(model_code = model_code,
     data = standata)

mcmc_areas(fishersExact, pars = paste0("theta_1[", 11:1, "]"), prob = 0.89) + scale_y_discrete(expand = c(0,0.1)) +
  xlim(c(-1,1))

mcmc_areas(fishersExact, pars = paste0("theta_2[", 11:1, "]"), prob = 0.89) + scale_y_discrete(expand = c(0,0.1)) +
  xlim(c(-1,1))

mcmc_areas(fishersExact, pars = paste0("diff[", 11:1, "]"), prob = 0.89) + scale_y_discrete(expand = c(0,0.1)) +
  xlim(c(-1,1))