24 Count Predicted Variable

Consider a situation in which we observe two nominal values for every item measured…. Across the whole sample, the result is a table of counts for each combination of values of the nominal variables. The counts are what we are trying to predict and the nominal variables are the predictors. This is the type of situation addressed in this chapter….

In the context of the generalized linear model (GLM) introduced in Chapter 15, this chapter’s situation involves a predicted value that is a count, for which we will use an inverse-link function that is exponential along with a Poisson distribution for describing noise in the data (Kruschke, 2015, pp. 703–704)

24.1 Poisson exponential model

Following Kruschke, we will “refer to the model that will be explained in this section as Poisson exponential because, as we will see, the noise distribution is a Poisson distribution and the inverse-link function is exponential” (p. 704).

24.1.1 Data structure.

Kruschke has the Snee (1974) data for Table 24.1 saved as the HairEyeColor.csv file.

library(tidyverse)
library(janitor)

my_data <- read_csv("data.R/HairEyeColor.csv")

glimpse(my_data)
## Rows: 16
## Columns: 3
## $ Hair  <chr> "Black", "Black", "Black", "Black", "Blond", "Blond", "Blond", "…
## $ Eye   <chr> "Blue", "Brown", "Green", "Hazel", "Blue", "Brown", "Green", "Ha…
## $ Count <dbl> 20, 68, 5, 15, 94, 7, 16, 10, 84, 119, 29, 54, 17, 26, 14, 14

In order to retain some of the ordering in Table 24.1, we’ll want to make Hair a factor and recode Brown as Brunette.

my_data <-
  my_data %>% 
  mutate(Hair = if_else(Hair == "Brown", "Brunette", Hair) %>% 
           factor(., levels = c("Black", "Brunette", "Red", "Blond")))

Here’s a quick way to use pivot_wider() to make most of the table.

my_data %>% 
  pivot_wider(names_from = Hair,
              values_from = Count)
## # A tibble: 4 x 5
##   Eye   Black Blond Brunette   Red
##   <chr> <dbl> <dbl>    <dbl> <dbl>
## 1 Blue     20    94       84    17
## 2 Brown    68     7      119    26
## 3 Green     5    16       29    14
## 4 Hazel    15    10       54    14

However, that didn’t get us the marginal totals. For those, we’ll uncount() the cells in the data and then make the full table with janitor::tabyl() and janitor::adorn_totals().

my_data %>%
  uncount(weights = Count, .remove = F) %>% 
  tabyl(Eye, Hair) %>%
  adorn_totals(c("row", "col")) %>% 
  knitr::kable()
Eye Black Brunette Red Blond Total
Blue 20 84 17 94 215
Brown 68 119 26 7 220
Green 5 29 14 16 64
Hazel 15 54 14 10 93
Total 108 286 71 127 592

That last knitr::kable() line just formatted the output a bit.

24.1.3 Poisson noise distribution.

Simon-Denis Poisson’s distribution follows the form

\[p(y | \lambda) = \frac{\lambda^y \exp (-\lambda)}{y!},\]

where \(y\) is a non-negative integer and \(\lambda\) is a non-negative real number. The mean of the Poisson distribution is \(\lambda\). Importantly, the variance of the Poisson distribution is also \(\lambda\) (i.e., the standard deviation is \(\sqrt \lambda\)). Thus, in the Poisson distribution, the variance is completely yoked to the mean. (p. 707)

We can work with that expression directly in base R. Here we use \(\lambda = 5.5\) and \(y = 2\).

lambda <- 5.5
y      <- 2

lambda^y * exp(-lambda) / factorial(y)
## [1] 0.06181242

If we’d like to simulate from the Poisson distribution, we’d use the rpois() function. It takes two arguments, n and lambda. Let’s generate 1,000 draws based on \(\lambda = 5\).

set.seed(24)

d <- tibble(y = rpois(n = 1000, lambda = 5))

Here are the mean and variance.

d %>% 
  summarise(mean     = mean(y),
            variance = var(y))
## # A tibble: 1 x 2
##    mean variance
##   <dbl>    <dbl>
## 1  4.98     5.17

They’re not exactly the same because of simulation variance, but they get that way real quick with a larger sample.

set.seed(24)

tibble(y = rpois(n = 100000, lambda = 5)) %>% 
  summarise(mean     = mean(y),
            variance = var(y))
## # A tibble: 1 x 2
##    mean variance
##   <dbl>    <dbl>
## 1  4.99     4.98

Let’s put rpois() to work and make Figure 24.1.

Before we go any further, let’s talk color and theme. For this chapter, we’ll take our color palette from the ochRe package (Allan et al., 2017), which provides several color palettes inspired by the art, landscapes, and wildlife of Australia. For the plots to follow, we will use the "healthy_reef" palette.

library(ochRe)

scales::show_col(ochre_palettes[["healthy_reef"]])

hr <- ochre_palettes[["healthy_reef"]]

Our overall plot theme will be based on cowplot::theme_minimal_grid() with a few adjustments.

library(cowplot)

theme_set(
  theme_minimal_grid() +
    theme(legend.key = element_rect(fill = hr[2], color = hr[2]),
          panel.background = element_rect(fill = hr[2], color = hr[2]),
          panel.border = element_rect(color = hr[2]),
          panel.grid = element_blank())
)

Now we have our palette and theme, we’re ready to make our version of the Poisson histograms in Figure 24.1.

set.seed(24)

tibble(lambda = c(1.8, 8.3, 32.1)) %>% 
  mutate(y = map(lambda, rpois, n = 1e5)) %>% 
  unnest(y) %>%
  
  ggplot(aes(x = y)) +
  geom_histogram(aes(y = stat(density)),
                 binwidth = 1, fill = hr[3]) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.05))) +
  scale_y_continuous("p(y)", expand = expansion(mult = c(0, 0.05))) +
  facet_wrap(~ lambda, ncol = 1,
             labeller = label_bquote(dpois(y*"|"*lambda == .(lambda))))

For more on our labeller = label_bquote syntax, check out this.

But anyway, given \(\lambda\), the Poisson distribution gives the probabilities of specific non-negative integers. And instead of using our hand-coded function from above, we can also use dpois() to get precise density values.

dpois(2, lambda = 5.5)
## [1] 0.06181242

24.1.4 The complete model and implementation in JAGS brms.

To get a sense of Kruschke’s hierarchical Bayesian alternative to the conventional ANOVA approach, let’s make our version of the model diagram in Figure 24.2.

library(patchwork)

# bracket
p1 <-
  tibble(x = .95,
         y = .5,
         label = "{_}") %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 10, hjust = 1, color = hr[4], family = "Times") +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) +
  ylim(0, 1) +
  theme_void()

##  plain arrow
# save our custom arrow settings
my_arrow <- arrow(angle = 20, length = unit(0.35, "cm"), type = "closed")
p2 <-
  tibble(x    = .68,
         y    = 1,
         xend = .68,
         yend = .25) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = hr[8]) +
  xlim(0, 1) +
  theme_void()

