22.Rmd

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

# Nominal Predicted Variable

> This chapter considers data structures that have a nominal predicted variable. When the nominal predicted variable has only two possible values, this reduces to the case of the dichotomous predicted variable considered in the previous chapter. In the present chapter, we generalize to cases in which the predicted variable has three or more categorical values...
>
> The traditional treatment of this sort of data structure is called multinomial logistic regression or conditional logistic regression. We will consider Bayesian approaches to these methods. As usual, in Bayesian software it is easy to generalize the traditional models so they are robust to outliers, allow different variances within levels of a nominal predictor, and have hierarchical structure to share information across levels or factors as appropriate. [@kruschkeDoingBayesianData2015, p. 649]

## Softmax regression

"The key descriptor of the [models in this chapter is their] inverse-link function, which is the softmax function (which will be defined below). Therefore, [Kruschke] refer[ed] to the method as softmax regression instead of multinomial logistic regression" (p. 650)

Say we have a metric predictor $x$ and a multinomial criterion $y$ with $k$ categories. We can express the basic linear model as

$$\lambda_k = \beta_{0, k} + \beta_{1, k} x,$$

for which the subscripts $k$ indicate there's a linear model for each of the $k$ categories. We call the possible set of $k$ outcomes $S$. Taking the case where $k = 3$, we'd have

\begin{align*}
\lambda_{[1]} & = \beta_{0, [1]} + \beta_{1, [1]} x, \\
\lambda_{[2]} & = \beta_{0, [2]} + \beta_{1, [2]} x, \text{and} \\
\lambda_{[3]} & = \beta_{0, [3]} + \beta_{1, [3]} x.
\end{align*}

In this scenario, what we want to know is the probability of $\lambda_{[1]}$, $\lambda_{[2]}$, and $\lambda_{[3]}$. The probability of a given outcome $k$ follows the formula

$$\phi_k = \operatorname{softmax}_S (\{\lambda_k\}) = \frac{\exp (\lambda_k)}{\sum_{c \in S} \exp  (\lambda_c)}.$$

> In words, [the equation] says that the probability of outcome $k$ is the exponentiated linear propensity of outcome $k$ relative to the sum of exponentiated linear propensities across all outcomes in the set $S$. You may be wondering, Why exponentiate? Intuitively, we have to go from propensities that can have negative values to probabilities that can only have non-negative values, and we have to preserve order. The exponential function satisfies that need. (p. 650)

You may be wondering what happened to $y$ and where all those $\lambda$s came from. Here we're using $\lambda$ to describe the propensity of outcome $k$, as indexed within our criterion $y$. So, the output of these models, $\phi_k$, is the relative probability we'll see each of our $k$ categories within our criterion $y$. What we want is $\phi_k$. The way we parameterize that with the softmax function is with $\lambda_k$.

There are are indeterminacies in the system of equations Kruschke covered in this section, the upshot of which is we'll end up making one of the $k$ categories the reference category, which we term $r$. Continuing on with our univariable model, we choose convenient constants for our parameters for $r$: $\beta_{0, r} = 0$ and $\beta_{1, r} = 0$. As such, *the regression coefficients for the remaining categories are relative to those for* $r$.

Kruschke saved the data for Figure 22.1 in the `SoftmaxRegData1.csv` and `SoftmaxRegData2.csv` files.

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

d1 <- read_csv("data.R/SoftmaxRegData1.csv")
d2 <- read_csv("data.R/SoftmaxRegData2.csv")

glimpse(d1)
glimpse(d2)
```

Before we explore these data in a plot, let's talk color and theme. For this chapter, we'll carry forward our practice from [Chapter 21][Dichotomous Predicted Variable] and take our color palette from the [**PNWColors** package](https://CRAN.R-project.org/package=PNWColors). This time, our color palette will be `"Lake"`.

```{r, warning = F, message = F, fig.width = 4, fig.height = 1}
library(PNWColors)

pl <- pnw_palette(name = "Lake")

pl
```

We'll base our overall plot theme on `cowplot::theme_minimal_grid()`, with many color adjustments from `PNWColors::pnw_palette(name = "Lake")`.

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

theme_set(
  theme_minimal_grid() +
    theme(text = element_text(color = pl[1]),
          axis.text = element_text(color = pl[1]),
          axis.ticks = element_line(color = pl[1]),
          legend.background = element_blank(),
          legend.box.background = element_blank(),
          legend.key = element_rect(fill = pl[8]),
          panel.background = element_rect(fill = pl[8], color = pl[8]),
          panel.grid = element_blank(),
          strip.background = element_rect(fill = pl[7], color = pl[7]),
          strip.text = element_text(color = pl[1]))
)
```

