# 19 Metric Predicted Variable with One Nominal Predictor

This chapter considers data structures that consist of a metric predicted variable and a nominal predictor…. This type of data structure can arise from experiments or from observational studies. In experiments, the researcher assigns the categories (at random) to the experimental subjects. In observational studies, both the nominal predictor value and the metric predicted value are generated by processes outside the direct control of the researcher. In either case, the same mathematical description can be applied to the data (although causality is best inferred from experimental intervention).

The traditional treatment of this sort of data structure is called single-factor analysis of variance (ANOVA), or sometimes one-way ANOVA. Our Bayesian approach will be a hierarchical generalization of the traditional ANOVA model. The chapter will also consider the situation in which there is also a metric predictor that accompanies the primary nominal predictor. The metric predictor is sometimes called a covariate, and the traditional treatment of this data structure is called analysis of covariance (ANCOVA). The chapter also considers generalizations of the traditional models, because it is straight forward in Bayesian software to implement heavy-tailed distributions to accommodate outliers, along with hierarchical structure to accommodate heterogeneous variances in the different groups, etc. (pp. 553–554)

## 19.1 Describing multiple groups of metric data

Figure 19.1 illustrates the conventional description of grouped metric data. Each group is represented as a position on the horizontal axis. The vertical axis represents the variable to be predicted by group membership. The data are assumed to be normally distributed within groups, with equal standard deviation in all groups. The group means are deflections from overall baseline, such that the deflections sum to zero. Figure 19.1 provides a specific numerical example, with data that were randomly generated from the model. (p. 554)

We’ll want a custom data-generating function for our primary group data.

```
library(tidyverse)
generate_data <- function(seed, mean) {
set.seed(seed)
rnorm(n, mean = grand_mean + mean, sd = 2)
}
n <- 20
grand_mean <- 101
d <-
tibble(group = 1:5,
deviation = c(4, -5, -2, 6, -3)) %>%
mutate(d = purrr::map2(group, deviation, generate_data)) %>%
unnest(d) %>%
mutate(iteration = rep(1:n, times = 5))
glimpse(d)
```

```
## Observations: 100
## Variables: 4
## $ group <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2…
## $ deviation <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, -5, -5, -5, -5, -5,…
## $ d <dbl> 103.74709, 105.36729, 103.32874, 108.19056, 105.65902, 103.35906, 105.97486, 10…
## $ iteration <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3,…
```

Here we’ll make a tibble containing the necessary data for the rotated Gaussians. As far as I can tell, Kruschke’s Gaussians only span to the bounds of percentile-based 98% intervals. We partition off those bounds for each `group`

by the `ll`

and `ul`

columns in the first `mutate()`

function. In the second `mutate()`

, we expand the dataset to include a sequence of 100 values between those lower- and upper-limit points. In the third `mutate()`

, we feed those points into the `dnorm()`

function, with group-specific means and a common `sd`

.

```
densities <-
d %>%
distinct(group, deviation) %>%
mutate(ll = qnorm(.01, mean = grand_mean + deviation, sd = 2),
ul = qnorm(.99, mean = grand_mean + deviation, sd = 2)) %>%
mutate(d = map2(ll, ul, seq, length.out = 100)) %>%
mutate(density = map2(d, grand_mean + deviation, dnorm, sd = 2)) %>%
unnest(c(d, density))
head(densities)
```

```
## # A tibble: 6 x 6
## group deviation ll ul d density
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 4 100. 110. 100. 0.0133
## 2 1 4 100. 110. 100. 0.0148
## 3 1 4 100. 110. 101. 0.0165
## 4 1 4 100. 110. 101. 0.0183
## 5 1 4 100. 110. 101. 0.0203
## 6 1 4 100. 110. 101. 0.0224
```

We’ll need two more supplementary tibbles to add the flourishes to the plot. The `arrow`

tibble will specify our light-gray arrows. The `text`

tibble will contain our annotation information.

```
arrow <-
tibble(group = 1:5,
d = grand_mean,
deviation = c(4, -5, -2, 6, -3),
offset = .1)
head(arrow)
```

```
## # A tibble: 5 x 4
## group d deviation offset
## <int> <dbl> <dbl> <dbl>
## 1 1 101 4 0.1
## 2 2 101 -5 0.1
## 3 3 101 -2 0.1
## 4 4 101 6 0.1
## 5 5 101 -3 0.1
```

```
text <-
tibble(group = c(0:5, 0),
d = grand_mean,
deviation = c(0, 4, -5, -2, 6, -3, 10),
offset = rep(c(1/4, 0), times = c(6, 1)),
angle = rep(c(90, 0), times = c(6, 1)),
label = c("beta[0] == 101", "beta['[1]'] == 4", "beta['[2]'] == -5", "beta['[3]'] == -2", "beta['[4]'] == 6", "beta['[5]'] == 3", "sigma['all'] == 2"))
head(text)
```

```
## # A tibble: 6 x 6
## group d deviation offset angle label
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 0 101 0 0.25 90 beta[0] == 101
## 2 1 101 4 0.25 90 beta['[1]'] == 4
## 3 2 101 -5 0.25 90 beta['[2]'] == -5
## 4 3 101 -2 0.25 90 beta['[3]'] == -2
## 5 4 101 6 0.25 90 beta['[4]'] == 6
## 6 5 101 -3 0.25 90 beta['[5]'] == 3
```

Now we’re ready to plot.

```
library(ggridges)
d %>%
ggplot(aes(x = d, y = group, group = group)) +
geom_vline(xintercept = grand_mean, color = "white") +
geom_jitter(height = .05, alpha = 1/4, shape = 1) +
# the Gausians
geom_ridgeline(data = densities,
aes(height = -density),
min_height = NA, scale = 3/2, size = 3/4,
fill = "transparent", color = "grey50") +
# the small arrows
geom_segment(data = arrow,
aes(xend = d + deviation,
y = group + offset, yend = group + offset),
color = "grey50", size = 1,
arrow = arrow(length = unit(.2, "cm"))) +
# the large arrow on the left
geom_segment(aes(x = 80, xend = grand_mean,
y = 0, yend = 0),
color = "grey50", size = 3/4,
arrow = arrow(length = unit(.2, "cm"))) +
# the text
geom_text(data = text,
aes(x = grand_mean + deviation, y = group + offset,
label = label, angle = angle),
size = 4, parse = T) +
scale_y_continuous(NULL, breaks = 1:5,
labels = c("<1,0,0,0,0>", "<0,1,0,0,0>", "<0,0,1,0,0>", "<0,0,0,1,0>", "<0,0,0,0,1>")) +
coord_flip(ylim = c(-0.2, 5.5), xlim = 90:112) +
xlab(NULL) +
theme(panel.grid = element_blank())
```

The descriptive model presented in Figure 19.1 is the traditional one used by classical ANOVA (which is described a bit more in the next section). More general models are straight forward to implement in Bayesian software. For example, outliers could be accommodated by using heavy-tailed noise distributions (such as a \(t\) distribution) instead of a normal distribution, and different groups could be given different standard deviations. (p. 556)

## 19.2 Traditional analysis of variance

The terminology, “analysis of variance,” comes from a decomposition of overall data variance into within-group variance and between-group variance (Fisher, 1925). Algebraically, the sum of squared deviations of the scores from their overall mean equals the sum of squared deviations of the scores from their respective group means plus the sum of squared deviations of the group means from the overall mean. In other words, the total variance can be partitioned into within-group variance plus between-group variance. Because one definition of the word “analysis” is separation into constituent parts, the term ANOVA accurately describes the underlying algebra in the traditional methods. That algebraic relation is not used in the hierarchical Bayesian approach presented here. The Bayesian method can estimate component variances, however. Therefore, the Bayesian approach is not ANOVA, but is analogous to ANOVA. (p. 556)

## 19.3 Hierarchical Bayesian approach

“Our goal is to estimate its parameters in a Bayesian framework. Therefore, all the parameters need to be given a meaningfully structured prior distribution” (p. 557). However, our approach will depart a little from the one in the text. All our parameters will **not** “have generic noncommittal prior distributions” (p. 557). Most importantly, we will not follow the example in Gelman (2006) of putting a broad uniform prior on \(\sigma_y\). Rather, we will continue using the half-Gaussian prior, as recommended by the Stan team. However, we will follow Kruschke’s lead for the overall intercept and use a Gaussian prior “made broad on the scale of the data” (p. 557). And like Kruschke, we will estimate \(\sigma_\beta\) from the data.

Later on, Kruschke opined

A crucial pre-requisite for estimating \(\sigma_\beta\) from all the groups is an assumption that all the groups are representative and informative for the estimate. It only makes sense to influence the estimate of one group with data from the other groups if the groups can be meaningfully described as representative of a shared higher-level distribution. (p. 559)

Although I agree with him in spirit, this doesn’t appear to strictly be the case. As odd and paradoxical as this sounds, partial pooling can be of use even when the some of the cases are of a different kind. For more on the topic, see Efron and Morris’s classic paper and my blog post walking out one of their examples in **brms**.

### 19.3.1 Implementation in ~~JAGS~~ brms.

The **brms** setup, of course, differs a bit from JAGS.

```
fit <-
brm(data = my_data,
family = gaussian,
y ~ 1 + (1 | categirical_variable),
prior = c(prior(normal(0, x), class = Intercept),
prior(normal(0, x), class = b),
prior(cauchy(0, x), class = sd),
prior(cauchy(0, x), class = sigma)))
```

The noise standard deviation \(\sigma_y\) is depicted in the prior statement including the argument `class = sigma`

. The grand mean is depicted by the first `1`

in the model formula and its prior is indicated by the `class = Intercept`

argument. We indicate we’d like group-based deviations from the grand mean with the `(1 | categirical_variable)`

syntax, where the `1`