# normal density
p3 <-
  tibble(x = seq(from = -3, to = 3, by = .1)) %>% 
  ggplot(aes(x = x, y = (dnorm(x)) / max(dnorm(x)))) +
  geom_area(fill = hr[2]) +
  annotate(geom = "text",
           x = 0, y = .2,
           label = "normal",
           size = 7, color = hr[8]) +
  annotate(geom = "text",
           x = c(0, 1.45), y = .6,
           hjust = c(.5, 0),
           label = c("italic(M)[0]", "italic(S)[0]"), 
           size = 7, color = hr[8], family = "Times", parse = T) +
  scale_x_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = hr[8]))

# second normal density
p4 <-
  tibble(x = seq(from = -3, to = 3, by = .1)) %>% 
  ggplot(aes(x = x, y = (dnorm(x)) / max(dnorm(x)))) +
  geom_area(fill = hr[2]) +
  annotate(geom = "text",
           x = 0, y = .2,
           label = "normal",
           size = 7, color = hr[8]) +
  annotate(geom = "text",
           x = c(0, 1.45), y = .6,
           hjust = c(.5, 0),
           label = c("0", "sigma[beta][1]"), 
           size = 7, color = hr[8], family = "Times", parse = T) +
  scale_x_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = hr[8]))

# third normal density
p5 <-
  tibble(x = seq(from = -3, to = 3, by = .1)) %>% 
  ggplot(aes(x = x, y = (dnorm(x)) / max(dnorm(x)))) +
  geom_area(fill = hr[2]) +
  annotate(geom = "text",
           x = 0, y = .2,
           label = "normal",
           size = 7, color = hr[8]) +
  annotate(geom = "text",
           x = c(0, 1.45), y = .6,
           hjust = c(.5, 0),
           label = c("0", "sigma[beta][2]"), 
           size = 7, color = hr[8], family = "Times", parse = T) +
  scale_x_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = hr[8]))

# fourth normal density
p6 <-
  tibble(x = seq(from = -3, to = 3, by = .1)) %>% 
  ggplot(aes(x = x, y = (dnorm(x)) / max(dnorm(x)))) +
  geom_area(fill = hr[2]) +
  annotate(geom = "text",
           x = 0, y = .2,
           label = "normal",
           size = 7, color = hr[8]) +
  annotate(geom = "text",
           x = 0, y = .6,
           hjust = .5,
           label = "0", 
           size = 7, color = hr[8], family = "Times", parse = T) +
  scale_x_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = hr[8]))

# likelihood formula
p7 <-
  tibble(x = .5,
         y = .25,
         label = "exp(beta[0]+sum()[italic(j)]*beta[1]['['*italic(j)*']']*italic(x)[1]['['*italic(j)*']'](italic(i))+sum()[italic(k)]*beta[2]['['*italic(k)*']']*italic(x)[2]['['*italic(k)*']'](italic(i))+sum()[italic(jk)]*beta[1%*%2]['['*italic(jk)*']']*italic(x)[1%*%2]['['*italic(jk)*']'](italic(i)))") %>% 
 
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 6.75, color = hr[8], parse = T, family = "Times") +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) +
  ylim(0, 1) +
  theme_void()

# sigma
p8 <-
  tibble(x = .13,
         y = .6,
         label = "sigma[beta][1%*%2]") %>% 
 
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(hjust = 0, size = 7, color = hr[8], parse = T, family = "Times") +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 2)) +
  ylim(0, 1) +
  theme_void()

# four annotated arrows
p9 <-
  tibble(x    = c(.05, .34, .64, .945),
         y    = c(1, 1, 1, 1),
         xend = c(.06, .2, .47, .77),
         yend = c(0, 0, 0, 0)) %>%
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = hr[8]) +
  annotate(geom = "text",
           x = c(.025, .24, .30, .53, .60, .83, .91), y = .5,
           label = c("'~'", "'~'", "italic(j)", "'~'", "italic(k)", "'~'", "italic(jk)"),
           size = c(10, 10, 7, 10, 7, 10, 7), 
           color = hr[8], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

# Poisson density
p10 <-
  tibble(x = 0:9) %>% 
  ggplot(aes(x = x, y = .87 * (dpois(x, lambda = 3.5)) / max(dpois(x, lambda = 3.5))) ) +
  geom_col(color = hr[2], fill = hr[2], width = .45) +
  annotate(geom = "text",
           x = 4.5, y = .2,
           label = "Poisson",
           size = 7, color = hr[8]) +
  annotate(geom = "text",
           x = 4.5, y = .89,
           label = "lambda[italic(i)]",
           size = 7, color = hr[8], family = "Times", parse = TRUE) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = hr[8]))

# one annotated arrow
p11 <-
  tibble(x    = .5,
         y    = 1,
         xend = .5,
         yend = .35) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = hr[8]) +
  annotate(geom = "text",
           x = .4, y = .725,
           label = "'='",
           size = 10, color = hr[8], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

# the final annotated arrow
p12 <-
  tibble(x     = c(.375, .625),
         y     = c(1/3, 1/3),
         label = c("'~'", "italic(i)")) %>%

  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = c(10, 7),
            color = hr[8], parse = T, family = "Times") +
  geom_segment(x = .5, xend = .5,
               y = 1, yend = 0,
               arrow = my_arrow, color = hr[8]) +
  xlim(0, 1) +
  theme_void()

# some text
p13 <-
  tibble(x     = .5,
         y     = .5,
         label = "italic(y[i])") %>%

  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 7, color = hr[8], parse = T, family = "Times") +
  xlim(0, 1) +
  theme_void()

# define the layout
layout <- c(
  area(t = 1, b = 1, l = 6, r = 7),
  area(t = 1, b = 1, l = 10, r = 11),
  area(t = 1, b = 1, l = 14, r = 15),
  area(t = 2, b = 3, l = 6, r = 7),
  area(t = 2, b = 3, l = 10, r = 11),
  area(t = 2, b = 3, l = 14, r = 15),
  area(t = 3, b = 4, l = 1, r = 3),
  area(t = 3, b = 4, l = 5, r = 7),
  area(t = 3, b = 4, l = 9, r = 11),
  area(t = 3, b = 4, l = 13, r = 15),
  area(t = 6, b = 7, l = 1, r = 15),
  area(t = 3, b = 4, l = 15, r = 16),
  area(t = 5, b = 6, l = 1, r = 15),
  area(t = 10, b = 11, l = 7, r = 9),
  area(t = 8, b = 10, l = 7, r = 9),
  area(t = 12, b = 12, l = 7, r = 9),
  area(t = 13, b = 13, l = 7, r = 9)
)

# combine and plot!
(p1 + p1 + p1 + p2 + p2 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + p11 + p12 + p13) + 
  plot_layout(design = layout) &
  ylim(0, 1) &
  theme(plot.margin = margin(0, 5.5, 0, 5.5))

Using Kruschke’s method,

the prior is supposed to be broad on the scale of the data, but we must be careful about what scale is being modeled by the baseline and deflections. The counts are being directly described by \(\lambda\), but it is \(\log (\lambda)\) being described by the baseline and deflections. Thus, the prior on the baseline and deflections should be broad on the scale of the logarithm of the data. To establish a generic baseline, consider that if the data points were distributed equally among the cells, the mean count would be the total count divided by the number of cells. The biggest possible standard deviation across cells would occur when all the counts were loaded into a single cell and all the other cells were zero. (pp. 709–710, emphasis added)

