ch08_interactions.qmd

# Conditional Manatees {#interactions}

```{r}
#| include: false
library(dplyr)
library(tidyr)
library(tidybayes)
library(rethinking)
library(brms)
library(modelr)
library(skimr)
library(loo)
library(ggplot2)
library(dagitty)
library(ggdag)
library(ggdist)
library(patchwork)
library(paletteer)
```

Some options to facilitate the computations

```{r}
#  For execution on a local, multicore CPU with excess RAM
options(mc.cores = parallel::detectCores())
#  To avoid recompilation of unchanged Stan programs
rstan::rstan_options(auto_write = TRUE)
```

The default theme used by `ggplot2`

```{r}
# The default theme used by ggplot2
ggplot2::theme_set(ggthemes::theme_pander())
ggplot2::theme_update(title = element_text(color = "midnightblue"))
```

## Building an interaction

Load the data, log transform the gdp measure, remove incomplete cases and create a character column for Africa or Not Africa.

```{r ch08_dataRugged}
data(rugged)
dataRugged <- rugged |>
  filter(complete.cases(rgdppc_2000)) |>
  mutate(log_gdp = log(rgdppc_2000),
         is_africa = if_else(cont_africa == 1, "Africa", "Not Africa"),
         is_africa = as.factor(is_africa))
rm(rugged)
# NOTE: Make sure as.vector() is outside of scale().
#       Otherwise it keeps the vector as an array and causes all sort of little
#       problems.  In particular, a very obscure, fine error message
#       in doing brms fit for b8.2.
dataRugged_nona <- dataRugged |>
  drop_na(rgdppc_2000) |>
  mutate(log_gdp_s = log_gdp / mean(log_gdp),
         rugged_s = scales::rescale(rugged),
         rugged_sc = as.vector(scale(rugged_s, center = TRUE, scale = FALSE)))
dataRugged_nona |>
  select(log_gdp, log_gdp_s, rugged, rugged_s, rugged_sc) |>
  skim() |>
  select(-n_missing, -complete_rate) |>
  mutate(across(.cols = where(is.numeric), .fns = round, digits = 2))
```

and we use the following DAG, see overthinking box in introduction of section 8.1 for another possible DAG.

```{r dagRugged}
dagRugged <- list()
dagRugged <- within(dagRugged, {
  coords <- tibble(name = c("C", "G", "R", "U"),
                     x = c(3, 2, 1, 2),
                     y = c(2, 2, 2, 1))
  dag <- dagify(G ~ C + R + U,
                R ~ U,
                latent = "U",
                outcome = "G",
                coords = coords)
  p <- dag |>
    ggdag_status(aes(color = status), as_factor = TRUE, node_size = 14,
                 text_size = 4, text_col = "midnightblue") +
    scale_color_paletteer_d("khroma::light", 
                          na.value = "honeydew3",
                          direction = 1) +
    theme_dag() +
    theme(legend.position = c(0.8, 0.2)) +
    labs(title = "African nations", subtitle = "Section 8.1")
})
dagRugged$p

```

```{r plotRugged}
plotRugged <- list()
plotRugged <- within(plotRugged, {
  Africa <- dataRugged_nona |>
    filter(grepl("^africa$", x = is_africa, ignore.case = TRUE)) |>
    ggplot(aes(x = rugged_s, y = log_gdp_s)) +
    geom_smooth(method = "lm", formula = y ~ x, fill = "lightblue", color = "royalblue") +
    geom_point(color = "burlywood4") +
    labs(title = "African nations", x = "ruggedness (rescale)", 
         y = "log GDP (prop of mean)")
  
  notAfrica <- dataRugged_nona |>
    filter(!grepl("^africa$", x = is_africa, ignore.case = TRUE)) |>
    ggplot(aes(x = rugged_s, y = log_gdp_s)) +
    geom_smooth(method = "lm", formula = y ~ x, fill = "burlywood1", color = "burlywood4") +
    geom_point(color = "royalblue") +
    labs(title = "Non-African nations", x = "ruggedness (rescale)",
         y = "log GDP (prop of mean)")
  title <- "Figure 8.2. Separate linear regressions inside and outside of Africa"
})
wrap_plots(plotRugged[c("Africa", "notAfrica")]) + 
  plot_annotation(title = plotRugged$title)
```