Now bind the two data frames together and plot our version of Figure 22.1.

```{r, fig.width = 5.5}
bind_rows(d1, d2) %>%
  mutate(data = rep(str_c("d", 1:2), each = n() / 2),
         Y = factor(Y)) %>% 
  
  ggplot(aes(x = X1, y = X2, label = Y, color = Y)) +
  geom_hline(yintercept = 0, color = "white") +
  geom_vline(xintercept = 0, color = "white") +
  geom_text(size = 3) +
  scale_color_manual(values = pl[2:5]) +
  labs(x = expression(x[1]),
       y = expression(x[2])) +
  coord_equal() +
  theme(legend.position = "none") +
  facet_wrap(~ data, ncol = 2)
```

### Softmax reduces to logistic for two outcomes.

"When there are only two outcomes, the softmax formulation reduces to the logistic regression of Chapter 21" (p. 653).

### Independence from irrelevant attributes.

> An important property of the softmax function of Equation 22.2 is known as independence from irrelevant attributes [@luceIndividualChoiceBehavior2012; @luceLuceChoiceAxiom2008]. The model implies that the ratio of probabilities of two outcomes is the same regardless of what other possible outcomes are included in the set. Let $S$ denote the set of possible outcomes. Then, from the definition of the softmax function, the ratio of outcomes $j$ and $k$ is
>
> $$\frac{\phi_j}{\phi_k} = \frac{\exp (\lambda_j) / \sum_{c \in S} \exp (\lambda_c)}{\exp (\lambda_k) / \sum_{c \in S} \exp (\lambda_c)}$$
>
> The summation in the denominators cancels and has no effect on the ratio of probabilities. Obviously if we changed the set of outcomes $S$ to any other set $S^*$ that still contains outcomes $j$ and $k$, the summation $\sum_{c \in S^*}$ would still cancel and have no effect on the ratio of probabilities. (p. 654)

Just to walk out that denominators-canceling business a little further,

\begin{align*}
\frac{\phi_j}{\phi_k} & = \frac{\exp (\lambda_j) / \sum_{c \in S} \exp (\lambda_c)}{\exp (\lambda_k) / \sum_{c \in S} \exp (\lambda_c)} \\
& = \frac{\exp (\lambda_j)}{\exp (\lambda_k)}.
\end{align*}

Thus even in the case of a very different set of possible outcomes $S^\text{very different}$, it remains that $\frac{\phi_j}{\phi_k} = \frac{\exp (\lambda_j)}{\exp (\lambda_k)}$.

Getting more applied, here's a tibble presentation of Kruschke's commute example with three modes of transportation.

```{r}
tibble(mode       = c("walking", "bicycling", "bussing"),
       preference = 3:1) %>% 
  mutate(`chance %` = (100 * preference / sum(preference)) %>% round(digits = 1))
```

Sticking with the example, if we take bicycling out of the picture, the `preference` values remain, but the `chance %` values change.

```{r}
tibble(mode       = c("walking", "bussing"),
       preference = c(3, 1)) %>% 
  mutate(`chance %` = 100 * preference / sum(preference))
```

Though we retain the same walking/bussing ratio, we end up with a different model of relative probabilities.

## Conditional logistic regression

> Softmax regression conceives of each outcome as an independent change in log odds from the reference outcome, and a special case of that is dichotomous logistic regression. But we can generalize logistic regression another way, which may better capture some patterns of data. The idea of this generalization is that we divide the set of outcomes into a hierarchy of two-set divisions, and use a logistic to describe the probability of each branch of the two-set divisions. (p. 655)

The model follows the generic equation

\begin{align*}
\phi_{S^* | S} = \operatorname{logistic}(\lambda_{S^* | S}) \\
\lambda_{S^* | S} = \beta_{0, S^* | S} + \beta_{1, {S^* | S}} x,
\end{align*}

where the conditional response probability (i.e., the goal of the analysis) is $\phi_{S^* | S}$. $S^*$ and $S$ denote the subset of outcomes and larger set of outcomes, respectively, and $\lambda_{S^* | S}$ is the propensity based on some linear model.

Kruschke saved the data for Figure 22.3 in the `CondLogistRegData1.csv` and `CondLogistRegData2.csv` files.

```{r, warning = F, message = F}
d3 <- read_csv("data.R/CondLogistRegData1.csv")
d4 <- read_csv("data.R/CondLogistRegData2.csv")

glimpse(d3)
glimpse(d4)
```