on the left side of the bar indicates we’d like our intercepts to vary by group and the `categirical_variable`

part simply represents the name of a given categorical variable we’d like those intercepts to vary by. The **brms** default is to do this with deviance scores, the mean for which will be zero. Although it’s not obvious in the formula syntax, the model presumes the group-based deviations are normally distributed with a mean of zero and a standard deviation, which Kruschke termed \(\sigma_\beta\). There is no prior for the mean. It’s set at zero. But there is a prior for \(\sigma_\beta\), which is denoted by the argument `class = sd`

. We, of course, are not using a uniform prior on any of our variance parameters. But in order to be weakly informative, we will use the half-Cauchy. Recall that since the **brms** default is to set the lower bound for any variance parameter to 0, there’s no need to worry about doing so ourselves. So even though the syntax only indicates `cauchy`

, it’s understood to mean Cauchy with a lower bound at zero; since the mean is usually 0, that makes is a half-Cauchy.

Kruschke set the upper bound for his \(\sigma_y\) to 10 times the standard deviation of the criterion variable. The tails of the half-Cauchy are sufficiently fat that, in practice, I’ve found it doesn’t matter much what you set the \(SD\) of its prior to. One is often a sensible default for reasonably-scaled data. But if we want to take a more principled approach, we can set it to the size of the criterion’s \(SD\) or perhaps even 10 times that.

Kruschke suggested using a gamma on \(\sigma_\beta\), which is a sensible alternative to half-Cauchy often used within the Stan universe. Especially in situations in which you would like to (a) keep the variance parameter above zero, but (b) still allow it to be arbitrarily close to zero, and also (c) let the likelihood dominate the posterior, the Stan team recommends the gamma(2, 0) prior, based on the paper by Chung and colleagues. But you should note that I don’t mean a literal 0 for the second parameter in the gamma distribution, but rather some small value like 0.1 or so. This is all clarified in Chung et al. Here’s what gamma(2, 0.1) looks like.

```
tibble(x = seq(from = 0, to = 110, by = .1)) %>%
ggplot(aes(x = x,
ymin = 0,
ymax = dgamma(x, 2, .1))) +
geom_ribbon(size = 0, fill = "grey67") +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = 0:100) +
theme(panel.grid = element_blank())
```

And if you’d like that prior be even less informative, just reduce it to like gamma(2, 0.01) or so. Kruschke goes further to recommend “the shape and rate parameters of the gamma distribution are set so its mode is `sd(y)/2`

and its standard deviation is `2*sd(y)`

, using the function `gammaShRaFromModeSD`