Before we show how to fit the model, we need the old gamma_a_b_from_omega_sigma() function.

gamma_a_b_from_omega_sigma <- function(mode, sd) {
  
  if (mode <= 0) stop("mode must be > 0")
  if (sd   <= 0) stop("sd must be > 0")
  rate <- (mode + sqrt(mode^2 + 4 * sd^2)) / (2 * sd^2)
  shape <- 1 + mode * rate
  
  return(list(shape = shape, rate = rate))
  
}

Here are a few intermediate values before we set the stanvars.

n_x1_level <- length(unique(my_data$x1))
n_x2_level <- length(unique(my_data$x2))
n_cell     <- nrow(my_data)

Now we’re ready to define the stanvars.

y_log_mean <-
  log(sum(my_data$y) / (n_x1_level * n_x2_level))

y_log_sd <-
  c(rep(0, n_cell - 1), sum(my_data$y)) %>% 
  sd() %>% 
  log()

s_r <- gamma_a_b_from_omega_sigma(mode = y_log_sd, sd = 2 * y_log_sd)

stanvars <- 
  stanvar(y_log_mean, name = "y_log_mean") + 
  stanvar(y_log_sd,   name = "y_log_sd") + 
  stanvar(s_r$shape,  name = "alpha") +
  stanvar(s_r$rate,   name = "beta")

You’d then fit a Poisson model with two nominal predictors using Kruschke’s hierarchical-shrinkage method like this.

fit <-
  brm(data = my_data,
      family = poisson,
      y ~ 1 + (1 | x1) + (1 | x2) + (1 | x1:x2),
      prior = c(prior(normal(y_log_mean, y_log_sd * 2), class = Intercept),
                prior(gamma(alpha, beta), class = sd)),
      stanvars = stanvars)

By brms default, family = poisson uses the log link. Thus family = poisson(link = "log") should return the same results. Notice the right-hand side of the model formula. We have three hierarchical variance parameters. This hierarchical-shrinkage approach to frequency-table data has its origins in Gelman (2005), Analysis of variance–why it is more important than ever.

24.2 Example: Hair eye go again

We’ll be using the same data, from above. As an alternative to Table 24.1, it might be handy to take a more colorful approach to wading into the data.

# wrangle
my_data %>%
  uncount(weights = Count, .remove = F) %>% 
  tabyl(Eye, Hair) %>%
  adorn_totals(c("row", "col")) %>% 
  data.frame() %>% 
  pivot_longer(-Eye, names_to = "Hair") %>% 
  mutate(Eye = fct_rev(Eye)) %>% 
  
  # plot
  ggplot(aes(x = Hair, y = Eye, label = value)) +
  geom_raster(aes(fill = value)) +
  geom_text(aes(color = value < 250)) +
  scale_fill_gradient(low = hr[9], high = hr[2]) +
  scale_color_manual(values = hr[c(9, 2)]) +
  scale_x_discrete(expand = c(0, 0), position = "top") +
  scale_y_discrete(expand = c(0, 0)) +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks = element_blank(),
        legend.position = "none")

Load the brms and tidybayes packages.

library(brms)
library(tidybayes)

Now we’ll save the preparatory values necessary for the stanvars.

n_x1_level <- length(unique(my_data$Hair))
n_x2_level <- length(unique(my_data$Eye))
n_cell     <- nrow(my_data)

n_x1_level
## [1] 4
n_x2_level
## [1] 4
n_cell
## [1] 16

Here are the values we’ll save as stanvars.

y_log_mean <-
  log(sum(my_data$Count) / (n_x1_level * n_x2_level))

y_log_sd <-
  c(rep(0, n_cell - 1), sum(my_data$Count)) %>% 
  sd() %>% 
  log()

s_r <- gamma_a_b_from_omega_sigma(mode = y_log_sd, sd = 2 * y_log_sd)

y_log_mean
## [1] 3.610918
y_log_sd
## [1] 4.997212
s_r$shape
## [1] 1.640388
s_r$rate
## [1] 0.1281491

As a quick detour, it might be interesting to see what the kind of gamma distribution is entailed by those last two values.

tibble(x = seq(from = 0, to = 70, length.out = 1e3)) %>% 
  mutate(density = dgamma(x, s_r$shape, s_r$rate)) %>% 
  
  ggplot(aes(x = x, y = density)) +
  geom_area(fill = hr[4]) +
  scale_x_continuous(NULL, expand = expansion(mult = c(0, -0.05))) +
  scale_y_continuous(NULL, breaks = NULL, expand = expansion(mult = c(0, 0.05))) +
  ggtitle(expression("Kruschke's wide prior for "*sigma[beta*x]))

Save the stanvars.

stanvars <- 
  stanvar(y_log_mean, name = "y_log_mean") + 
  stanvar(y_log_sd,   name = "y_log_sd") + 
  stanvar(s_r$shape,  name = "alpha") +
  stanvar(s_r$rate,   name = "beta")

Fit Kruschke’s model with brms.

