Exercises

First load the environment and install the required packages:

renv::restore(prompt = FALSE)

library(posterior)
library(priorsense)
library(bayesplot)

Part 1: Posterior draws objects

For this set of exercises we will use draws from the classic eight schools model. These draws are included in both posterior and ArviZ.

eight_schools_draws <- example_draws("eight_schools")

Understanding draws objects

Inspect the posterior object. How many chains, iterations and variables does it contain?

You can print the draws object, or you can use the nchains(), niterations() and nvariables() functions.

# Your code here

Subsetting draws

Extract only the first chain.

You can use the subset_draws() function, specifying the chain argument.

# Your code here

Extract only the first 10 iterations, from all the chains.

You can use the subset_draws() function, specifying the iteration argument.

# Your code here

Thinning draws

Thin the draws so that only half of the draws are included. Then try automatic thinning.

You can use the function thin_draws(), optionally specifying the thin argument.

# Your code here

Summarising draws

Extract the variable mu and summarise it via mean and sd.

You can use the functions subset_draws(), specifying the variable argument and summarise_draws(), specifying the functions.

# Your code here

Creating new variables

Create a new variable that is the difference between school 1 and school 2 means. Call it diff_1_2. Do the same with school 3 and school 4. Then summarise these new variables with median, 0.3 and 0.7 quantiles.

You can use the functions mutate_variables(), summarise_draws() and quantile2()

# Your code here

Marginal posteriors

Plot the marginal posteriors for the school means.

You can use the functions mcmc_intervals() or mcmc_areas().

# Your code here

Pairs plot

Plot the variables mu and tau in an pairs plot.

You can use the functions mcmc_pairs()

# Your code here

Part 2: Convergence diagnostics and uncertainty

R-hat

Calculate the R-hat for all the variables in the model. Which variables have high R-hat (> 1.01)?

You can use the function summarise_draws() with rhat

# Your code here

Effective sample size (ESS)

Calculate the bulk and tail ESS for all the variables in the model.

You can use the function summarise_draws() with ess_tail and ess_bulk

# Your code here

Monte Carlo standard error

Calculate the mean of each variable, and also the Monte Carlo standard error of the mean.

You can use the functions summarise_draws() with mean and mcse_mean

# Your code here

Then do the same for the 0.05 and 0.95 quantiles. Think about how the Monte Carlo standard error might influence how you report the quantiles.

You can use the functions summarise_draws() with quantile2 and mcse_quantile

# Your code here

Pareto diagnostics

Calculate the minimum sample size for stable estimates for each variable in the model. Which has the highest minimum sample size?

You can use the function summarise_draws() with pareto_min_ss

# Your code here

Part 3: Model evaluation and critique

We can generate prior predictions from the eight schools model, using the following function.

eight_schools_prior <- function(ndraws, mu_prior_sd = 1, tau_prior_sd = 1, sigma_obs = c(15, 10, 16, 11, 9, 11, 10, 18)) {

  J <- length(sigma_obs)

  mu <- rnorm(ndraws, 0, mu_prior_sd)
  tau <- abs(rnorm(ndraws, 0, tau_prior_sd))

  theta <- sapply(1:J, \(j) rnorm(ndraws, mu, tau))
  yrep <- sapply(1:J, \(j) rnorm(ndraws, theta[, j], sigma_obs[j])) 
  
  prior_draws <- as_draws(yrep)
  variables(prior_draws) <- paste0("yrep[", 1:J, "]")
  
  return(prior_draws)
}

Generate 1000 prior predictive draws, and plot the distributions for each school. Try with different mu_prior_sd and tau_prior_sd values (e.g. 1, 10, 100).

You can use the ppd_intervals function.

prior_predictive_draws <- eight_schools_prior(
  ndraws = 1000,
  mu_prior_sd = 1,
  tau_prior_sd = 1
)

# Your code here

Posterior predictive checks

We can create posterior predictive draws from our posterior draws and plot against our actual observations.

sigma_obs <- c(15, 10, 16, 11, 9, 11, 10, 18)

