19.Rmd

---
title: "Chapter 19. Metric Predicted Variable with One Nominal Predictor"
author: "A Solomon Kurz"
date: "`r format(Sys.Date())`"
output:
  github_document
---

```{r, echo = FALSE, cache = FALSE}
knitr::opts_chunk$set(fig.retina = 2.5)
knitr::opts_chunk$set(fig.align = "center")
options(width = 100)
```

# 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. [@kruschkeDoingBayesianData2015, pp. 553--554]

## 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.

```{r, warning = F, message = F}
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 = map2(group, deviation, generate_data)) %>% 
  unnest(d) %>% 
  mutate(iteration = rep(1:n, times = 5))

glimpse(d)
```

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`.

```{r}
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)
```

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.

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

head(arrow)

text <-
  tibble(d         = grand_mean,
         group     = c(0:5, 0),
         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)
```

Now we're ready to plot.

```{r, fig.width = 6, fig.height = 2.75, warning = F, message = F}
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>")) +
  xlab(NULL) +
  coord_flip(xlim = c(90, 112),
             ylim = c(-0.2, 5.5)) +
  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)

## 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 [@fisherStatisticalMethodsResearch1925]. 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)

## 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 [@gelmanPriorDistributionsVariance2006] of putting a broad uniform prior on $\sigma_y$. Rather, we will continue using the half-Gaussian prior, as [recommended by the Stan team](https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations). 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 [-@efronSteinParadoxStatistics1977] paper, [*Stein's paradox in statistics*](http://statweb.stanford.edu/~ckirby/brad/other/Article1977.pdf), and [my blog post](https://solomonkurz.netlify.com/post/stein-s-paradox-and-what-partial-pooling-can-do-for-you/) walking out one of their examples in **brms**.

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

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

```{r, eval = F}
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 [-@chungNondegeneratePenalizedLikelihood2013, click [here](http://www.stat.columbia.edu/~gelman/research/published/chung_etal_Pmetrika2013.pdf)]. 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 @chungNondegeneratePenalizedLikelihood2013. Here's what $\operatorname{gamma} (2, 0.1)$ looks like. 

```{r, fig.width = 6, fig.height = 2}
tibble(x = seq(from = 0, to = 110, by = .1)) %>% 
  
  ggplot(aes(x = x, ymin = 0, ymax = dgamma(x, 2, 0.1))) +
  geom_ribbon(size = 0, fill = "grey67") +
  annotate(geom = "text", x = 14.25, y = 0.015, label = "'gamma'~(2*', '*0.1)", 
           parse = T, color = "grey92", size = 4.25) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(0, 100)) +
  theme(panel.grid = element_blank())
```

If you'd like that prior be even less informative, just reduce it to like $\operatorname{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.

```{r}
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.

```{r}
sd_y  <- 1 

omega <- sd_y / 2
sigma <- 2 * sd_y

(s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma))
```

And that produces the following gamma distribution.

```{r, fig.width = 6, fig.height = 2}
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 = c(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.

### Example: Sex and death.

Let's load and `glimpse()` at Hanley and Shapiro's [-@hanleySexualActivityLifespan1994] fruit-fly data.

```{r, message = F}
my_data <- read_csv("data.R/FruitflyDataReduced.csv")

glimpse(my_data)
```

We can use `geom_density_ridges()` to help get a sense of how our criterion `Longevity` is distributed across groups of `CompanionNumber`.

```{r, fig.width = 6, fig.height = 2.5, message = F}
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(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        legend.position = "none",
        panel.grid = element_blank())
```

Let's fire up **brms**.

```{r, warning = F, message = F}
library(brms)
```

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

```{r}
(mean_y <- mean(my_data$Longevity))
(sd_y <- sd(my_data$Longevity))

omega <- sd_y / 2
sigma <- 2 * sd_y

(s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma))
```

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

```{r}
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.

```{r fit19.1}
fit19.1 <-
  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,
      file = "fits/fit19.01")
``` 

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

```{r, fig.width = 8, fig.height = 4}
plot(fit19.1)
```

Also like Kruschke, our chains appear moderately autocorrelated.

```{r, message = F, warning = F, fig.width = 7, fig.height = 4}
post <- posterior_samples(fit19.1, add_chain = T)

library(bayesplot)

theme_set(theme_gray() +
            theme(panel.grid = element_blank()))

mcmc_acf(post, pars = c("b_Intercept", "sd_CompanionNumber__Intercept", "sigma"), lags = 10)
```

Here's the model summary.

```{r}
print(fit19.1)
```

With the `ranef()` function, we can get the summaries of the group-specific deflections.

```{r}
ranef(fit19.1)
```

And with the `coef()` function, we can get those same group-level summaries in a non-deflection metric.

```{r}
coef(fit19.1)
```

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

```{r}
posterior_summary(fit19.1)["sigma", ]
```

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.

```{r}
# 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)
```

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()`.

```{r, fig.height = 2.5, fig.width = 4}
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) +
  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.") +
  coord_cartesian(xlim = c(0, 110))
```

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**.

```{r, warning = F, message = F}
library(tidybayes)
```

For the first part, we'll take `tidybayes::add_fitted_draws()` for a whirl.

```{r}
densities <-
  my_data %>% 
  distinct(CompanionNumber) %>% 
  add_fitted_draws(fit19.1, n = 20, seed = 19, dpar = c("mu", "sigma"))

glimpse(densities)
```

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** [see the [Posterior fits](https://mjskay.github.io/tidybayes/articles/tidy-brms.html#posterior-fits) section of @kayExtractingVisualizingTidy2020a]. 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, `fit19.1`. 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.

```{r}
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)
```

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.

```{r, fig.height = 4.5, fig.width = 4}
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)) +
  labs(title = "Data with Posterior Predictive Distrib.", 
       y = NULL) +
  coord_cartesian(xlim = c(0, 110),
                  ylim = c(1.25, 5.25)) +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())
```

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](https://sites.google.com/site/doingbayesiandataanalysis/corrigenda). 

### 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.

```{r, fig.width = 3, fig.height = 2.5}
post %>% 
  transmute(c = (`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`) %>% 
  
  ggplot(aes(x = c, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T) +
  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.

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

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

```{r, fig.width = 3, fig.height = 2.5}
post %>% 
  transmute(es = ((`r_CompanionNumber[Pregnant1,Intercept]` + `r_CompanionNumber[Pregnant1,Intercept]`) / 2 - `r_CompanionNumber[None0,Intercept]`) / sigma) %>% 
  
  ggplot(aes(x = es, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T) +
  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.

```{r, fig.width = 8, fig.height = 2.5}
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, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("Difference") +
  facet_wrap(~key, scales = "free")
```

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.

```{r, fig.width = 8, fig.height = 2.5}
differences %>% 
  mutate_all(.funs = ~ . / post$sigma) %>% 
  gather() %>% 
  
  ggplot(aes(x = value, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  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 [see the [Comparing levels of a factor](https://mjskay.github.io/tidybayes/articles/tidy-brms.html#comparing-levels-of-a-factor) section of @kayExtractingVisualizingTidy2020a]. Here we'll demonstrate with `geom_halfeyeh()`.

```{r, fig.height = 4.5, fig.width = 6}
fit19.1 %>%
  # 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.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())
```

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.

```{r fit19.2}
fit19.2 <-
  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,
      file = "fits/fit19.02")
``` 

Here's the model summary.

```{r}
print(fit19.2)
```

Here are their LOO values and their difference score.

```{r, message = F}
fit19.1 <- add_criterion(fit19.1, criterion = "loo")
fit19.2 <- add_criterion(fit19.2, criterion = "loo")

loo_compare(fit19.1, fit19.2) %>% 
  print(simplify = F)
```

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

```{r}
(mw <- model_weights(fit19.1, fit19.2))
```

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

```{r}
mw %>% 
  round(digits = 3)
```

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.

```{r, fig.width = 4, fig.height = 2.5}
post %>% 
  ggplot(aes(x = sd_CompanionNumber__Intercept, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 100, slab_size = .25, outline_bars = T) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = c(0, 50)) +
  labs(title = expression("Behold the fit19.1 posterior for "*sigma[beta]*"."),
       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.

### 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.

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

For our hierarchical model `fit19.1`, the posterior means are rank ordered in the same way as the empirical data.

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

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.

```{r}
coef(fit19.1)$CompanionNumber[, , "Intercept"] %>% 
  as_tibble() %>% 
  summarise(range = max(Estimate) - min(Estimate))
```

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).

```{r}
coef(fit19.1, summary = F)$CompanionNumber[, , "Intercept"] %>% 
  as_tibble() %>% 
  transmute(range = Pregnant1 - Virgin8) %>% 
  mode_hdi(range)
```

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.

### 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()`.

```{r fit19.3}
fit19.3 <-
  update(fit19.1,
         newdata = my_data %>% 
           filter(CompanionNumber %in% c("Pregnant1", "Virgin8")),
         control = list(adapt_delta = 0.999,
                        max_treedepth = 12),
         seed = 19,
         file = "fits/fit19.03")
```

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

```{r}
print(fit19.3)
```

If you compare the posteriors for $\sigma_\beta$ across the two models, you'll see how the one for `fit19.3` is substantially larger.

```{r}
posterior_summary(fit19.1)["sd_CompanionNumber__Intercept", ]
posterior_summary(fit19.3)["sd_CompanionNumber__Intercept", ]
```

Here that is in a coefficient plot using `tidybayes::stat_intervalh()`.

```{r, fig.width = 6, fig.height = 1.25}
bind_rows(posterior_samples(fit19.1) %>% select(sd_CompanionNumber__Intercept),
          posterior_samples(fit19.3) %>% select(sd_CompanionNumber__Intercept)) %>% 
  mutate(fit = rep(c("fit19.1", "fit19.3"), each = n() / 2)) %>% 
  
  ggplot(aes(x = sd_CompanionNumber__Intercept, y = fit)) +
  stat_intervalh(point_interval = mode_hdi, .width = c(.5, .8, .95)) +
  scale_color_grey("HDI", start = 2/3, end = 0, 
                   labels = c("95%", "80%", "50%")) +
  labs(x = expression(sigma[beta]),
       y = NULL) +
  coord_cartesian(xlim = c(0, 80)) +
  theme(legend.key.size = unit(0.45, "cm"))
```

This all implies less shrinkage and a larger range. 

```{r, warning = F, message = F}
coef(fit19.3, summary = F)$CompanionNumber[, , "Intercept"] %>% 
  as_tibble() %>% 
  transmute(range = Pregnant1 - Virgin8) %>% 
  mode_hdi(range)
```

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.

```{r, eval = F}
fit16.3 <-
  brm(data = my_data,
      family = student,
      bf(Score ~ 0 + Group, 
         sigma ~ 0 + Group),
      prior = c(prior(normal(mean_y, sd_y * 100), 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,
      seed = 16,
      file = "fits/fit16.03")
```

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.

```{r}
# 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))
(sd_y <- sd(my_small_data$Longevity))

omega <- sd_y / 2
sigma <- 2 * sd_y

(s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma))
```

Here we update `stanvars`.

```{r}
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, `fit19.1`.

```{r, fit19.4}
fit19.4 <-
  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,
      file = "fits/fit19.04")
```

Here's the model summary.

```{r}
print(fit19.4)
```

Man, look at those `Bulk_ESS` values! As it turns out, they can be [greater than the number of post-warmup samples](https://andrewgelman.com/2018/01/18/measuring-speed-stan-incorrectly-faster-thought-cases-due-antithetical-sampling/). And here's the range in posterior means.

```{r, warning = F, message = F}
fixef(fit19.4, summary = F) %>% 
  as_tibble() %>% 
  transmute(range = CompanionNumberPregnant1 - CompanionNumberVirgin8) %>% 
  mode_hdi(range)
```

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.

## 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.

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

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)

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

Kruschke recalled `fit19.1`'s estimate for $\sigma_y$ had a posterior mode around 14.8. Let's confirm with a plot.

```{r, fig.width = 4, fig.height = 2}
posterior_samples(fit19.1) %>% 
  ggplot(aes(x = sigma, y = 0)) +
  geom_halfeyeh(point_interval = 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.

```{r}
my_data <-
  my_data %>% 
  mutate(thorax_c = Thorax - mean(Thorax))

head(my_data)
```

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`.

```{r}
(mean_y <- mean(my_data$Longevity))
(sd_y <- sd(my_data$Longevity))
(sd_thorax_c <- sd(my_data$thorax_c))

omega <- sd_y / 2
sigma <- 2 * sd_y

(s_r <- gamma_a_b_from_omega_sigma(mode = omega, sd = sigma))

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.

```{r fit19.5}
fit19.5 <-
  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,
      file = "fits/fit19.05")
``` 

Here's the model summary.

```{r}
print(fit19.5)
```

Let's see if that $\sigma_y$ posterior shrank.

```{r, fig.width = 4, fig.height = 2}
posterior_samples(fit19.5) %>% 
  ggplot(aes(x = sigma, y = 0)) +
  geom_halfeyeh(point_interval = 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.

```{r, fig.width = 4, fig.height = 2.25}
bind_cols(posterior_samples(fit19.1) %>% select(sigma),
          posterior_samples(fit19.5) %>% select(sigma)) %>% 
  transmute(dif = sigma - sigma1) %>% 
  
  ggplot(aes(x = dif, y = 0)) +
  geom_halfeyeh(point_interval = mode_hdi, .width = c(.5, .95)) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "This is a difference distribution",
       x = expression(sigma[y][" | fit19.1"]-sigma[y][" | fit19.5"]))
```

If you want a quick and dirty plot of the relation between `thorax_c` and `Longevity`, you might employ `brms::conditional_effects()`.

```{r, fig.width = 3, fig.height = 3}
conditional_effects(fit19.5)
```

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()`.

```{r}
(r <- range(my_data$Thorax))
mean(r)
```

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.

```{r}
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(fit19.5, 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)
```

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. 

```{r}
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(fit19.5, n = n_draws, seed = 19, value = "Longevity")

glimpse(f)
```

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

```{r, fig.width = 8, fig.height = 2.5}
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 = c(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 we make and save them in the simple difference metric.

```{r}
post <- posterior_samples(fit19.5)

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)

p1 <-
  differences %>% 
  gather() %>%   
  
  ggplot(aes(x = value, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  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 leg work above, it's really easy to just convert the differences in bulk with `mutate_all()`. After the conversion, we'll bind the two rows of subplots together with a little **patchwork** and display the results.

```{r, fig.width = 8, fig.height = 5}
p2 <-
  differences %>% 
  mutate_all(.funs = ~. / post$sigma) %>% 
  gather() %>%   
  
  ggplot(aes(x = value, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.95, .5),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  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)

library(patchwork)
p1 / p2
```

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

### 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)

### 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.

```{r, eval = F}
fit17.4 <-
  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),
                # the next line is new
                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"),
      seed = 17,
      file = "fits/fit17.04")
```

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

```{r eval = F}
fit19.5 <-
  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,
      file = "fits/fit19.05")
``` 

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 model 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)

## 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. If you recall, @Bürkner2020Distributional calls these [distributional models](https://CRAN.R-project.org/package=brms/vignettes/brms_distreg.html) 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.

```{r, fig.width = 6, fig.height = 2}
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`.

```{r, fig.width = 3, fig.height = 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 = c(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.

```{r, eval = F}
# 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)
``` 

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

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

```{r, message = F}
my_data <- read_csv("data.R/NonhomogVarData.csv")

head(my_data)
```

Here are the means and $SD$s for each `Group`.

```{r}
my_data %>% 
  group_by(Group) %>% 
  summarise(mean = mean(Y),
            sd   = sd(Y))
```

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`.

```{r}
(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")
```

Now fit the ANOVA-like homogeneous-variances model.

```{r fit19.6}
fit19.6 <-
  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,
      file = "fits/fit19.06")
``` 

Here's the model summary.

```{r}
print(fit19.6)
```

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

```{r}
# how many model-implied Gaussians would you like?
n_draws <- 20

densities <-
  my_data %>% 
  distinct(Group) %>% 
  add_fitted_draws(fit19.6, 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)
```

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.

```{r, fig.height = 3.5, fig.width = 3.75}
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)) +
  labs(title = "Data with Posterior Predictive Distrib.", 
       y = NULL) +
  coord_cartesian(xlim = c(75, 125),
                  ylim = c(1.25, 4.5)) +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())
```

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

```{r, fig.width = 4, fig.height = 2.5}
post <- posterior_samples(fit19.6)

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, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.5, .95),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  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.

```{r, fig.width = 4, fig.height = 2.5}
differences %>% 
  mutate_all(.funs = ~. / post$sigma) %>% 
  gather() %>%
  mutate(key = factor(key, levels = c("D vs A", "C vs B"))) %>% 
  
  ggplot(aes(x = value, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.5, .95),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 30, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  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.

```{r, fig.height = 4, fig.width = 6}
fit19.6 %>%
  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.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())
```

But to get back on track, here are the `stanvars` for the robust hierarchical variances model.

```{r}
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.

```{r fit19.7}
fit19.7 <-
  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,
      file = "fits/fit19.07")
``` 

The chains look good.

```{r, fig.width = 8, fig.height = 5}
plot(fit19.7)
```

Here's the parameter summary.

```{r}
print(fit19.7)
```

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

```{r}
densities <-
  my_data %>% 
  distinct(Group) %>% 
  add_fitted_draws(fit19.7, 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)
```

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 and save our version of the top panel of Figure 19.7.

```{r}
p1 <-
  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)) +
  labs(title = "Data with Posterior Predictive Distrib.", 
       y = NULL) +
  coord_cartesian(xlim = c(75, 125),
                  ylim = c(1.25, 4.5)) +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank())