fit24.1 <-
  brm(data = my_data,
      family = poisson,
      Count ~ 1 + (1 | Hair) + (1 | Eye) + (1 | Hair:Eye),
      prior = c(prior(normal(y_log_mean, y_log_sd * 2), class = Intercept),
                prior(gamma(alpha, beta), class = sd)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 24,
      stanvars = stanvars,
      file = "fits/fit24.01")

As it turns out, if you try to fit Kruschke’s model with brms as is, you’ll run into difficulties with divergent transitions and the like. One approach is to try tuning the adapt_delta and max_treedepth parameters. I had no luck with that approach. E.g., cranking adapt_delta up past 0.9999 still returned a divergent transition or two.

Another approach is to step back and assess the model. We’re trying to fit a multilevel model with two grouping variables and their interaction with a total of 16 data points. That’s not a lot of data for fitting such a model. If you take a close look at our priors, you’ll notice they’re really quite weak. If you’re willing to tighten them up just a bit, the model can fit more smoothly. That will be our approach.

For the \(\sigma\) hyperparameter of the overall intercept’s Gaussian prior, Kruschke would have us multiply y_log_sd by 2. Here we’ll tighten up that \(\sigma\) hyperparameter by simply setting it to y_log_sd. The gamma priors for the upper-level variance parameters were based on a mode of y_log_sd and a standard deviation of the same but multiplied by 2 (i.e., 2 * y_log_sd). We’ll tighten that up a bit by simply defining those gammas by a standard deviation of y_log_sd. When you make those adjustments, the model fits with less fuss. In case you’re curious, here is what those priors look like.

# redifine our shape and rate
s_r <- gamma_a_b_from_omega_sigma(mode = y_log_sd, sd = y_log_sd)

# wrangle
bind_rows(
  # define beta[0]
  tibble(x = seq(from = y_log_mean - (4 * y_log_sd), to = y_log_mean + (4 * y_log_sd), length.out = 1e3)) %>%
    mutate(density = dnorm(x, y_log_mean, y_log_sd)),
  # define sigma[beta[x]]
  tibble(x = seq(from = 0, to = 40, length.out = 1e3)) %>% 
    mutate(density = dgamma(x, s_r$shape, s_r$rate))
) %>%
  mutate(prior = rep(c("beta[0]", "sigma[beta*x]"), each = n() / 2)) %>% 
  
  # plot
  ggplot(aes(x = x, y = density)) +
  geom_area(fill = hr[6]) +
  scale_y_continuous(NULL, breaks = NULL, expand = expansion(mult = c(0, 0.05))) +
  labs(title = "Priors",
       x = NULL) +
  facet_wrap(~ prior, scales = "free", labeller = label_parsed)

Update the stanvars.

stanvars <- 
  stanvar(y_log_mean, name = "y_log_mean") + 
  stanvar(y_log_sd,   name = "y_log_sd") + 
  stanvar(s_r$shape,  name = "alpha") +
  stanvar(s_r$rate,   name = "beta")

Now we’ve updated our stanvars, we’ll fit the modified model. We should note that even this version required some adjustment to the adapt_delta and max_treedepth parameters. But it wasn’t nearly the exercise in frustration entailed in the version, above.

fit24.1 <-
  brm(data = my_data,
      family = poisson,
      Count ~ 1 + (1 | Hair) + (1 | Eye) + (1 | Hair:Eye),
      prior = c(prior(normal(y_log_mean, y_log_sd), class = Intercept),
                prior(gamma(alpha, beta), class = sd)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999,
                     max_treedepth = 11.5),
      stanvars = stanvars,
      seed = 24,
      file = "fits/fit24.01")

Take a look at the parameter summary.

print(fit24.1)
##  Family: poisson 
##   Links: mu = log 
## Formula: Count ~ 1 + (1 | Hair) + (1 | Eye) + (1 | Hair:Eye) 
##    Data: my_data (Number of observations: 16) 
## Samples: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~Eye (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.29      1.47     0.13     5.16 1.00     2492     3672
## 
## ~Hair (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.36      1.37     0.20     5.04 1.00     2806     4091
## 
## ~Hair:Eye (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.89      0.25     0.52     1.49 1.00     3617     4791
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     3.18      1.27     0.58     5.72 1.00     3800     3274
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

You’ll notice that even though we tightened up the priors, the parameter estimates are still quite small relative to the values they allowed for. So even our tightened priors were quite permissive.

Let’s post process in preparation for Figure 24.3.

nd <-
  my_data %>% 
  arrange(Eye, Hair) %>% 
  # make the titles for the facet strips
  mutate(strip = str_c("E:", Eye, " H:", Hair, "\nN = ", Count))

f <-
  fitted(fit24.1,
         newdata = nd,
         summary = F) %>% 
  data.frame() %>%
  set_names(pull(nd, strip))

glimpse(f)
## Rows: 8,000
## Columns: 16
## $ `E:Blue H:Black\nN = 20`      <dbl> 16.12780, 19.34630, 22.69872, 26.77390, …
## $ `E:Blue H:Brunette\nN = 84`   <dbl> 91.02143, 90.52450, 91.02996, 71.96672, …
## $ `E:Blue H:Red\nN = 17`        <dbl> 24.02637, 17.39775, 21.90677, 18.44875, …
## $ `E:Blue H:Blond\nN = 94`      <dbl> 83.86624, 98.59621, 96.58789, 89.52365, …
## $ `E:Brown H:Black\nN = 68`     <dbl> 78.77468, 64.90828, 49.26054, 70.48308, …
## $ `E:Brown H:Brunette\nN = 119` <dbl> 127.0663, 120.1429, 121.2332, 121.8616, …
## $ `E:Brown H:Red\nN = 26`       <dbl> 29.81537, 17.18848, 25.97654, 20.65013, …
## $ `E:Brown H:Blond\nN = 7`      <dbl> 8.447764, 18.973227, 15.766449, 6.808252…
## $ `E:Green H:Black\nN = 5`      <dbl> 11.187876, 2.042147, 9.896912, 4.305856,…
## $ `E:Green H:Brunette\nN = 29`  <dbl> 27.07728, 27.94305, 33.49773, 15.13343, …
## $ `E:Green H:Red\nN = 14`       <dbl> 8.540725, 15.668955, 6.798632, 23.700314…
## $ `E:Green H:Blond\nN = 16`     <dbl> 13.671059, 12.225308, 12.223708, 23.9272…
## $ `E:Hazel H:Black\nN = 15`     <dbl> 11.43878, 15.80869, 11.21282, 32.82198, …
## $ `E:Hazel H:Brunette\nN = 54`  <dbl> 42.31558, 41.24963, 44.52099, 51.63807, …
## $ `E:Hazel H:Red\nN = 14`       <dbl> 18.21927, 13.09931, 15.40381, 15.26440, …
## $ `E:Hazel H:Blond\nN = 10`     <dbl> 12.707261, 15.630540, 15.094075, 8.44273…

Notice that when working with a Poisson model, fitted() defaults to returning estimates in the \(\lambda\) metric. If we want proportions/probabilities, we’ll have to compute them by dividing by the total \(N\). In this case \(N = 592\), which we get with sum(my_data$Count). Here we convert the data to the long format, compute the proportions, and plot to make the top portion of Figure 24.3.

f %>% 
  gather(key, count) %>% 
  mutate(proportion = count / 592) %>% 
  
  ggplot(aes(x = proportion, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[9], color = hr[5],
                     normalize = "panels") +
  scale_x_continuous(breaks = c(0, .1, .2)) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(0, .25)) +
  facet_wrap(~ key, scales = "free_y")

We’ll have to work a bit to get the deflection differences. If this was a simple multilevel model with a single random grouping variable, we could just use the ranef() function to return the deflections. Like fitted(), it’ll return summaries by default. But you can get the posterior draws with the summary = F argument. But since we used two grouping variables and their interaction, it’d be a bit of a pain to work that way. Happily, we do have a handy alternative. First, if we use the scale = "linear" argument, fitted() will return the draws in the \(\lambda\) scale rather than the original count metric. With the group-level draws in the \(\lambda\) metric, all we need to do is subtract the fixed effect (i.e., the grand mean, the population estimate) from each to convert them to the deflection metric. So below, we’ll

  1. make a custom make_deflection() function to do the conversions,
  2. redefine our nd data to make our naming conventions a little more streamlined,
  3. use fitted() and its scale = "linear" argument to get the draws in the \(\lambda\) metric,
  4. wrangle a touch, and
  5. use our handy make_deflection() function to convert the results to the deflection metric.

I know; that’s a lot. If you get lost, just go step by step and examine the results along the way.

# a. make a custom function
make_deflection <- function(x) {
  x - posterior_samples(fit24.1)$b_Intercept
}

# b. streamline `nd`
nd <-
  my_data %>% 
  arrange(Eye, Hair) %>% 
  mutate(strip = str_c("E:", Eye, " H:", Hair))

# c. use `fitted()`
deflections <-
  fitted(fit24.1,
         newdata = nd,
         summary = F,
         scale = "linear") %>% 
  # d. wrangle
  data.frame() %>%
  set_names(pull(nd, strip)) %>% 
  # e. use the `make_deflection()` function
  mutate_all(.funs = make_deflection)

# what have we done?
glimpse(deflections)
## Rows: 8,000
## Columns: 16
## $ `E:Blue H:Black`     <dbl> -0.19039053, 0.02047397, -0.44515915, 0.11303471,…
## $ `E:Blue H:Brunette`  <dbl> 1.5401600, 1.5635936, 0.9437209, 1.1018109, 1.825…
## $ `E:Blue H:Red`       <dbl> 0.20821721, -0.08568613, -0.48067188, -0.25939634…
## $ `E:Blue H:Blond`     <dbl> 1.4582882, 1.6490058, 1.0029856, 1.3201099, 1.988…
## $ `E:Brown H:Black`    <dbl> 1.3956567, 1.2309481, 0.3296556, 1.0809798, 1.489…
## $ `E:Brown H:Brunette` <dbl> 1.87377396, 1.84665461, 1.23024816, 1.62849272, 2…
## $ `E:Brown H:Red`      <dbl> 0.42408911, -0.09778779, -0.31027403, -0.14667113…
## $ `E:Brown H:Blond`    <dbl> -0.837033216, 0.001001884, -0.809583559, -1.25625…
## $ `E:Green H:Black`    <dbl> -0.55610429, -2.22802520, -1.27524499, -1.7144169…
## $ `E:Green H:Brunette` <dbl> 0.327760120, 0.388141435, -0.055989959, -0.457486…
## $ `E:Green H:Red`      <dbl> -0.826089091, -0.190345635, -1.650746385, -0.0089…
## $ `E:Green H:Blond`    <dbl> -0.3556538377, -0.4385187786, -1.0640904068, 0.00…
## $ `E:Hazel H:Black`    <dbl> -0.5339256068, -0.1814672463, -1.1504103112, 0.31…
## $ `E:Hazel H:Brunette` <dbl> 0.774220352, 0.777615105, 0.228493070, 0.76986626…
## $ `E:Hazel H:Red`      <dbl> -0.0684549414, -0.3694673118, -0.8328527281, -0.4…
## $ `E:Hazel H:Blond`    <dbl> -0.42876136, -0.19280029, -0.85316548, -1.0410869…

Now we’re ready to define our difference columns and plot our version of the lower panels in Figure 24.3.

deflections %>% 
  transmute(`Blue − Brown @ Black` = `E:Blue H:Black` - `E:Brown H:Black`,
            `Blue − Brown @ Blond` = `E:Blue H:Blond` - `E:Brown H:Blond`) %>% 
  mutate(`Blue.v.Brown\n(x)\nBlack.v.Blond` = `Blue − Brown @ Black` - `Blue − Brown @ Blond`) %>% 
  pivot_longer(everything(), values_to = "difference") %>% 

  ggplot(aes(x = difference, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[9], color = hr[5],
                    normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~ name, scales = "free")

If you’re curious, here are the precise summary values.

deflections %>% 
  mutate(`Blue − Brown @ Black` = `E:Blue H:Black` - `E:Brown H:Black`,
         `Blue − Brown @ Blond` = `E:Blue H:Blond` - `E:Brown H:Blond`) %>% 
  mutate(`Blue.v.Brown\n(x)\nBlack.v.Blond` = `Blue − Brown @ Black` - `Blue − Brown @ Blond`) %>% 
  pivot_longer(`Blue − Brown @ Black`:`Blue.v.Brown\n(x)\nBlack.v.Blond`) %>% 
  group_by(name) %>% 
  mode_hdi(value)
## # A tibble: 3 x 7
##   name                               value .lower .upper .width .point .interval
##   <chr>                              <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 "Blue − Brown @ Black"             -1.23  -1.68 -0.699   0.95 mode   hdi      
## 2 "Blue − Brown @ Blond"              2.36   1.74  3.15    0.95 mode   hdi      
## 3 "Blue.v.Brown\n(x)\nBlack.v.Blond" -3.54  -4.47 -2.74    0.95 mode   hdi

24.3 Example: Interaction contrasts, shrinkage, and omnibus test

“In this section, we consider some contrived data to illustrate aspects of interaction contrasts. Like the eye and hair data, the fictitious data have two attributes with four levels each” (p. 713). Let’s make the data.

my_data <-
  crossing(a = str_c("a", 1:4),
           b = str_c("b", 1:4)) %>% 
  mutate(count = c(rep(c(22, 11), each = 2) %>% rep(., times = 2),
                   rep(c(11, 22), each = 2) %>% rep(., times = 2)))

head(my_data)
## # A tibble: 6 x 3
##   a     b     count
##   <chr> <chr> <dbl>
## 1 a1    b1       22
## 2 a1    b2       22
## 3 a1    b3       11
## 4 a1    b4       11
## 5 a2    b1       22
## 6 a2    b2       22

In the last section, we covered how Kruschke’s broad priors can make fitting these kinds of models difficult when using HMC, particularly with so few cells. Our solution was to reign in the \(\sigma\) hyperparameter for the level-one intercept and to compute the gamma prior for the hierarchical deflections based on a standard deviation of the log of the maximum standard deviation for the data rather than two times that value. Let’s explore more options.

This data set has 16 cells. With so few cells, one might argue for a more conservative prior on the hierarchical deflections. Why not ditch the gamma altogether for a half normal centered on zero and with a \(\sigma\) hyperparameter of 1? Even though this is much tighter than Kruschke’s gamma prior approach, it’s still permissive on the \(\log\) scale. As for our intercept, we’ll continue with the same approach from last time.

With that in mind, make the stanvars.

n_x1_level <- length(unique(my_data$a))
n_x2_level <- length(unique(my_data$b))
n_cell     <- nrow(my_data)

y_log_mean <-
  log(sum(my_data$count) / (n_x1_level * n_x2_level))

y_log_sd <-
  c(rep(0, n_cell - 1), sum(my_data$count)) %>% 
  sd() %>% 
  log()

stanvars <- 
  stanvar(y_log_mean, name = "y_log_mean") + 
  stanvar(y_log_sd,   name = "y_log_sd")

Just for kicks, let’s take a quick look at our priors.

bind_rows(
  # define beta[0]
  tibble(x = seq(from = y_log_mean - (4 * y_log_sd), to = y_log_mean + (4 * y_log_sd), length.out = 1e3)) %>%
    mutate(density = dnorm(x, y_log_mean, y_log_sd)),
  # define sigma[beta[x]]
  tibble(x = seq(from = 0, to = 5, length.out = 1e3)) %>% 
    mutate(density = dnorm(x, 0, 1))
) %>%
  mutate(prior = rep(c("beta[0]", "sigma[beta*x]"), each = n() / 2)) %>% 
  
  ggplot(aes(x = x, y = density)) +
  geom_area(fill = hr[6]) +
  scale_y_continuous(NULL, breaks = NULL, expand = expansion(mult = c(0, 0.05))) +
  labs(title = "Priors",
       x = NULL) +
  facet_wrap(~ prior, scales = "free", labeller = label_parsed)

Fit the model.

fit24.2 <-
  brm(data = my_data,
      family = poisson,
      count ~ 1 + (1 | a) + (1 | b) + (1 | a:b),
      prior = c(prior(normal(y_log_mean, y_log_sd), class = Intercept),
                prior(normal(0, 1), class = sd)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999),
      stanvars = stanvars,
      seed = 24,
      file = "fits/fit24.02")

Review the summary.

print(fit24.2)
##  Family: poisson 
##   Links: mu = log 
## Formula: count ~ 1 + (1 | a) + (1 | b) + (1 | a:b) 
##    Data: my_data (Number of observations: 16) 
## Samples: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~a (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.16      0.18     0.00     0.67 1.00     2743     4003
## 
## ~a:b (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.31      0.12     0.10     0.57 1.00     2397     2214
## 
## ~b (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.16      0.18     0.00     0.64 1.00     2921     4134
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     2.77      0.20     2.36     3.18 1.00     2966     2456
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

We might plot our \(\sigma_\text{<group>}\) posteriors against our prior to get a sense of how strong it was.

posterior_samples(fit24.2) %>% 
  dplyr::select(starts_with("sd")) %>% 
  # set_names(str_c("expression(sigma", c("*a", "*ab", "*b"), ")")) %>% 
  pivot_longer(everything(), names_to = "sd") %>% 
  
  ggplot(aes(x = value)) +
  # prior
  geom_area(data = tibble(value = seq(from = 0, to = 3.5, by  =.01)),
            aes(y = dnorm(value, 0, 1)),
            fill = hr[8]) +
  # posterior
  geom_density(aes(fill = sd),
               size = 0, alpha = 2/3) +
  scale_fill_manual(values = hr[c(4, 6, 3)],
                    labels = c(expression(sigma[a]),
                                  expression(sigma[ab]),
                                  expression(sigma[b])),
                    guide = guide_legend(label.hjust = 0)) +
  scale_x_continuous(NULL, expand = expansion(mult = c(0, 0.05))) +
  scale_y_continuous(NULL, breaks = NULL, expand = expansion(mult = c(0, 0.05))) +
  theme(legend.background = element_rect(fill = "transparent"),
        legend.position = c(.9, .8))

Our \(\operatorname{Normal}^+ (0, 1)\) prior is that short dark shape in the background. The posteriors are the taller and more colorful mounds in the foreground. Here’s the top part of Figure 24.4.

nd <-
  my_data %>% 
  mutate(strip = str_c("a:", a, " b:", b, "\nN = ", count))

fitted(fit24.2,
       newdata = nd,
       summary = F) %>% 
  data.frame() %>%
  set_names(pull(nd, strip)) %>% 
  pivot_longer(everything(), values_to = "count") %>% 
  mutate(proportion = count / sum(my_data$count)) %>% 
  
  # plot!
  ggplot(aes(x = proportion, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[7], color = hr[6], normalize = "panels") +
  scale_x_continuous(breaks = c(.05, .1, .15)) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(0, .15)) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~ name, scales = "free_y")

Like before, we’ll have to work a bit to get the deflection differences.

# streamline `nd`
nd <-
  my_data %>% 
  mutate(strip = str_c("a:", a, " b:", b))

# use `fitted()`
deflections <-
  fitted(fit24.2,
         newdata = nd,
         summary = F,
         scale = "linear") %>% 
  # wrangle
  data.frame() %>%
  set_names(pull(nd, strip)) %>% 
  # use the `make_deflection()` function
  mutate_all(.funs = make_deflection)

# what have we done?
glimpse(deflections)
## Rows: 8,000
## Columns: 16
## $ `a:a1 b:b1` <dbl> 0.03981493, -0.07178858, -0.36990910, -0.23738049, 0.55729…
## $ `a:a1 b:b2` <dbl> -0.001506227, 0.365939888, -0.436705581, -0.277089109, 0.4…
## $ `a:a1 b:b3` <dbl> -0.4070548, -0.4810801, -0.6599968, -0.7367982, 0.1333308,…
## $ `a:a1 b:b4` <dbl> -0.3936501527, -0.3538763218, -0.9500306448, -0.5856615350…
## $ `a:a2 b:b1` <dbl> -0.171216075, 0.042794837, -0.644010150, -0.395430742, 0.6…
## $ `a:a2 b:b2` <dbl> -0.03820470, 0.43781558, -0.61775787, -0.39377575, 0.30727…
## $ `a:a2 b:b3` <dbl> -0.8794338, -0.4889285, -0.9245606, -0.7630403, 0.1623093,…
## $ `a:a2 b:b4` <dbl> -0.21752585, -0.32058580, -1.03354020, -0.57156581, -0.064…
## $ `a:a3 b:b1` <dbl> -0.53327070, -0.65578508, -0.84520194, -0.79649806, 0.2845…
## $ `a:a3 b:b2` <dbl> -0.44126816, -0.80726314, -0.79396536, -0.77117806, 0.0737…
## $ `a:a3 b:b3` <dbl> 0.17445861, 0.10526627, -0.77478472, -0.24744978, 0.397445…
## $ `a:a3 b:b4` <dbl> 0.17117765, -0.08145822, -0.80838316, -0.20029842, 0.40089…
## $ `a:a4 b:b1` <dbl> -0.7034799, -0.9251032, -0.6594339, -0.3938122, -0.1090155…
## $ `a:a4 b:b2` <dbl> -0.36079857, -0.32747118, -0.59256365, -0.71921276, -0.074…
## $ `a:a4 b:b3` <dbl> 0.009070538, -0.076615980, -0.550780205, -0.161385123, 0.3…
## $ `a:a4 b:b4` <dbl> 0.16546774, -0.01999708, -0.57196644, -0.22746664, 0.50365…

Now we’re ready to define some of the difference columns and plot our version of the leftmost lower panel in Figure 24.4.

deflections %>% 
  transmute(`a2 - a3 @ b2` = `a:a2 b:b2` - `a:a3 b:b2`,
            `a2 - a3 @ b3` = `a:a2 b:b3` - `a:a3 b:b3`) %>% 
  mutate(`a2.v.a3\n(x)\nb2.v.b3` = `a2 - a3 @ b2` - `a2 - a3 @ b3`) %>% 
  pivot_longer(`a2.v.a3\n(x)\nb2.v.b3`, values_to = "difference") %>% 

  ggplot(aes(x = difference, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[7], color = hr[6]) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(-.5, 2.5)) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~ name, scales = "free")