explained in Section 9.2.2.” (pp. 560–561). Let’s make that 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))
}
```

So in the case of standardized data where `sd(1)`

= 1, we’d use our `gamma_a_b_from_omega_sigma()`

function like so.

```
sd_y <- 1
omega <- sd_y / 2
sigma <- 2 * sd_y
(s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma))
```

```
## $shape
## [1] 1.283196
##
## $rate
## [1] 0.5663911
```

And that produces the following gamma distribution.

```
tibble(x = seq(from = 0, to = 50, by = .01)) %>%
ggplot(aes(x = x,
ymin = 0,
ymax = dgamma(x, s_r$shape, s_r$rate))) +
geom_ribbon(size = 0, fill = "grey67") +
scale_x_continuous(breaks = c(0, 1, 5, 10, 20)) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = 0:20) +
theme(panel.grid = element_blank())
```

In the parameter space that matters, from zero to one, that gamma is pretty noninformative. It peaks between the two, slopes very gently rightward, but has the nice steep slope on the left keeping the estimates off the zero boundary. And even though that right slope is very gentle given the scale of the data, it’s aggressive enough that it should keep the MCMC chains from spending a lot of time in ridiculous parts of the parameter space. I.e., when working with finite numbers of iterations, we want our MCMC chains wasting exactly zero iterations investigating what the density might be for \(\sigma_\beta \approx 1e10\) for standardized data.

### 19.3.2 Example: Sex and death.

Let’s load and `glimpse()`

at Hanley and Shapiro’s (1994) fruit fly data.

```
## Observations: 125
## Variables: 3
## $ Longevity <dbl> 35, 37, 49, 46, 63, 39, 46, 56, 63, 65, 56, 65, 70, 63, 65, 70, 77, 81, 8…
## $ CompanionNumber <chr> "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnan…
## $ Thorax <dbl> 0.64, 0.68, 0.68, 0.72, 0.72, 0.76, 0.76, 0.76, 0.76, 0.76, 0.80, 0.80, 0…
```

We can use `geom_density_ridges()`

to help get a sense of how our criterion `Longevity`

is distributed across groups of `CompanionNumber`

.

```
my_data %>%
group_by(CompanionNumber) %>%
mutate(group_mean = mean(Longevity)) %>%
ungroup() %>%
mutate(CompanionNumber = fct_reorder(CompanionNumber, group_mean)) %>%
ggplot(aes(x = Longevity, y = CompanionNumber, fill = group_mean)) +
geom_density_ridges(scale = 3/2, size = .2, color = "grey92") +
scale_fill_viridis_c(option = "A", end = .92) +
ylab(NULL) +
theme(panel.grid = element_blank(),
legend.position = "none",
axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

Let’s fire up **brms**.

We’ll want to do the preparatory work to define our `stanvars`

.

`## [1] 57.44`

`## [1] 17.56389`

```
## $shape
## [1] 1.283196
##
## $rate
## [1] 0.03224747
```

With the prep work is done, here are our `stanvars`

.

```
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta")
```

Now fit the model, our hierarchical Bayesian alternative to an ANOVA.

```
fit1 <-
brm(data = my_data,
family = gaussian,
Longevity ~ 1 + (1 | CompanionNumber),
prior = c(prior(normal(mean_y, sd_y * 5), class = Intercept),
prior(gamma(alpha, beta), class = sd),
prior(cauchy(0, sd_y), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
control = list(adapt_delta = 0.99),
stanvars = stanvars)
```

Much like Kruschke’s JAGS chains, our **brms** chains are well behaved.

Also like Kruschke, our chains appear moderately autocorrelated.

```
post <- posterior_samples(fit1, add_chain = T)
library(bayesplot)
theme_set(theme_grey() +
theme(panel.grid = element_blank()))
mcmc_acf(post, pars = c("b_Intercept", "sd_CompanionNumber__Intercept", "sigma"), lags = 10)
```

Here’s the model summary.

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Longevity ~ 1 + (1 | CompanionNumber)
## Data: my_data (Number of observations: 125)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~CompanionNumber (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 15.03 8.04 6.24 36.42 1.00 2311 3554
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 57.59 7.51 43.03 73.39 1.00 2505 3147
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 14.93 0.97 13.22 16.98 1.00 5824 5872
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

With the `ranef()`

function, we can get the summaries of the group-specific deflections.

```
## $CompanionNumber
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## None0 5.5794036 7.818598 -10.238748 20.800280
## Pregnant1 6.7890585 7.816206 -9.059062 22.125098
## Pregnant8 5.3766714 7.839145 -10.687290 20.717351
## Virgin1 -0.7754446 7.856840 -17.380637 14.139118
## Virgin8 -17.7177658 7.939150 -34.940259 -3.181382
```

And with the `coef()`

function, we can get those same group-level summaries in a non-deflection metric.

```
## $CompanionNumber
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## None0 63.17168 2.921267 57.43048 68.80815
## Pregnant1 64.38133 2.927249 58.70815 70.16666
## Pregnant8 62.96894 2.920661 57.21388 68.65320
## Virgin1 56.81683 2.934214 51.12637 62.52715
## Virgin8 39.87451 3.091971 33.85827 46.00826
```

Those are all estimates of the group-specific means. Since it wasn’t modeled, all have the same parameter estimates for \(\sigma_y\).

```
## Estimate Est.Error Q2.5 Q97.5
## 14.9346100 0.9662639 13.2174322 16.9832791
```

To prepare for our version of the top panel of Figure 19.3, we’ll use `sample_n()`

to randomly sample from the posterior draws.

```
# how many random draws from the posterior would you like?
n_draws <- 20
set.seed(19)
post_draws <-
post %>%
sample_n(size = n_draws, replace = F)
glimpse(post_draws)
```

```
## Observations: 20
## Variables: 11
## $ b_Intercept <dbl> 51.37055, 42.46025, 58.43885, 56.01002, 48.52338…
## $ sd_CompanionNumber__Intercept <dbl> 16.262790, 17.931260, 10.294159, 14.104043, 27.4…
## $ sigma <dbl> 13.71128, 13.91828, 16.24448, 14.03107, 16.30758…
## $ `r_CompanionNumber[None0,Intercept]` <dbl> 14.6642542, 18.4782479, 0.3908729, 6.5999786, 14…
## $ `r_CompanionNumber[Pregnant1,Intercept]` <dbl> 12.207309, 23.184843, 9.985389, 15.058426, 19.47…
## $ `r_CompanionNumber[Pregnant8,Intercept]` <dbl> 9.7663362, 15.5679928, 7.3901457, 7.7165676, 15.…
## $ `r_CompanionNumber[Virgin1,Intercept]` <dbl> 6.04861233, 15.33460132, -3.33245392, 0.58641261…
## $ `r_CompanionNumber[Virgin8,Intercept]` <dbl> -12.425826, -3.181759, -14.673815, -15.790020, -…
## $ lp__ <dbl> -526.6386, -528.8238, -531.0766, -528.5373, -526…
## $ chain <fct> 3, 4, 2, 3, 1, 2, 2, 4, 1, 4, 4, 4, 2, 4, 2, 2, …
## $ iter <dbl> 2677, 1910, 2483, 2557, 3745, 2803, 3677, 2255, …
```

Before we make our version of the top panel, let’s make a corresponding plot of the fixed intercept, the grand mean. The most important lines in the code, below are the ones where we used `stat_function()`

within `mapply()`

.

```
tibble(x = c(0, 150)) %>%
ggplot(aes(x = x)) +
mapply(function(mean, sd) {
stat_function(fun = dnorm,
args = list(mean = mean, sd = sd),
alpha = 2/3,
size = 1/3,
color = "grey50")
},
# enter means and standard deviations here
mean = post_draws[, "b_Intercept"],
sd = post_draws[, "sigma"]
) +
geom_jitter(data = my_data, aes(x = Longevity, y = -0.001),
height = .001,
alpha = 1/2) +
scale_x_continuous("Longevity", breaks = seq(from = 0, to = 100, by = 25)) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = 0:110) +
labs(title = "Posterior Predictive Distribution",
subtitle = "The jittered dots are the ungrouped Longevity data. The\nGaussians are posterior draws depicting the overall\ndistribution, the grand mean.")
```

Unfortunately, we can’t extend our `mapply(stat_function())`

method to the group-level estimates. To my knowledge, there isn’t a way to show the group estimates at different spots along the y-axis. And our `mapply(stat_function())`

approach has other limitations, too. Happily, we have some great alternatives. To use them, we’ll need a little help from **tidybayes**.

For the first part, we’ll take `tidybayes::add_fitted_draws()`

for a whirl.

```
densities <-
my_data %>%
distinct(CompanionNumber) %>%
add_fitted_draws(fit1, n = 20, seed = 19, dpar = c("mu", "sigma"))
glimpse(densities)
```

```
## Observations: 100
## Variables: 8
## Groups: CompanionNumber, .row [5]
## $ CompanionNumber <chr> "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnan…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2…
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .draw <int> 1115, 2133, 2685, 2745, 4031, 4436, 4483, 4803, 5677, 5784, 7173, 7557, 7…
## $ .value <dbl> 60.74737, 63.17496, 59.68698, 63.73562, 60.19609, 59.72371, 65.82900, 61.…
## $ mu <dbl> 60.74737, 63.17496, 59.68698, 63.73562, 60.19609, 59.72371, 65.82900, 61.…
## $ sigma <dbl> 15.11354, 14.27563, 16.01779, 16.30758, 14.23716, 15.40663, 16.24448, 14.…
```

With the first two lines, we made a \(5 \times 1\) tibble containing the five levels of the experimental grouping variable, `CompanionNumber`

. The `add_fitted_draws()`

function comes from **tidybayes**. The first argument of the `add_fitted_draws()`

is `newdata`

, which works much like it does in `brms::fitted()`

; it took our \(5 \times 1\) tibble. The next argument took our **brms** model fit, `fit1`

. With the `n`

argument, we indicated we just wanted 20 random draws from the posterior. The `seed`

argument makes those random draws reproducible. With `dpar`

, we requested distributional regression parameters in the output. In our case, those were the \(\mu\) and \(\sigma\) values for each level of `CompanionNumber`

. Since we took 20 draws across 5 groups, we ended up with a 100-row tibble.

The next steps are a direct extension of the method we used to make our Gaussians for our version of Figure 19.1.

```
densities <-
densities %>%
mutate(ll = qnorm(.025, mean = mu, sd = sigma),
ul = qnorm(.975, mean = mu, sd = sigma)) %>%
mutate(Longevity = map2(ll, ul, seq, length.out = 100)) %>%
unnest(Longevity) %>%
mutate(density = dnorm(Longevity, mu, sigma))
glimpse(densities)
```

```
## Observations: 10,000
## Variables: 12
## Groups: CompanionNumber, .row [5]
## $ CompanionNumber <chr> "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnant8", "Pregnan…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .draw <int> 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1…
## $ .value <dbl> 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.…
## $ mu <dbl> 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.74737, 60.…
## $ sigma <dbl> 15.11354, 15.11354, 15.11354, 15.11354, 15.11354, 15.11354, 15.11354, 15.…
## $ ll <dbl> 31.12538, 31.12538, 31.12538, 31.12538, 31.12538, 31.12538, 31.12538, 31.…
## $ ul <dbl> 90.36936, 90.36936, 90.36936, 90.36936, 90.36936, 90.36936, 90.36936, 90.…
## $ Longevity <dbl> 31.12538, 31.72381, 32.32223, 32.92066, 33.51908, 34.11750, 34.71593, 35.…
## $ density <dbl> 0.003867068, 0.004175850, 0.004502224, 0.004846502, 0.005208934, 0.005589…
```

If you look at the code we used to make `ll`

and `ul`

, you’ll see we used 95% intervals, this time. Our second `mutate()`

function is basically the same. After unnesting the tibble, we just needed to plug in the `Longevity`

, `mu`

, and `sigma`

values into the `dnorm()`

function to compute the corresponding density values.

```
densities %>%
ggplot(aes(x = Longevity, y = CompanionNumber)) +
# here we make our density lines
geom_ridgeline(aes(height = density, group = interaction(CompanionNumber, .draw)),
fill = NA, color = adjustcolor("grey50", alpha.f = 2/3),
size = 1/3, scale = 25) +
# the original data with little jitter thrown in
geom_jitter(data = my_data,
height = .04, alpha = 1/2) +
# pretty much everything below this line is aesthetic fluff
scale_x_continuous(breaks = seq(from = 0, to = 100, by = 25)) +
coord_cartesian(xlim = 0:110,
ylim = c(1.25, 5.25)) +
labs(title = "Data with Posterior Predictive Distrib.",
y = NULL) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

Do be aware that when you use this method, you may have to fiddle around with the `geom_ridgeline()`

`scale`

argument to get the Gaussian’s heights on reasonable-looking relative heights. Stick in different numbers to get a sense of what I mean. I also find that I’m often not a fan of the way the spacing on the y axis ends up with default `geom_ridgeline()`

. It’s easy to overcome this with a little `ylim`

fiddling.

To return to the more substantive interpretation, the top panel of

Figure 19.3 suggests that the normal distributions with homogeneous variances appear to be reasonable descriptions of the data. There are no dramatic outliers relative to the posterior predicted curves, and the spread of the data within each group appears to be reasonably matched by the width of the posterior normal curves. (Be careful when making visual assessments of homogeneity of variance because the visual spread of the data depends on the sample size; for a reminder see the [see the right panel of Figure 17.1, p. 478].) The range of credible group means, indicated by the peaks of the normal curves, suggests that the group Virgin8 is clearly lower than the others, and the group Virgin1 might be lower than the controls. To find out for sure, we need to examine the differences of group means, which we do in the next section. (p. 564)

For clarity, the “see the right panel of Figure 17.1, p. 478” part was changed following Kruschke’s Corrigenda.

### 19.3.3 Contrasts.

It is straight forward to examine the posterior distribution of credible differences. Every step in the MCMC chain provides a combination of group means that are jointly credible, given the data. Therefore, every step in the MCMC chain provides a credible difference between groups…

To construct the credible differences of group 1 and group 2, at every step in the MCMC chain we compute

\[\begin{align*} \mu_1 - \mu_2 & = (\beta_0 + \beta_1) - (\beta_0 + \beta_2) \\ & = (+1) \cdot \beta_1 + (-1) \cdot \beta_2 \end{align*}\]

In other words, the baseline cancels out of the calculation, and the difference is a sum of weighted group deflections. Notice that the weights sum to zero. To construct the credible differences of the average of groups 1-3 and the average of groups 4-5, at every step in the MCMC chain we compute

\[\begin{align*} (\mu_1 + \mu_2 + \mu_3) / 3 - (\mu_4 + \mu_5) / 2 & = ((\beta_0 + \beta_1) + (\beta_0 + \beta_2) + (\beta_0 + \beta_3) ) / 3 - ((\beta_0 + \beta_4) + (\beta_0 + \beta_5) ) / 2 \\ & = (\beta_1 + \beta_2 + \beta_3) / 3 - (\beta_4 + \beta_5) / 2 \\ & = (+ 1/3) \cdot \beta_1 + (+ 1/3) \cdot \beta_2 + (+ 1/3) \cdot \beta_3 + (- 1/2) \cdot \beta_4 + (- 1/2) \cdot \beta_5 \end{align*}\]

Again, the difference is a sum of weighted group deflections. The coefficients on the group deflections have the properties that they sum to zero, with the positive coefficients summing to +1 and the negative coefficients summing to −1. Such a combination is called a contrast. The differences can also be expressed in terms of effect size, by dividing the difference by \(\sigma_y\) at each step in the chain. (pp. 565–566)

To warm up, here’s how to compute the first contrast shown in the lower portion of Kruschke’s Figure 19.3–the contrast between the two pregnant conditions and the none-control condition.

```
post %>%
transmute(c = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`) %>%
ggplot(aes(x = c)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "Pregnant1.Pregnant8 vs None0",
x = "Difference")
```

In case you were curious, here are the HMC-based posterior mode and 95% HDIs.

```
post %>%
transmute(difference = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`) %>%
mode_hdi(difference)
```

```
## # A tibble: 1 x 6
## difference .lower .upper .width .point .interval
## <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 1.67 -6.83 9.07 0.95 mode hdi
```

Little difference, there. Now let’s quantify the same contrast as an effect size.

```
post %>%
transmute(es = ((`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`) / sigma) %>%
ggplot(aes(x = es)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "Pregnant1.Pregnant8 vs None0",
x = "Effect Size")
```

Tiny.

Okay, now let’s do the rest in bulk. First we’ll do the difference scores.

```
differences <-
post %>%
transmute(`Pregnant1.Pregnant8.None0 vs Virgin1` = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[None0,Intercept]`) / 3 - `r_CompanionNumber[Virgin1,Intercept]`,
`Virgin1 vs Virgin8` = `r_CompanionNumber[Virgin1,Intercept]` - `r_CompanionNumber[Virgin8,Intercept]`,
`Pregnant1.Pregnant8.None0 vs Virgin1.Virgin8` = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[None0,Intercept]`) / 3 - (`r_CompanionNumber[Virgin1,Intercept]` + `r_CompanionNumber[Virgin8,Intercept]`) / 2)
differences %>%
gather() %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Difference") +
facet_wrap(~key, scales = "free_x")
```

Because we save our data wrangling labor from above as `differences`

, it won’t take much more effort to compute and plot the corresponding effect sizes as displayed in the bottom row of Figure 19.3.

```
differences %>%
mutate_all(.funs = ~ . / post$sigma) %>%
gather() %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Effect Size") +
facet_wrap(~key, scales = "free_x")
```

In traditional ANOVA, analysts often perform a so-called omnibus test that asks whether it is plausible that all the groups are simultaneously exactly equal. I find that the omnibus test is rarely meaningful, however…. In the hierarchical Bayesian estimation used here, there is no direct equivalent to an omnibus test in ANOVA, and the emphasis is on examining all the meaningful contrasts. (p. 567)

Speaking of all meaningful contrasts, if you’d like to make all pairwise comparisons in a hierarchical model of this form, **tidybayes** offers a convenient way to do so. Here we’ll demonstrate with `geom_halfeyeh()`

.

```
fit1 %>%
# these two lines are where the magic is at
spread_draws(r_CompanionNumber[CompanionNumber,]) %>%
compare_levels(r_CompanionNumber, by = CompanionNumber) %>%
ggplot(aes(x = r_CompanionNumber, y = CompanionNumber)) +
geom_vline(xintercept = 0, color = "white") +
geom_halfeyeh(point_interval = mode_hdi, .width = c(.95, .5)) +
labs(x = "Contrast",
y = NULL) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

But back to that omnibus test notion. If you really wanted to, I suppose one rough analogue would be to use information criteria to compare the hierarchical model to one that includes a single intercept with no group-level deflections. Here’s what the simpler model would look like.

```
fit1_without_deflections <-
brm(data = my_data,
family = gaussian,
Longevity ~ 1,
prior = c(prior(normal(mean_y, sd_y * 5), class = Intercept),
prior(cauchy(0, sd_y), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
stanvars = stanvars)
```

Here’s the model summary.

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Longevity ~ 1
## Data: my_data (Number of observations: 125)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 57.42 1.59 54.32 60.56 1.00 9746 7638
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 17.69 1.16 15.56 20.12 1.00 7742 7149
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Here are their LOO values and their difference score.

```
fit1 <- add_criterion(fit1, "loo")
fit1_without_deflections <- add_criterion(fit1_without_deflections, "loo")
loo_compare(fit1, fit1_without_deflections) %>%
print(simplify = F)
```

```
## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
## fit1 0.0 0.0 -517.7 6.9 5.5 0.6 1035.4 13.8
## fit1_without_deflections -19.4 5.3 -537.1 7.1 1.8 0.3 1074.1 14.1
```

The hierarchical model has a better LOO. Here are the stacking-based model weights.

```
## fit1 fit1_without_deflections
## 9.999994e-01 6.437616e-07
```

If you don’t like scientific notation, just `round()`

.

```
## fit1 fit1_without_deflections
## 1 0
```

Yep, in complimenting the LOO difference, virtually all the stacking weight went to the hierarchical model. You might think of this another way. The conceptual question we’re asking is *Does it make sense to say that the* \(\sigma_\beta\) *parameter is zero? Is zero a credible value?* We’ll, I suppose we could just look at the posterior to assess for that.

```
post %>%
ggplot(aes(x = sd_CompanionNumber__Intercept)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, binwidth = 1) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = 0:50) +
labs(title = expression(paste("Behold the fit1 posterior for the ", sigma[beta], " parameter.")),
subtitle = "This parameter's many things, but zero isn't one of them.",
x = NULL)
```

Yeah, zero and other values close to zero don’t look credible for that parameter. 95% of the mass is between 5 and 30, with the bulk hovering around 10. We don’t need an \(F\)-test or even a LOO model comparison to see the writing on wall.

### 19.3.4 Multiple comparisons and shrinkage.

The previous section suggested that an analyst should investigate all contrasts of interest. This recommendation can be thought to conflict with traditional advice in the context on null hypothesis significance testing, which instead recommends that a minimal number of comparisons should be conducted in order to maximize the power of each test while keeping the overall false alarm rate capped at 5% (or whatever maximum is desired)…. Instead, a Bayesian analysis can mitigate false alarms by incorporating prior knowledge into the model. In particular, hierarchical structure (which is an expression of prior knowledge) produces shrinkage of estimates, and shrinkage can help rein in estimates of spurious outlying data. For example, in the posterior distribution from the fruit fly data, the modal values of the posterior group means have a range of 23.2. The sample means of the groups have a range of 26.1. Thus, there is some shrinkage in the estimated means. The amount of shrinkage is dictated only by the data and by the prior structure, not by the intended tests. (p. 568)

We may as well compute those ranges by hand. Here’s the range of the observed data.

```
my_data %>%
group_by(CompanionNumber) %>%
summarise(mean = mean(Longevity)) %>%
summarise(range = max(mean) - min(mean))
```

```
## # A tibble: 1 x 1
## range
## <dbl>
## 1 26.1
```

For our hierarchical model `fit1`

, the posterior means are rank ordered in the same way as the empirical data.

```
coef(fit1)$CompanionNumber[, , "Intercept"] %>%
data.frame() %>%
rownames_to_column(var = "companion_number") %>%
arrange(Estimate) %>%
mutate_if(is.double, round, digits = 1)
```

```
## companion_number Estimate Est.Error Q2.5 Q97.5
## 1 Virgin8 39.9 3.1 33.9 46.0
## 2 Virgin1 56.8 2.9 51.1 62.5
## 3 Pregnant8 63.0 2.9 57.2 68.7
## 4 None0 63.2 2.9 57.4 68.8
## 5 Pregnant1 64.4 2.9 58.7 70.2
```

If we compute the range by a difference of the point estimates of the highest and lowest posterior means, we can get a quick number.

```
coef(fit1)$CompanionNumber[, , "Intercept"] %>%
as_tibble() %>%
summarise(range = max(Estimate) - min(Estimate))
```

```
## # A tibble: 1 x 1
## range
## <dbl>
## 1 24.5
```

Note that wasn’t fully Bayesian of us. Those means and their difference carry uncertainty with them and that uncertainty can be fully expressed if we use all the posterior draws (i.e., use `summary = F`

and wrangle).

```
coef(fit1, summary = F)$CompanionNumber[, , "Intercept"] %>%
as_tibble() %>%
transmute(range = Pregnant1 - Virgin8) %>%
mode_hdi(range)
```

```
## # A tibble: 1 x 6
## range .lower .upper .width .point .interval
## <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 24.7 15.7 32.9 0.95 mode hdi
```

Happily, the central tendency of the range is near equivalent with both methods, but now we have 95% intervals, too. Do note how wide they are. This is why we work with the full set of posterior draws.

### 19.3.5 The two-group case.

A special case of our current scenario is when there are only two groups. The model of the present section could, in principle, be applied to the two-group case, but the hierarchical structure would do little good because there is virtually no shrinkage when there are so few groups (and the top-level prior on \(\sigma_\beta\) is broad as assumed here). (p. 568)

For kicks and giggles, let’s practice. Since `Pregnant1`

and `Virgin8`

had the highest and lowest empirical means—making them the groups best suited to define our range, we’ll use them to fit the 2-group hierarchical model. To fit it with haste, just use `update()`

.

```
fit2 <-
update(fit1,
newdata = my_data %>%
filter(CompanionNumber %in% c("Pregnant1", "Virgin8")),
control = list(adapt_delta = 0.999,
max_treedepth = 12),
seed = 19)
```

Even with just two groups, there were no gross issues with fitting the model.

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Longevity ~ 1 + (1 | CompanionNumber)
## Data: my_data %>% filter(CompanionNumber %in% c("Pregnan (Number of observations: 50)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~CompanionNumber (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 33.07 22.62 8.86 93.59 1.00 3269 4606
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 52.23 24.82 1.13 105.83 1.00 3447 3274
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 14.24 1.49 11.66 17.51 1.00 5500 4602
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

If you compare the posteriors for \(\sigma_\beta\) across the two models, you’ll see how the one for `fit2`

is substantially larger.

```
## Estimate Est.Error Q2.5 Q97.5
## 15.031865 8.041133 6.238274 36.416484
```

```
## Estimate Est.Error Q2.5 Q97.5
## 33.065688 22.615175 8.861131 93.587013
```

This implies less shrinkage and a larger range.

```
coef(fit2, summary = F)$CompanionNumber[, , "Intercept"] %>%
as_tibble() %>%
transmute(range = Pregnant1 - Virgin8) %>%
mode_hdi(range)
```

```
## # A tibble: 1 x 6
## range .lower .upper .width .point .interval
## <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 25.9 17.6 33.4 0.95 mode hdi
```

And indeed, the range between the two groups is larger. Now the posterior mode for their difference has almost converged to that of the raw data. Kruschke then went on to recommend using a single-level model in such situations, instead.

That is why the two-group model in Section 16.3 did not use hierarchical structure, as illustrated in Figure 16.11 (p. 468). That model also used a \(t\) distribution to accommodate outliers in the data, and that model allowed for heterogeneous variances across groups. Thus, for two groups, it is more appropriate to use the model of Section 16.3. The hierarchical multi-group model is generalized to accommodate outliers and heterogeneous variances in Section 19.5. (p. 568)

As a refresher, here’s what the **brms** code for that Chapter 16 model looked like.

```
fit3 <-
brm(data = my_data,
family = student,
bf(Score ~ 0 + Group, sigma ~ 0 + Group),
prior = c(prior(normal(mean_y, sd_y * 10), class = b),
prior(normal(0, log(sd_y)), class = b, dpar = sigma),
prior(exponential(one_over_twentynine), class = nu)),
chains = 4, cores = 4,
stanvars = stanvars)
```

Let’s adjust it for our data. Since we have a reduced data set, we’ll need to re-compute our `stanvars`

values, which were based on the raw data.

```
# it's easier to just make a reduced data set
my_small_data <-
my_data %>%
filter(CompanionNumber %in% c("Pregnant1", "Virgin8"))
(mean_y <- mean(my_small_data$Longevity))
```

`## [1] 51.76`

`## [1] 19.11145`

```
## $shape
## [1] 1.283196
##
## $rate
## [1] 0.02963623
```

Here we update `stanvars`

.

```
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta") +
stanvar(1/29, name = "one_over_twentynine")
```

Note that our priors, here, are something of a blend of those from Chapter 16 and those from our hierarchical model, `fit1`

.

```
fit3 <-
brm(data = my_small_data,
family = student,
bf(Longevity ~ 0 + CompanionNumber, sigma ~ 0 + CompanionNumber),
prior = c(prior(normal(mean_y, sd_y * 10), class = b),
prior(normal(0, log(sd_y)), class = b, dpar = sigma),
prior(exponential(one_over_twentynine), class = nu)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
stanvars = stanvars)
```

Here’s the model summary.

```
## Family: student
## Links: mu = identity; sigma = log; nu = identity
## Formula: Longevity ~ 0 + CompanionNumber
## sigma ~ 0 + CompanionNumber
## Data: my_small_data (Number of observations: 50)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## CompanionNumberPregnant1 64.65 3.26 58.29 71.19 1.00 12438 8155
## CompanionNumberVirgin8 38.80 2.53 33.88 43.78 1.00 13341 8746
## sigma_CompanionNumberPregnant1 2.73 0.15 2.45 3.05 1.00 13345 9114
## sigma_CompanionNumberVirgin8 2.48 0.16 2.18 2.80 1.00 13031 8729
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu 39.61 31.02 6.01 120.48 1.00 12182 8961
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Man, look at those effective samples! As it turns out, they can be greater than the number of post-warmup samples. And here’s the range in posterior means.

```
fixef(fit3, summary = F) %>%
as_tibble() %>%
transmute(range = CompanionNumberPregnant1 - CompanionNumberVirgin8) %>%
mode_hdi(range)
```

```
## # A tibble: 1 x 6
## range .lower .upper .width .point .interval
## <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 25.7 18.1 34.4 0.95 mode hdi
```

The results are pretty much the same as that of the two-group hierarchical model, maybe a touch larger. Yep, Kruschke was right. Hierarchical models with two groups and permissive priors on \(\sigma_\beta\) don’t shrink the estimates to the grand mean all that much.

## 19.4 Including a metric predictor

“In Figure 19.3, the data within each group have a large standard deviation. For example, longevities in the Virgin8 group range from 20 to 60 days” (p. 568). Turns out Kruschke’s slightly wrong on this. Probably just a typo.

```
my_data %>%
group_by(CompanionNumber) %>%
summarise(min = min(Longevity),
max = max(Longevity),
range = max(Longevity) - min(Longevity))
```

```
## # A tibble: 5 x 4
## CompanionNumber min max range
## <chr> <dbl> <dbl> <dbl>
## 1 None0 37 96 59
## 2 Pregnant1 42 97 55
## 3 Pregnant8 35 86 51
## 4 Virgin1 21 81 60
## 5 Virgin8 16 60 44
```

But you get the point. For each group, there was quite a range. We might add predictors to the model to help account for those ranges.

The additional metric predictor is sometimes called a covariate. In the experimental setting, the focus of interest is usually on the nominal predictor (i.e., the experimental treatments), and the covariate is typically thought of as an ancillary predictor to help isolate the effect of the nominal predictor. But mathematically the nominal and metric predictors have equal status in the model. Let’s denote the value of the metric covariate for subject \(i\) as \(x_\text{cov}(i)\). Then the expected value of the predicted variable for subject \(i\) is

\[\mu (i) = \beta_0 + \sum_j \beta_{[j]} x_{[j]} (i) + \beta_\text{cov} x_\text{cov}(i)\]

with the usual sum-to-zero constraint on the deflections of the nominal predictor stated in Equation 19.2. In words, Equation 19.5 says that the predicted value for subject \(i\) is a baseline plus a deflection due to the group of \(i\) plus a shift due to the value of \(i\) on the covariate. (p. 569)

And the \(j\) subscript, recall, denotes group membership. In this context, it often

makes sense to set the intercept as the mean of predicted values if the covariate is re-centered at its mean value, which is denoted \(\overline x_\text{cov}\). Therefore Equation 19.5 is algebraically reformulated to make the baseline respect those constraints…. The first equation below is simply Equation 19.5 with \(x_\text{cov}\) recentered on its mean, \(\overline x_\text{cov}\). The second line below merely algebraically rearranges the terms so that the nominal deflections sum to zero and the constants are combined into the overall baseline:

\[\begin{align*} \mu & = \alpha_0 + \sum_j \alpha_{[j]} x_{[j]} + \alpha_\text{cov} (x_\text{cov} - \overline{x}_\text{cov}) \\ & = \underbrace{\alpha_0 + \overline{\alpha} - \alpha_\text{cov} \overline{x}_\text{cov}}_{\beta_0} + \sum_j \underbrace{(\alpha_{[j]} - \overline{\alpha})}_{\beta_[j]} x_{[j]} + \underbrace{\alpha_\text{cov}}_{\beta_{\text{cov}}} x_\text{cov} \\ & \text{where } \overline{\alpha} = \frac{1}{J} \sum^J_{j = 1} \alpha_{[j]} \end{align*}\] (pp. 569–570)

### 19.4.1 Example: Sex, death, and size.

Kruschke recalled `fit1`

’s estimate for \(\sigma_y\) had a posterior mode around 14.8. Let’s confirm with a plot.

```
posterior_samples(fit1) %>%
ggplot(aes(x = sigma, y = 0)) +
geom_halfeyeh(point_range = mode_hdi, .width = c(.5, .95)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(expression(sigma[y])) +
theme(panel.grid = element_blank())
```

Yep, that looks about right. That large of a difference in days would indeed make it difficult to detect between-group differences if those differences were typically on the scale of just a few days. Since `Thorax`

is moderately correlated with `Longevity`

, including `Thorax`

in the statistical model should help shrink that \(\sigma_y\) estimate, making it easier to compare group means. Following the sensibilities from the equations just above, here we’ll mean-center our covariate, first.

```
## # A tibble: 6 x 4
## Longevity CompanionNumber Thorax thorax_c
## <dbl> <chr> <dbl> <dbl>
## 1 35 Pregnant8 0.64 -0.181
## 2 37 Pregnant8 0.68 -0.141
## 3 49 Pregnant8 0.68 -0.141
## 4 46 Pregnant8 0.72 -0.101
## 5 63 Pregnant8 0.72 -0.101
## 6 39 Pregnant8 0.76 -0.0610
```

Our model code follows the structure of that in Kruschke’s `Jags-Ymet-Xnom1met1-MnormalHom-Example.R`

and `Jags-Ymet-Xnom1met1-MnormalHom.R`

files. As a preparatory step, we redefine the values necessary for `stanvars`

.

`## [1] 57.44`

`## [1] 17.56389`

`## [1] 0.07745367`

```
## $shape
## [1] 1.283196
##
## $rate
## [1] 0.03224747
```

```
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(sd_thorax_c, name = "sd_thorax_c") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta")
```

Now we’re ready to fit the `brm()`

model, our hierarchical alternative to ANCOVA.

```
fit4 <-
brm(data = my_data,
family = gaussian,
Longevity ~ 1 + thorax_c + (1 | CompanionNumber),
prior = c(prior(normal(mean_y, sd_y * 5), class = Intercept),
prior(normal(0, 2 * sd_y / sd_thorax_c), class = b),
prior(gamma(alpha, beta), class = sd),
prior(cauchy(0, sd_y), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
control = list(adapt_delta = 0.99),
stanvars = stanvars)
```

Here’s the model summary.

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Longevity ~ 1 + thorax_c + (1 | CompanionNumber)
## Data: my_data (Number of observations: 125)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~CompanionNumber (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 14.07 7.53 5.92 33.51 1.00 2575 4157
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 57.35 6.82 44.01 70.95 1.00 2429 3081
## thorax_c 136.07 12.48 111.30 160.49 1.00 7984 7348
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 10.61 0.69 9.34 12.05 1.00 7505 7602
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Let’s see if that \(\sigma_y\) posterior shrank.

```
posterior_samples(fit4) %>%
ggplot(aes(x = sigma, y = 0)) +
geom_halfeyeh(point_range = mode_hdi, .width = c(.5, .95)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(expression(sigma[y]))
```

Yep, sure did! Now our between-group comparisons should be more precise. Heck, if we wanted to we could even make a difference plot.

```
dif <- posterior_samples(fit1) %>% select(sigma) - posterior_samples(fit4) %>% select(sigma)
dif %>%
ggplot(aes(x = sigma, y = 0)) +
geom_halfeyeh(point_range = mode_hdi, .width = c(.5, .95)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = "This is a difference distribution",
x = expression(sigma[y]))
```

If you want a quick and dirty plot of the relation between `thorax_c`

and `Longevity`

, you can use `brms::conditional_effects()`

.

But to make plots like the ones at the top of Figure 19.5, we’ll have to work a little harder. First, we need some intermediary values marking off the three values along the `Thorax`

-axis Kruschke singled out in his top panel plots. As far as I can tell, they were the `min()`

, the `max()`

, and their `mean()`

.

`## [1] 0.64 0.94`

`## [1] 0.79`

Next, we’ll make the data necessary for our side-tipped Gaussians. For kicks and giggles, we’ll choose 80 draws instead of 20. But do note how we used our `r`

values, from above, to specify both `Thorax`

and `thorax_c`

values in addition to the `CompanionNumber`

categories for the `newdata`

argument. Otherwise, this workflow is very much the same as in previous plots.

```
n_draws <- 80
densities <-
my_data %>%
distinct(CompanionNumber) %>%
expand(CompanionNumber, Thorax = c(r[1], mean(r), r[2])) %>%
mutate(thorax_c = Thorax - mean(my_data$Thorax)) %>%
add_fitted_draws(fit4, n = n_draws, seed = 19, dpar = c("mu", "sigma")) %>%
mutate(ll = qnorm(.025, mean = mu, sd = sigma),
ul = qnorm(.975, mean = mu, sd = sigma)) %>%
mutate(Longevity = map2(ll, ul, seq, length.out = 100)) %>%
unnest(Longevity) %>%
mutate(density = dnorm(Longevity, mu, sigma))
glimpse(densities)
```

```
## Observations: 120,000
## Variables: 14
## Groups: CompanionNumber, Thorax, thorax_c, .row [15]
## $ CompanionNumber <chr> "None0", "None0", "None0", "None0", "None0", "None0", "None0", "None0", "…
## $ Thorax <dbl> 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0…
## $ thorax_c <dbl> -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .draw <int> 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 3…
## $ .value <dbl> 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.…
## $ mu <dbl> 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.56409, 38.…
## $ sigma <dbl> 12.08948, 12.08948, 12.08948, 12.08948, 12.08948, 12.08948, 12.08948, 12.…
## $ ll <dbl> 14.86915, 14.86915, 14.86915, 14.86915, 14.86915, 14.86915, 14.86915, 14.…
## $ ul <dbl> 62.25903, 62.25903, 62.25903, 62.25903, 62.25903, 62.25903, 62.25903, 62.…
## $ Longevity <dbl> 14.86915, 15.34783, 15.82652, 16.30520, 16.78389, 17.26258, 17.74126, 18.…
## $ density <dbl> 0.004834375, 0.005220395, 0.005628408, 0.006058804, 0.006511894, 0.006987…
```

Here, we’ll use a simplified workflow to extract the `fitted()`

values in order to make the regression lines. Since these are straight lines, all we need are two values for each draw, one at the extremes of the `Thorax`

axis.

```
f <-
my_data %>%
distinct(CompanionNumber) %>%
expand(CompanionNumber, Thorax = c(r[1], mean(r), r[2])) %>%
mutate(thorax_c = Thorax - mean(my_data$Thorax)) %>%
add_fitted_draws(fit4, n = n_draws, seed = 19, value = "Longevity")
glimpse(f)
```

```
## Observations: 1,200
## Variables: 8
## Groups: CompanionNumber, Thorax, thorax_c, .row [15]
## $ CompanionNumber <chr> "None0", "None0", "None0", "None0", "None0", "None0", "None0", "None0", "…
## $ Thorax <dbl> 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0.64, 0…
## $ thorax_c <dbl> -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.18096, -0.…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ .draw <int> 31, 208, 586, 800, 1099, 1115, 1353, 1800, 2120, 2133, 2195, 2220, 2492, …
## $ Longevity <dbl> 38.56409, 30.73483, 41.54440, 39.10625, 36.16563, 31.61745, 39.23632, 36.…
```

Now we’re ready to make our plots for the top row of Figure 19.3.

```
densities %>%
ggplot(aes(x = Longevity, y = Thorax)) +
# the Gaussians
geom_ridgeline(aes(height = -density, group = interaction(Thorax, .draw)),
fill = NA, size = 1/5, scale = 5/3,
color = adjustcolor("grey50", alpha.f = 1/5),
min_height = NA) +
# the vertical lines below the Gaussians
geom_line(aes(group = interaction(Thorax, .draw)),
color = "grey50", alpha = 1/5, size = 1/5) +
# the regression lines
geom_line(data = f,
aes(group = .draw),
alpha = 1/5, size = 1/5, color = "grey50") +
# the data
geom_point(data = my_data,
alpha = 1/2) +
coord_flip(xlim = 0:110,
ylim = c(.58, 1)) +
facet_wrap(~CompanionNumber, ncol = 5)
```

Now we have a covariate in the model, we have to decide on which of its values we want to base our group comparisons. Unless there’s a substantive reason for another value, the mean is a good standard choice. And since the covariate `thorax_c`

is already mean centered, that means we can effectively leave it out of the equation. Here they are in the simple difference metric.

```
post <- posterior_samples(fit4)
differences <-
post %>%
transmute(`Pregnant1.Pregnant8 vs None0` = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`,
`Pregnant1.Pregnant8.None0 vs Virgin1` = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[None0,Intercept]`) / 3 - `r_CompanionNumber[Virgin1,Intercept]`,
`Virgin1 vs Virgin8` = `r_CompanionNumber[Virgin1,Intercept]` - `r_CompanionNumber[Virgin8,Intercept]`,
`Pregnant1.Pregnant8.None0 vs Virgin1.Virgin8` = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[None0,Intercept]`) / 3 - (`r_CompanionNumber[Virgin1,Intercept]` + `r_CompanionNumber[Virgin8,Intercept]`) / 2)
differences %>%
gather() %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Difference") +
theme(strip.text = element_text(size = 6.4)) +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