### Making a rugged model

and split the data into countries from Africa and not.

```{r}
# lst <- d |>
#   split(d$is_africa)
# str(lst)
```

and now creating a simple univariate model

$$
\begin{align*}
\log{(log\_gdp\_s_i)} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= \alpha + \beta \cdot rugged\_sc_i \\
\alpha &\sim \mathcal{N}(1, 1) \\
\beta &\sim \mathcal{N}(0, 1) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

Now fit the model. Get the **prior samples** by using `sample_prior = TRUE`.

```{r ch08_fit08_01a}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
fit08_01a <- xfun::cache_rds({
  out <- brm(
    data = dataRugged_nona,
    family = gaussian,
    log_gdp_s ~ 1 + rugged_sc,
    prior = c(
      prior(normal(1, 1), class = Intercept),
      prior(normal(0, 1), class = b),
      prior(exponential(1), class = sigma)),
    sample_prior = TRUE,
    iter = 1000, warmup = 500, chains = 4, cores = detectCores(),
    seed = 809)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_01a")
tictoc::toc()
```

```{r ch08_b08_01a}
# tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
# b8.1a <- xfun::cache_rds({
#   out <- brm(
#     data = dd,
#     family = gaussian,
#     log_gdp_s ~ 1 + rugged_sc,
#     prior = c(
#       prior(normal(1, 1), class = Intercept),
#       prior(normal(0, 1), class = b),
#       prior(exponential(1), class = sigma)),
#     sample_prior = TRUE,
#     iter = 2000, warmup = 1000, chains = 4, cores = detectCores(),
#     seed = 8)
#   out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
#   file = "ch08_b08_01a")
# tictoc::toc()
```

```{r}
posterior_summary(fit08_01a) |>
  round(digits = 3)
```

the estimates are described in section 8.1.1 of McElreath but he seems to have

```{r}
prior08_01a <- list()
prior08_01a <- within(prior08_01a, {
  draws <- prior_draws(fit08_01a)
  df <- draws |>
    slice_sample(n = 50) |>
    tibble::rownames_to_column(var = "id") |>
    expand(nesting(id, Intercept, b), rugged_sc = c(-2, 2)) |>
    mutate(log_gdp_s = Intercept + b * rugged_sc,
         rugged_s  = rugged_sc + mean(dataRugged_nona$rugged_s))
  est <- c("fixed" = min(dataRugged_nona$log_gdp_s), "b" = diff(range(dataRugged_nona$log_gdp_s)))
  
  p <- ggplot(df, aes(x = rugged_s, y = log_gdp_s, group = id)) +
    geom_line(color = "orchid") +
    geom_hline(yintercept = range(dataRugged_nona$log_gdp_s), 
               size = 1, linetype = 2, color = "royalblue") +
    geom_abline(intercept = est["fixed"], slope = est["b"], color = "purple", size = 1) +
    coord_cartesian(xlim = c(0, 1), ylim = c(0.5, 1.5)) +
    labs(
      subtitle = "Intercept ~ dnorm(1, 1)\nb ~ dnorm(0, 1)",
      x = "ruggedness (rescaled)",
      y = "log GDP (prop of mean)")
})
# prior08_01a$p
```

Now using the prior where we want the intercept to be around 1 with extremes from 0.8 to 1.2 (i.e. a mean of 1 and sd of 0.1) and the slope to have extremes about $\pm 0.6$, that is a mean of 0 with sd of 0.3 (i.e. 2 sd with sd = 3 from a mean of 0).

```{r ch08_fit08_01b}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
fit08_01b <- xfun::cache_rds({
  out <- update(
    fit08_01a,
    newdata = dataRugged_nona,
    prior = c(
      prior(normal(1, 0.1), class = Intercept),
      prior(normal(0, 0.3), class = b),
      prior(exponential(1), class = sigma)),
    sample_prior = TRUE,
    seed = 809)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_01b")
tictoc::toc()
```