For Figure 24.4, bottom right, we average across the four cells in each quadrant and then compute the contrast.

deflections %>% 
  # in this intermediate step, we compute the quadrant averages
  # `tl` = top left, `br` = bottom right, and so on
  transmute(tl = (`a:a1 b:b1` + `a:a1 b:b2` + `a:a2 b:b1` + `a:a2 b:b2`) / 4,
            tr = (`a:a1 b:b3` + `a:a1 b:b4` + `a:a2 b:b3` + `a:a2 b:b4`) / 4,
            bl = (`a:a3 b:b1` + `a:a3 b:b2` + `a:a4 b:b1` + `a:a4 b:b2`) / 4,
            br = (`a:a3 b:b3` + `a:a3 b:b4` + `a:a4 b:b3` + `a:a4 b:b4`) / 4) %>%
  # compute the contrast of interest
  mutate(`A1.A2.v.A3.A4\n(x)\nB1.B2.v.B3.B4` = (tl - bl) - (tr - br)) %>% 
  pivot_longer(`A1.A2.v.A3.A4\n(x)\nB1.B2.v.B3.B4`, values_to = "difference") %>%
  
  # plot
  ggplot(aes(x = difference, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[7], color = hr[6]) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(-.5, 2.5)) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~ name, scales = "free")

The model presented here has no way to conduct an “ominbus” test of interaction. However, like the ANOVA-style models presented in Chapters 19 and 20, it is easy to extend the model so it has an inclusion coefficient on the interaction deflections. The inclusion coefficient can have values of 0 or 1, and is given a Bernoulli prior. (p. 716)

