As Kruschke put it, “There are many realistic situations that involve meaningful hierarchical structure. Bayesian modeling software makes it straightforward to specify and analyze complex hierarchical models” (2015, p. 221). IMO, brms makes it even easier than JAGS. Further down, we read:
The parameters at different levels in a hierarchical model are all merely parameters that coexist in a joint parameter space. We simply apply Bayes’ rule to the joint parameter space, as we did for example when estimating two coin biases back in Figure 7.5, p. 167. To say it a little more formally with our parameters \(\theta\) and \(\omega\), Bayes’ rule applies to the joint parameter space: \(p(\theta, \omega \mid D) \propto p(D \mid \theta, \omega) p(\theta, \omega)\). What is special to hierarchical models is that the terms on the right-hand side can be factored into a chain of dependencies, like this:
The refactoring in the second line means that the data depend only on the value of \(\theta\), in the sense that when the value \(\theta\) is set then the data are independent of all other parameter values. Moreover, the value of \(\theta\) depends on the value of \(\omega\) and the value of \(\theta\) is conditionally independent of all other parameters. Any model that can be factored into a chain of dependencies like [this] is a hierarchical model. (pp. 222–223)
9.1 A single coin from a single mint
Recall from the last chapter that our likelihood is the Bernoulli distribution,
\[y_i \sim \operatorname{Bernoulli}(\theta).\]
We’ll use the beta density for our prior distribution for \(\theta\),
On page 224, Kruschke wrote: “The value of \(\kappa\) governs how near \(\theta\) is to \(\omega\), with larger values of \(\kappa\) generating values of \(\theta\) more concentrated near \(\omega\).” To give a sense of that, we’ll simulate 20 beta distributions, all with \(\omega = 0.25\) but with \(\theta\) increasing from 10 to 200, by 10. We’ll then plot them with a little help from the ggridges package(Wilke, 2021).
# Load packageslibrary(tidyverse)library(cowplot)library(ggridges)library(ggforce)library(patchwork)library(viridis)library(brms)library(bayesplot)library(tidybayes)# Define functionbeta_by_k <-function(k) { w <-0.25tibble(x =seq(from =0, to =1, length.out =1000)) |>mutate(theta =dbeta(x = x,shape1 = w * (k -2) +1,shape2 = (1- w) * (k -2) +1))}# Datatibble(k =seq(from =10, to =200, by =10)) |>mutate(theta =map(k, beta_by_k)) |>unnest(theta) |># Plotggplot(aes(x = x, y = k,height = theta,group = k, fill = k)) +geom_vline(xintercept =0.25, color ="grey85", linewidth =1/2) +geom_ridgeline(color ="grey92", scale =2, size =1/5) +scale_y_continuous(expression(kappa), breaks =1:20*10) +scale_fill_viridis_c(expression(kappa), option ="A") +xlab(expression(theta)) +theme_minimal_hgrid()
Holding \(\omega\) constant, the density gets more concentrated around \(\omega\) as \(\kappa\) increases. But back to the text: “Now we make the essential expansion of our scenario into the realm of hierarchical models. Instead of thinking of \(\omega\) as fixed by prior knowledge, we think of it as another parameter to be estimated” (p. 224). In the hierarchical model diagram of Figure 9.1, Kruschke depicted how we might treat \(\omega\) as a parameter controlled by the prior distribution, \(\operatorname{Beta}(A_\omega, B_\omega)\). Here’s our version of the diagram.
p1 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill ="grey67") +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +annotate(geom ="text",x =0.5, y =0.6,label ="italic(A[omega])*', '*italic(B[omega])", family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))p2 <-tibble(x =c(0.5, 0.475, 0.26, 0.08, 0.06, 0.5, 0.55, 0.85, 1.15, 1.2),y =c(1, 0.7, 0.6, 0.5, 0.2, 1, 0.7, 0.6, 0.5, 0.2),line =rep(letters[2:1], each =5)) |>ggplot(aes(x = x, y = y)) +geom_bspline(aes(color = line),linewidth =2/3, show.legend =FALSE) +annotate(geom ="text",x =0, y =0.125,label ="omega(italic(K)-2)+1*', '*(1-omega)(italic(K)-2)+1",family ="Times", hjust =0, parse =TRUE, size =7) +annotate(geom ="text",x =1/3, y =0.7,label ="'~'",family ="Times", parse =TRUE, size =10) +scale_x_continuous(expand =c(0, 0), limits =c(0, 2)) +scale_color_manual(values =c("grey75", "black")) +ylim(0, 1) +theme_void()p3 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill ="grey67") +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))## An annotated arrow# Save our custom arrow settingsmy_arrow <-arrow(angle =20, length =unit(0.35, "cm"), type ="closed")p4 <-tibble(x =0.5,y =1,xend =0.5,yend =0) |>ggplot(aes(x = x, xend = xend,y = y, yend = yend)) +geom_segment(arrow = my_arrow) +annotate(geom ="text",x =0.375, y =1/3,label ="'~'",family ="Times", parse =TRUE, size =10) +xlim(0, 1) +theme_void()# Bar plot of Bernoulli datap5 <-tibble(x =0:1,d = (dbinom(x, size =1, prob =0.6)) /max(dbinom(x, size =1, prob =0.6))) |>ggplot(aes(x = x, y = d)) +geom_col(fill ="grey67", width =0.4) +annotate(geom ="text",x =0.5, y =0.2,label ="Bernoulli",size =7) +annotate(geom ="text",x =0.5, y =0.94,label ="theta", family ="Times", parse =TRUE, size =7) +xlim(-0.75, 1.75) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# another annotated arrowp6 <-tibble(x =c(0.375, 0.625),y =c(1/3, 1/3),label =c("'~'", "italic(i)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Some textp7 <-tibble(x =0.5,y =0.5,label ="italic(y[i])") |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =7) +xlim(0, 1) +theme_void()# Define the layoutlayout <-c(area(t =1, b =2, l =1, r =1),area(t =4, b =5, l =1, r =1),area(t =3, b =4, l =1, r =2),area(t =6, b =6, l =1, r =1),area(t =7, b =8, l =1, r =1),area(t =9, b =9, l =1, r =1),area(t =10, b =10, l =1, r =1))# Combine, adjust, and plot!(p1 + p3 + p2 + p4 + p5 + p6 + p7) +plot_layout(design = layout) &ylim(0, 1) &theme(plot.margin =margin(0, 5.5, 0, 5.5))
Note that whereas this model includes a hierarchical prior for \(\omega\), the hyperparameter \(K\) is fixed across cases.
9.1.1 Posterior via grid approximation.
When the parameters extend over a finite domain, and there are not too many of them, then we can approximate the posterior via grid approximation. In our present situation, we have the parameters \(\theta\) and \(\omega\) that both have finite domains, namely the interval \([0, 1]\). Therefore, a grid approximation is tractable and the distributions can be readily graphed. (p. 226)
Given \(\alpha\) and \(\beta\), we can compute the corresponding mode \(\omega\). To foreshadow, consider \(\text{beta}(2, 2)\).
alpha <-2beta <-2(alpha -1) / (alpha + beta -2)
[1] 0.5
That is, the mode of \(\operatorname{Beta}(2, 2)\) is \(.5\).
We won’t be able to make the wireframe plots on the left of Figure 9.2, but we can make the others. We’ll make the initial data following Kruschke’s (p. 226) formulas.
Next we’ll define the parameter space as a tightly-spaced sequence of values ranging from 0 to 1.
parameter_space <-seq(from =0, to =1, by =0.01)
Now we’ll use parameter_space to define the ranges for the two variables, theta and omega, which we’ll save in a tibble. We’ll then sequentially feed those theta and omega values into our make_prior() while manually specifying the desired values for alpha, beta, and kappa.
# Define the grid for our grid approximationd <-crossing(theta = parameter_space,omega = parameter_space) |># Compute the joint priormutate(prior =make_prior(theta, omega, alpha =2, beta =2, kappa =100)) |># Convert the prior from the density metric to the probability mass metricmutate(prior = prior /sum(prior))head(d)
You could also make this with geom_tile(), but geom_raster() with interpolate = TRUE smooths the color transitions. Since we are going to be making a lot of plots like this in this chapter, we should consider streamlining our plotting code. In Chapter 19 of Wichkam’s (2016)ggplot2: Elegant graphics for data analysis, we learn how to make a custom geom. Here we’ll use those skills to wrap the bulk of the plot code from above into a single geom we’ll call geom_2dd(), for 2D-density plots.
d |>ggplot(aes(x = theta, y = omega, fill = prior)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
If we collapse “the joint prior across \(\theta\)” (i.e., group_by(omega) and then sum(prior)), we plot the marginal distribution for \(p(\omega)\) as seen in the top right panel.
Note how we loaded the viridis package(Garnier, 2021). That gave us access to the viridis_pal() function, which will allow us to discretize the viridis palettes and save the color names as objects. In our case, we discretized the "A" palette into nine colors and saved the fourth as a_purple. Here’s the color name.
a_purple
[1] "#822681FF"
We’ll use that color in many of the plots to follow. It’ll be something of a signature color for this chapter.
Anyway, since we are going to be making a lot of plots like this in this chapter, we’ll make another custom geom called geom_marginal().
d |>group_by(omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = omega, y = prior)) +geom_marginal(ul =0.035) +labs(x =expression(omega),y =expression(Marginal~p(omega))) +coord_flip()
We’ll follow a similar procedure to get the marginal probability distribution for theta.
d |>group_by(theta) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta, y = prior)) +geom_marginal(ul =0.035) +labs(x =expression(theta),y =expression(Marginal~p(theta)))
We’ll use the filter() function to take the two slices from the posterior grid. Since we’re taking slices, we’re no longer working with the joint probability distribution. As such, our two marginal prior distributions for theta no longer sum to 1, which means they’re no longer in a probability metric. No worries. After we group by omega, we can simply divide prior by the sum() of prior which renormalizes the two slices “so that they are individually proper probability densities that sum to \(1.0\) over \(\theta\)” (p. 226).
As Kruschke pointed out at the top of page 228, these are indeed beta densities. Here’s proof.
# We'll want this for the annotationd_text <-tibble(theta =c(0.75, 0.25),y =10,label =c("Beta(74.5, 25.5)", "Beta(25.5, 74.5)"),omega = letters[1:2])# Here's the primary data for the plottibble(theta =rep(parameter_space, times =2),alpha =rep(c(74.5, 25.5), each =101),beta =rep(c(25.5, 74.5), each =101),omega =rep(letters[1:2], each =101)) |># The plotggplot(aes(x = theta, fill = omega)) +geom_area(aes(y =dbeta(x = theta, shape1 = alpha, shape2 = beta))) +geom_text(data = d_text,aes(y = y, label = label, color = omega)) +scale_x_continuous(expression(theta), breaks =0:5/5, expand =c(0, 0), limits =c(0, 1)) +scale_y_continuous("density", expand =c(0, 0), limits =c(0, 11)) +scale_fill_viridis_d(option ="A", begin =2/9, end =6/9) +scale_color_viridis_d(option ="A", begin =2/9, end =6/9) +theme_cowplot() +panel_border() +theme(axis.line =element_blank(),legend.position ="none")
But back on track, we need the dbern() function for the lower three rows of Figure 9.2.
dbern <-function(x, z =NULL, n =NULL) { x^z * (1- x)^(n - z)}
Time to feed theta and our data into the dbern() function, which will allow us to make the 2-dimensional density plot in the middle of Figure 9.2.
# Define the datan <-12z <-9# Compute the likelihoodd <- d |>mutate(likelihood =dbern(x = theta, z = z, n = n))# Plotd |>ggplot(aes(x = theta, y = omega, fill = likelihood)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
Note how this plot demonstrates how the likelihood is solely dependent on \(\theta\); it’s orthogonal to \(\omega\). This is the visual consequence of Kruschke’s Formula 9.6,
That is, in the second line of the equation, the probability of \(y\) was only conditional on \(\theta\). But the reason we call this a hierarchical model is because the probability of \(\theta\) itself is conditioned on \(\omega\). The prior itself had a prior.
From Formula 9.1, the posterior \(p(\theta, \omega \mid D)\) is proportional to \(p(D \mid \theta) \; p(\theta \mid \omega) \; p(\omega)\). Divide that by the normalizing constant and we’ll have it in a proper probability metric. Recall that we’ve already saved the results of \(p(\theta \mid \omega) \; p(\omega)\) in the prior column. So we just need to multiply prior by likelihood and divide by their sum.
Our first depiction will be the middle panel of the second row from the bottom.
d <- d |>mutate(posterior = (likelihood * prior) /sum(likelihood * prior)) d |>ggplot(aes(x = theta, y = omega, fill = posterior)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
Although the likelihood was orthogonal to \(\omega\), conditioning the prior for \(\theta\) on \(\omega\) resulted in a posterior that was conditioned on both \(\theta\) and \(\omega\).
Making the marginal plots for posterior is much like when making them for prior, above.
# For omegad |>group_by(omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = omega, y = posterior)) +geom_marginal(ul =0.035) +labs(x =expression(omega),y =expression(Marginal~p(omega*"|"*D))) +coord_flip()
# For thetad |>group_by(theta) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta, y = posterior)) +geom_marginal(ul =0.035) +labs(x =expression(theta),y =expression(Marginal~p(theta*"|"*D))) +coord_cartesian()
Note that after we slice with filter(), the next two wrangling lines renormalize those posterior slices into probability metrics. That is, when we take a slice through the joint posterior at a particular value of \(\omega\), and renormalize by dividing the sum of discrete probability masses in that slice, we get the conditional distribution \(p(\theta \mid \omega, D)\).
In the next example depicted in Figure 9.3, we consider what happens when we combine the same data of 9 heads out of 12 trials to the same Bernoulli likelihood \(p(y \mid \theta)\), but his time with a much lower \(K\) values expressing greater uncertainty in the \(\operatorname{Beta} \big (\theta \mid \omega (6 - 2) + 1, (1 - \omega) (6 - 2) + 1 \big )\) portion of the joint prior and with a more certain hyperprior for \(\omega\), \(\operatorname{Beta}(\omega \mid 20, 20)\).
To repeat the process for Figure 9.3, we’ll first compute the new joint prior.
Here’s the initial data and the 2-dimensional density plot for the prior, the middle plot in the top row of Figure 9.3.
d |>ggplot(aes(x = theta, y = omega, fill = prior)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
That higher certainty in \(\omega\) resulted in a two-dimensional density plot where the values on the \(y\)-axis were concentrated near 0.5. This will have down-the-road consequences for the posterior. But before we get there, we’ll average over omega and theta to plot their marginal prior distributions.
# For omegad |>group_by(omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = omega, y = prior)) +geom_marginal(ul =0.052) +labs(x =expression(omega),y =expression(Marginal~p(omega))) +coord_flip()
# For thetad |>group_by(theta) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta, y = prior)) +geom_marginal(ul =0.039) +labs(x =expression(theta),y =expression(Marginal~p(theta)))
Here are the two short plots in the right panel of the second row from the top of Figure 9.3.
# Computed <- d |>mutate(likelihood =dbern(x = theta, z = z, n = n))# Plotd |>ggplot(aes(x = theta, y = omega, fill = likelihood)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
Now on to the posterior. Our first depiction will be the middle panel of the second row from the bottom of Figure 9.3. This will be \(p(\theta, \omega \mid y)\).
# Compute the posteriord <- d |>mutate(posterior = (likelihood * prior) /sum(likelihood * prior)) # Plotd |>ggplot(aes(x = theta, y = omega, fill = posterior)) +geom_2dd() +labs(x =expression(theta),y =expression(omega))
Here are the marginal plots for the two dimensions in our posterior.
# For omegad |>group_by(omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = omega, y = posterior)) +geom_marginal(ul =0.052) +labs(x =expression(omega),y =expression(Marginal~p(omega*"|"*D))) +coord_flip()
# For thetad |>group_by(theta) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta, y = posterior)) +geom_marginal(ul =0.039) +labs(x =expression(theta),y =expression(Marginal~p(theta*"|"*D)))
And we’ll finish off with the plots of Figure 9.3’s lower right panel.
In summary, Bayesian inference in a hierarchical model is merely Bayesian inference on a joint parameter space, but we look at the joint distribution (e.g., \(p(\theta, \omega)\)) in terms of its marginal on a subset of parameters (e.g., \(p(\omega)\)) and its conditional distribution for other parameters (e.g., \(p(\theta \mid \omega)\)). We do this primarily because it is meaningful in the context of particular models. (p. 230)
9.2 Multiple coins from a single mint
What if we collect data from more than one coin created by the mint? If each coin has its own distinct bias \(\theta_s\), then we are estimating a distinct parameter value for each coin, and using all the data to estimate \(\omega\). (p. 230)
Kruschke broke down a model of this form with his diagram in Figure 9.4. Here’s our version of that figure.
p1 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill = a_purple) +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +annotate(geom ="text",x =0.5, y =0.6,label ="italic(A[omega])*', '*italic(B[omega])", family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))p2 <-tibble(x =c(0.5, 0.475, 0.26, 0.08, 0.06, 0.5, 0.55, 0.85, 1.15, 1.2),y =c(1, 0.7, 0.6, 0.5, 0.2, 1, 0.7, 0.6, 0.5, 0.2),line =rep(letters[2:1], each =5)) |>ggplot(aes(x = x, y = y)) +geom_bspline(aes(color = line),linewidth =2/3, show.legend =FALSE) +annotate(geom ="text",x =0, y =0.125,label ="omega(italic(K)-2)+1*', '*(1-omega)(italic(K)-2)+1",family ="Times", hjust =0, parse =TRUE, size =7) +annotate(geom ="text",x =1/3, y =0.7,label ="'~'",family ="Times", parse =TRUE, size =10) +scale_x_continuous(expand =c(0, 0), limits =c(0, 2)) +scale_color_manual(values =c("grey75", "black")) +ylim(0, 1) +theme_void()p3 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill = a_purple) +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# An annotated arrowp4 <-tibble(x =c(0.35, 0.65),y =c(1/3, 1/3),label =c("'~'", "italic(s)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Bar plot of Bernoulli datap5 <-tibble(x =0:1,d = (dbinom(x, size =1, prob =0.6)) /max(dbinom(x, size =1, prob =0.6))) |>ggplot(aes(x = x, y = d)) +geom_col(fill = a_purple, width =0.4) +annotate(geom ="text",x =0.5, y =0.2,label ="Bernoulli",size =7) +annotate(geom ="text",x =0.5, y =0.92,label ="theta[italic(s)]", family ="Times", parse =TRUE, size =7) +xlim(-0.75, 1.75) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Another annotated arrowp6 <-tibble(x =c(0.35, 0.65),y =c(1/3, 1/3),label =c("'~'", "italic(i)*'|'*italic(s)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Some textp7 <-tibble(x =1,y =0.5,label ="italic(y)[italic(i)*'|'*italic(s)]") |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =7) +xlim(0, 2) +theme_void()# Define the layoutlayout <-c(area(t =1, b =2, l =1, r =1),area(t =4, b =5, l =1, r =1),area(t =3, b =4, l =1, r =2),area(t =6, b =6, l =1, r =1),area(t =7, b =8, l =1, r =1),area(t =9, b =9, l =1, r =1),area(t =10, b =10, l =1, r =1))# Plot!(p1 + p3 + p2 + p4 + p5 + p6 + p7) +plot_layout(design = layout) &ylim(0, 1) &theme(plot.margin =margin(0, 5.5, 0, 5.5))
The diagram accounts for multiple coins with the \(s\) index.
9.2.1 Posterior via grid approximation.
Now we have two coins,
the full prior distribution is a joint distribution over three parameters: \(\omega\), \(\theta_1\), and \(\theta_2\). In a grid approximation, the prior is specified as a three-dimensional (3D) array that holds the prior probability at various grid points in the 3D space. (p. 233)
The biases for both coins, \(\theta_1\), and \(\theta_2\), have the same prior \(\operatorname{Beta} \big(\theta_j \mid \omega (5 - 2) + 1, (1 - \omega)(5 - 2) + 1 \big)\), which, if it’s not apparent, is marked by the rather uncertain \(K = 5\). As in our first example depicted in Figure 9.2, we have a gentle hyperprior \(\operatorname{Beta}(\omega \mid 2, 2)\), which centers the posterior mode for \(\omega\) at .5. To express this in plots, we’re going to have to update our make_prior() function. It was originally designed to handle two dimensions, \(\theta\) and \(\omega\). But now we have to update it to handle our three dimensions.
Unlike what Kruschke said in the text (p. 233), we’re not using a 3D data array. Rather, we’re just using a tibble with which prior has been expanded across all possible dimensions of the three indexing variables: theta_1, theta_2, and omega. As you can see from the ‘Rows’ count, above, this makes for a very long tibble.
“Because the parameter space is 3D, a distribution on it cannot easily be displayed on a 2D page. Instead, Figure 9.5 shows various marginal distributions” (p. 234). The consequence of that is when we marginalize, we’ll have to group by the two variables we’d like to retain for the plot. For example, the plots in the left and middle columns of the top row are the same save for their indices. So let’s just do the plot for theta_1. In order to marginalize over theta_2, we’ll need to group_by(theta_1, omega) and then summarise(prior = sum(prior)).
d |>group_by(theta_1, omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_1, y = omega, fill = prior)) +geom_2dd() +labs(x =expression(theta[1]),y =expression(omega))
But we just have to average over omega and theta_1 to plot their marginal prior distributions.
# For omegad |>group_by(omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = omega, y = prior)) +geom_marginal(ul =0.041) +labs(x =expression(omega),y =expression(p(omega))) +coord_flip()
# For thetad |>group_by(theta_1) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_1, y = prior)) +geom_marginal(ul =0.041) +labs(x =expression(theta[1]),y =expression(p(theta[1])))
Before we make the plots in the middle row of Figure 9.5, we need to add the likelihoods. Recall that we’re presuming the coin flips contained in \(D_1\) and \(D_2\) are independent. Kruschke explained in Section 7.4.1, that
independence of the data across the two coins means that the data from coin 1 depend only on the bias in coin 1, and the data from coin 2 depend only on the bias in coin 2, which can be expressed formally as \(p(y_1 \mid \theta_1, \theta_2) = p(y_1 \mid \theta_1)\) and \(p(y_2 \mid \theta_1, \theta_2) = p(y_2 \mid \theta_2)\). (p. 164)
The likelihood function for our two series of coin flips is then
Now after a little group_by() followed by summarise() we can plot the two marginal likelihoods, the two plots in the middle row of Figure 9.5.
# likelihood_1d |>group_by(theta_1, omega) |>summarise(likelihood =sum(likelihood)) |>ggplot(aes(x = theta_1, y = omega, fill = likelihood)) +geom_2dd() +labs(x =expression(theta[1]),y =expression(omega))
# likelihood_2d |>group_by(theta_2, omega) |>summarise(likelihood =sum(likelihood)) |>ggplot(aes(x = theta_2, y = omega, fill = likelihood)) +geom_2dd() +labs(x =expression(theta[2]),y =expression(omega))
The likelihoods look good. Next we compute the posterior in the same way we’ve done before: multiply the prior and the likelihood and then divide by their sum in order to convert the results to a probability metric.
# Computed <- d |>mutate(posterior = (prior * likelihood) /sum(prior * likelihood))# posterior_1d |>group_by(theta_1, omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_1, y = omega, fill = posterior)) +geom_2dd() +labs(x =expression(theta[1]),y =expression(omega))
# posterior_2d |>group_by(theta_2, omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_2, y = omega, fill = posterior)) +geom_2dd() +labs(x =expression(theta[2]),y =expression(omega))
Here’s the right plot on the second row from the bottom, the posterior distribution for \(\omega\).
# For omegad |>group_by(omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = omega, y = posterior)) +geom_marginal(ul =0.041) +labs(x =expression(omega),y =expression(p(omega*"|"*D))) +coord_flip()
Now here are the marginal posterior plots on the bottom row of Figure 9.5.
# For theta_1d |>group_by(theta_1) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_1, y = posterior)) +geom_marginal(ul =0.041) +labs(x =expression(theta[1]),y =expression(p(theta[1]*"|"*D)))
# For theta_2d |>group_by(theta_2) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_2, y = posterior)) +geom_marginal(ul =0.041) +labs(x =expression(theta[2]),y =expression(p(theta[2]*"|"*D)))
We’ll do this dog and pony one more time for Figure 9.6. Keeping the data and likelihood constant, we now set \(K = 75\) and but retain both \(A_\omega = 2\) and \(B_\omega = 2\). First, we make our new data object, d.
Again, note how the only thing we changed from the last time was increasing kappa to 75. Also like last time, the plots in the left and middle columns of the top row are the same save for their indices. But unlike last time, we’ll make both in preparation for a grand plotting finale. You’ll see.
One more step: for the 2D density plots in this section, we’ll omit the coord_equal() line from our custom geom_2dd() geom. This will help us with the formatting of our final plot.
p11 <- d |>group_by(theta_1, omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_1, y = omega, fill = prior)) +geom_2dd(font_size =10) +annotate(geom ="text", x =0.05, y =0.925, label ="p(list(theta[1], omega))", color ="white", hjust =0, parse =TRUE) +labs(x =expression(theta[1]),y =expression(omega))p12 <- d |>group_by(theta_2, omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_2, y = omega, fill = prior)) +geom_2dd(font_size =10) +annotate(geom ="text", x =0.05, y =0.925, label ="p(list(theta[2], omega))", color ="white", hjust =0, parse =TRUE) +labs(x =expression(theta[2]),y =expression(omega))p11
p12
Now we’ll average over omega and theta to plot their marginal prior distributions.
# For omegap13 <- d |>group_by(omega) |>summarise(prior =sum(prior)) |>ggplot(aes(x = omega, y = prior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(omega),y =expression(p(omega))) +coord_flip()# For theta_1p21 <- d |>group_by(theta_1) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_1, y = prior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(theta[1]),y =expression(p(theta[1])))# For theta_2p22 <- d |>group_by(theta_2) |>summarise(prior =sum(prior)) |>ggplot(aes(x = theta_2, y = prior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(theta[2]),y =expression(p(theta[2])))p13
p21
p22
Let’s get those likelihoods in there and plot.
# D1: 3 heads, 12 tailsn1 <-15z1 <-3# D2: 4 heads, 1 tailn2 <-5z2 <-4# Compute the likelihoodsd <- d |>mutate(likelihood_1 =dbern(x = theta_1, z = z1, n = n1),likelihood_2 =dbern(x = theta_2, z = z2, n = n2)) |>mutate(likelihood = likelihood_1 * likelihood_2)# Plot likelihood_1p31 <- d |>group_by(theta_1, omega) |>summarise(likelihood =sum(likelihood)) |>ggplot(aes(x = theta_1, y = omega, fill = likelihood)) +geom_2dd(font_size =10) +labs(x =expression(theta[1]),y =expression(omega))# Plot likelihood_2p32 <- d |>group_by(theta_2, omega) |>summarise(likelihood =sum(likelihood)) |>ggplot(aes(x = theta_2, y = omega, fill = likelihood)) +geom_2dd(font_size =10) +labs(x =expression(theta[2]),y =expression(omega))p31
p32
Compute the posterior and make the left and middle plots of the second row to the bottom of Figure 9.6.
d <- d |>mutate(posterior = (prior * likelihood) /sum(prior * likelihood))# posterior_1p41 <- d |>group_by(theta_1, omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_1, y = omega, fill = posterior)) +geom_2dd(font_size =10) +annotate(geom ="text", x =0.05, y =0.925, label =expression(p(list(theta[1], omega)*"|"*D)), color ="white", hjust =0, parse =TRUE) +labs(x =expression(theta[1]),y =expression(omega))# posterior_2p42 <- d |>group_by(theta_2, omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_2, y = omega, fill = posterior)) +geom_2dd(font_size =10) +annotate(geom ="text", x =0.05, y =0.925, label =expression(p(list(theta[2], omega)*"|"*D)), color ="white", hjust =0, parse =TRUE) +labs(x =expression(theta[2]),y =expression(omega))p41
p42
Here’s the right plot on the same row, the posterior distribution for \(\omega\).
# For omegap43 <- d |>group_by(omega) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = omega, y = posterior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(omega),y =expression(p(omega*"|"*D))) +coord_flip()p43
Finally, here are the marginal posterior plots on the bottom row of Figure 9.6.
# For theta_1p51 <- d |>group_by(theta_1) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_1, y = posterior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(theta[1]),y =expression(p(theta[1]*"|"*D)))# For theta_2p52 <- d |>group_by(theta_2) |>summarise(posterior =sum(posterior)) |>ggplot(aes(x = theta_2, y = posterior)) +geom_marginal(ul =0.041, font_size =10) +labs(x =expression(theta[2]),y =expression(p(theta[2]*"|"*D)))p51
p52
Did you notice how we saved each of plot from this last batch as objects? For the grand finale of this subsection, we’ll be stitching all those subplots together using syntax from the patchwork package. But before we do, we need to define three more subplots: the subplots with the annotation.
The grid approximation displayed in Figures 9.5 and 9.6 used combs of only \([101]\) points on each parameter (\(\omega\), \(\theta_1\), and \(\theta_2\)). This means that the 3D grid had \([101^3 = 1{,}030{,}301]\) points, which is a size that can be handled easily on an ordinary desktop computer of the early \(21\)st century. It is interesting to remind ourselves that the grid approximation displayed in Figures 9.5 and 9.6 would have been on the edge of computability \(50\) years ago, and would have been impossible \(100\) years ago. The number of points in a grid approximation can get hefty in a hurry. If we were to expand the example by including a third coin, with its parameter \(\theta_3\), then the grid would have \([101^4 = 104{,}060{,}401]\) points, which already strains small computers. Include a fourth coin, and the grid contains over \([10\) billion\(]\) points. Grid approximation is not a viable approach to even modestly large problems, which we encounter next. (p. 235)
In case you didn’t catch it, we used different numbers of points to evaluate each parameter. Whereas Kruschke indicated in the |> he only used 50, we used 101. That value of 101 came from how we defined our parameter_space with the code seq(from = 0, to = 1, by = 0.01). The reason we used a more densely-packed parameter space was to get smoother-looking 2D density plots.
9.2.2 A realistic model with MCMC.
In this section, Kruschke freed up the previously fixed value of \(K\), now letting \(\kappa\) vary hierarchically with a gamma prior. He depicted the model in Figure 9.7. Here’s our version of the figure.
# A beta densityp1 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill = a_purple) +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +annotate(geom ="text",x =0.5, y =0.6,label ="italic(A[omega])*', '*italic(B[omega])", family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# A gamma densityp2 <-tibble(x =seq(from =0, to =5, by =0.01),d = (dgamma(x, 1.75, 0.85) /max(dgamma(x, 1.75, 0.85)))) |>ggplot(aes(x = x, y = d)) +geom_area(fill = a_purple) +annotate(geom ="text",x =2.5, y =0.2,label ="gamma",size =7) +annotate(geom ="text",x =2.5, y =0.6,label ="list(italic(S)[kappa], italic(R)[kappa])",family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))p3 <-tibble(x =c(0.5, 0.475, 0.26, 0.08, 0.06, 0.5, 0.55, 0.85, 1.15, 1.175,1.5, 1.4, 1, 0.25, 0.2, 1.5, 1.49, 1.445, 1.4, 1.39),y =c(1, 0.7, 0.6, 0.5, 0.2, 1, 0.7, 0.6, 0.5, 0.2,1, 0.7, 0.6, 0.5, 0.2, 1, 0.75, 0.6, 0.45, 0.2),line =rep(letters[2:1], each =5) |>rep(times =2),plot =rep(1:2, each =10)) |>ggplot(aes(x = x, y = y, group =interaction(plot, line))) +geom_bspline(aes(color = line),linewidth =2/3, show.legend =FALSE) +annotate(geom ="text",x =0, y =0.1,label ="omega(kappa-2)+1*', '*(1-omega)(kappa-2)+1",family ="Times", hjust =0, parse =TRUE, size =7) +annotate(geom ="text",x =c(1/3, 1.15), y =0.7,label ="'~'",family ="Times", parse =TRUE, size =10) +scale_color_manual(values =c("grey75", "black")) +scale_x_continuous(expand =c(0, 0), limits =c(0, 2)) +ylim(0, 1) +theme_void()# Another beta densityp4 <-tibble(x =seq(from =0.01, to =0.99, by =0.01),d = (dbeta(x, 2, 2)) /max(dbeta(x, 2, 2))) |>ggplot(aes(x = x, y = d)) +geom_area(fill = a_purple) +annotate(geom ="text",x =0.5, y =0.2,label ="beta",size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# An annotated arrowp5 <-tibble(x =c(0.375, 0.625),y =c(1/3, 1/3),label =c("'~'", "italic(s)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Bar plot of Bernoulli datap6 <-tibble(x =0:1,d = (dbinom(x, size =1, prob =0.6)) /max(dbinom(x, size =1, prob =0.6))) |>ggplot(aes(x = x, y = d)) +geom_col(fill = a_purple, width =0.4) +annotate(geom ="text",x =0.5, y =0.2,label ="Bernoulli",size =7) +annotate(geom ="text",x =0.5, y =0.94,label ="theta", family ="Times", parse =TRUE, size =7) +xlim(-0.75, 1.75) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Another annotated arrowp7 <-tibble(x =c(0.35, 0.65),y =c(1/3, 1/3),label =c("'~'", "italic(i)*'|'*italic(s)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Some textp8 <-tibble(x =0.5,y =0.5,label ="italic(y[i])['|'][italic(s)]") |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =7) +xlim(0, 1) +theme_void()# Define the layoutlayout <-c(area(t =1, b =2, l =1, r =1),area(t =1, b =2, l =2, r =2),area(t =4, b =5, l =1, r =1),area(t =3, b =4, l =1, r =2),area(t =6, b =6, l =1, r =1),area(t =7, b =8, l =1, r =1),area(t =9, b =9, l =1, r =1),area(t =10, b =10, l =1, r =1))# Plot!(p1 + p2 + p4 + p3 + p5 + p6 + p7 + p8) +plot_layout(design = layout) &ylim(0, 1) &theme(plot.margin =margin(0, 5.5, 0, 5.5))
“Because the value of \(\kappa − 2\) must be non-negative, the prior distribution on \(\kappa − 2\) must not allow negative values” (p. 237). Gamma is one of the distributions with that property. The gamma distribution is defined by two parameters, its shape and rate. To get a sense of how those play out, here’ a look at the gamma densities of Figure 9.8.
# How many points do you want in your sequence of x values?length <-150# Wrangletibble(shape =c(0.01, 1.56, 1, 6.25),rate =c(0.01, 0.0312, 0.02, 0.125)) |>expand_grid(x =seq(from =0, to =200, length.out = length)) |>mutate(mean = shape *1/ rate,sd =sqrt(shape * (1/ rate)^2)) |>mutate(label =str_c("shape = ", shape, ", rate = ", rate, "\nmean = ", mean, ", sd = ", round(sd, 4))) |># Plotggplot(aes(x = x, y =dgamma(x = x, shape = shape, rate = rate))) +geom_area(aes(fill = label)) +scale_x_continuous(expression(kappa), expand =expansion(mult =c(0, 0.05))) +scale_y_continuous(expression(p(kappa*"|"*s*","*r)), breaks =c(0, 0.01, 0.02),expand =expansion(mult =c(0, 0.05))) +scale_fill_viridis_d(option ="A", end =0.9, breaks =NULL) +coord_cartesian(xlim =c(0, 150)) +theme_cowplot(line_size =0) +panel_border() +facet_wrap(~ label)
You can find the formulas for the mean and \(\textit{SD}\) for a given gamma distribution here. We used those formulas in the second mutate() statement for the data-prep stage of that last figure.
Using \(s\) for shape and \(r\) for rate, Kruschke’s Equations 9.7 and 9.8 are as follows:
\[
s = \frac{\mu^2}{\sigma^2} \;\;\; \text{and} \;\;\; r = \frac{\mu}{\sigma^2} \;\;\; \text{for mean} \;\;\; \mu > 0 \\
s = 1 + \omega r \;\;\; \text{where} \;\;\; r = \frac{\omega + \sqrt{\omega^2 + 4\sigma^2}}{2\sigma^2} \;\;\; \text{for mode} \;\;\; \omega > 0.
\]
With those in hand, we can follow Kruschke’s DBDA2E-utilities.R file to make a couple convenience functions.
gamma_s_and_r_from_mean_sd <-function(mean, sd) {if (mean <=0) stop("`mean` must be > 0")if (sd <=0) stop("`sd` must be > 0") shape <- mean^2/ sd^2 rate <- mean / sd^2list(shape = shape, rate = rate)}gamma_s_and_r_from_mode_sd <-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 * ratelist(shape = shape, rate = rate)}
They’re easy to put to use:
gamma_s_and_r_from_mean_sd(mean =10, sd =100)
$shape
[1] 0.01
$rate
[1] 0.001
gamma_s_and_r_from_mode_sd(mode =10, sd =100)
$shape
[1] 1.105125
$rate
[1] 0.01051249
Here’s a more detailed look at the structure of their output.
In applied statistics, the typical way to model a Bernoulli variable is with logistic regression. Instead of going through the pain of setting up a model in brms that mirrors the one in the text, I’m going to set up a hierarchical logistic regression model, instead.
Note the family = bernoulli(link = logit) argument. In work-a-day regression with vanilla Gaussian variables, the prediction space is unbounded. But when we want to model the probability of a success for a Bernoulli variable (i.e., \(\theta\)), we need to constrain the model to only produce predictions between 0 and 1. With logistic regression, we use a link function to do just that. The consequence is that instead of modeling the probability, \(\theta\), we’re modeling the logit probability.
In case you’re curious, the logit of \(\theta\) follows the formula
But anyway, we’ll be doing logistic regression using the logit link. Kruschke covered this in detail in Chapter 21.
The next new part of our syntax is (1 | s). As in the popular frequentist lme4 package(Bates et al., 2015, 2022), you specify random effects or group-level parameters with the (|) syntax in brms. On the left side of the |, you tell brms what parameters you’d like to make random (i.e., vary by group). On the right side of the |, you tell brms what variable you want to group the parameters by. In our case, we want the intercepts to vary over the grouping variable s.
fit9.1<-brm(data = my_data,family =bernoulli(link = logit), y ~1+ (1| s),prior =c(prior(normal(0, 1.5), class = Intercept),prior(normal(0, 1), class = sd)),iter =20000, warmup =1000, thin =10, chains =4, cores =4,seed =9,file ="fits/fit09.01")
As it turns out, the \(N(0, 1.5)\) prior is flat in the probability space for the intercept in a logistic regression model. We’ll explore that a little further down. The \(N(0, 1)\) prior for the random effect is actually a half Normal. That’s because brms defaults to bound \(\textit{SD}\) parameters to zero and above. The half Normal prior for a hierarchical \(\textit{SD}\) parameter in a logistic regression model is weakly regularizing and is conservative in the sense that it presumes some pooling is preferable to no pooling. If you wanted to take a lighter approach, you might use something like a cauchy(0, 5), instead. See the prior wiki by the Stan team for more ideas on priors.
Here are the trace plots and posterior densities of the main parameters.
plot(fit9.1, widths =2:3)
The trace plots indicate no problems with convergence. We’ll need to extract the posterior draws with as_draws_df().
draws <-as_draws_df(fit9.1)
One of the nice things about bayesplot is it returns ggplot2 objects. As such, we can amend their theme settings to be consistent with our other ggplot2 plots. Here we’ll amend bayesplot::mcmc_acf() to the theme_cowplot() theme.
draws |>mutate(chain = .chain) |>mcmc_acf(pars =vars(b_Intercept, sd_s__Intercept), lags =10) +theme_cowplot()
It appears fit9.1 had very low autocorrelations. Here we’ll examine the \(N_\textit{eff}/N\) ratio.
The \(N_\textit{eff}/N\) ratio values for our model parameters were excellent. Here’s a numeric summary of the model.
print(fit9.1)
Family: bernoulli
Links: mu = logit
Formula: y ~ 1 + (1 | s)
Data: my_data (Number of observations: 280)
Draws: 4 chains, each with iter = 20000; warmup = 1000; thin = 10;
total post-warmup draws = 7600
Multilevel Hyperparameters:
~s (Number of levels: 28)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.28 0.18 0.01 0.68 1.00 7103 7350
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.25 0.14 -0.54 0.02 1.00 7885 7200
Draws were sampled 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’ll need the base plogis() function to convert the model parameters to predict \(\theta\) rather than \(\operatorname{logit}(\theta)\). After the conversions, we’ll be ready to make the histograms in the lower portion of Figure 9.10.
# Wrangledraws_small <- draws |># Convert the linear model parameters to the probability space with `plogis()`mutate(`theta[1]`= (b_Intercept +`r_s[S01,Intercept]`) |>plogis(),`theta[14]`= (b_Intercept +`r_s[S14,Intercept]`) |>plogis(),`theta[28]`= (b_Intercept +`r_s[S28,Intercept]`) |>plogis()) |># Make the difference distributionsmutate(`theta[1] - theta[14]`=`theta[1]`-`theta[14]`,`theta[1] - theta[28]`=`theta[1]`-`theta[28]`,`theta[14] - theta[28]`=`theta[14]`-`theta[28]`) |>select(starts_with("theta"))draws_small |>pivot_longer(cols =everything()) |># This line is unnecessary, but will help order the plots mutate(name =factor(name, levels =c("theta[1]", "theta[14]", "theta[28]", "theta[1] - theta[14]", "theta[1] - theta[28]", "theta[14] - theta[28]"))) |>ggplot(aes(x = value, y =0)) +stat_histinterval(point_interval = mode_hdi, .width =0.95,breaks =40, fill = a_purple, normalize ="panels") +scale_y_continuous(NULL, breaks =NULL) +xlab(NULL) +theme_minimal_hgrid() +facet_wrap(~ name, ncol =3, scales ="free")
If you wanted the specific values of the posterior modes and 95% HDIs, you could execute this.
p1 <- draws_small |>ggplot(aes(x =`theta[1]`, y =`theta[14]`)) +geom_abline(linetype =2) +geom_point(alpha =1/8, color = a_purple, size =1/8)p2 <- draws_small |>ggplot(aes(x =`theta[1]`, y =`theta[28]`)) +geom_abline(linetype =2) +geom_point(alpha =1/8, color = a_purple, size =1/8)p3 <- draws_small |>ggplot(aes(x =`theta[14]`, y =`theta[28]`)) +geom_abline(linetype =2) +geom_point(alpha =1/8, color = a_purple, size =1/8)(p1 + p2 + p3) &coord_cartesian(xlim =c(0, 1),ylim =c(0, 1)) &theme_minimal_grid()
This is posterior distribution for the population estimate for \(\theta\), which roughly corresponds to the upper right histogram of \(\omega\) in Figure 9.10.
# This part makes it easier to set the break points in `scale_x_continuous()` labels <- draws |>transmute(theta = b_Intercept |>plogis()) |>mode_hdi() |>pivot_longer(cols = theta:.upper) |>mutate(label = value |>round(digits =3) |>as.character())draws |>mutate(theta = b_Intercept |>plogis()) |>ggplot(aes(x = theta, y =0)) +stat_histinterval(point_interval = mode_hdi, .width =0.95,breaks =40, fill = a_purple) +scale_x_continuous(expression(theta), breaks = labels$value,labels = labels$label) +scale_y_continuous(NULL, breaks =NULL) +theme_minimal_hgrid()
I’m not aware there’s a straight conversion to get \(\sigma\) in a probability metric. As far as I can tell, you have to first use coef() to “extract [the] model coefficients, which are the sum of population-level effects and corresponding group-level effects” (Bürkner, 2022a, p. 58). With the model coefficient draws in hand, you can index them by posterior iteration, group them by that index, compute the iteration-level \(\textit{SD}\)’s, and then plot the distribution of the \(\textit{SD}\)’s.
# The tibble of the primary datasigmas <-coef(fit9.1, summary =FALSE)$s |>as_tibble() |>mutate(iter =row_number()) |>group_by(iter) |>pivot_longer(cols =-iter) |>mutate(theta =plogis(value)) |>summarise(sd =sd(theta))# This, again, is just to customize `scale_x_continuous()`labels <- sigmas |>mode_hdi(sd) |>pivot_longer(cols = sd:.upper) |>mutate(label = value |>round(3) |>as.character())# The plotsigmas |>ggplot(aes(x = sd, y =0)) +stat_histinterval(point_interval = mode_hdi, .width =0.95,breaks =40, fill = a_purple) +scale_x_continuous(expression(paste(sigma, " of ", theta, " in a probability metric")),breaks = labels$value,labels = labels$label) +scale_y_continuous(NULL, breaks =NULL) +theme_minimal_hgrid()
And now you have a sense of how to do all those by hand, bayesplot::mcmc_pairs() offers a fairly quick way to get a good portion of Figure 9.10.
Did you see how we slipped in the color_scheme_set() and bayesplot_theme_set() lines at the top? Usually, the plots made with bayesplot are easy to modify with ggplot2 syntax. Plots made with mcmc_pairs() function are one notable exception. On the back end, these made by combining multiple ggplot into a grid, a down-the-line result of which is they are difficult to modify. Happily, one can make some modifications beforehand by altering the global settings with the color_scheme_set() and bayesplot_theme_set() functions. You can learn more in the discussion on issue #128 on the bayesplot GitHub repo.
Kruschke used a \(\operatorname{Beta}(1, 1)\) prior for \(\omega\). If you randomly draw from that prior and plot a histogram, you’ll see it was flat.
You’ll note that plot corresponds to the upper right panel of Figure 9.11.
Recall that we used a logistic regression model with a normal(0, 1.5) prior on the intercept. If you sample from normal(0, 1.5) and then convert the draws using plogis(), you’ll discover that our normal(0, 1.5) prior was virtually flat on the probability scale. Here we’ll show the consequence of a variety of zero-mean Gaussian priors for the intercept of a logistic regression model.
# Define a functionr_norm <-function(i, n =1e4) {set.seed(1)rnorm(n = n, mean =0, sd = i) |>plogis()}# Simulate and wrangletibble(sd =seq(from =0.25, to =3, by =0.25)) |>group_by(sd) |>mutate(prior =map(sd, r_norm)) |>unnest(prior) |>ungroup() |>mutate(sd =str_c("sd = ", sd)) |># Plot!ggplot(aes(x = prior)) +geom_histogram(binwidth =0.05, boundary =0, color ="white",fill = a_purple, linewidth =0.2) +scale_x_continuous(labels =c(0, ".25", ".5", ".75", 1)) +scale_y_continuous(NULL, breaks =NULL) +coord_cartesian(xlim =c(0, 1)) +theme_minimal_hgrid() +panel_border() +facet_wrap(~ sd)
It appears that as \(\sigma\) goes lower than 1.25, the prior becomes increasingly regularizing, pulling the estimate for \(\theta\) to a neutral .5. However, as the prior’s \(\sigma\) gets larger than 1.25, more and more of the probability mass ends up at extreme values.
Next, Kruschke examined the prior distribution. There are a few ways to do this. Like in the last chapter, the one we’ll explore involves adding the sample_prior = "only" argument to the brm() function. When you do so, the results of the model are just the prior. That is, brm() leaves out the likelihood. This returns a bunch of draws from the prior predictive distribution.
If we feed fit9.1_prior into the as_draws_df() function, we’ll get back a data frame of draws from the prior, but with the same parameter names we’d get from the posterior.
Those plots clarify our hierarchical logistic regression model was a little more regularizing than Kruschke’s. The consequence of our priors was more aggressive regularization, greater shrinkage toward zero. The prose in the next section of the text clarifies this isn’t necessarily a bad thing.
Finally, here are the plots for the lower triangle in Figure 9.11.
p1 <- prior_draws |>ggplot(aes(x =`theta[1]`, y =`theta[14]`)) +geom_point(alpha =1/8, color = a_purple, size =1/8)p2 <- prior_draws |>ggplot(aes(x =`theta[1]`, y =`theta[28]`)) +geom_point(alpha =1/8, color = a_purple, size =1/8)p3 <- prior_draws |>ggplot(aes(x =`theta[14]`, y =`theta[28]`)) +geom_point(alpha =1/8, color = a_purple, size =1/8)(p1 + p2 + p3) &geom_abline(linetype =2) &coord_cartesian(xlim =c(0, 1),ylim =c(0, 1)) &theme_minimal_grid()
In case you were curious, here are the Pearson’s correlation coefficients among the priors.
“In typical hierarchical models, the estimates of low-level parameters are pulled closer together than they would be if there were not a higher-level distribution. This pulling together is called shrinkage of the estimates” (p. 245, emphasis in the original)
Further,
shrinkage is a rational implication of hierarchical model structure, and is (usually) desired by the analyst because the shrunken parameter estimates are less affected by random sampling noise than estimates derived without hierarchical structure. Intuitively, shrinkage occurs because the estimate of each low-level parameter is influenced from two sources: (1) the subset of data that are directly dependent on the low-level parameter, and (2) the higher-level parameters on which the low-level parameter depends. The higher- level parameters are affected by all the data, and therefore the estimate of a low-level parameter is affected indirectly by all the data, via their influence on the higher-level parameters. (p. 247)
This isn’t in the text, but it might help if we gave a sense of multilevel shrinkage by plotting the phenomena using the results from our model fit9.1.
my_data |>group_by(s) |>summarise(p =mean(y)) |>mutate(theta =coef(fit9.1)$s[, 1, "Intercept"] |>plogis()) |>pivot_longer(cols =-s) |># Add a little jitter to reduce the overplottingmutate(value = value +runif(n =n(), min =-0.02, max =0.02),name =if_else(name =="p", "italic(z/N)", "theta")) |>ggplot(aes(x = value, y = name, group = s)) +geom_point(color =alpha(a_purple, 1/2)) +geom_line(linewidth =1/3, alpha =1/3) +scale_x_continuous(breaks =0:5/5, expand =c(0.01, 0.01), limits =0:1) +scale_y_discrete(NULL, labels = ggplot2:::parse_safe) +labs(x ="data proportion or theta value",title ="Multilevel shrinkage in fit9.1") +theme_minimal_hgrid() +panel_border()
The dots in the \(z/N\) row are the sample statistics. The dots in the \(\theta\) row are the posterior means for each of the levels of s, the grouping variable in the my_data data. You’ll note that we jittered the values for both within the second mutate() line to help reduce the overplotting. If you don’t understand what that means, run the code without that line or set the values within runif() closer to zero. You’ll see. Anyway, for more on multilevel shrinkage and for plots of this kind, check out Efron and Morris’s classic (1977) paper, Stein’s paradox in statistics, and my blog post walking out one of their examples in brms.
9.4 Speeding up JAGS brms
Here we’ll compare the time it takes to fit fit1 as either bernoulli(link = logit) or binomial(link = logit).
See how we’re using proc.time() to record when we began and finished evaluating our brm() code? The last time we covered that was way back in Section 3.7.5. In that section, we also learned how subtracting the former from the latter yields the total elapsed time.
stop_time_bernoulli - start_time_bernoulli
user system elapsed
20.863 1.269 31.873
stop_time_binomial - start_time_binomial
user system elapsed
20.602 1.138 32.010
These times are based on my current laptop (a 2023 MacBook Pro). Your mileage may vary. If you wanted to be rigorous about this, you could do this multiple times in a mini simulation.
As to the issue of parallel processing, we’ve been doing this all along. Note our chains = 4, cores = 4 arguments.
Since Kruschke wrote his text, we have other options for speeding up your brms models related to within-chain parallelization and the backend = "cmdstanr" option. For all the details, see Weber & Bürkner’s (2022) vignette, Running brms models with within-chain parallelization.
9.5 Extending the hierarchy: Subjects within categories
Many data structures invite hierarchical descriptions that may have multiple levels. Software such as JAGS [brms] makes it easy to implement hierarchical models, and Bayesian inference makes it easy to interpret the parameter estimates, even for complex nonlinear hierarchical models. Here, we take a look at one type of extended hierarchical model. (p. 251)
As we will address below, our version of Figure 9.13 will look rather different from Kruschke’s. It’s something of a combination of the sensibilities from Figures 20.2 and 21.10. Even still, the diagram is of three-level model that shares many similarities to Kruschke’s and, as we will see, yields very similar results.
# Half-normal densityp1 <-tibble(x =seq(from =0, to =3, by =0.01)) |>ggplot(aes(x = x, y = (dnorm(x)) /max(dnorm(x)))) +geom_area(fill = a_purple) +annotate(geom ="text",x =1.5, y =0.2,label ="half-normal",size =7) +annotate(geom ="text",x =1.5, y =0.6,label ="0*','*~italic(S)[italic(s)*'|'*italic(c)]", family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Second half-normal densityp2 <-tibble(x =seq(from =0, to =3, by =0.01)) |>ggplot(aes(x = x, y = (dnorm(x)) /max(dnorm(x)))) +geom_area(fill = a_purple) +annotate(geom ="text",x =1.5, y =0.2,label ="half-normal",size =7) +annotate(geom ="text",x =1.5, y =0.6,label ="0*','*~italic(S)[italic(c)]", family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Annotated arrowp3 <-tibble(x =0.85,y =1,xend =0.5,yend =0.25) |>ggplot(aes(x = x, xend = xend,y = y, yend = yend)) +geom_segment(arrow = my_arrow) +annotate(geom ="text",x =0.54, y =0.6, label ="'~'", family ="Times", parse =TRUE, size =10) +xlim(0, 1) +theme_void()# Normal densityp4 <-tibble(x =seq(from =-3, to =3, by =0.1)) |>ggplot(aes(x = x, y = (dnorm(x)) /max(dnorm(x)))) +geom_area(fill = a_purple) +annotate(geom ="text",x =0, y =0.2,label ="normal",size =7) +annotate(geom ="text",x =c(0, 1.45), y =0.6,hjust =c(0.5, 0),label =c("italic(M)[0]", "italic(S)[0]"), family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Second normal densityp5 <-tibble(x =seq(from =-3, to =3, by =0.1)) |>ggplot(aes(x = x, y = (dnorm(x)) /max(dnorm(x)))) +geom_area(fill = a_purple,) +annotate(geom ="text",x =0, y =0.2,label ="normal",size =7) +annotate(geom ="text",x =c(0, 1.45), y =0.6,hjust =c(.5, 0),label =c("0", "sigma[italic(s)*'|'*italic(c)]"), family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Third normal densityp6 <-tibble(x =seq(from =-3, to =3, by =0.1)) |>ggplot(aes(x = x, y = (dnorm(x)) /max(dnorm(x)))) +geom_area(fill = a_purple) +annotate(geom ="text",x =0, y =0.2,label ="normal",size =7) +annotate(geom ="text",x =c(0, 1.45), y =0.6,hjust =c(0.5, 0),label =c("0", "sigma[italic(c)]"), family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Three annotated arrowsp7 <-tibble(x =c(0.09, 0.48, 0.9),y =1,xend =c(0.2, 0.425, 0.775),yend =0.2) |>ggplot(aes(x = x, xend = xend,y = y, yend = yend)) +geom_segment(arrow = my_arrow) +annotate(geom ="text",x =c(0.10, 0.42, 0.49, 0.81, 0.87), y =0.6,label =c("'~'", "'~'", "italic(s)*'|'*italic(c)", "'~'", "italic(c)"), family ="Times", parse =TRUE, size =c(10, 10, 7, 10, 7)) +xlim(0, 1) +theme_void()# Likelihood formulap8 <-tibble(x =0.5,y =0.5,label ="logistic(beta[0]+sum()[italic(s)*'|'*italic(c)]*beta['['*italic(s)*'|'*italic(c)*']']*italic(x)['['*italic(s)*'|'*italic(c)*']'](italic(i))+sum()[italic(c)]*beta['['*italic(c)*']']*italic(x)['['*italic(c)*']'](italic(i)))") |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =7) +scale_x_continuous(expand =c(0, 0), limits =c(0, 1)) +ylim(0, 1) +theme_void()# A second annotated arrowp9 <-tibble(x =c(0.375, 0.5),y =c(0.75, 0.3),label =c("'='", "mu[italic(i)*'|'*italic(sc)]")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0.4,arrow = my_arrow) +xlim(0, 1) +theme_void()# Binomial densityp10 <-tibble(x =0:7) |>ggplot(aes(x = x, y = (dbinom(x, size =7, prob =0.625)) /max(dbinom(x, size =7, prob =0.625)))) +geom_col(fill = a_purple, width =0.4) +annotate(geom ="text",x =3.5, y =0.2,label ="binomial",size =7) +annotate(geom ="text",x =7, y =0.85,label ="italic(N)[italic(i)*'|'*italic(s)]",family ="Times", parse =TRUE, size =7) +coord_cartesian(xlim =c(-1, 8),ylim =c(0, 1.2)) +theme_void() +theme(axis.line.x =element_line(linewidth =0.5))# Another annotated arrowp11 <-tibble(x =c(0.375, 0.7),y =c(1/3, 1/3),label =c("'~'", "italic(i)*'|'*italic(sc)")) |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =c(10, 7)) +geom_segment(x =0.5, xend =0.5,y =1, yend =0,arrow = my_arrow) +xlim(0, 1) +theme_void()# Some textp12 <-tibble(x =1,y =0.5,label ="italic(y)[italic(i)*'|'*italic(sc)]") |>ggplot(aes(x = x, y = y, label = label)) +geom_text(family ="Times", parse =TRUE, size =7) +xlim(0, 2) +theme_void()# Define the layoutlayout <-c(area(t =1, b =2, l =7, r =9),area(t =1, b =2, l =11, r =13),area(t =4, b =5, l =1, r =3),area(t =4, b =5, l =5, r =7),area(t =4, b =5, l =9, r =11),area(t =3, b =4, l =6, r =8),area(t =3, b =4, l =10, r =12),area(t =7, b =8, l =1, r =11),area(t =6, b =7, l =1, r =11),area(t =11, b =12, l =5, r =7),area(t =9, b =11, l =5, r =7),area(t =13, b =14, l =5, r =7),area(t =15, b =15, l =5, r =7))# Combine and plot!(p1 + p2 + p4 + p5 + p6 + p3 + p3 + p8 + p7 + p10 + p9 + p11 + p12) +plot_layout(design = layout) &ylim(0, 1) &theme(plot.margin =margin(0, 5.5, 0, 5.5))
Though we will be fitting a hierarchical model with subjects \(s\) within categories \(c\), the higher-level parameters will not be \(\omega\) and \(\kappa\). As we’ll go over, below, we will use the binomial distribution within a more conventional hierarchical logistic regression paradigm. In this paradigm, we have an overall intercept, often called \(\alpha\) or \(\beta_0\), which will be our analogue to Kruschke’s overall \(\omega\). For the two grouping categories, \(s\) and \(c\), we will have \(\sigma\) estimates, which express the variability within those grouping. You’ll see when we get there.
9.5.1 Example: Baseball batting abilities by position.
In his footnote #6, Kruschke indicated he retrieved the data from http://www.baseball-reference.com/leagues/MLB/2012-standard-batting.shtml on December of 2012. To give a sense of the data, here are the number of occasions by primary position, PriPos, with their median at bat, AtBats, values.
# A tibble: 9 × 3
PriPos n median
<chr> <int> <dbl>
1 Pitcher 324 4
2 Catcher 103 170
3 Left Field 103 164
4 1st Base 81 265
5 3rd Base 75 267
6 2nd Base 72 228.
7 Center Field 67 259
8 Shortstop 63 205
9 Right Field 60 340.
As these data are aggregated, we’ll fit with an aggregated binomial model. This is still logistic regression. The Bernoulli distribution is a special case of the binomial distribution when the number of trials in each data point is 1 (see Bürkner, 2022b for details). Since our data are aggregated, the information encoded in Hits is a combination of multiple trials, which requires us to jump up to the more general binomial likelihood. Note the Hits | trials(AtBats) syntax. With that bit, we instructed brms that our criterion, Hits, is an aggregate of multiple trials and the number of trials is encoded in AtBats.
Also note the (1 | PriPos) + (1 | PriPos:Player) syntax. In this model, we have two grouping factors, PriPos and Player. Thus we have two (|) arguments. But since players are themselves nested within positions, we have encoded that nesting with the (1 | PriPos:Player) syntax. For more on this style of syntax, see Kristoffer Magnusson’s handy blog post, Using R and lme/lmer to fit different two- and three-level longitudinal models. Since brms syntax is based on that from the earlier lme4 package, the basic syntax rules apply. Bürkner (2022a), of course, also covered these topics in the brmsformula subsection of the brms reference manual.
Happily, most have a very favorable ratio. Here’s a numeric summary of the primary model parameters.
print(fit9.2)
Family: binomial
Links: mu = logit
Formula: Hits | trials(AtBats) ~ 1 + (1 | PriPos) + (1 | PriPos:Player)
Data: my_data (Number of observations: 948)
Draws: 3 chains, each with iter = 3500; warmup = 500; thin = 1;
total post-warmup draws = 9000
Multilevel Hyperparameters:
~PriPos (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.33 0.10 0.19 0.59 1.00 2786 4735
~PriPos:Player (Number of levels: 948)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.14 0.01 0.12 0.15 1.00 4243 6049
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.16 0.12 -1.39 -0.93 1.00 783 1842
Draws were sampled 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).
As far as I’m aware, brms offers three major ways to get the group-level parameters for a hierarchical model: using one of the as_draws functions, coef(), or fitted(). We’ll cover each, beginning with as_draws. In order to look at the autocorrelation plots, above, we already saved the results from as_draws_df(fit9.2) as draws. Let’s look at its structure with head().
In the text, Kruschke described the model as having 968 parameters. Our draws data frame has one vector for each, with a couple others tacked onto the end. In the hierarchical logistic regression model, the group-specific parameters for the levels of PriPos are additive combinations of the global intercept vector, b_Intercept and each position-specific vector, r_PriPos[i.Base,Intercept], where i is a fill-in for the position of interest. And recall that since the linear model is of the logit of the criterion, we’ll need to use plogis() to convert that to the probability space.
If you take a glance at Figures 9.14 through 9.16 in the text, we’ll be making a lot of histograms of the same basic structure. To streamline our code a bit, we can make a custom histogram plotting function.
We could follow the same procedure to make the right portion of Figure 9.14. But instead, let’s switch gears and explore the second way brms affords us for plotting group-level parameters. This time, we’ll use coef().
Up in Section 9.2.4, we learned that we can use coef() to “extract [the] model coefficients, which are the sum of population-level effects and corresponding group-level effects” (Bürkner, 2022a, p. 58). The grouping level we’re interested in is PriPos, so we’ll use that to index the information returned by coef(). Since coef() returns a matrix, we’ll use as_tibble() to convert it to a tibble.
tibble [9,000 × 9] (S3: tbl_df/tbl/data.frame)
$ 1st Base.Intercept : num [1:9000] -1.07 -1.1 -1.08 -1.08 -1.07 ...
$ 2nd Base.Intercept : num [1:9000] -1.1 -1.1 -1.1 -1.1 -1.08 ...
$ 3rd Base.Intercept : num [1:9000] -1.02 -1.03 -1.03 -1.06 -1.08 ...
$ Catcher.Intercept : num [1:9000] -1.19 -1.15 -1.14 -1.14 -1.14 ...
$ Center Field.Intercept: num [1:9000] -1.04 -1.04 -1.01 -1.03 -1.07 ...
$ Left Field.Intercept : num [1:9000] -1.08 -1.08 -1.02 -1.04 -1.08 ...
$ Pitcher.Intercept : num [1:9000] -1.92 -1.84 -1.96 -1.9 -1.84 ...
$ Right Field.Intercept : num [1:9000] -1.05 -1.07 -1.04 -1.01 -1.05 ...
$ Shortstop.Intercept : num [1:9000] -1.13 -1.16 -1.09 -1.12 -1.14 ...
Keep in mind that coef() returns the values in the logit scale when used for logistic regression models. So we’ll have to use plogis() to convert the estimates to the probability metric. After we’re done converting the estimates, we’ll then make the difference distributions.
coef_small <- coef_primary_position |>select(`1st Base.Intercept`, Catcher.Intercept, Pitcher.Intercept) |>transmute(`1st Base`=`1st Base.Intercept`, Catcher = Catcher.Intercept, Pitcher = Pitcher.Intercept) |>mutate_all(.funs = plogis) |># Here we make the difference distributionsmutate(`Pitcher - Catcher`= Pitcher - Catcher,`Catcher - 1st Base`= Catcher -`1st Base`)head(coef_small)
While we’re at it, we should capitalize on the opportunity to show how these results are the same as those derived from our as_draws_df() approach, above.
For Figures 9.15 and 9.16, Kruschke drilled down further into the posterior. To drill along with him, we’ll take the opportunity to showcase fitted(), the third way brms affords us for plotting group-level parameters.
# This will make life easier. Just go with itname_list <-c("Kyle Blanks", "Bruce Chen", "ShinSoo Choo", "Ichiro Suzuki", "Mike Leake", "Wandy Rodriguez", "Andrew McCutchen", "Brett Jackson")# We'll define the data we'd like to feed into `fitted()`, herend <- my_data |>filter(Player %in% name_list) |># These last two lines aren't typically necessary, but they allow us to# arrange the rows in the same order we find the names in Figures 9.15 and 9.16mutate(Player =factor(Player, levels = name_list)) |>arrange(Player)fitted_players <-fitted( fit9.2, newdata = nd,scale ="linear",summary =FALSE) |>as_tibble() |># Rename the values as returned by `as_tibble()`set_names(name_list) |># Convert the values from the logit scale to the probability scalemutate_all(.funs = plogis) |># In this last section, we make our difference distributions mutate(`Kyle Blanks - Bruce Chen`=`Kyle Blanks`-`Bruce Chen`,`ShinSoo Choo - Ichiro Suzuki`=`ShinSoo Choo`-`Ichiro Suzuki`,`Mike Leake - Wandy Rodriguez`=`Mike Leake`-`Wandy Rodriguez`,`Andrew McCutchen - Brett Jackson`=`Andrew McCutchen`-`Brett Jackson`)glimpse(fitted_players)
Note our use of the scale = "linear" argument in the fitted() function. By default, fitted() returns predictions on the scale of the criterion. But we don’t want a list of successes and failures; we want player-level parameters. When you specify scale = "linear", you request fitted() return the values in the parameter scale.
To make our version of Figure 9.7, we’ll have to switch gears from player-specific effects to those specific to positions averaged over individual players. The fitted() approach will probably make this the easiest. To do this, we’ll need to specify values for AtBats, which will need to be a positive integer. However, since we’re asking for fitted values of the linear predictors by setting scale = "linear", any value meeting the positive-integer criterion will return the same results. We’ll keep things simple and set AtBats = 1.
Another consideration is if we’d like to use fitted() to average across one of the hierarchical grouping parameters (i.e., (1 | PriPos:Player)), we’ll need to employ the re_formula argument. With the line re_formula = ~ (1 | PriPos), we’ll instruct fitted() to return the PriPos-specific effects after averaging across levels of Player. The rest is quire similar to our method from above.
Now we make and save the nine position-specific subplots for Figure 9.17.
p1 <- fitted_positions |>pivot_longer(cols =everything(), values_to ="theta") |># Though technically not needed, this line reorders the panels to match the textmutate(name =factor(name, levels =c("1st Base", "Catcher", "Pitcher", "2nd Base", "Center Field", "Right Field", "3rd Base", "Left Field", "Shortstop"))) |>ggplot(aes(x = theta)) +geom_histogram(binwidth =0.0025, color ="white", fill = a_purple, linewidth =0.1) +stat_pointinterval(aes(y =0), point_interval = mode_hdi, .width =0.95, size =1) +scale_y_continuous(NULL, breaks =NULL) +xlab(expression(theta)) +coord_cartesian(xlim =c(0.1, 0.28)) +theme_minimal_hgrid() +panel_border() +facet_wrap(~ name, nrow =3, scales ="free_y")
In this code block, we’ll make the subplot for the overall batting average. Given the size of the model, it’s perhaps easiest to pull that information from the model with the fixef() function.
Do note that, unlike Kruschke’s Figure 9.17, our subplots are all based on \(\theta\) rather than \(\omega\). This is because of our use of a hierarchical aggregated binomial, rather than the approach Kruschke took. Even so, look how similar the results are.
Finally, we have only looked at a tiny fraction of the relations among the \(968\) parameters. We could investigate many more comparisons among parameters if we were specifically interested. In traditional statistical testing based on \(p\)-values (which will be discussed in Chapter 11), we would pay a penalty for even intending to make more comparisons. This is because a \(p\) value depends on the space of counter-factual possibilities created from the testing intentions. In a Bayesian analysis, however, decisions are based on the posterior distribution, which is based only on the data (and the prior), not on the testing intention. More discussion of multiple comparisons can be found in Section 11.4. (pp. 259–260)
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Bates, D., Maechler, M., Bolker, B., & Steven Walker. (2022). lme4: Linear mixed-effects models using Eigen’ and S4. https://CRAN.R-project.org/package=lme4
Rosa, L., Rosa, E., Sarner, L., & Barrett, S. (1998). A close look at therapeutic touch. JAMA : The Journal of the American Medical Association, 279(13), 1005–1010. https://doi.org/10.1001/jama.279.13.1005