Let's bind the two data frames together and plot in bulk.

```{r, fig.width = 5}
bind_rows(d3, d4) %>%
  mutate(data = rep(str_c("d", 3:4), each = n() / 2),
         Y = factor(Y)) %>% 
  
  ggplot(aes(x = X1, y = X2, label = Y, color = Y)) +
  geom_hline(yintercept = 0, color = "white") +
  geom_vline(xintercept = 0, color = "white") +
  geom_text(size = 3) +
  scale_color_manual(values = pl[2:5]) +
  coord_equal() +
  labs(x = expression(x[1]),
       y = expression(x[2])) +
  theme(legend.position = "none") +
  facet_wrap(~ data, ncol = 2)
```

## Implementation in ~~JAGS~~ brms

### Softmax model.

Kruschke pointed out in his Figure 22.4 and the surrounding prose that we speak of the *categorical distribution* when fitting softmax models. Our **brms** paradigm will be much the same. To fit a softmax model with the `brm()` function, you specify `family = categorical`. The default is to use the logit link. In his [-@Bürkner2021Parameterization] [*Parameterization of response distributions in brms*](https://CRAN.R-project.org/package=brms/vignettes/brms_families.html#ordinal-and-categorical-models) vignette, Bürkner clarified:

> The **categorical** family is currently only implemented with the multivariate logit link function and has density
>
> $$f(y) = \mu_y = \frac{\exp (\eta_y)}{\sum_{k = 1}^K \exp (\eta_k)}$$
>
> Note that $\eta$ does also depend on the category $k$. For reasons of identifiability, $\eta_1$ is set to $0$.

Though there's no explicit softmax talk in that vignette, you can find it documented in his code [here](https://github.com/paul-buerkner/brms/blob/bc550ff3a2d41656a6711737faf1049207657800/R/distributions.R), starting in line 1891.

Now onto our **ggplot2** + **patchwork** version of the model diagram in Figure 22.4. I'm not gonna lie. The requisite code is a slog. We'll take the task in bits. First, we make and save the elements for the diagram on the left.

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

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

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

## an annotated arrow
# save our custom arrow settings
my_arrow <- arrow(angle = 20, length = unit(0.35, "cm"), type = "closed")
p3 <-
  tibble(x    = .5,
         y    = 1,
         xend = .73,
         yend = 0) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = pl[1]) +
  annotate(geom = "text",
           x = c(.48, .72), y = .5,
           label = c("'~'", "italic(k)"),
           size = c(10, 7), color = pl[1], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

## another annotated arrow
p4 <-
  tibble(x    = .5,
         y    = 1,
         xend = .4,
         yend = 0) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = pl[1]) +
  annotate(geom = "text",
           x = c(.34, .6), y = .5,
           label = c("'~'", "italic(jk)"),
           size = c(10, 7), color = pl[1], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

# likelihood formula
p5 <-
  tibble(x = .5,
         y = .5,
         label = "softmax(beta[0]['['*italic(k)*']']+sum()[italic(j)]~beta[italic(j)]['['*italic(k)*']']~italic(x)[italic(ji)])") %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 7, color = pl[1], parse = T, family = "Times") +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) +
  ylim(0, 1) +
  theme_void()

# a third annotated arrow
p6 <-
  tibble(x     = c(.375, .6),
         y     = c(1/2, 1/2),
         label = c("'='", "italic(i)")) %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = c(10, 7), color = pl[1], parse = T, family = "Times") +
  geom_segment(x = .5, xend = .5,
               y = 1, yend = 0,
               arrow = my_arrow, color = pl[1]) +
  xlim(0, 1) +
  theme_void()

# bar plot of categorical data
p7 <-
  tibble(x = 0:3,
         d = c(.5, .85, .5, .85)) %>% 
  ggplot(aes(x = x, y = d)) +
  geom_col(fill = pl[6], color = pl[5], width = .45) +
  annotate(geom = "text",
           x = 1.5, y = .2,
           label = "categorical",
           size = 7, color = pl[1]) +
  annotate(geom = "text",
           x = 1.25, y = .9, hjust = 0,
           label = "mu[italic(i)*'['*italic(k)*']']",
           size = 7, color = pl[1], family = "Times", parse = TRUE) +
  coord_cartesian(xlim = c(-.5, 3.5),
                  ylim = 0:1) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = pl[1]),
        plot.background = element_rect(fill = pl[8], color = "white", size = 1))