Like we discussed in earlier chapters, this isn’t a feasible approach for brms. However, we can compare this model with a simpler one that omits the interaction. First, fit the model.

fit24.3 <-
  brm(data = my_data,
      family = poisson,
      count ~ 1 + (1 | a) + (1 | b),
      prior = c(prior(normal(y_log_mean, y_log_sd), class = Intercept),
                prior(normal(0, 1), class = sd)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9999),
      stanvars = stanvars,
      seed = 24,
      file = "fits/fit24.03")

Now we can compare them by their stacking weights.

model_weights(fit24.2, fit24.3) %>% 
  round(digits = 3)
## fit24.2 fit24.3 
##       1       0

Virtually all the weight went to the interaction model. Also, if we step back and ask ourselves what the purpose of an omnibus text of an interaction is for in this context, we’d conclude such a test is asking the question Is \(\sigma_{a \times b}\) the same as zero? Let’s look again at that posterior from fit24.2.

posterior_samples(fit24.2) %>% 
  
  ggplot(aes(x = `sd_a:b__Intercept`, y = 0)) +
  stat_halfeye(.width = .95, fill = hr[7], color = hr[6]) +
  scale_x_continuous(expression(sigma[a%*%b]), 
                     expand = expansion(mult = c(0, 0.05)), limits = c(0, NA)) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(subtitle = "Does this look the same as zero, to you?")