```{r}
posterior_summary(fit08_01b) |>
  round(digits = 3)
```

```{r}
prior08_01b <- list()
prior08_01b <- within(prior08_01b, {
  draws <- prior_draws(fit08_01b) |>
    mutate(.draw = seq_len(n())) |>
    relocate(.draw)
  df <- draws |>
    slice_sample(n = 50) |>
    expand(nesting(.draw, Intercept, b), rugged_sc = c(-2, 2)) |>
    mutate(log_gdp_s = Intercept + b * rugged_sc,
         rugged_s  = rugged_sc + mean(dataRugged_nona$rugged_s))
  est <- c("fixed" = min(dataRugged_nona$log_gdp_s), "b" = diff(range(dataRugged_nona$log_gdp_s)))
  
  p <- ggplot(df, aes(x = rugged_s, y = log_gdp_s, group = .draw)) +
    geom_line(color = "orchid") +
    geom_hline(yintercept = range(dataRugged_nona$log_gdp_s), 
               size = 1, linetype = 2, color = "royalblue") +
    geom_abline(intercept = est["fixed"], slope = est["b"], color = "purple", size = 1) +
    coord_cartesian(xlim = c(0, 1), ylim = c(0.5, 1.5)) +
    labs(
      subtitle = "Intercept ~ dnorm(1, 0.1)\nb ~ dnorm(0, 0.3)",
      x = "ruggedness (rescaled)",
      y = "log GDP (prop of mean)")
})
# prior08_01b$p
```

```{r}
prior08_01a$p + prior08_01b$p + 
  plot_annotation(
    title = "Figure 8.3. Simulating different priors to evaluate their fit")
```

### Adding an indicator variable isn't enough

We add the `cid` variable to identify the continent.

```{r}
dataRugged_nona <- dataRugged_nona |>
  mutate(cid = as.factor(if_else(cont_africa == 1, "1", "2")))
```

and fitting the data to the following model

$$
\begin{align*}
\log{(log\_gdp\_s_i)} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= \alpha[cid] + \beta \cdot rugged\_sc_i \\
\alpha &\sim \mathcal{N}(1, 0.1) \\
\beta &\sim \mathcal{N}(0, 0.3) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

```{r ch08_fit08_02}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit08_02 <- xfun::cache_rds({
  out <- brm(
    data = dataRugged_nona,
    family = gaussian,
    log_gdp_s ~ 0 + cid + rugged_sc,
    prior = c(
      prior(normal(1, 0.1), class = b, coef = cid1),
      prior(normal(1, 0.1), class = b, coef = cid2),
      prior(normal(0, 0.3), class = b, coef = rugged_sc),
      prior(exponential(1), class = sigma)),
    sample_prior = TRUE,
    iter = 1000, warmup = 500, chains = 4, cores = detectCores(),
    seed = 811)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_02")
tictoc::toc()
```

```{r}
posterior_summary(fit08_02) |>
  round(digits = 3)
```

and measuring the models' performance

```{r}
loo::loo_compare(fit08_01b, fit08_02, criterion = "waic") |> 
  print(simplify = FALSE)
```

with thew model weights

```{r}
brms::model_weights(fit08_01b, fit08_02) |>
  round(digits = 2)
```

and create the fitted data used for the plot

```{r}
tidybayes::get_variables(fit08_02)
```

```{r ch08_lpred08_02}
# get the fitted values
lpred08_02 <- crossing(
  cid = as.factor(1:2),
  rugged_sc = seq(from = -0.2, to = 1.2, length.out = 30)) |>
  mutate(rugged_sc = as.vector(scale(rugged_sc))) |>
  add_linpred_draws(fit08_02, ndraws = 50) |>
  mean_qi() |>
  mutate(is_africa = if_else(cid == 1, "Africa", "Not Africa")) |>
  mutate(is_africa = as.factor(is_africa))
# glimpse(lpred08_02)
```