# the final annotated arrow
p8 <-
  tibble(x     = c(.375, .625),
         y     = c(1/3, 1/3),
         label = c("'~'", "italic(i)")) %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = c(10, 7), color = pl[1], parse = T, family = "Times") +
  geom_segment(x = .5, xend = .5,
               y = 1, yend = 0, 
               color = pl[1], arrow = my_arrow) +
  xlim(0, 1) +
  theme_void()

# some text
p9 <-
  tibble(x     = 1,
         y     = .5,
         label = "italic(y[i])") %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 7, color = pl[1], parse = T, family = "Times") +
  xlim(0, 2) +
  theme_void()

# define the layout
layout <- c(
  area(t = 1, b = 2, l = 1, r = 2),
  area(t = 1, b = 2, l = 3, r = 4),
  area(t = 3, b = 3, l = 1, r = 2),
  area(t = 3, b = 3, l = 3, r = 4),
  area(t = 4, b = 4, l = 1, r = 4),
  area(t = 5, b = 5, l = 2, r = 3),
  area(t = 6, b = 7, l = 2, r = 3),
  area(t = 8, b = 8, l = 2, r = 3),
  area(t = 9, b = 9, l = 2, r = 3)
)

# combine and plot!
a <-
  (
    (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) + 
      plot_layout(design = layout) &
      ylim(0, 1) &
      theme(plot.margin = margin(0, 5.5, 0, 5.5))
  )
```

Now we make and save the elements for the diagram on the right.

```{r, warning = F, message = F}
# third normal density
p2 <-
  tibble(x = seq(from = -3, to = 3, by = .1)) %>% 
  ggplot(aes(x = x, y = (dnorm(x)) / max(dnorm(x)))) +
  geom_area(fill = pl[6], color = pl[5]) +
  annotate(geom = "text",
           x = 0, y = .2,
           label = "normal",
           size = 7, color = pl[1]) +
  annotate(geom = "text",
           x = c(0, 1.5), y = .6,
           label = c("italic(M[j])", "italic(S[j])"), 
           size = 7, color = pl[1], hjust = 0, family = "Times", parse = T) +
  scale_x_continuous(expand = c(0, 0)) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = pl[1]),
        plot.background = element_rect(fill = pl[8], color = "white", size = 1))

## an annotated arrow
p3 <-
  tibble(x    = .5,
         y    = 1,
         xend = .85,
         yend = 0) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = pl[1]) +
  annotate(geom = "text",
           x = .49, y = .5,
           label = "'~'",
           size = 10, color = pl[1], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

## another annotated arrow
p4 <-
  tibble(x    = .5,
         y    = 1,
         xend = .4,
         yend = 0) %>%
  
  ggplot(aes(x = x, xend = xend,
             y = y, yend = yend)) +
  geom_segment(arrow = my_arrow, color = pl[1]) +
  annotate(geom = "text",
           x = c(.35, .57), y = .5,
           label = c("'~'", "italic(j)"),
           size = c(10, 7), color = pl[1], family = "Times", parse = T) +
  xlim(0, 1) +
  theme_void()

# likelihood formula
p5 <-
  tibble(x = .5,
         y = .5,
         label = "logistic(beta[0]+sum()[italic(j)]~beta[italic(j)]~italic(x)[italic(ji)])") %>% 
  
  ggplot(aes(x = x, y = y, label = label)) +
  geom_text(size = 7, color = pl[1], parse = T, family = "Times") +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 1)) +
  ylim(0, 1) +
  theme_void()

# bar plot of Bernoulli data
p7 <-
  tibble(x = 0:1,
         d = (dbinom(x, size = 1, prob = .6)) / max(dbinom(x, size = 1, prob = .6))) %>% 
  
  ggplot(aes(x = x, y = d)) +
  geom_col(fill = pl[6], color = pl[5], width = .4) +
  annotate(geom = "text",
           x = .5, y = .2,
           label = "Bernoulli",
           size = 7, color = pl[1]) +
  annotate(geom = "text",
           x = .5, y = .9,
           label = "mu[italic(i)]", 
           size = 7, color = pl[1], family = "Times", parse = T) +
  xlim(-.75, 1.75) +
  theme_void() +
  theme(axis.line.x = element_line(size = 0.5, color = pl[1]),
        plot.background = element_rect(fill = pl[8], color = "white", size = 1))

# combine and plot!
c <-
  (
    (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9) + 
      plot_layout(design = layout) &
      ylim(0, 1) &
      theme(plot.margin = margin(0, 5.5, 0, 5.5))
  )