Sure, there’s some uncertainty in that posterior. But that is not zero. We didn’t need an omnibus test or even model comparison via stacking weights to figure that one out.

If you wanted to get fancy with it, we might even do a hierarchical variance decomposition. Here the question is what percentage of the hierarchical variance is attributed to a, b and their interaction? Recall that brms returns those variance parameters in the \(\sigma\) metric. So before we can compare them in terms of percentages of the total variance, we have to first have to square them.

post <-
  posterior_samples(fit24.2) %>% 
  transmute(`sigma[a]^2`  = sd_a__Intercept^2,
            `sigma[b]^2`  = sd_b__Intercept^2,
            `sigma[ab]^2` = `sd_a:b__Intercept`^2) %>% 
  mutate(`sigma[total]^2` = `sigma[a]^2` + `sigma[b]^2` + `sigma[ab]^2`)

head(post)
##     sigma[a]^2  sigma[b]^2 sigma[ab]^2 sigma[total]^2
## 1 1.182039e-01 1.198310832   0.1024887     1.41900343
## 2 5.228487e-03 0.140569145   0.1430617     0.28885930
## 3 6.035757e-02 0.023973823   0.1116029     0.19593425
## 4 3.257825e-03 0.006660097   0.0818224     0.09174032
## 5 8.411625e-05 0.149055015   0.1264936     0.27563273
## 6 8.240073e-04 0.036842948   0.1836877     0.22135465

Now we just need to divide the individual variance parameters by their total and multiply by 100 to get the percent of variance for each. We’ll look at the results in a plot.

post %>% 
  pivot_longer(-`sigma[total]^2`) %>% 
  mutate(`% hierarchical variance` = 100 * value / `sigma[total]^2`) %>% 
  
  ggplot(aes(x = `% hierarchical variance`, y = name)) +
  stat_halfeye(.width = .95, fill = hr[7], color = hr[6]) +
  scale_x_continuous(limits = c(0, 100), expand = c(0, 0)) +
  scale_y_discrete(NULL, labels = ggplot2:::parse_safe) +
  coord_cartesian(ylim = c(1.5, 3.5)) +
  theme(plot.margin = margin(5.5, 10, 5.5, 5.5))

Just as each of the variance parameters was estimated with uncertainty, all that uncertainty got propagated into their transformations. Even in the midst of all this uncertainty, it’s clear that a good healthy portion of the hierarchical variance is from the interaction. Again, whatever you might think about \(a \times b\), it’s definitely not zero.

24.4 Log-linear models for contingency tables Bonus: Alternative parameterization

The Poisson distribution is widely used for count data. But notice how in our figures, we converted the results to the proportion metric. Once you’re talking about proportions, it’s not hard to further adjust your approach to thinking in terms of probabilities. So instead of thinking about the \(n\) within each cell of our contingency table, we might also think about the probability of a given condition. To approach the data this way, we could use a multilevel aggregated binomial model. McElreath covered this in Chapter 10 of his Statistical rethinking. See my (2020) translation of the text into brms code, too.

Here’s how to fit that model.

fit24.4 <-
  brm(data = my_data,
      family = binomial,
      count | trials(264) ~ 1 + (1 | a) + (1 | b) + (1 | a:b),
      prior = c(prior(normal(0, 2), class = Intercept),
                prior(normal(0, 1), class = sd)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999),
      seed = 24,
      file = "fits/fit24.04")

A few things about the syntax: The aggregated binomial model uses the logit link, just like with typical logistic regression. So when you specify family = binomial, you’re requesting the logit link. The left side of the formula argument, count | trials(264) indicates a few things. First, our criterion is count. The bar | that follows on its right indicates we’d like add additional information about the criterion. In the case of binomial regression, brms requires we specify how many trials the value in each cell of the data is referring to. When we coded trials(264), we indicated each cell was a total count of 264 trials. In case it isn’t clear, here is where the value 264 came from.

my_data %>% 
  summarise(total_trials = sum(count))
## # A tibble: 1 x 1
##   total_trials
##          <dbl>
## 1          264

Now look over the summary.

print(fit24.4)
##  Family: binomial 
##   Links: mu = logit 
## Formula: count | trials(264) ~ 1 + (1 | a) + (1 | b) + (1 | a:b) 
##    Data: my_data (Number of observations: 16) 
## Samples: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
##          total post-warmup samples = 8000
## 
## Group-Level Effects: 
## ~a (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.18      0.19     0.00     0.72 1.00     2613     3238
## 
## ~a:b (Number of levels: 16) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.34      0.12     0.12     0.61 1.00     2449     2437
## 
## ~b (Number of levels: 4) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.19     0.00     0.70 1.00     2297     3314
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -2.72      0.21    -3.09    -2.24 1.00     3640     3282
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