```

Here we make the difference distributions in the middle of Figure 19.8.

```{r, fig.width = 4, fig.height = 2.5}
post <- posterior_samples(fit19.7)

p2 <-
  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, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.5, .95),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 40, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("Difference") +
  facet_wrap(~key, scales = "free_x", ncol = 4)
```

And here we make the plots of the corresponding effect sizes at the bottom of the Figure 19.8.

```{r, fig.width = 4, fig.height = 2.5}
p3 <-
  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, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.5, .95),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 30, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("Effect Size") +
  facet_wrap(~key, scales = "free_x", ncol = 4)
```

Combine them all and plot!

```{r, fig.width = 4.5, fig.height = 8}
p1 / p2 / p3 + plot_layout(heights = c(2, 1, 1))
```

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.

```{r, fig.width = 8, fig.height = 2.5}
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, y = 0)) +
  stat_histintervalh(point_interval = mode_hdi, .width = c(.5, .95),
                     fill = "grey67", slab_color = "grey92",
                     breaks = 30, slab_size = .25, outline_bars = T,
                     normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(expression(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](https://twitter.com/wdonald_1985) and colleagues' work on what they call Mixed Effect Location and Scale Models [MELSM; e.g., @williamsBayesianMultivariateMixedeffects2019; @williamsBayesianNonlinearMixedeffects2019]. They're quite handy and I've begun using them in my applied work (e.g., [here](https://osf.io/vekpf/)).

## ~~Exercises~~ Walk out an effect size

We computed a lot of effect sizes in this chapter. They were all standardized mean differences. @cohenStatisticalPowerAnalysis1988 discussed these kinds of effect sizes in this way:

> We need a "pure" number, one free of our original measurement unit, with which to index what can be alternately called the degree of departure from the null hypothesis of the alternate hypothesis, or the ES (effect size) we wish to detect. This is accomplished by standardizing the raw effect size as expressed in the measurement unit of the dependent variable by dividing it by the (common) standard deviation of the measures in their respective populations, the latter also in the original measurement unit. (p. 20)

Though Cohen framed his discussion in terms of null-hypothesis significance testing, we can just as easily apply it to our Bayesian modeling framework. The main thing is we can use his definitions form above to define a particular kind of effect size—the standardized mean difference between two groups. This is commonly referred to as a Cohen's $d$, which follows the formula

$$d = \frac{\bar y_A - \bar y_B}{s_y},$$

where the unstandardized means of the variable of interest $y$ are compared between two groups, A and B. From the raw data, we compute their two means, $\bar y_A$ and $\bar y_B$, and divide their difference by the common (i.e., pooled) standard deviation $s_y$. As is the typical case, the empirically-derived means and standard deviations are stand-ins (i.e., estimates) of the population parameters. If we're willing to ignore uncertainty, we can do this all by hand.

Let's walk this out with the fruit-fly data from section 19.3.2.

```{r, message = F}
my_data <- read_csv("data.R/FruitflyDataReduced.csv")

glimpse(my_data)
```

Recall we have five groups indexed by `CompanionNumber`, each with $n = 25$.

```{r}
my_data %>% 
  count(CompanionNumber)
```

Let's focus on just two groups, the male fruit flies for which individual males were supplied access to one or with virgin female fruit flies per day. In the data, these are `CompanionNumber == Virgin1` and `CompanionNumber == Virgin8`, respectively. Here's a look at their mean `Longevity` values.

```{r}
my_data %>% 
  filter(str_detect(CompanionNumber, "Virgin")) %>% 
  group_by(CompanionNumber) %>% 
  summarise(mean = mean(Longevity))
```

If we're willing to treat the males in the `Virgin1` group as group "a" and those in the `Virgin8` group as group "b", we can save those mean values like so.

```{r}
y_bar_a <- 56.76
y_bar_b <- 38.72
```

Now we'll compute their pooled standard deviation.

```{r}
my_data %>% 
  filter(str_detect(CompanionNumber, "Virgin")) %>% 
  summarise(s = sd(Longevity))
```

Save that value.

```{r}
s_y <- 16.24533
```

Now computing Cohen's $d$ is just simple arithmetic.

```{r}
(y_bar_a - y_bar_b) / s_y
```

Though I'm not up on contemporary standards in fruit fly research, a Cohen’s $d$ of that size would be considered [conspicuously] large in most areas of my field (psychology). If we'd like to compute the $d$ estimates for any other combination of experimental conditions, we'd just follow the corresponding arithmetic. 

As I hinted at earlier, the problem with this approach is it ignores uncertainty. Frequentists use various formulas to express this in terms of 95% confidence intervals. Our approach will be to express it with the posterior distribution of a Bayesian model. We've already accomplished this with our `fit19.1` from above. Here we'll use three other approaches. 

Instead of the Bayesian hierarchical alternative to the frequentist ANOVA, we can use a single-level model where we predict a metric variable with separate intercepts for the two levels of `CompanionNumber`. First, we subset the data and define our `stanvars`.

```{r}
my_data <-
  my_data %>% 
  filter(str_detect(CompanionNumber, "Virgin"))

mean_y <- (y_bar_a + y_bar_b) / 2

stanvars <- 
  stanvar(mean_y, name = "mean_y") + 
  stanvar(s_y,    name = "sd_y")
```

Fit the model with `brm()`.

```{r fit19.8}
fit19.8 <-
  brm(data = my_data,
      family = gaussian,
      Longevity ~ 0 + CompanionNumber,
      prior = c(prior(normal(mean_y, sd_y * 5), class = b),
                prior(cauchy(0, sd_y), class = sigma)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 19,
      stanvars = stanvars,
      file = "fits/fit19.08")
```

Check the summary.

```{r}
print(fit19.8)
```

Extract the posterior draws.

```{r}
post <- posterior_samples(fit19.8)
```

Here we'll plot the three dimensions of the posterior, each with the corresponding value from the Cohen's $d$ formula marked off as a vertical line in the foreground.

```{r, fig.width = 8, fig.height = 2.5}
lines <-
  tibble(name  = c("b_CompanionNumberVirgin1", "b_CompanionNumberVirgin8", "sigma"),
         value = c(y_bar_a, y_bar_b, s_y))

post %>% 
  pivot_longer(-lp__) %>% 
  
  ggplot(aes(x = value, y = 0)) +
  stat_halfeyeh(point_interval = mode_hdi, .width = .95,
                fill = "grey67", slab_color = "grey92",
                normalize = "panels") +
  geom_vline(data = lines,
             aes(xintercept = value),
             color = "grey33", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("posterior") +
  facet_wrap(~name, scales = "free")
```

The model did a great job capturing the group means. Notice how the posterior `sigma` is a bit lower than the value $s_y$. This is because $\sigma$ in the model is conditioned on the mean. Because we have two means, one for each group, that has 'explained some of the variation' in the criterion variable `Longevity`. Thus, our posterior for $\sigma$ might be more fully expressed as $\sigma | \text{CompanionNumber}$. If we would like to compute our Cohen's $d$ using the posterior iterations from `fit19.8`, we'd execute something like this.

```{r, fig.width = 3, fig.height = 2.5}
post %>% 
  mutate(d = (b_CompanionNumberVirgin1 - b_CompanionNumberVirgin8) / s_y) %>% 
  
  ggplot(aes(x = d, y = 0)) +
  stat_halfeyeh(point_interval = mode_hdi, .width = .95,
                fill = "grey67", slab_color = "grey92") +
  geom_vline(xintercept = (y_bar_a - y_bar_b) / s_y,
             color = "grey33", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(expression(paste("Cohen's ", italic(d), " expressed as a posterior")))
```

Similar to the previous plots, this time we superimposed the density with the estimate for $d$ we computed above, `(y_bar_a - y_bar_b) / s_y`. Happily, the hand-calculated estimate coheres nicely with the central tendency of our posterior distribution. But now we get a full measure of uncertainty. Notice how wide those 95% HDIs are. Hopefully this isn't a surprise given our noncommittal priors and only $n = 25$ for both groups. There's a lot of uncertainty in that posterior.

A second way we might use a single-level model to compute a Cohen's $d$ effect size is using a dummy variable. We'll convert our nominal variable `CompanionNumber` into a binary variable `Virgin1` for which 1 corresponds to `CompanionNumber == Virgin1` and 0 corresponds to `CompanionNumber == Virgin8`. Compute the dummy.

```{r}
my_data <-
  my_data %>% 
  mutate(Virgin1 = if_else(CompanionNumber == "Virgin1", 1, 0))
```

Now fit the dummy-predictor model with `brm()`.

```{r fit19.9}
fit19.9 <-
  brm(data = my_data,
      family = gaussian,
      Longevity ~ 1 + Virgin1,
      prior = c(prior(normal(mean_y, sd_y * 5), class = Intercept),
                prior(normal(0, sd_y * 5), class = b),
                prior(cauchy(0, sd_y), class = sigma)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 19,
      stanvars = stanvars,
      file = "fits/fit19.09")
```

```{r}
print(fit19.9)
```

With this parameterization, our posterior for `Intercept` is the same, within simulation variation, as `CompanionNumberVirgin8` from `fit7`. The posterior for `sigma` is about the same for both models, too. But focus on `Virgin1`. This is the unstandardized mean difference, what we called $\bar y_A - \bar y_B$ in our formula for Cohen's $d$. Here's a look at its posterior distribution with its empirical estimate superimposed with a vertical line.

```{r, fig.width = 3, fig.height = 2.5}
post <- posterior_samples(fit19.9)

post %>% 
  ggplot(aes(x = b_Virgin1, y = 0)) +
  stat_halfeyeh(point_interval = mode_hdi, .width = .95,
                fill = "grey67", slab_color = "grey92") +
  geom_vline(xintercept = y_bar_a - y_bar_b,
             color = "grey33", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("Unstandardized mean difference")
```

Here's how to standardize that unstandardized effect size into a Cohen's-$d$ metric.

```{r, fig.width = 3, fig.height = 2.5}
post %>% 
  mutate(d = b_Virgin1 / s_y) %>% 
  
  ggplot(aes(x = d, y = 0)) +
  stat_halfeyeh(point_interval = mode_hdi, .width = .95,
                fill = "grey67", slab_color = "grey92") +
  geom_vline(xintercept = (y_bar_a - y_bar_b) / s_y,
             color = "grey33", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(expression(paste("Cohen's ", italic(d), " expressed as a posterior")))
```

Let's work this one more way. By simple algebra, a standardized mean difference is the same as the difference between two standardized means. If we standardize the criterion `Longevity` before fitting the model and continue using the dummy variable approach, the `Virgin1` posterior will be the same as a Cohen's $d$.

Standardize the criterion.

```{r}
my_data <-
  my_data %>% 
  mutate(Longevity_s = (Longevity - mean(Longevity)) / sd(Longevity))
```

Because our criterion in a standardized metric, we no longer need our `stanvars`.

```{r fit19.10}
fit19.10 <-
  brm(data = my_data,
      family = gaussian,
      Longevity_s ~ 1 + Virgin1,
      prior = c(prior(normal(0, 1 * 5), class = Intercept),
                prior(normal(0, 1 * 5), class = b),
                prior(cauchy(0, 1), class = sigma)),
      iter = 3000, warmup = 1000, chains = 4, cores = 4,
      seed = 19,
      file = "fits/fit19.10")
```

Behold our out-of-the-box Bayesian Cohen's $d$.

```{r, fig.width = 3, fig.height = 2.5}
# no transformation necessary
posterior_samples(fit19.10) %>% 
  
  ggplot(aes(x = b_Virgin1, y = 0)) +
  stat_halfeyeh(point_interval = mode_hdi, .width = .95,
                fill = "grey67", slab_color = "grey92") +
  geom_vline(xintercept = (y_bar_a - y_bar_b) / s_y,
             color = "grey33", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(expression(paste("Cohen's ", italic(d), " expressed as a posterior")))
```

If you work them through, the approaches we took for `fit19.8` and `fit19.9` can be generalized to models with more than two groups. You just need to be careful how to compute the $s_y$ for each comparison. The same cannot be said for our `fit19.10` approach. If you were to, say, fit a model that produced estimates for all five groups, you could not just standardize the data beforehand using the overall standard deviation. If you did, the pooled standard deviation you'd be basing your not-quite-Cohen's $d$ posteriors on would be pooled across all groups, not just the two groups entailed in a given two-group comparison.

It's also the case the that standardized mean differences we computed for `fit19.1`, above, are not quite Cohen's $d$ effect sizes in the same way these have been. This is because the hierarchical approach we used partially pooled the estimates for each group toward the grand mean. You might say they were hierarchically-regularized Cohen's $d$s. But then again, Cohen's formula for his $d$ statistic did not account for Bayesian priors, either. So perhaps a purist would deny that any of the standardized mean differences we've computed in this chapter were proper Cohen's $d$ effect sizes. To be on the safe side, tell your readers exactly how you computed your models and what formulas you used to compute your effect sizes.

### Populations and samples.

You may have noticed that in our equation for $d$, above, we defined it in terms of sample statistics. Here it is again:

$$d = \frac{\bar y_A - \bar y_B}{s_y}$$

Sometimes we speak of the population parameter $\delta$, which is correspondingly defined as

$$\delta = \frac{\mu_A - \mu_B}{\sigma},$$

where $\mu_A$ and $\mu_B$ are the population means for the two groups under consideration and $\sigma$ is the pooled standard deviation in the population. Often times we don't have access to these values, which is why we run experiments and fit statistical models. But sometimes we do have access to the population parameters. In those cases, we can just plug them into the formula rather than estimate them in our models or with our sample statistics. 

Back in section 16.1.2, we saw an example of this with the `TwoGroupIQ` data. Let's load them, again.

```{r, message = F}
my_data <- read_csv("data.R/TwoGroupIQ.csv")

glimpse(my_data)
```

The data are IQ scores for participants in two groups. They look like this.

```{r, fig.width = 5, fig.height = 2.75}
my_data %>% 
  ggplot(aes(x = Score, Group)) +
  geom_halfeyeh(point_interval = mode_hdi, .width = c(.95, .5)) +
  xlab("IQ score") +
  coord_cartesian(ylim = c(1.5, 2.25))
```

If we wanted to compute the point estimate for Cohen's $d$ by hand, we might extract the necessary sample statistics and go.

```{r}
# save the sample means for the groups
y_bar_a <- filter(my_data, Group == "Smart Drug") %>% summarise(mean = mean(Score)) %>% pull()
y_bar_b <- filter(my_data, Group == "Placebo")    %>% summarise(mean = mean(Score)) %>% pull()

# save the pooled standard deviation in the sample
s_y <- sd(my_data$Score)

# compute Cohen's d
(y_bar_a - y_bar_b) / s_y
```

To get full posteriors for $d$, you might use any of the three methods we practiced with, above. Because we're only working with sample statistics, these are all only approximations of $\delta$. The thing about IQ scores is we actually know the population mean and standard deviation for IQ. They are 100 and 15, respectively. We know this because the people who make IQ tests design them that way. Let's see how well out sample statistics approximate the population parameters.

```{r}
my_data %>% 
  group_by(Group) %>% 
  summarise(mean = mean(Score),
            sd = sd(Score))
```

Unsurprisingly, the values for the `Smart Drug` group are notably different from the population parameters. But notice how close the values from the `Placebo` group are to the population parameters. If they weren't, we’d be concerned the `Placebo` condition was not a valid control. Looks like it was.

However, notice that the mean and standard deviation for the `Placebo` group are not the exact values of 100 and 15 the way they are in the population. If we wanted to compute a standardized mean difference between our `Smart Drug` group and folks in the population, we could just plug those values directly into our effect size equation. Here's what that would look like if we plug in the population mean for the control group.

```{r}
(y_bar_a - 100) / s_y
```

The result is very close to the one above. But this time our equation for $d$ was

$$d = \frac{\bar y_A - \mu_B}{s_y},$$

where we used the population mean $\mu_B$, but the other two terms were based on values from the sample. As long as you are defining the `Placebo` control as a stand-in for the population, this is a more precice way to compute $d$. Going further, we can also replace our $s_y$ with $\sigma$ (i.e., 15). 

```{r}
(y_bar_a - 100) / 15
```

Now our estimate is quite different. *Why*? Recall that sample standard deviations for both groups were larger than 15. That resulted in a $s_y$ value that was larger than 15, which meant we divided the unstandardized difference by a larger denominator. That larger denominator returned a smaller product. 

We have more than simple algebra to contend with, though. Replacing $s_y$ with 15 actually changed the meaning of the effect size. We were no longer using a pooled standard deviation. Rather, we were using the unconditional population standard deviation. Again, as long as we're equating our `Placebo` group with the population, we might restate that 15 as $\sigma_B$. Sometimes you might see this expressed as $\sigma_C$, where $C$ indicates "Control" (i.e., the reference category). That changed our formula to this:

$$d = \frac{\bar y_A - \mu_B}{\sigma_B}$$

So yes, that's right. We don't have to use pooled standard deviations to compute our standardized mean differences. If we have a good reference category, we could just use the sample estimate or population parameter for that, instead. Confusingly, you might see all these variants referred to as Cohen's $d$ within the literature. As with all the other decisions you make with experimental design and data analysis, use carefully reasoning to decide on how you'd like to compute your effect sizes, report them transparently to your audience, and prepare yourself for a collegial debate with Reviewer #2.

## Session info {-}

```{r}
sessionInfo()
```

```{r, echo = F}
# remove our objects
rm(generate_data, n, grand_mean, d, densities, arrow, text, gamma_a_b_from_omega_sigma, sd_y, omega, sigma, s_r, my_data, mean_y, stanvars, fit19.1, post, n_draws, post_draws, fit19.2, mw, fit19.3, fit19.4, my_small_data, sd_thorax_c, fit19.5, r, f, differences, p1, p2, fit19.6, fit19.7, p3, y_bar_a, y_bar_b, s_y, fit19.8, lines, fit19.9, fit19.10)
```

```{r, echo = F, message = F, warning = F, results = "hide"}
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)

ggplot2::theme_set(ggplot2::theme_grey())

bayesplot::color_scheme_set("blue")
```