```

Here we combine the two model diagrams and plot!

```{r, fig.width = 8, fig.height = 5, warning = F, message = F}
b <- plot_spacer()

(a | b | c) + plot_layout(widths = c(4, 1, 4))
```

### Conditional logistic model.

The conditional logistic regression models are not natively supported in **brms** at this time. However, if you follow [issue #560](https://github.com/paul-buerkner/brms/issues/560), you'll see there are ways to fit them using the nonlinear syntax. If you compare the syntax Bürkner used in that thread on January 30^th^ to the JAGS syntax Kruschke showed on pages 661 and 662, you'll see they appear to follow contrasting parameterizations.

I think that's about as far as I'm going with this model type at this time. If you work through the solution, please share your code in my [GitHub issue #22](https://github.com/ASKurz/Doing-Bayesian-Data-Analysis-in-brms-and-the-tidyverse/issues/22).

### Results: Interpreting the regression coefficients.

#### Softmax model.

Load **brms**.

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

Along with Kruschke, we'll be modeling the `d1` data. In case it's not clear, the `X1` and `X2` variables are already in a standardized metric.

```{r, message = F}
d1 %>% 
  pivot_longer(-Y) %>% 
  group_by(name) %>% 
  summarise(mean = mean(value), 
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 2)
```

This will make it easier to set the priors. Here we'll just use the rather wide priors Kruschke indicated on page 662. Note our use of the `dpar` argument within the `prior()` function.

```{r fit22.1}
fit22.1 <-
  brm(data = d1, 
      family = categorical(link = logit),
      Y ~ 0 + Intercept + X1 + X2,
      prior = c(prior(normal(0, 20), class = b, dpar = mu2),
                prior(normal(0, 20), class = b, dpar = mu3),
                prior(normal(0, 20), class = b, dpar = mu4)),
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      seed = 22,
      file = "fits/fit22.01")
```

Since it's the default, we didn't have to include the `(link = logit)` bit in the `family` argument. I'm just being explicit for the sake of pedagogy. Take a look at the parameter summary.

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

As indicated in the formulas, above, we get posteriors for each level of `Y`, except for `Y == 1`. That serves as the reference category. The values for $\beta_{i, k = 1}$ are all fixed at $0$.

Here's how we might make the histograms in Figure 22.5.

```{r, fig.width = 8, fig.height = 4.5, warning = F, message = F}
library(tidybayes)

# extract the posterior draws
post <- posterior_samples(fit22.1)

# wrangle
post %>% 
  pivot_longer(-lp__) %>% 
  mutate(name = str_remove(name, "b_")) %>% 
  mutate(lambda    = str_extract(name, "[2-4]+") %>% str_c("lambda==", .),
         parameter = if_else(str_detect(name, "Intercept"), "beta[0]",
                             if_else(str_detect(name, "X1"), "beta[1]", "beta[2]"))) %>% 
  
  # plot
  ggplot(aes(x = value, y = 0)) +
  stat_histinterval(point_interval = mode_hdi, .width = .95,
                    fill = pl[4], slab_color = pl[3], color = pl[1], point_color = pl[2],
                    normalize = "panels") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("marginal posterior") +
  facet_grid(lambda ~ parameter, labeller = label_parsed, scales = "free_x")
```

Because the $\beta$ values for when $\lambda = 1$ are all fixed to 0, we left those plots out of our version of the figure. If you really wanted them, you'd have to enter the corresponding cells into the data before plotting. If you summarize each parameter by it's posterior mean, `round()`, and wrangle a little, you can arrange the results in a similar way that the equations for $\lambda_2$ through $\lambda_4$ are displayed on the left side of Figure 22.5.

```{r, message = F}
post %>% 
  pivot_longer(-lp__) %>% 
  mutate(name = str_remove(name, "b_")) %>% 
  mutate(lambda    = str_extract(name, "[2-4]+") %>% str_c("lambda[", ., "]"),
         parameter = if_else(str_detect(name, "Intercept"), "beta[0]",
                             if_else(str_detect(name, "X1"), "beta[1]", "beta[2]"))) %>% 
  group_by(lambda, parameter) %>% 
  summarise(mean = mean(value) %>% round(digits = 1)) %>% 
  pivot_wider(names_from = parameter,
              values_from = mean)