```{r ch08_plot08_02}
ggplot(dataRugged_nona, aes(x = rugged_sc, y = log_gdp_s, fill = is_africa, color = is_africa)) +
  geom_smooth(data = lpred08_02, aes(x = rugged_sc, y = .linpred, ymin = .lower, ymax = .upper),
              stat = "identity",
              alpha = 1/4, size = 1/2) +
  geom_point(size = 1) +
  scale_fill_manual(values = c("Africa" = "springgreen2", "Not Africa" = "violet")) +
  scale_color_manual(values = c("Africa" = "springgreen2", "Not Africa" = "violet")) +
  coord_cartesian(xlim = c(0, 1)) +
  theme(legend.position = c(.80, .90),
        legend.title = element_blank()) +
  labs(title = "Figure 8.4",
       subtitle = "model b8.2",
       x = "ruggedness (standardized)",
       y = "log GDP (as proportion of mean)")
```

### Adding an interaction does work

$$
\begin{align*}
\log{(log\_gdp\_s_i)} &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= \alpha_{[cid]} + \beta_{[cid]} \cdot rugged\_sc_i \\
\alpha &\sim \mathcal{N}(1, 0.1) \\
\beta &\sim \mathcal{N}(0, 0.3) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

```{r ch08_fit08_03}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "80 secs."))
fit08_03 <- xfun::cache_rds({
  out <- brm(data = dataRugged_nona,
      family = gaussian,
      formula = bf(log_gdp_s ~ 0 + a + b * rugged_sc,
         a ~ 0 + cid,
         b ~ 0 + cid,
         nl = TRUE),
      prior = c(prior(normal(1, 0.1), class = b, coef = cid1, nlpar = a),
                prior(normal(1, 0.1), class = b, coef = cid2, nlpar = a),
                prior(normal(0, 0.3), class = b, coef = cid1, nlpar = b),
                prior(normal(0, 0.3), class = b, coef = cid2, nlpar = b),
                prior(exponential(1), class = sigma)),
      iter = 1000, warmup = 500, chains = 4, cores = detectCores(),
      seed = 821)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_03")
tictoc::toc()
```

```{r}
posterior_summary(fit08_03) |>
  round(digits = 3)
```

```{r}
loo::loo_compare(fit08_01b, fit08_02, fit08_03, criterion = "waic") |> 
  print(simplify = FALSE)
```

with thew model weights

```{r}
brms::model_weights(fit08_01b, fit08_02, fit08_03) |>
  round(digits = 2)
```

### Plotting the interaction

```{r ch08_lpred08_03}
lpred08_03 <- crossing(
  cid = as.factor(1:2),
  rugged_sc = seq(from = -0.2, to = 1.2, length.out = 30)) |>
  mutate(rugged_sc = as.vector(scale(rugged_sc))) |>
  add_linpred_draws(fit08_03, ndraws = 50) |>
  mean_qi() |>
  mutate(is_africa = if_else(cid == 1, "Africa", "Not Africa")) |>
  mutate(is_africa = as.factor(is_africa))
# glimpse(lpred08_03)
```

```{r ch08_plot08_03}
ggplot(dataRugged_nona, aes(x = rugged_sc, y = log_gdp_s, fill = is_africa, color = is_africa)) +
  geom_smooth(data = lpred08_03, aes(x = rugged_sc, y = .linpred, ymin = .lower, ymax = .upper),
              stat = "identity",
              alpha = 1/4, size = 1/2) +
  geom_point(size = 1) +
  scale_fill_manual(values = c("Africa" = "springgreen2", "Not Africa" = "violet")) +
  scale_color_manual(values = c("Africa" = "springgreen2", "Not Africa" = "violet")) +
  coord_cartesian(xlim = c(0, 1)) +
  theme(legend.position = "none") +
  labs(title = "Figure 8.5",
       subtitle = "model b8.3",
       x = "ruggedness (standardized)",
       y = "log GDP (as proportion of mean)") +
  facet_wrap(~ is_africa)
```

## Symmetry of interactions

## Continuous interactions

### A winter flower

```{r ch08_dataTulips}
data(tulips, package = "rethinking")
dataTulips <- tulips |>
  mutate(blooms_r = scales::rescale(blooms),
         water_c = as.vector(scale(water, scale = FALSE)),
         shade_c = as.vector(scale(shade, scale = FALSE)))
rm(tulips)
```