Now we’ll look at the differences in the effect size metric. Since we saved our work above, it’s really easy to just convert the differences in bulk with `mutate_all()`

.

```
differences %>%
mutate_all(.funs = ~. / post$sigma) %>%
gather() %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Effect Size") +
theme(strip.text = element_text(size = 6.4)) +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

“The HDI widths of all the contrasts have gotten smaller by virtue of including the covariate in the analysis” (p. 571).

### 19.4.2 Analogous to traditional ANCOVA.

In contrast with ANCOVA,

Bayesian methods do not partition the least-squares variance to make estimates, and therefore the Bayesian method is analogous to ANCOVA but is not ANCOVA. Frequentist practitioners are urged to test (with \(p\) values) whether the assumptions of (a) equal slope in all groups, (b) equal standard deviation in all groups, and (c) normally distributed noise can be rejected. In a Bayesian approach, the descriptive model is generalized to address these concerns, as will be discussed in Section 19.5. (p. 572)

### 19.4.3 Relation to hierarchical linear regression.

Here Kruschke contrasts our last model with the one from way back in Chapter 17, section 3. As a refresher, here’s what that code looked like.

```
fit4 <-
brm(data = my_data,
family = student,
y_z ~ 1 + x_z + (1 + x_z || Subj),
prior = c(prior(normal(0, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(normal(0, 1), class = sigma),
prior(normal(0, 1), class = sd),
prior(exponential(one_over_twentynine) + 1, class = nu)),
chains = 4, cores = 4,
stanvars = stanvar(1/29, name = "one_over_twentynine"))
```

And for convenience, here’s the code from the model we just fit.

```
fit5 <-
brm(data = my_data,
family = gaussian,
Longevity ~ 1 + thorax_c + (1 | CompanionNumber),
prior = c(prior(normal(mean_y, sd_y * 5), class = Intercept),
prior(normal(0, 2 * sd_y / sd_Thorax_c), class = b),
prior(gamma(alpha, beta), class = sd),
prior(cauchy(0, sd_y), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
control = list(adapt_delta = 0.99),
stanvars = stanvars)
```

It’s easy to get lost in the differences in the priors and the technical details with the model chains and such. The main thing to notice, here, is the differences in the model formulas (i.e., the likelihoods). Both models had intercepts and slopes. But whereas the model from 17.3 set both parameters to random, only the intercept in our last modes was random. The covariate `thorax_c`

was fixed–it did not vary by group. Had we wanted it to, our `formula`

syntax would have been something like `Longevity ~ 1 + thorax_c + (1 + thorax_c || CompanionNumber)`

. And again, as noted in Chapter 17, the `||`

portion of the syntax set the random intercepts and slopes to be orthogonal (i.e., correlate exactly at zero). As we’ll see, this will often not be the case. But let’s not get ahead of ourselves.

Conceptually, the main difference between the models is merely the focus of attention. In the hierarchical linear regression model, the focus was on the slope coefficient. In that case, we were trying to estimate the magnitude of the slope, simultaneously for individuals and overall. The intercepts, which describe the levels of the nominal predictor, were of ancillary interest. In the present section, on the other hand, the focus of attention is reversed. We are most interested in the intercepts and their differences between groups, with the slopes on the covariate being of ancillary interest. (p. 573)

## 19.5 Heterogeneous variances and robustness against outliers

On page 574, Kruschke laid out the schematic for a hierarchical Student’s-\(t\) model in for which both the \(\mu\) and \(\sigma\) parameters are random. Bürkner calls these distributional models and they are indeed available within the **brms** framework. But there’s a catch. Though we can model \(\sigma\) all day long and we can even make it hierarchical, **brms** limits us to modeling the hierarchical \(\sigma\) parameters within the typical Gaussian framework. That is, we will depart from Kruschke’s schematic in that we will be

- modeling the log of \(\sigma\),
- indicating its grand mean with the
`sigma ~ 1`

syntax, - modeling the group-level deflections as Gaussian with a mean of 0 and standard deviation \(\sigma_\sigma\) estimated from the data,
- and choosing a sensible prior for \(\sigma_\sigma\) that is left-bound at 0 and gently slopes to the right (i.e., a folded \(t\) or gamma distribution).

Since we’re modeling \(\log(\sigma)\), we might use Gaussian prior centered on `sd(my_data$y) %>% log()`

and a reasonable spread like 1. We can simulate a little to get a sense of what those distributions look like.

```
n_draws <- 1e3
set.seed(19)
tibble(prior = rnorm(n_draws, mean = log(1), sd = 1)) %>%
mutate(prior_exp = exp(prior)) %>%
gather(key, value) %>%
ggplot(aes(x = value)) +
geom_density(fill = "grey50", color = "transparent") +
facet_wrap(~key, scales = "free")
```

Here’s what is looks like with `sd = 2`

.

```
set.seed(19)
tibble(prior = rnorm(n_draws, mean = log(1), sd = 2)) %>%
mutate(prior_exp = exp(prior)) %>%
ggplot(aes(x = prior_exp)) +
geom_density(fill = "grey50", color = "transparent") +
coord_cartesian(xlim = 0:17)
```

Though we’re still peaking around 1, there’s more mass in the tail, making it easier for the likelihood to pull away from the prior mode.

But all this is the prior on the fixed effect, the grand mean of \(\log (\sigma)\). Keep in mind we’re also estimating group-level deflections using a hierarchical model. The good old folded \(t\) on the unit scale is already pretty permissive for an estimate that is itself on the log scale. To make it more conservative, set \(\nu\) to infinity and go with a folded Gaussian. Or keep your regularization loose and go with a low-\(\nu\) folded \(t\) or even a folded Cauchy. And, of course, one could even go with a gamma.

Consider we have data `my_data`

for which our primary variable of interest is `y`

. Starting from preparing our `stanvars`

values, here’s what the model code might look like.

```
# get ready for `stanvars`
mean_y <- mean(my_data$y)
sd_y <- sd(my_data$y)
omega <- sd_y / 2
sigma <- 2 * sd_y
s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma)
# define `stanvars`
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta") +
stanvar(1/29, name = "one_over_twentynine")
# fit the model
fit <-
brm(data = my_data,
family = student,
bf(Longevity ~ 1 + (1 | CompanionNumber),
sigma ~ 1 + (1 | CompanionNumber)),
prior = c(# grand means
prior(normal(mean_y, sd_y * 5), class = Intercept),
prior(normal(log(sd_y), 1), class = Intercept, dpar = sigma),
# the priors controlling the spread for our hierarchical deflections
prior(gamma(alpha, beta), class = sd),
prior(normal(0, 1), class = sd, dpar = sigma),
# don't forget our student-t nu
prior(exponential(one_over_twentynine), class = nu)),
stanvars = stanvars)
```

### 19.5.1 Example: Contrast of means with different variances.

Let’s load and take a look at Kruschke’s simulated group data.

```
## # A tibble: 6 x 2
## Group Y
## <chr> <dbl>
## 1 A 97.8
## 2 A 99.9
## 3 A 92.4
## 4 A 96.9
## 5 A 101.
## 6 A 80.7
```

Here are the means and \(SD\)s for each `Group`

.

```
## # A tibble: 4 x 3
## Group mean sd
## <chr> <dbl> <dbl>
## 1 A 97.0 8.00
## 2 B 99.0 1.000
## 3 C 102. 1.
## 4 D 104. 8.00
```

First we’ll fit the model with homogeneous variances. To keep things simple, here we’ll fit a conventional model following the form of our original `fit1`

. Here are our `stanvars`

.

`## [1] 100.5`

`## [1] 6.228965`

```
## $shape
## [1] 1.283196
##
## $rate
## [1] 0.09092861
```

```
# define `stanvars`
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta")
```

Now fit the ANOVA-like homogeneous-variances model.

```
fit5 <-
brm(data = my_data,
family = gaussian,
Y ~ 1 + (1 | Group),
prior = c(prior(normal(mean_y, sd_y * 10), class = Intercept),
prior(gamma(alpha, beta), class = sd),
prior(cauchy(0, sd_y), class = sigma)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
control = list(adapt_delta = 0.999),
stanvars = stanvars)
```

Here’s the model summary.

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Y ~ 1 + (1 | Group)
## Data: my_data (Number of observations: 96)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~Group (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.96 3.56 1.46 14.67 1.00 1814 2546
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 100.47 3.19 94.43 106.42 1.00 2255 1885
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.76 0.43 4.99 6.68 1.00 5581 5873
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Let’s get ready to make our version of the top of Figure 19.7. First we wrangle.

```
# how many model-implied Gaussians would you like?
n_draws <- 20
densities <-
my_data %>%
distinct(Group) %>%
add_fitted_draws(fit5, n = n_draws, seed = 19, dpar = c("mu", "sigma")) %>%
mutate(ll = qnorm(.025, mean = mu, sd = sigma),
ul = qnorm(.975, mean = mu, sd = sigma)) %>%
mutate(Y = map2(ll, ul, seq, length.out = 100)) %>%
unnest(Y) %>%
mutate(density = dnorm(Y, mu, sigma)) %>%
group_by(.draw) %>%
mutate(density = density / max(density))
glimpse(densities)
```

```
## Observations: 8,000
## Variables: 12
## Groups: .draw [20]
## $ Group <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A"…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ .draw <int> 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, …
## $ .value <dbl> 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235…
## $ mu <dbl> 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235, 99.42235…
## $ sigma <dbl> 6.640072, 6.640072, 6.640072, 6.640072, 6.640072, 6.640072, 6.640072, 6.640072…
## $ ll <dbl> 86.40805, 86.40805, 86.40805, 86.40805, 86.40805, 86.40805, 86.40805, 86.40805…
## $ ul <dbl> 112.4367, 112.4367, 112.4367, 112.4367, 112.4367, 112.4367, 112.4367, 112.4367…
## $ Y <dbl> 86.40805, 86.67096, 86.93388, 87.19679, 87.45971, 87.72262, 87.98554, 88.24846…
## $ density <dbl> 0.1465288, 0.1582290, 0.1705958, 0.1836410, 0.1973740, 0.2118018, 0.2269281, 0…
```

In our wrangling code, the main thing to notice is those last two lines. If you look closely to Kruschke’s Gaussians, you’ll notice they all have the same maximum height. Up to this point, ours haven’t. This has to do with technicalities on how densities are scaled. In brief, the wider densities have been shorter. So those last two lines scaled all the densities within the same group to the same metric. Otherwise the code was business as usual.

Anyway, here’s our version of the top panel of Figure 19.7.

```
densities %>%
ggplot(aes(x = Y, y = Group)) +
geom_ridgeline(aes(height = density, group = interaction(Group, .draw)),
fill = NA, color = adjustcolor("grey50", alpha.f = 2/3),
size = 1/3, scale = 3/4) +
geom_jitter(data = my_data,
height = .04, alpha = 1/2) +
scale_x_continuous(breaks = seq(from = 80, to = 120, by = 10)) +
coord_cartesian(xlim = 75:125,
ylim = c(1.25, 4.5)) +
labs(title = "Data with Posterior Predictive Distrib.",
y = NULL) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

Here are the difference distributions in the middle of Figure 19.7.

```
post <- posterior_samples(fit5)
differences <-
post %>%
transmute(`D vs A` = `r_Group[D,Intercept]` - `r_Group[A,Intercept]`,
`C vs B` = `r_Group[C,Intercept]` - `r_Group[B,Intercept]`)
differences %>%
gather() %>%
mutate(key = factor(key, levels = c("D vs A", "C vs B"))) %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(x = "Difference") +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

Now here are the effect sizes at the bottom of the figure.

```
differences %>%
mutate_all(.funs = ~. / post$sigma) %>%
gather() %>%
mutate(key = factor(key, levels = c("D vs A", "C vs B"))) %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(x = "Effect Size") +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

Oh and remember, if you’d like to get all the possible contrasts in bulk, **tidybayes** has got your back.

```
fit5 %>%
spread_draws(r_Group[Group,]) %>%
compare_levels(r_Group, by = Group) %>%
# these next two lines allow us to reorder the contrasts along the y
ungroup() %>%
mutate(Group = reorder(Group, r_Group)) %>%
ggplot(aes(x = r_Group, y = Group)) +
geom_vline(xintercept = 0, color = "white") +
geom_halfeyeh(point_interval = mode_hdi, .width = c(.95, .5)) +
labs(x = "Contrast",
y = NULL) +
coord_cartesian(ylim = c(1.33, 6.5)) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

But to get back on track, here are the `stanvars`

for the robust hierarchical variances model.

```
stanvars <-
stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y") +
stanvar(s_r$shape, name = "alpha") +
stanvar(s_r$rate, name = "beta") +
stanvar(1/29, name = "one_over_twentynine")
```

Now fit that robust better-than-ANOVA model.

```
fit6 <-
brm(data = my_data,
family = student,
bf(Y ~ 1 + (1 | Group),
sigma ~ 1 + (1 | Group)),
prior = c(# grand means
prior(normal(mean_y, sd_y * 10), class = Intercept),
prior(normal(log(sd_y), 1), class = Intercept, dpar = sigma),
# the priors controlling the spread for our hierarchical deflections
prior(gamma(alpha, beta), class = sd),
prior(normal(0, 1), class = sd, dpar = sigma),
# don't forget our student-t nu
prior(exponential(one_over_twentynine), class = nu)),
iter = 4000, warmup = 1000, chains = 4, cores = 4,
seed = 19,
control = list(adapt_delta = 0.99,
max_treedepth = 12),
stanvars = stanvars)
```

The chains look good.

Here’s the summary.

```
## Family: student
## Links: mu = identity; sigma = log; nu = identity
## Formula: Y ~ 1 + (1 | Group)
## sigma ~ 1 + (1 | Group)
## Data: my_data (Number of observations: 96)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Group-Level Effects:
## ~Group (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 4.55 3.14 1.32 13.02 1.00 3148 4776
## sd(sigma_Intercept) 1.23 0.40 0.65 2.18 1.00 4350 5529
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 100.48 2.73 94.73 105.97 1.00 4094 4128
## sigma_Intercept 1.23 0.54 0.21 2.37 1.00 4568 5951
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## nu 32.67 28.51 4.65 109.84 1.00 10110 8384
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Let’s get ready to make our version of the top of Figure 19.7. First we wrangle.

```
densities <-
my_data %>%
distinct(Group) %>%
add_fitted_draws(fit6, n = n_draws, seed = 19, dpar = c("mu", "sigma", "nu")) %>%
mutate(ll = qt(.025, df = nu),
ul = qt(.975, df = nu)) %>%
mutate(Y = map2(ll, ul, seq, length.out = 100)) %>%
unnest(Y) %>%
mutate(density = dt(Y, nu)) %>%
# notice the conversion
mutate(Y = mu + Y * sigma) %>%
group_by(.draw) %>%
mutate(density = density / max(density))
glimpse(densities)
```

```
## Observations: 8,000
## Variables: 13
## Groups: .draw [20]
## $ Group <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A"…
## $ .row <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ .chain <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ .iteration <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA…
## $ .draw <int> 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, 1115, …
## $ .value <dbl> 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906…
## $ mu <dbl> 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906, 97.33906…
## $ sigma <dbl> 6.611823, 6.611823, 6.611823, 6.611823, 6.611823, 6.611823, 6.611823, 6.611823…
## $ nu <dbl> 24.33915, 24.33915, 24.33915, 24.33915, 24.33915, 24.33915, 24.33915, 24.33915…
## $ ll <dbl> -2.062378, -2.062378, -2.062378, -2.062378, -2.062378, -2.062378, -2.062378, -…
## $ ul <dbl> 2.062378, 2.062378, 2.062378, 2.062378, 2.062378, 2.062378, 2.062378, 2.062378…
## $ Y <dbl> 83.70299, 83.97846, 84.25394, 84.52941, 84.80489, 85.08037, 85.35584, 85.63132…
## $ density <dbl> 0.1299849, 0.1401936, 0.1510374, 0.1625369, 0.1747115, 0.1875784, 0.2011531, 0…
```

If you look closely at our code, above, you’ll note switching from the Gaussian to the Student \(t\) required changes in our flow. Most obviously, we switched from `qnorm()`

and `dnorm()`

to `qt()`

and `dt()`

, respectively. The base **R** Student \(t\) functions don’t take arguments for \(\mu\) and \(\sigma\). Rather, they’re presumed to be 0 and 1, respectively. That means that for our first three `mutate()`

functions, the computations were all based on the standard Student \(t\), with only the \(\nu\) parameter varying according to the posterior. The way we corrected for that was with the fourth `mutate()`

.

Now we’re ready to make our version of the top panel of Figure 19.7.

```
densities %>%
ggplot(aes(x = Y, y = Group)) +
geom_ridgeline(aes(height = density, group = interaction(Group, .draw)),
fill = NA, color = adjustcolor("grey50", alpha.f = 2/3),
size = 1/3, scale = 3/4) +
geom_jitter(data = my_data,
height = .04, alpha = 1/2) +
scale_x_continuous(breaks = seq(from = 80, to = 120, by = 10)) +
coord_cartesian(xlim = 75:125,
ylim = c(1.25, 4.5)) +
labs(title = "Data with Posterior Predictive Distrib.",
y = NULL) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
```

Here are the difference distributions in the middle of Figure 19.8.

```
post <- posterior_samples(fit6)
post %>%
transmute(`D vs A` = `r_Group[D,Intercept]` - `r_Group[A,Intercept]`,
`C vs B` = `r_Group[C,Intercept]` - `r_Group[B,Intercept]`) %>%
gather() %>%
mutate(key = factor(key, levels = c("D vs A", "C vs B"))) %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Difference") +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

And here are the corresponding effect sizes at the bottom of the Figure 19.8.

```
post %>%
transmute(`D vs A` = (`r_Group[D,Intercept]` - `r_Group[A,Intercept]`) /
exp((b_sigma_Intercept + `r_Group__sigma[D,Intercept]` + b_sigma_Intercept + `r_Group__sigma[A,Intercept]`) / 2),
`C vs B` = (`r_Group[C,Intercept]` - `r_Group[B,Intercept]`) /
exp((b_sigma_Intercept + `r_Group__sigma[C,Intercept]` + b_sigma_Intercept + `r_Group__sigma[B,Intercept]`) / 2)) %>%
gather() %>%
mutate(key = factor(key, levels = c("D vs A", "C vs B"))) %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("Effect Size") +
facet_wrap(~key, scales = "free_x", ncol = 4)
```

Notice that because (a) the sigma parameters were heterogeneous and (b) they were estimated on the log scale, we had to do quite a bit more data processing before they effect size estimates were ready.

“Finally, because each group has its own estimated scale (i.e., \(\sigma_j\)), we can investigate differences in scales across groups” (p. 578). That’s not a bad idea. Even though Kruschke didn’t show this in the text, we may as well give it a go.

```
post %>%
transmute(`D vs A` = exp(b_sigma_Intercept + `r_Group__sigma[D,Intercept]`) - exp(b_sigma_Intercept + `r_Group__sigma[A,Intercept]`),
`C vs B` = exp(b_sigma_Intercept + `r_Group__sigma[C,Intercept]`) - exp(b_sigma_Intercept + `r_Group__sigma[B,Intercept]`),
`D vs C` = exp(b_sigma_Intercept + `r_Group__sigma[D,Intercept]`) - exp(b_sigma_Intercept + `r_Group__sigma[C,Intercept]`),
`B vs A` = exp(b_sigma_Intercept + `r_Group__sigma[B,Intercept]`) - exp(b_sigma_Intercept + `r_Group__sigma[A,Intercept]`)) %>%
gather() %>%
mutate(key = factor(key, levels = c("D vs A", "C vs B", "D vs C", "B vs A"))) %>%
ggplot(aes(x = value)) +
geom_histogram(color = "grey92", fill = "grey67",
size = .2, bins = 40) +
stat_pointintervalh(aes(y = 0),
point_interval = mode_hdi, .width = c(.95, .5)) +
scale_y_continuous(NULL, breaks = NULL) +
xlab(expression(paste("Differences in ", sigma))) +
facet_wrap(~key, scales = "free", ncol = 4)
```

For more on models including a hierarchical structure on both the mean and scale structures, check out Donald Williams and colleagues’ work on what they call Mixed Effect Location and Scale Models (MELSM; e.g., here or here). They’re quite handy and I’ve begun using them in my applied work (e.g., here).

## Session info

```
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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_1.1.0 bayesplot_1.7.0 brms_2.10.3 Rcpp_1.0.2 ggridges_0.5.1 forcats_0.4.0
## [7] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3 readr_1.3.1 tidyr_1.0.0 tibble_2.1.3
## [13] ggplot2_3.2.1 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.4-1 ellipsis_0.3.0 rsconnect_0.8.15
## [4] ggstance_0.3.2 markdown_1.1 base64enc_0.1-3
## [7] rstudioapi_0.10 rstan_2.19.2 svUnit_0.7-12
## [10] DT_0.9 fansi_0.4.0 lubridate_1.7.4
## [13] xml2_1.2.0 bridgesampling_0.7-2 knitr_1.23
## [16] shinythemes_1.1.2 zeallot_0.1.0 jsonlite_1.6
## [19] broom_0.5.2 shiny_1.3.2 compiler_3.6.0
## [22] httr_1.4.0 backports_1.1.5 assertthat_0.2.1
## [25] Matrix_1.2-17 lazyeval_0.2.2 cli_1.1.0
## [28] later_1.0.0 htmltools_0.4.0 prettyunits_1.0.2
## [31] tools_3.6.0 igraph_1.2.4.1 coda_0.19-3
## [34] gtable_0.3.0 glue_1.3.1.9000 reshape2_1.4.3
## [37] cellranger_1.1.0 vctrs_0.2.0 nlme_3.1-139
## [40] crosstalk_1.0.0 xfun_0.10 ps_1.3.0
## [43] rvest_0.3.4 mime_0.7 miniUI_0.1.1.1
## [46] lifecycle_0.1.0 gtools_3.8.1 zoo_1.8-6
## [49] scales_1.0.0 colourpicker_1.0 hms_0.4.2
## [52] promises_1.1.0 Brobdingnag_1.2-6 parallel_3.6.0
## [55] inline_0.3.15 shinystan_2.5.0 yaml_2.2.0
## [58] gridExtra_2.3 loo_2.1.0 StanHeaders_2.19.0
## [61] stringi_1.4.3 dygraphs_1.1.1.6 pkgbuild_1.0.5
## [64] rlang_0.4.1 pkgconfig_2.0.3 matrixStats_0.55.0
## [67] HDInterval_0.2.0 evaluate_0.14 lattice_0.20-38
## [70] rstantools_2.0.0 htmlwidgets_1.5 labeling_0.3
## [73] tidyselect_0.2.5 processx_3.4.1 plyr_1.8.4
## [76] magrittr_1.5 R6_2.4.0 generics_0.0.2
## [79] pillar_1.4.2 haven_2.1.0 withr_2.1.2
## [82] xts_0.11-2 abind_1.4-5 modelr_0.1.4
## [85] crayon_1.3.4 arrayhelpers_1.0-20160527 utf8_1.1.4
## [88] rmarkdown_1.13 grid_3.6.0 readxl_1.3.1
## [91] callr_3.3.2 threejs_0.3.1 digest_0.6.21
## [94] xtable_1.8-4 httpuv_1.5.2 stats4_3.6.0
## [97] munsell_0.5.0 viridisLite_0.3.0 shinyjs_1.0
```