```

As Kruschke mentioned in the text, "the estimated parameter values should be near the generating values, but not exactly the same because the data are merely a finite random sample" (pp. 662--663). Furthermore,

>  interpreting the parameters is always contextualized relative to the model. For the softmax model, the regression coefficient for outcome $k$ on predictor $x_j$ indicates that rate at which the log odds of that outcome increase relative to the reference outcome for a one unit increase in $x_j$, assuming that a softmax model is a reasonable description of the data. (p. 663)

Unfortunately, this makes the parameters difficult to interpret directly. Kruschke didn't show a plot like this, but it might be helpful to further understand what this model means in terms of probabilities for a given `y` value. Here we'll use the `fitted()` function to return the conditional probabilities for all four response options for `Y` based on various combinations of `X1` and `X2`.

```{r, fig.width = 6, fig.height = 5, warning = F, message = F}
nd <- crossing(X1 = seq(from = -2, to = 2, length.out = 50),
               X2 = seq(from = -2, to = 2, length.out = 50))

fitted(fit22.1,
       newdata = nd) %>% 
  as_tibble() %>% 
  select(contains("Estimate")) %>% 
  set_names(str_c("Y==", 1:4)) %>% 
  bind_cols(nd) %>% 
  pivot_longer(contains("Y"),
               values_to = "phi") %>% 
  
  ggplot(aes(x = X1, y = X2, fill = phi)) +
  geom_raster(interpolate = T) +
  scale_fill_gradientn(expression(phi[italic(k)*"|"*italic(S)]),
                       colours = pnw_palette(name = "Lake", n = 101),
                       limits = c(0, 1)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  facet_wrap(~ name, labeller = label_parsed)
```

Now use that plot while you walk through the final paragraph in this subsection.

> It is easy to transform the estimated parameter values to a different reference category. Recall from Equation 22.3 (p. 651) that arbitrary constants can be added to all the regression coefficients without changing the model prediction. Therefore, to change the parameters estimates so they are relative to outcome $R$, we simply subtract $\beta_{j, R}$ from $\beta_{j, k}$ for all predictors $j$ and all outcomes $k$. We do this at every step in the MCMC chain. For example, in Figure 22.5, consider the regression coefficient on $x_1$ for outcome 2. Relative to reference outcome 1, this coefficient is positive, meaning that the probability of outcome 2 increases relative to outcome 1 when $x_1$ increases. You can see this in the data graph, as the region of 2's falls to right side (positive $x_1$ direction) of the region of 1's. But if the reference outcome is changed to outcome 4, then the coefficient on $x_1$ for outcome 2 changes to a negative value. Algebraically this happens because the coefficient on $x_1$ for outcome 4 is larger than for outcome 2, so when the coefficient for outcome 4 is subtracted, the result is a negative value for the coefficient on outcome 2. Visually, you can see this in the data graph, as the region of 2's falls to the left side (negative $x_1$ direction) of the region of 4's. Thus, *interpreting regression coefficients in a softmax model is rather different than in linear regression. In linear regression, a positive regression coefficient implies that* $y$ *increases when the predictor increases. But not in softmax regression, where a positive regression coefficient is only positive with respect to a particular reference outcome*. (p. 664, *emphasis* added)

##### Bonus: Consider the interceps-only softmax model.

Models like `fit22.1`, above, are great when you want to explore predictors for your nominal variables. However, these models are also really useful when you aren't interested in predictor variables. In these cases, the intercepts-only model will help you compute high-quality Bayesian intervals around the category percentages. Let's walk through an example to see what I mean. Fit an intercepts-only version of the model, above.

```{r fit22.2}
fit22.2 <-
  brm(data = d1, 
      family = categorical(link = logit),
      Y ~ 1,
      prior = c(prior(normal(0, 20), class = Intercept, dpar = mu2),
                prior(normal(0, 20), class = Intercept, dpar = mu3),
                prior(normal(0, 20), class = Intercept, dpar = mu4)),
      iter = 2000, warmup = 1000, cores = 4, chains = 4,
      seed = 22,
      file = "fits/fit22.02")
```

Review the model summary.

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

Even with an intercepts-only model, the parameters in a softmax model are difficulty to interpret directly. Though you might be tempted to presume `mu2_Intercept` through `mu4_Intercept` were probabilities on the log-odds scale, they're not. This, recall, is due to their connection to the reference category. If we return to the equation from [Section 22.1][Softmax regression],

$$\phi_k = \operatorname{softmax}_S (\{\lambda_k\}) = \frac{\exp (\lambda_k)}{\sum_{c \in S} \exp  (\lambda_c)},$$

we can get a sense of how to convert these parameters to relative probabilities. Here's how to do so by hand with the posterior draws.

```{r}
posterior_samples(fit22.2) %>% 
  transmute(iter        = 1:n(),
            `lambda[1]` = 0,  # recall this is the default
            `lambda[2]` = b_mu2_Intercept,
            `lambda[3]` = b_mu3_Intercept,
            `lambda[4]` = b_mu4_Intercept) %>% 
  pivot_longer(-iter) %>% 
  # the next two rows are where the magic happens
  group_by(iter) %>% 
  mutate(phi = exp(value) / sum(exp(value))) %>%
  mutate(Y = str_extract(name, "\\d")) %>% 
  group_by(Y) %>% 
  mean_qi(phi) %>% 
  mutate_if(is.double, round, digits = 2)
```

We can compute these relative probability values ($\phi_k$) much easier with `fitted()`.

```{r}
f <- fitted(fit22.2)

f[1, , ] %>% 
  t() %>% 
  round(digits = 2)
```

Anyway, the reason you might want to go through the trouble of fitting an intercepts-only softmax model is to improve on the kinds of par plots people often report in their manuscripts. Consider these two:

```{r, fig.width = 8, fig.height = 4}
# descriptive statistics
p1 <-
  d1 %>% 
  count(Y) %>%  
  mutate(p = n / sum(n)) %>% 
  
  ggplot(aes(x = Y, y = p)) +
  geom_col(fill = pl[6]) +
  scale_y_continuous(NULL, labels = scales::percent, breaks = 0:4 / 10,
                     expand = expansion(mult = c(0, 0.05)), limits = c(0, .4)) +
  labs(subtitle = "sample statistics")

# population percentages
p2 <-
  f[1, , ] %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column("level") %>% 
  mutate(Y = str_extract(level, "\\d")) %>% 
  
  ggplot(aes(x = Y, y = Estimate)) +
  geom_col(fill = pl[6]) +
  geom_linerange(aes(ymin = Q2.5, ymax = Q97.5),
                 color = pl[2], size = 1) +
  scale_y_continuous(NULL, labels = scales::percent, breaks = 0:4 / 10,
                     expand = expansion(mult = c(0, 0.05)), limits = c(0, .4)) +
  labs(subtitle = "model-based population percentages\n(with 95% interval bars)")

p1 + p2 + plot_annotation(title = "The softmax model adds information to the conventional sample-based\nbar plot.")
```

The plot on the left is the kind of sample data summary you'll see in countless articles and data presentations. Though it's a great way to quickly summarize the relative percentages of each category, it does nothing to express how (un)certain we are those sample statistics will describe the population. The intercepts-only softmax model returns the posterior distributions for the population probabilities. In the plot on the right, the bars mark off the posterior means and the vertical lines mark off the 95% intervals. This is why we want the intercepts-only softmax model.

Okay, now just for kicks and giggles, I'd like to go on a plotting tangent. If data analysts broadly replaced the typical sample-based plots (on the left) for the model-based plots (on the right), it would be a great improvement. At a personal level, though, I think simple bar plots are over used. Let's explore four alternatives. In this block, we'll make and save the first three.

```{r}
# rotated bar plot
p1 <-
  f[1, , ] %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column("level") %>% 
  mutate(Y = str_extract(level, "\\d")) %>% 
  
  ggplot(aes(x = Estimate, y = Y)) +
  geom_col(fill = pl[6]) +
  geom_linerange(aes(xmin = Q2.5, xmax = Q97.5),
                 color = pl[2], size = 1) +
  scale_x_continuous(NULL, labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05))) +
  labs(subtitle = "rotated bar plot")