### The models

#### Calibrating the priors

Our preliminary model, as a first jest in terms of prior is

$$
\begin{align*}
blooms\_r_i = blooms - \max(blooms)  &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= \alpha + \beta_W \cdot (water - \overline{water})+ \beta_S \cdot (shade - \overline{shade}) \\
&= \alpha + \beta_W \cdot water\_c_i+ \beta_S \cdot shade\_c_i
\alpha &\sim \mathcal{N}(0.5, 1) \\
\beta_W &\sim \mathcal{N}(0, 1) \\
\beta_S &\sim \mathcal{N}(0, 1) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

When looking at the data with `skim()` to evaluate the priors we obtain

```{r}
skim(dataTulips) |>
  select(-n_missing, -complete_rate) |>
  mutate(across(.cols = where(is.numeric), .fns = round, digits = 2))
```

We see that `blooms_r` must be between 0 and 1. The prior used assign most probability outside of that range

```{r}
param <- c("mean" = 0.5, "sd" = 1)
pnorm(q = -1, mean = param["mean"], sd = param["sd"]) + 
  pnorm(q = 1, mean = param["mean"], sd = param["sd"], lower.tail = FALSE)
```

lets say that we we want only 5% of the values outside the range (2.5% on each side) then, going with trial an error, the boundaries would be about

```{r}
param <- c("mean" = 0.5, "sd" = 0.25)
pnorm(q = -1, mean = param["mean"], sd = param["sd"]) + 
  pnorm(q = 1, mean = param["mean"], sd = param["sd"], lower.tail = FALSE)
```

Therefore we will use

$$
\alpha \sim \mathcal{N}(0.5, 0.25)
$$ and since the range for `water_c` and `shade_c` is -1 to 1 then we can use the same logic for both as follows

```{r}
param <- c("mean" = 0, "sd" = 0.25)
pnorm(q = -1, mean = param["mean"], sd = param["sd"]) + 
  pnorm(q = 1, mean = param["mean"], sd = param["sd"], lower.tail = FALSE)
```

which means virtually almost all values will be between -1 and 1. When looking at the `skim()` summary we see that there are many extreme values so this prior covers this situation well.

Therefore, our model with a little more informative priors is

$$
\begin{align*}
blooms\_r_i &\sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i &= \alpha + \beta_W \cdot water\_c_i+ \beta_S \cdot shade\_c_i \\
\alpha &\sim \mathcal{N}(0.5, 0.25) \\
\beta_W &\sim \mathcal{N}(0, 0.25) \\
\beta_S &\sim \mathcal{N}(0, 0.25) \\
\sigma &\sim \mathcal{Exp}(1)
\end{align*}
$$

and we fit that model to the data

```{r ch08_fit08_04}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
fit08_04 <- xfun::cache_rds({
  out <- brm(data = dataTulips,
      family = gaussian,
      formula = blooms_r ~ 1 + water_c + shade_c,
      prior = c(prior(normal(0.5, 0.25), class = Intercept),
                prior(normal(0, 0.25), class = b, coef = water_c),
                prior(normal(0, 0.25), class = b, coef = shade_c),
                prior(exponential(1), class = sigma)),
      iter = 1000, warmup = 500, chains = 2, cores = detectCores(),
      seed = 823)
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_04")
tictoc::toc()
```

```{r}
posterior_summary(fit08_04) |>
  round(digits = 3)
```

#### Adding an interaction

Using the notation $\gamma_{W, i} = \beta_W+\beta_{WS} \cdot shade\_c_i$ we get the new model with interactions

$$
\begin{align*}
blooms\_r_i & \sim \mathcal{N}(\mu_i, \sigma) \\
\mu_i & = 
 \alpha + \gamma_{W, i} \cdot water\_c_i+ \beta_S \cdot shade\_c_i \\
 &= \alpha + \beta_W \cdot water\_c_i + \beta_S \cdot shade\_c_i + \beta_{WS} \cdot shade\_c_i \cdot water\_c_i
 \\