eight_schools_post_pred <- eight_schools_draws |>
  mutate_variables(
    `yrep[1]` = rnorm(length(`theta[1]`), `theta[1]`, sigma_obs[1]),
    `yrep[2]` = rnorm(length(`theta[2]`), `theta[2]`, sigma_obs[2]),
    `yrep[3]` = rnorm(length(`theta[3]`), `theta[3]`, sigma_obs[3]),
    `yrep[4]` = rnorm(length(`theta[4]`), `theta[4]`, sigma_obs[4]),
    `yrep[5]` = rnorm(length(`theta[5]`), `theta[5]`, sigma_obs[5]),
    `yrep[6]` = rnorm(length(`theta[6]`), `theta[6]`, sigma_obs[6]),
    `yrep[7]` = rnorm(length(`theta[7]`), `theta[7]`, sigma_obs[7]),
    `yrep[8]` = rnorm(length(`theta[8]`), `theta[8]`, sigma_obs[8])
  ) |>
  subset_draws(variable = "yrep") |>
  as_draws_matrix()

# Or alternatively in a nicer way, using rvars

eight_schools_post_pred <- eight_schools_draws |>
  as_draws_rvars() |>
  mutate_variables(
    yrep = rvar_rng(
      rnorm,
      n = length(sigma_obs),
      mean = theta,
      sd = sigma_obs
    )
  ) |>
  subset_draws(variable = "yrep") |>
  as_draws_matrix()

Plot the posterior predictions on top of the actual observations. Then use the PIT-ECDF plot.

You can use the ppc_intervals() and ppc_pit_ecdf() functions.

y <- c(28,  8, -3,  7, -1,  1, 18, 12)
# Your code here

Sensitivity checks

Check for prior and likelihood sensitivity in the eight schools model. First check by power-scaling all priors jointly, then select only the “mu” and only the “tau” prior separately.

You can use powerscale_sensitivity() with the prior_selection argument.

eight_schools_draws_ps <- example_powerscale_model("eight_schools")$draws
# Your code here

Next plot sensitivity using as density plots. Plot only the mu and tau variables.

# Your code here

Part 4: Bringing it all together

We have provided four sets of posterior draws from hierarchical models of observed migratory bird counts recorded between 2000 and 2020 at the Hanko Bird Observatory (Halias).

For species \(j\),

\(y \sim \mathrm{Poisson}(\lambda_j)\)

or

\(y \sim \mathrm{NegativeBinomial}(\lambda_j, \phi).\)

The species-specific abundances are linked through a hierarchical prior,

\(\log(\lambda_j) \sim \mathrm{Normal}(\mu, \sigma)\)

where \(\mu\) represents the average abundance across species and \(\sigma\) controls the amount of pooling between species.

Your task is to explore the posterior draws and diagnostic outputs for the four fitted models.

As you work through the diagnostics, try to identify which model corresponds to each of the following situations:

• Convergence issues caused by poor chain mixing (for example, insufficient warmup).
• Inadequate fit to the data caused by the choice of observation model, shown by posterior predictive checks.
• Issues caused by priors in conflict with the likelihood.
• No major issues, although there is still room for model improvement.

Use posterior summaries, convergence diagnostics, posterior predictive checks, and sensitivity analyses to guide your investigation.

You can also look at the Stan model (birds_per_year.stan) and consider how you might improve it.

# Use the following as a starting point

# Load the draws
draws1 <- readRDS("data/draws1.rds")
draws2 <- readRDS("data/draws2.rds")
draws3 <- readRDS("data/draws3.rds")
draws4 <- readRDS("data/draws4.rds")

# Load the data fed to Stan
stan_data <- readRDS("data/stan_data.rds")

# Get the actual observations and the corresponding species
y <- stan_data$count
species <- stan_data$species # useful for grouped plots

# Create yrep matrices for posterior predictive checks
yrep1 <-  as_draws_matrix(subset_draws(draws1, "y_rep"))
yrep2 <-  as_draws_matrix(subset_draws(draws2, "y_rep"))
yrep3 <-  as_draws_matrix(subset_draws(draws3, "y_rep"))
yrep4 <-  as_draws_matrix(subset_draws(draws4, "y_rep"))