# coefficient plot
p2 <-
  f[1, , ] %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column("level") %>% 
  mutate(Y = str_extract(level, "\\d")) %>% 
  
  ggplot(aes(x = Y, y = Estimate)) +
  geom_pointrange(aes(ymin = Q2.5, ymax = Q97.5),
                 color = pl[2], size = 1, fatten = 2) +
  scale_y_continuous(NULL, labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05)), limits = c(0, NA)) +
  labs(subtitle = "coefficient plot")

# CCDF bar plots
p3 <-
  fitted(fit22.2, summary = F)[, 1, ] %>% 
  data.frame() %>% 
  set_names(1:4) %>% 
  pivot_longer(everything(), values_to = "p") %>% 
  mutate(Y = factor(name)) %>% 
  
  ggplot(aes(x = Y, y = p)) +
  stat_ccdfinterval(.width = .95, fill = pl[6], color = pl[2],
                    size = 1.5, point_size = 2) +
  scale_y_continuous(NULL, labels = scales::percent, breaks = 0:4 / 10,
                     expand = expansion(mult = c(0, 0.05)), limits = c(0, .4)) +
  labs(subtitle = "CCDF bar plot")
```

The fourth alternative is a little weird AND it's going to take a bit more work than the first three.

```{r}
# for annotation
text <-
 fitted(fit22.2)[1, , ] %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column("level") %>% 
  mutate(y       = (Estimate / 2) + lag(cumsum(Estimate), default = 0),
         summary = str_c(round(100 * Estimate, 1), "*'% ± '*", round(200 * Est.Error, 1)),
         label_y = str_c("Y = ", 1:4))