\alpha & \sim \mathcal{N}(0.5, 0.25) \\
\beta_W & \sim \mathcal{N}(0, 0.25) \\
\beta_S & \sim \mathcal{N}(0, 0.25) \\
\beta_{WS} & \sim \mathcal{N}(0, 0.25) \\
\sigma & \sim \mathcal{Exp}(1)
\end{align*}
$$

and now fitting the model with interaction

```{r ch08_fit08_05}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "60 secs."))
fit08_05 <- xfun::cache_rds({
  out <- update(fit08_04,
                newdata = dataTulips,
                formula. = blooms_r ~ 1 + water_c + shade_c + water_c:shade_c,
                prior = c(prior(normal(0.5, 0.25), class = Intercept),
                  prior(normal(0, 0.25), class = b, coef = water_c),
                  prior(normal(0, 0.25), class = b, coef = shade_c),
                  prior(normal(0, 0.25), class = b, coef = water_c:shade_c),
                  prior(exponential(1), class = sigma)))
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_05")
tictoc::toc()
```

```{r}
posterior_summary(fit08_05) |>
  round(digits = 3)
```

### Plotting continuous interactions

```{r ch08_plot08_04post}
plot08_04post <- list()
plot08_04post <- within(plot08_04post, {
  draws <- data.frame(water_c = -1:1, shade_c = -1:1) |>
    expand(water_c, shade_c) |>
    add_linpred_draws(fit08_04, ndraws = 20) |>
    mutate(label = paste("shade_c =", shade_c)) |>
    identity()
  
  intrvl <- draws |>
    select(-label) |>
    mean_qi() |>
    mutate(.draw = 0,
           label = paste("shade_c =", shade_c))
  
  p <- ggplot(draws, aes(x = water_c, y = .linpred, group = .draw)) +
    geom_line(color = "orchid") +
    geom_line(data = intrvl, aes(x = water_c, y = .linpred, group = .draw), 
              color = "darkgreen", linewidth = 1) +
    geom_hline(
      data = data.frame(y = c(0, 1)), aes(yintercept = y),
      size = 1, linetype = 2, color = "royalblue") +
  coord_cartesian(xlim = c(-1, 1), ylim = c(-0.5, 1.5)) +
    facet_wrap(. ~ label, nrow = 1) +
  labs(
    subtitle = "Model 8.4 NO INTERACTIONS\nIntercept ~ dnorm(0.5, 0.25)\nb ~ dnorm(0, 0.25)",
    x = "water (centered)",
    y = "bloom (rescaled)")
})
plot08_04post$p
```

```{r ch08_plot08_05post}
plot08_05post <- list()
plot08_05post <- within(plot08_05post, {
  draws <- data.frame(water_c = -1:1, shade_c = -1:1) |>
    expand(water_c, shade_c) |>
    add_linpred_draws(fit08_05, ndraws = 20) |>
    mutate(label = paste("shade_c =", shade_c))
  
  intrvl <- draws |>
    select(-label) |>
    mean_qi() |>
    mutate(.draw = 0,
           label = paste("shade_c =", shade_c))
  
  p <- ggplot(draws, aes(x = water_c, y = .linpred, group = .draw)) +
    geom_line(color = "orchid") +
    geom_line(data = intrvl, aes(x = water_c, y = .linpred, group = .draw), 
              color = "darkgreen", linewidth = 1) +
    geom_hline(
      data = data.frame(y = c(0, 1)), aes(yintercept = y),
      size = 1, linetype = 2, color = "royalblue") +
  # geom_abline(intercept = est["fixed"], slope = est["b"], color = "purple", size = 1) +
  coord_cartesian(xlim = c(-1, 1), ylim = c(-0.5, 1.5)) +
    facet_wrap(. ~ label, nrow = 1) +
  labs(
    subtitle = "Model 8.5 WITH INTERACTIONS\nIntercept ~ dnorm(0.5, 0.25)\nb ~ dnorm(0, 0.25)",
    x = "water (centered)",
    y = "bloom (rescaled)")
})
# plot08_05post$intrvl
plot08_05post$p
```

```{r}
plot08_04post$p / plot08_05post$p +
  plot_annotation(title = "Tryptich plot of predicted bloom by level of shade (-1, 0, 1)")
```