See that mu = logit part in the second line of the summary? Yep, that’s our link function. Since we used a different likelihood and link function from earlier models, it shouldn’t be surprising the parameters look different. But notice how the aggregated binomial model yields virtually the same results for the top portion of Figure 24.4.

nd <-
  my_data %>% 
  mutate(strip = str_c("a:", a, " b:", b, "\nN = ", count))

fitted(fit24.4,
       newdata = nd,
       summary = F) %>% 
  data.frame() %>%
  set_names(pull(nd, strip)) %>% 
  pivot_longer(everything(), values_to = "count") %>% 
  mutate(proportion = count / sum(my_data$count)) %>% 
  
  # plot!
  ggplot(aes(x = proportion, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = hr[5], color = hr[3]) +
  scale_x_continuous(breaks = c(.05, .1, .15)) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(0, .15)) +
  theme(panel.grid = element_blank()) +
  facet_wrap(~ name, scales = "free_y")

To further demonstrate the similarity of this approach to Kruschke’s multilevel Poisson approach, let’s compare the model-based cell estimates for each of the combinations of a and b, by both fit24.2 and fit24.4.

# compute the fitted summary statistics
rbind(fitted(fit24.2),
      fitted(fit24.4)) %>% 
  data.frame() %>% 
  # add an augmented version of the data
  bind_cols(expand(my_data, 
                   fit = c("fit2 (Poisson likelihood)", "fit4 (binomial likelihood)"),
                   nesting(a, b, count))) %>% 
  mutate(cell = str_c(a, "\n", b)) %>% 
  
  # plot
  ggplot(aes(x = cell)) +
  geom_hline(yintercept = c(11, 22), color = hr[3], linetype = 2) +
  geom_pointrange(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5, color = fit),
                  fatten = 1.5, position = position_dodge(width = 0.5)) +
  geom_point(aes(y = count),
             size = 2, color = hr[9]) +
  scale_color_manual(NULL, values = hr[c(7, 5)]) +
  scale_y_continuous("count", breaks = c(0, 11, 22, 33), limits = c(0, 33)) +
  theme(legend.position = "top")

The black points are the raw data. The colored point-ranges to the left and right of each data point are the posterior means and percentile-based 95% intervals for each of the cells. The results are virtually the same between the two models. Also note how both models partially pooled towards the grand mean. That’s one of the distinctive features of using the hierarchical approach.

Wrapping up, this chapter focused on how one might use the Poisson likelihood to model contingency-table data from a multilevel modeling framework. The Poisson likelihood is also handy for count data within a single-level structure, with metric predictors, and with various combinations of metric and nominal predictors. For more practice along those lines, check out Section 10.2 in my project recoding McElreath’s Statistical rethinking.

Session info

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidybayes_2.3.1 brms_2.15.0     Rcpp_1.0.6      patchwork_1.1.1
##  [5] cowplot_1.1.1   ochRe_1.0.0     janitor_2.1.0   forcats_0.5.1  
##  [9] stringr_1.4.0   dplyr_1.0.5     purrr_0.3.4     readr_1.4.0    
## [13] tidyr_1.1.3     tibble_3.1.1    ggplot2_3.3.3   tidyverse_1.3.0
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1         backports_1.2.1      plyr_1.8.6          
##   [4] igraph_1.2.6         sp_1.4-4             svUnit_1.0.3        
##   [7] splines_4.0.4        crosstalk_1.1.0.1    TH.data_1.0-10      
##  [10] rstantools_2.1.1     inline_0.3.17        digest_0.6.27       
##  [13] htmltools_0.5.1.1    rsconnect_0.8.16     fansi_0.4.2         
##  [16] magrittr_2.0.1       modelr_0.1.8         RcppParallel_5.0.2  
##  [19] matrixStats_0.57.0   xts_0.12.1           sandwich_3.0-0      
##  [22] prettyunits_1.1.1    colorspace_2.0-0     rvest_0.3.6         
##  [25] ggdist_2.4.0.9000    haven_2.3.1          xfun_0.22           
##  [28] callr_3.7.0          crayon_1.4.1         jsonlite_1.7.2      
##  [31] lme4_1.1-25          survival_3.2-10      zoo_1.8-8           
##  [34] glue_1.4.2           gtable_0.3.0         emmeans_1.5.2-1     
##  [37] V8_3.4.0             distributional_0.2.2 pkgbuild_1.2.0      
##  [40] rstan_2.21.2         abind_1.4-5          scales_1.1.1        
##  [43] mvtnorm_1.1-1        DBI_1.1.0            miniUI_0.1.1.1      
##  [46] xtable_1.8-4         HDInterval_0.2.2     stats4_4.0.4        
##  [49] StanHeaders_2.21.0-7 DT_0.16              htmlwidgets_1.5.2   
##  [52] httr_1.4.2           threejs_0.3.3        arrayhelpers_1.1-0  
##  [55] ellipsis_0.3.2       pkgconfig_2.0.3      loo_2.4.1           
##  [58] farver_2.1.0         sass_0.3.1           dbplyr_2.0.0        
##  [61] utf8_1.2.1           tidyselect_1.1.0     labeling_0.4.2      
##  [64] rlang_0.4.11         reshape2_1.4.4       later_1.2.0         
##  [67] munsell_0.5.0        cellranger_1.1.0     tools_4.0.4         
##  [70] cli_2.5.0            generics_0.1.0       broom_0.7.5         
##  [73] ggridges_0.5.3       evaluate_0.14        fastmap_1.1.0       
##  [76] processx_3.5.2       knitr_1.31           fs_1.5.0            
##  [79] nlme_3.1-152         mime_0.10            projpred_2.0.2      
##  [82] xml2_1.3.2           compiler_4.0.4       bayesplot_1.8.0     
##  [85] shinythemes_1.1.2    rstudioapi_0.13      gamm4_0.2-6         
##  [88] curl_4.3             reprex_0.3.0         statmod_1.4.35      
##  [91] bslib_0.2.4          stringi_1.5.3        highr_0.8           
##  [94] ps_1.6.0             Brobdingnag_1.2-6    lattice_0.20-41     
##  [97] Matrix_1.3-2         nloptr_1.2.2.2       markdown_1.1        
## [100] shinyjs_2.0.0        vctrs_0.3.8          pillar_1.6.0        
## [103] lifecycle_1.0.0      jquerylib_0.1.4      bridgesampling_1.0-0
## [106] estimability_1.3     raster_3.4-5         httpuv_1.6.0        
## [109] R6_2.5.0             bookdown_0.21        promises_1.2.0.1    
## [112] gridExtra_2.3        codetools_0.2-18     boot_1.3-26         
## [115] colourpicker_1.1.0   MASS_7.3-53          gtools_3.8.2        
## [118] assertthat_0.2.1     withr_2.4.2          shinystan_2.5.0     
## [121] multcomp_1.4-16      mgcv_1.8-33          parallel_4.0.4      
## [124] hms_0.5.3            grid_4.0.4           coda_0.19-4         
## [127] minqa_1.2.4          rmarkdown_2.7        snakecase_0.11.0    
## [130] shiny_1.6.0          lubridate_1.7.9.2    base64enc_0.1-3     
## [133] dygraphs_1.1.1.6