p4 <-
  # wrangle
  fitted(fit22.2, summary = F)[, 1, ] %>% 
  data.frame() %>% 
  set_names(1:4) %>% 
  mutate(row = 1:n()) %>% 
  pivot_longer(-row) %>% 
  mutate(Y = fct_rev(name)) %>% 

  # plot
  ggplot(aes(x = row, y = value)) +
  geom_col(aes(fill = Y),
           position = "stack", size = 0, width = 1) +
  geom_text(data = text,
            aes(x = -2000, y = y, label = label_y),
            color = pl[1]) +
  geom_text(data = text,
            aes(x = 6000, y = y, label = summary),
            color = pl[1], parse = T) +
  scale_fill_manual(values = pl[7:4], breaks = NULL) +
  scale_x_continuous(NULL, breaks = NULL, limits = c(-4000, 8000)) +
  scale_y_continuous(NULL, labels = scales::percent, expand = c(0, 0)) +
  labs(subtitle = "stacked bar plot with uncertain boundaries")
```

Okay, now combine the four ggplots and behold!

```{r, fig.width = 8, fig.height = 7.5}
# combine
p1 + p2 + p3 + p4 + 
  plot_annotation(title = "Alternatives to the conventional bar plot",
                  tag_levels = "a", tag_prefix = "(", tag_suffix = ")")
```

The rotated bar plot (a) is most useful when the names of the levels are longer character strings. For example, imagine that instead of `1` through `4`, the four levels were countries or pharmaceutical drugs. You wouldn't want to mess with formatting those on an $x$-axis. Use a rotated bar plot, instead. For my taste, the simple coefficient plot (b) gets the job done nicely without the unnecessary clutter of the bars. Matthew Kay's CCDF bar plots (c) provide a fuller expression of the shape of the posterior uncertainty in each percentage. If desired, you could even omit the dot intervals from those and the visual expression of uncertainty would still remain. The final stacked bar plot with the fuzzy boundaries (d) is the result of a little back-and-forth on twitter ([original tweet](https://twitter.com/SolomonKurz/status/1372632774285348864)). It's the oddball of the group, but what it does uniquely well is show how the percentages of the groups all depend on one another and must, be definition, sum to 100.

#### Conditional logistic model.

I'm not pursuing this model type at this time. If you work through the solution, please share your code in my [GitHub issue #22](https://github.com/ASKurz/Doing-Bayesian-Data-Analysis-in-brms-and-the-tidyverse/issues/22) or chime in on the Stan forum thread called [Nominal data and Kruschke's "conditional logistic" approach](https://discourse.mc-stan.org/t/nominal-data-and-kruschkes-conditional-logistic-approach/21433).

## Generalizations and variations of the models

These models can be generalized to include different kinds of predictors, variants robust to outliers, and model comparison via information criteria and so forth. You can find a few more examples with softmax regression in Chapter 10 of the first edition of McElreath's [-@mcelreathStatisticalRethinkingBayesian2015] [*Statistical rethinking*](https://xcelab.net/rm/statistical-rethinking/) and Chapter 11 of his second edition [@mcelreathStatisticalRethinkingBayesian2020]. See also Kurz [-@kurzStatisticalRethinkingSecondEd2021, [Section 11.3](https://bookdown.org/content/4857/god-spiked-the-integers.html#multinomial-and-categorical-models)] and Kurz [-@kurzStatisticalRethinkingBrms2020, [Section 10.3.1](https://bookdown.org/content/3890/counting-and-classification.html#multinomial)] for walkthroughs with **brms**.

## Session info {-}

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

```{r, message = F, echo = F}
# remove the objects
rm(d1, d2, pl, d3, d4, p1, p2, my_arrow, p3, p4, p5, p6, p7, p8, p9, layout, a, c, b, fit22.1, post, nd, fit22.2, f, text)
```

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