### Plotting prior predictions

::: callout-warning
The plot in this section are different than what McElreath and Kurz have. Yet, McElreath"s are different than Kurz's, so there doesn't seem to be a consensus between the 2. I keep the plots below as I don't see anything wrong with them.
:::

```{r ch08_fit08_04prior}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "20 secs."))
fit08_04prior <- xfun::cache_rds({
  out <- update(fit08_04,
                sample_prior = "only")
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_04prior")
tictoc::toc()
```

```{r}
posterior_summary(fit08_04prior)|>
  round(digits = 3)
```

```{r ch08_plot08_04prior}
plot08_04prior <- list()
plot08_04prior <- within(plot08_04prior, {
  draws <- data.frame(water_c = -1:1, shade_c = -1:1) |>
    expand(water_c, shade_c) |>
    add_linpred_draws(fit08_04prior, ndraws = 20) |>
    mutate(label = paste("shade_c =", shade_c))
  
  intrvl <- draws |>
    select(-label) |>
    mean_qi() |>
    mutate(.draw = 0,
           label = paste("shade_c =", shade_c))
  
  p <- ggplot(draws, aes(x = water_c, y = .linpred, group = .draw)) +
    geom_line(color = "orchid") +
    geom_line(data = intrvl, aes(x = water_c, y = .linpred, group = .draw), 
              color = "darkgreen", linewidth = 1) +
    geom_hline(
      data = data.frame(y = c(0, 1)), aes(yintercept = y),
      size = 1, linetype = 2, color = "royalblue") +
  coord_cartesian(xlim = c(-1, 1), ylim = c(-0.5, 1.5)) +
    facet_wrap(. ~ label, nrow = 1) +
  labs(
    subtitle = "Model 8.4 WITHOUT INTERACTIONS\nIntercept ~ dnorm(0.5, 0.25)\nb ~ dnorm(0, 0.25)",
    x = "water (centered)",
    y = "bloom (rescaled)")
})
# plot08_04prior$intrvl
plot08_04prior$p
```

```{r ch08_fit08_05prior}
tictoc::tic(msg = sprintf("run time of %s, use the cache.", "20 secs."))
fit08_05prior <- xfun::cache_rds({
  out <- update(fit08_05,
                sample_prior = "only")
  out <- brms::add_criterion(out, criterion = c("waic", "loo"))},
  file = "ch08_fit08_05prior")
tictoc::toc()
```

```{r}
posterior_summary(fit08_05prior)|>
  round(digits = 3)
```

```{r ch08_plot08_05prior}
plot08_05prior <- list()
plot08_05prior <- within(plot08_05prior, {
  draws <- data.frame(water_c = -1:1, shade_c = -1:1) |>
    expand(water_c, shade_c) |>
    add_linpred_draws(fit08_05prior, ndraws = 20) |>
    mutate(label = paste("shade_c =", shade_c))
  
  intrvl <- draws |>
    select(-label) |>
    mean_qi() |>
    mutate(.draw = 0,
           label = paste("shade_c =", shade_c))
  
  p <- ggplot(draws, aes(x = water_c, y = .linpred, group = .draw)) +
    geom_line(color = "orchid") +
    geom_line(data = intrvl, aes(x = water_c, y = .linpred, group = .draw), 
              color = "darkgreen", linewidth = 1) +
    geom_hline(
      data = data.frame(y = c(0, 1)), aes(yintercept = y),
      size = 1, linetype = 2, color = "royalblue") +
  coord_cartesian(xlim = c(-1, 1), ylim = c(-0.5, 1.5)) +
    facet_wrap(. ~ label, nrow = 1) +
  labs(
    subtitle = "Model 8.5 WITH INTERACTIONS\nIntercept ~ dnorm(0.5, 0.25)\nb ~ dnorm(0, 0.25)",
    x = "water (centered)",
    y = "bloom (rescaled)")
})
# plot08_05prior$intrvl
plot08_05prior$p
```

```{r}
plot08_04prior$p / plot08_05prior$p +
  plot_annotation(title = "Tryptich plot of prior predicted bloom by level of shade (-1, 0, 1)")
```

## Summary