beta_binomial_example

Libraries

library(tidyverse)
library(tidybayes)
library(rstan)
library(posterior)
set.seed(236)
options(mc.cores = parallelly::availableCores())
rstan_options(auto_write = TRUE)

Data

rat_data = read.table("https://raw.githubusercontent.com/vnminin/sismid_mcmc_one/master/2023/code/rat_tumor.txt", header=TRUE)

num_obs = dim(rat_data)[1]

Fitting the model using rstan

The following code fits the beta_binom.stan model to the rat tumor data. The posterior_object contains both the data, as well as parameters specifying the model priors. These are all treated as “data” by stan.

posterior_object <- list(priors_only = FALSE,
                         num_obs = num_obs,
                         x = rat_data$x,
                         n = rat_data$n,
                         lambda_alpha = 0.1,
                         lambda_beta = 0.1)

chain = 4
iterations = 1000
seed = 236
thin = 1

posterior_fit <- stan(file = "https://raw.githubusercontent.com/vnminin/sismid_mcmc_one/main/2023/code/beta_binom.stan",
                  data = posterior_object,
                  seed = seed,
                  iter = iterations,
                  thin = thin,
                  chain = chain)

Assessing MCMC Diagnostics

We can assess the traceplot for the scaled likelihood.

traceplot(posterior_fit, pars = "lp__")

We can also summarise Rhat and effective sample size.

summary <- posterior_fit %>% 
           as_draws() %>% 
           summarise_draws() %>% 
           summarise(
             min_rhat = min(rhat),
             max_rhat = max(rhat),
             min_ess_bulk = min(ess_bulk),
             max_ess_bulk = max(ess_bulk),
             min_ess_tail = min(ess_tail),
             max_ess_tail = max(ess_tail)
           )

summary

## # A tibble: 1 × 6
##   min_rhat max_rhat min_ess_bulk max_ess_bulk min_ess_tail max_ess_tail
##      <dbl>    <dbl>        <dbl>        <dbl>        <dbl>        <dbl>
## 1    0.999     1.01         343.        3693.         572.        2098.

Visualizing the results

theta_draws <- posterior_fit %>% spread_draws(theta[i])

vis_theta <- theta_draws %>% filter(i %% 2 == 1) %>%
             arrange(i) %>%
             mutate(theta_lab = paste("Theta ", i)) 


theta_plot <- vis_theta %>% 
              ggplot(aes(x = fct_inorder(theta_lab), y = theta)) +
              geom_boxplot() +
              theme_bw() +
              xlab("Theta") +
              ylab("Value") + 
              theme(axis.text.x = element_text(angle = 90),
                    text = element_text(size = 20))

theta_plot

fixed_vars <- posterior_fit %>% spread_draws(alpha, beta) 

plot_vars <- fixed_vars %>% pivot_longer(-c(.chain, .iteration, .draw), names_to = "name")

fixed_plot <- plot_vars %>% 
              ggplot(aes(x=value)) +
              geom_histogram() +
              facet_wrap(vars(name)) +
              theme_bw() +
              xlab("Value") + 
              theme(text = element_text(size = 20))

fixed_plot