ch8_conditional.qmd

---
title: "ch8_conditional_manatees"
editor: visual
execute:
  echo: false
  warning: false
  message: false
  cache: true
  cache.lazy: false
  fig-align: center
---

```{r setup}
knitr::opts_chunk$set(fig.align = "center")

library(tidyverse)
library(brms)
library(bayesplot)
library(tidybayes)
library(patchwork)
library(GGally)
library(dagitty)
library(ggdag)
library(ggrepel)
library(ggthemes)
library(rcartocolor)

library(rethinking)
data(rugged)
data(tulips)

detach(package:rethinking, unload = T)

at <- c(-3, -2, -1, 0, 1, 2, 3)


```

The Armstrong Whitworth A.W.38 Whitley was a frontline Royal Air Force bomber. During the second World War, the A.W.38 carried bombs and pamphlets into German territory. The A.W.38 has fierce natural enemies: artillery and interceptor fire. Many planes never returned from their missions. And those that survived had the scars to prove it. Most observers intuit that helping bombers means reducing the kind of damage we see on them, perhaps by adding armor to the parts of the plane that show the most damage.

Up-armoring the damaged portions of returning bombers did little good. Instead, improving the A.W.38 bomber meant armoring the undamaged sections. The evidence from surviving bombers is misleading, because it is conditional on survival. Bombers that returned home conspicuously lacked damage to the cockpit and engines. They got lucky. Bombers that never returned home were less so. To get the right answer, in either context, we have to realize that the kind of damage seen is conditional on survival.

![](plots/airplane.png){fig-align="center" width=50%}

Conditioning is one of the most important principles of statistical inference. Data, like the bomber damage, are conditional on how they get into our sample. Posterior distributions are conditional on the data. All model-based inference is conditional on the model. Every inference is conditional on something, whether we notice it or not.

Simple linear models frequently fail to provide enough conditioning, however. Every model so far in this book has assumed that each predictor has an independent association with the mean of the outcome. What if we want to allow the association to be conditional? For example, in the primate milk data from the previous chapters, suppose the relationship between milk energy and brain size varies by taxonomic group (ape, monkey, prosimian). This is the same as suggesting that the influence of brain size on milk energy is conditional on taxonomic group. The linear models of previous chapters cannot address this question.

To model deeper conditionality—where the importance of one predictor depends upon another predictor—we need interaction (also known as moderation). Interaction is a kind of conditioning, a way of allowing parameters (really their posterior distributions) to be conditional on further aspects of the data.

More generally, interactions are central to most statistical models beyond the cozy world of Gaussian outcomes and linear models of the mean. In generalized linear models (GLMs), even when one does not explicitly define variables as interacting, they will always interact to some degree. Multilevel models induce similar effects. Common sorts of multilevel models are essentially massive interaction models, in which estimates (intercepts and slopes) are conditional on clusters (person, genus, village, city, galaxy) in the data. Multilevel interaction effects are complex. They’re not just allowing the impact of a predictor variable to change depending upon some other variable, but they are also estimat- ing aspects of the distribution of those changes. This may sound like genius, or madness, or both. Regardless, you can’t have the power of multilevel modeling without it.


##  Discrete ~Continuous Interactions


### Ruggedness vs GDP in Africa 🌍

In this vignette we look at the Log gdp per capita of different countries. The more rugged the terrain the lower the gdp is on average, which fits a nice causal story of the harder it is to cross the land the less trade and wealth is created. This story is reversed though when we look at countries in Africa, as seen below in the figure.

```{r}
rugged_clean <-
  rugged %>%
  mutate(log_gdp = log(rgdppc_2000)) %>% 
  filter(complete.cases(rgdppc_2000)) %>%  # countries with GDP data
  # re-scale variables
  mutate(log_gdp_std = log_gdp / mean(log_gdp), 
         rugged_std  = rugged / max(rugged),
         rugged_std_centered = rugged_std - mean(rugged_std))
```


```{r}
p1 <- rugged_clean %>% 
  filter(cont_africa == 1) %>% 
  ggplot(aes(x = rugged_std, y = log_gdp_std)) +
  geom_smooth(method = "lm", formula = y ~ x,
              fill = palette_pander(n = 2)[1],
              color = palette_pander(n = 2)[1]) +
  geom_point(color = palette_pander(n = 2)[1]) +
  geom_text_repel(data = . %>% 
                    filter(country %in% c("Lesotho", "Seychelles")),  
                  aes(label = country), 
                  size = 3, family = "Times", seed = 8) +
  labs(subtitle = "African nations",
       x = "ruggedness (standardized)",
       y = "log GDP (as proportion of mean)")+
  theme_minimal()

p2 <-
  rugged_clean %>% 
  filter(cont_africa == 0) %>% 
  ggplot(aes(x = rugged_std, y = log_gdp_std)) +
  geom_smooth(method = "lm", formula = y ~ x,
              fill = palette_pander(n = 2)[2],
              color = palette_pander(n = 2)[2]) +
  geom_point(color = palette_pander(n = 2)[2]) +
  geom_text_repel(data = . %>% 
                    filter(country %in% c("Switzerland", "Tajikistan")),  
                  aes(label = country), 
                  size = 3, family = "Times", seed = 8) +
  xlim(0, 1) +
  labs(subtitle = "Non-African nations",
       x = "ruggedness (standardized)",
       y = "log GDP (as proportion of mean)")+
  theme_minimal()

# combine
p1 + p2 + plot_annotation(title = "Figure 8.2. Separate linear regressions inside and outside of Africa")
```

What's to make of the story? Maybe historical slavery has a lasting effect on the economies of african countries, the flatter more accessible regions of africa had slaving?

```{r}
#| fig-width: 4
#| fig-height: 3
dag_coords <- tibble(name = c("R", "G", "U", "C"),
                     x = c(1, 2, 2, 2.5),
                     y = c(2, 2, 1, 2))

m6 <- dagify("G" ~ "R" + "U" + "C",
             "R" ~ "U",
             coords = dag_coords) %>% 
  ggplot(aes(x = x, y = y, xend = xend, yend = yend))+
  geom_dag_point(color = "firebrick", alpha = 1/4, size = 10)+
  geom_dag_point(data = . %>% filter(name == "U"), 
                 color = "firebrick4", fill = NA, size = 10, shape = 21, 
                 stroke = 2, linetype = "dashed")+
  geom_dag_text( color = "firebrick", parse = TRUE) +
  geom_dag_edges(edge_color = "firebrick") +
  scale_x_continuous(NULL, breaks = NULL, expand = c(0.1, 0.1)) +
  scale_y_continuous(NULL, breaks = NULL, expand = c(0.2, 0.2)) +
  theme_bw() +
  theme(panel.grid = element_blank())+
  labs(caption = "G = GDP, R = Ruggedness, C = Continent, U = Unobserved Effect")

m6
```

In the graph above we formalize this idea by saying that ruggedness $R$ influences the current GDP $G$. Both $R$ & $G$ are influenced by some set of unknown confounders $U$ like proximity to coastline, which we'll ignore for the moment. Finally $C$ the continent effects $G$ as well, crucially $R$ and $C$ could be independent or interact on their influence on $G$. The DAG does not display an interaction, instead we declare outside of the graph like this: $G = f(R, C)$.

How do we estimate that function $f$, we could split up the data and make separate models 1 for africa and 1 for all the other continents countries. But this lead to a poor estimate of other parameters like $\sigma$. Additionally if we wanted to compare models with an information criteria we'd need to use the same data, so splitting the data also hurts that process.

Our first model is:

$$ log(y_i) \ sim Normal(\mu_i, \sigma)$$

$$ \mu_i = \alpha + \beta(\text{rugged}_i - \overline{rugged})$$
```{r}
b8.0 <- brm(data = rugged_clean,
            family = gaussian(),
            log_gdp_std ~ 1 + rugged_std_centered,
            prior = c(
              prior(normal(1, 0.1), class = Intercept),
              prior(normal(0, 0.3), class = b),
              prior(exponential(1), class = sigma)
            ),
            iter = 2000, warmup = 500, cores = 4, seed = 5,
            backend = "cmdstanr", silent = 2, file = "fits/b08.0.1")
```

::: panel-tabset

##### 📈 μ heatmap

```{r}

min_rugged = min(rugged_clean$rugged_std_centered)
max_rugged = max(rugged_clean$rugged_std_centered)

simModel <- as_tibble(b8.0) %>% 
  mutate(simRuggedness = seq(from = min_rugged, to = max_rugged, length.out = n()),
         simGDPEst= Intercept + (b_rugged_std_centered * simRuggedness),
         simGDP = rnorm(n(), simGDPEst, sd = sigma))
         
# ymin <- min(c(simModel$simGDP, simModel$simGDPEst, rugged_clean$log_gdp_std))
# ymax <- max(c(simModel$simGDP, simModel$simGDPEst, rugged_clean$log_gdp_std))

at <- seq(from = min_rugged, to = max_rugged, length.out = 5)

modelEst_plot <- ggplot() +
  stat_density_2d(data = simModel, 
                  aes(x = simRuggedness, y = simGDPEst, fill = after_stat(ndensity)),
                  geom = "raster", contour = FALSE) +
  scale_fill_viridis_c(option = "magma") +
  geom_point(data = rugged_clean,
             aes(x = rugged_std_centered, y = log_gdp_std), 
             shape = 21, color = "white", fill = "black", lwd = 3, alpha = .8)+
  labs(y = "Log GDP (1.0 = average country)",
       title = "Log GDP ~ Ruggedness", subtitle = "Mu estimate")+
  theme_minimal()+
  scale_x_continuous("Ruggedness (0 = minimum ruggedness, 100 = maximum ruggedness)",
                     breaks = at,
                     labels = round(at + mean(rugged_clean$rugged_std), 1) * 100) +
  guides(fill = "none")

modelEst_plot
```
#####  🎛 Parameters

```{r}
#| fig-width: 10
#| fig-height: 2

as_tibble(b8.0) %>% 
  rename("Ruggedness" = b_rugged_std_centered) %>% 
  dplyr::select(c(`Ruggedness`, `Intercept`)) %>% 
  pivot_longer(cols = everything(),
               names_to = "Covariate",
               values_to = "Effect") %>% 
  ggplot(aes(x = `Effect`, y = reorder(Covariate, `Effect`))) +
  stat_halfeye(point_interval = median_qi, .width = .95,
               fill = "firebrick4") +
  labs(x = "Effect on Log GDP per capita",
       y = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank())

```



:::

#### Indicator Variable Solution

The first thing to realize is that just including an indicator variable for African nations, won't reveal the reversed slope. 

To build a model that alls nations inside nd outside Africa to have different intercepts, we need to modify the model for $\mu_i$ so that the mean is conditional on continent. The conventional way to do this would be to just add another term to the linear model:

$$ \mu_i = \alpha + \beta_1(\text{rugged}_i - \overline{rugged}) - \beta_2  \mathbf{I} (\text{Africa}_i)$$

Where $\text{Africa}_i$ is a 0/1 indicator variable for Africa or not Africa. But this model assumes that african countries have more uncertainty inherently built into their $mu$ estimate, which makes no sense.

Our simple solution is to create separate intercepts for the different categories, like so:

$$ \mu_i = \alpha_\text{Africa [i]} + \beta_1(\text{rugged}_i - \overline{rugged})$$
Where $\text{Africa [i]}$ is an index variable which takes the value 1 for African nations and 2 for all other nations.

```{r}
rugged_clean <- rugged_clean %>% 
  mutate(cont = ifelse(cont_africa == 1, "African", "Not African"))

b8.1 <- brm(data = rugged_clean,
            family = gaussian(),
            log_gdp_std ~ 0 + cont + rugged_std_centered,
            prior = c(
              prior(normal(0.9, 0.1), class = b, coef = "contAfrican"),
              prior(normal(1.1, 0.1), class = b, coef = "contNotAfrican"),
              prior(normal(0, 0.3), class = b, coef = "rugged_std_centered"),
              prior(exponential(1), class = sigma)
            ),
            iter = 2000, warmup = 1000, cores = 4, seed = 5,
            backend = "cmdstanr", silent = 2, file = "fits/b08.1.2")

# get_prior(data = rugged_clean,
#             family = gaussian(),
#             log_gdp_std ~ 0 + cont + rugged_std_centered)


```



::: panel-tabset

##### 📈 μ heatmap

::: panel-tabset

###### Plot 1

```{r}
simCov <- 
  crossing(cont = c("African", "Not African"),
           rugged_std = seq(from = -0.2, to = 1.2, length.out = 30)) %>% 
  mutate(rugged_std_centered = rugged_std - mean(rugged_clean$rugged_std))

fitted(b8.1, 
       newdata = simCov,
       summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  bind_cols(expand_grid(draws = 1:4000, simCov)) %>% 
  ggplot(aes(x = rugged_std_centered, y = value, fill = cont, color = cont))+
  stat_lineribbon(.width = seq(from = .03, to = .99, by =.03),
                  alpha = .1, size = 0) +
  geom_point(data = rugged_clean, aes(x = rugged_std_centered, y = log_gdp_std,
                                      fill = cont, color = cont))+
  theme_minimal() + 
  labs(y = "log GDP (as proportion of mean)") +
  scale_x_continuous("Ruggedness (0 = minimum ruggedness, 100 = maximum ruggedness)",
                   breaks = at,
                   labels = round(at + mean(rugged_clean$rugged_std), 1) * 100,
                   limits = c(floor(min(at) * 10) / 10, ceiling(max(at) * 10) / 10))+
  theme(legend.background = element_blank(),
      legend.direction = "horizontal",
      legend.position = c(.67, .93),
      legend.title = element_blank())

```

###### plot 2

```{r}
fitted(b8.1, 
       newdata = simCov,
       summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  bind_cols(expand_grid(draws = 1:4000, simCov)) %>% 
  ggplot(aes(x = rugged_std_centered, y = value)) +
  stat_density_2d(aes(fill = after_stat(ndensity)), 
                  geom = "raster", contour = FALSE) +
  scale_fill_viridis_c(option = "magma") +
  geom_point(data = rugged_clean, 
             aes(x = rugged_std_centered, y = log_gdp_std),
             shape = 21, color = "white", fill = "black", size = 2, alpha = 0.8) +
  facet_wrap(~ cont, labeller = label_both) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "black"),
    plot.background = element_rect(fill = "black"),
    strip.text = element_text(color = "white", size = 12),
    axis.text = element_text(color = "white"),
    axis.title.x = element_text(color = "white"),
    axis.title.y = element_text(color = "white", angle = 90, vjust = 0.5),  # Rotated y-axis title
    plot.title = element_text(color = "white")
  ) +
  labs(y = "log GDP (as proportion of mean)") +
  scale_x_continuous("Ruggedness (0 = minimum ruggedness, 100 = maximum ruggedness)",
                     breaks = at,
                     labels = round(at + mean(rugged_clean$rugged_std), 1) * 100) +
  guides(fill = "none")
```
:::

#####  🎛 Parameters

```{r}
#| fig-width: 10
#| fig-height: 2

as_tibble(b8.1) %>% 
  rename("Ruggedness" = b_rugged_std_centered,
         "African" = b_contAfrican,
         "Non-African" = b_contNotAfrican) %>% 
  dplyr::select(c(`Ruggedness`, `African`, `Non-African`)) %>% 
  pivot_longer(cols = everything(),
               names_to = "Covariate",
               values_to = "Effect") %>% 
  ggplot(aes(x = `Effect`, y = reorder(Covariate, `Effect`))) +
  stat_halfeye(point_interval = median_qi, .width = .95,
               fill = "firebrick4") +
  labs(x = "Effect on Log GDP per capita",
       y = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank())

```

:::

Adding Africa as an index covariate didn't get that reversal of slopes that we saw in our first plot. But when we compare this model to not having the Africa index covariate with WAIC, its clearly superior. So the Africa variable is clearly picking up some important association in the data. African countries have a lower average Log GDP regardless of their ruggedness and that's what our African index covariate supports.



```{r}
#| fig-width: 10
#| fig-height: 2

b8.0 <- add_criterion(b8.0, criterion = "waic")
b8.1 <- add_criterion(b8.1, criterion = "waic") 

w <- loo_compare(b8.0, b8.1, criterion = "waic")

w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column("model_name") %>% 
  mutate(model_name = fct_reorder(model_name, waic, .desc = T)) %>% 
  
  ggplot(aes(x = waic, y = model_name, 
             xmin = waic - se_waic, 
             xmax = waic + se_waic)) +
  geom_pointrange(color = carto_pal(7, "BurgYl")[7], 
                  fill = carto_pal(7, "BurgYl")[5], shape = 21) +
  labs(title = "WAIC model comparison plot",
       x = NULL, y = NULL) +
  theme(axis.ticks.y = element_blank())+
  theme_minimal()

b8.0 <- add_criterion(b8.0, criterion = "loo")
b8.1 <- add_criterion(b8.1, criterion = "loo")

loo_compare(b8.0, b8.1, criterion = "loo") %>% 
  print(simplify = F)
```

#### Adding an Interaction does work

How can we we recover the difference in slope that we saw at the beginning of this section? We need a proper interaction effect. This just means we make the slope conditional on whether it's part of Africa or not. 

Just above we modeled 

$$ \mu_i = \alpha_\text{Africa [i]} + \beta_1(\text{rugged}_i - \overline{rugged})$$

But now we're going to make an index for $\beta$ as well.

$$ \mu_i = \alpha_\text{Africa [i]} + \beta_\text{Africa[i]}(\text{rugged}_i - \overline{rugged})$$

```{r}
b8.2 <- brm(data = rugged_clean,
            family = gaussian(),
            bf(log_gdp_std ~ 0 + a + b * rugged_std_centered,
               a ~ 0 + cont,
               b  ~ 0 + cont,
               nl = TRUE),
            prior = c(
              prior(normal(1, 0.1), class = b, nlpar = a, coef = contAfrican),
              prior(normal(1, 0.1), class = b, nlpar = a, coef = contNotAfrican),
              prior(normal(0, 0.3), class = b, nlpar = b, coef = contAfrican),
              prior(normal(0, 0.3), class = b, nlpar = b, coef = contNotAfrican),
              prior(exponential(1), class = sigma)
            ),
            iter = 2000, warmup = 1000, cores = 4, seed = 5,
            backend = "cmdstanr", silent = 2, file = "fits/b08.2.1")

# get_prior(data = rugged_clean,
#             family = gaussian(),
#             bf(log_gdp_std ~ 0 + a + b * rugged_std_centered,
#                a ~ 0 + cont,
#                b  ~ 0 + cont,
#                nl = TRUE))
```

::: panel-tabset

##### 📈 μ heatmap

::: panel-tabset

###### Plot 1

```{r}
simCov <- 
  crossing(cont = c("African", "Not African"),
           rugged_std = seq(from = -0.2, to = 1.2, length.out = 30)) %>% 
  mutate(rugged_std_centered = rugged_std - mean(rugged_clean$rugged_std))

fitted(b8.2, 
       newdata = simCov,
       summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  bind_cols(expand_grid(draws = 1:4000, simCov)) %>% 
  ggplot(aes(x = rugged_std_centered, y = value, fill = cont, color = cont))+
  stat_lineribbon(.width = seq(from = .03, to = .99, by =.03),
                  alpha = .1, size = 0) +
  geom_point(data = rugged_clean, aes(x = rugged_std_centered, y = log_gdp_std,
                                      fill = cont, color = cont))+
  theme_minimal() + 
  labs(y = "log GDP (as proportion of mean)") +
  scale_x_continuous("Ruggedness (0 = minimum ruggedness, 100 = maximum ruggedness)",
                   breaks = at,
                   labels = round(at + mean(rugged_clean$rugged_std), 1) * 100,
                   limits = c(floor(min(at) * 10) / 10, ceiling(max(at) * 10) / 10))+
  theme(legend.background = element_blank(),
      legend.direction = "horizontal",
      legend.position = c(.67, .93),
      legend.title = element_blank())

```





###### plot 2

```{r}
fitted(b8.2, 
       newdata = simCov,
       summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  bind_cols(expand_grid(draws = 1:4000, simCov)) %>% 
  ggplot(aes(x = rugged_std_centered, y = value)) +
  stat_density_2d(aes(fill = after_stat(ndensity)), 
                  geom = "raster", contour = FALSE) +
  scale_fill_viridis_c(option = "magma") +
  geom_point(data = rugged_clean, 
             aes(x = rugged_std_centered, y = log_gdp_std),
             shape = 21, color = "white", fill = "black", size = 2, alpha = 0.8) +
  facet_wrap(~ cont, labeller = label_both) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "black"),
    plot.background = element_rect(fill = "black"),
    strip.text = element_text(color = "white", size = 12),
    axis.text = element_text(color = "white"),
    axis.title.x = element_text(color = "white"),
    axis.title.y = element_text(color = "white", angle = 90, vjust = 0.5),  # Rotated y-axis title
    plot.title = element_text(color = "white")
  ) +
  labs(y = "log GDP (as proportion of mean)") +
  scale_x_continuous("Ruggedness (0 = minimum ruggedness, 100 = maximum ruggedness)",
                     breaks = at,
                     labels = round(at + mean(rugged_clean$rugged_std), 1) * 100) +
  guides(fill = "none")
```

:::

##### 🎛 Parameters

```{r}
#| fig-width: 10
#| fig-height: 2

as_tibble(b8.2) %>% 
  rename("Intercept African" = b_a_contAfrican,
         "Intercept Non-African" = b_a_contNotAfrican,
         "Ruggedness African" = b_b_contAfrican,
         "Ruggedness Non-African" = b_b_contNotAfrican,
         ) %>% 
  dplyr::select(c(`Intercept African`, `Intercept Non-African`, `Ruggedness African`, `Ruggedness Non-African`)) %>% 
  pivot_longer(cols = everything(),
               names_to = "Covariate",
               values_to = "Effect") %>% 
  ggplot(aes(x = `Effect`, y = reorder(Covariate, `Effect`))) +
  stat_halfeye(point_interval = median_qi, .width = .95,
               fill = "firebrick4") +
  labs(x = "Effect on Log GDP per capita",
       y = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank())

```

:::

```{r}
#| fig-width: 10
#| fig-height: 2

b8.2 <- add_criterion(b8.2, criterion = "loo")
b8.2 <- add_criterion(b8.2, criterion = "waic")

loo_compare(b8.0, b8.1, b8.2, criterion = "loo") %>% 
  print(simplify = F)



w <- loo_compare(b8.0, b8.1, b8.2, criterion = "waic")

w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column("model_name") %>% 
  mutate(model_name = fct_reorder(model_name, waic, .desc = T)) %>% 
  
  ggplot(aes(x = waic, y = model_name, 
             xmin = waic - se_waic, 
             xmax = waic + se_waic)) +
  geom_pointrange(color = carto_pal(7, "BurgYl")[7], 
                  fill = carto_pal(7, "BurgYl")[5], shape = 21) +
  labs(title = "WAIC model comparison plot",
       x = NULL, y = NULL) +
  theme(axis.ticks.y = element_blank())+
  theme_minimal()
```

Out interaction model performs best yet, it receives more than 95% of the model weight. Though the remaining weight percentage placed on our 2nd model, does mean there is some reason to think that our model is a bit overfit on its slopes.

#####
The interaction there has two equally valid phrasings.
(1) How much does the association between ruggedness and log GDP depend upon whether the nation is in Africa?
(2) How much does the association of Africa with log GDP depend upon ruggedness?

While these two possibilities sound different to most humans, our statistical engine thinks they are identical.

```{r}
fitted(b8.2, 
       newdata = simCov,
       summary = F) %>%
  data.frame() %>%
  pivot_longer(everything()) %>% 
  bind_cols(expand_grid(draws = 1:4000, simCov)) %>% 
  select(-name) %>% 
  pivot_wider(names_from = cont, values_from = value) %>% 
  mutate(delta =   `African` - `Not African`) %>% 
  ggplot(aes(x = rugged_std, y = delta)) + 
  stat_lineribbon(.width = .95, fill = palette_pander(n = 8)[8], alpha = 3/4) +
  geom_hline(yintercept = 0, linetype = 2) +
  annotate(geom = "text",
           x = .2, y = 0,
           label = "Africa higher GDP\nAfrica lower GDP",
           family = "Times") +
  labs(x = "ruggedness (standardized)",
       y = "expected difference log GDP") +
  coord_cartesian(xlim = c(0, 1),
                  ylim = c(-0.3, 0.2))+
  theme_minimal()
  
```

This plot is counter-factual. There is no raw data here. Instead we are seeing through the model’s eyes and imagining comparisons between iden- tical nations inside and outside Africa, as if we could independently manipulate continent and also terrain ruggedness. Below the horizontal dashed line, African nations have lower expected GDP. This is the case for most terrain ruggedness values. But at the highest rugged- ness values, a nation is possibly better off inside Africa than outside it. Really it is hard to find any reliable difference inside and outside Africa, at high ruggedness values. It is only in smooth nations that being in Africa is a liability for the economy.

Simple interactions are symmetric. Within the model, there’s no basis to prefer one interpretation over the other, because in fact they are the same interpretation. But when we reason causally about models, our minds tend to prefer one interpretation over the other, because it’s usually easier to imagine manipulating one of the predictor variables instead of the other. In this case, it’s hard to imagine manipulating which continent a nation is on. But it’s easy to imagine manipulating terrain ruggedness, by flattening hills or blasting tunnels through mountains. Africa’s unusually positive relationship with terrain ruggedness is due to historical causes, not contemporary terrain, then tunnels might improve economies in the present. At the same time, continent is not really a cause of economic activity. Rather there are historical and political factors associated with continents, and we use the continent variable as a proxy for those factors. It is manipulation of those other factors that would matter.


### Continuous - Continuous Interactions

Richard McElreath writes:

>I want to convince the reader that interaction effects are difficult to interpret. They are nearly impossible to interpret, using only posterior means and standard deviations. Once interactions exist, multiple parameters are always in play at the same time. It is hard enough with the simple, categorical interactions from the terrain ruggedness example. Once we start modeling interactions among more than one continuous variables, it gets much harder. It’s one thing to make a slope conditional upon a category. In such a context, the model reduces to estimating a different slope for each category. But it’s quite a lot harder to understand that a slope varies in a continuous fashion with a continuous variable. Interpretation is much harder in this case, even though the mathematics of the model are essentially the same as in the categorical case.

#### Tulips Example

The data in this example are sizes of blooms from beds of tulips grown in greenhouses, under different soil and light conditions. The blooms column will be our outcome—what we wish to predict. The water and shade columns will be our predictor variables.

Since both light and water help plants grow and produce blooms, it stands to reason that the independent effect of each will be to produce bigger blooms. But we’ll also be interested in the interaction between these two variables. In the absence of light, for example, it’s hard to see how water will help a plant—photosynthesis depends upon both light and water. Like- wise, in the absence of water, sunlight does a plant little good. One way to model such an interdependency is to use an interaction effect. In the absence of a good mechanistic model of the interaction, one that uses a theory about the plant’s physiology to hypothesize the functional relationship between light and water, then a simple linear two-way interaction is a good start. But ultimately it’s not close to the best that we could do.

The causal DAG for this model is quite simple we have $\text{Water} \rightarrow \text{Bloom} \leftarrow \text{Shade}$. As before the DAG doesn't tell us the function through which Water and Shade jointly include Bloom. In principle, every unique combination of Water and Shade could have a different mean Bloom. But we'll start simpler.

##### Model with No Interaction



$$ \text{Bloom}_i \sim Normal(\mu_\text{Bloom}_i, \sigma)$$

$$ \mu_{\text{Bloom}_i} = \alpha + \beta_W(\text{Water}_i - \overline{Water}) + \beta_S(\text{Shade}_i - \overline{Shade})$$
```{r}

tulips_clean <- tulips %>% 
  mutate(water_z = (water - mean(water)),
         shade_z = (shade - mean(shade)),
         blooms_s = blooms / max(blooms))

ggplot(data = tulips_clean)+
  geom_jitter(aes(x = water_z, y = blooms_s, color = shade_z), height = 0, width = .03)


n = 1e4
prior_pred <- tibble(Intercept = rnorm(n, 0.5, 0.12),
                     shade_beta = rnorm(n, 0, 0.05),
                     water_beta = rnorm(n, 0, 0.05),
                     shade_water_interaction_beta = rnorm(n, 0, 0.05),
                     sigma = rexp(n, 2),
                     water_z = sample(c(-1, 0, 1), size = n, replace = TRUE),
                     shade_z = sample(c(-1, 0, 1), size = n, replace = TRUE),
                     mu = Intercept + (shade_beta * shade_z) + (water_beta * water_z) + 
                       (shade_water_interaction_beta * shade_z * water_z),
                     ppd = rnorm(n, mu, sigma)
                     )

hist(prior_pred$mu, br = 30)
hist(prior_pred$ppd, br = 100)



b8.3 <- brm(data = tulips_clean,
            family = gaussian(),
            blooms_s ~ 1 + water_z + shade_z,
            prior = c(
              prior(normal(0.5, 0.12), class = Intercept),
              prior(normal(0, 0.05), class = b, coef = shade_z),
              prior(normal(0, 0.05), class = b, coef = water_z),
              prior(exponential(2), class = sigma)),
            cores = 4, iter = 2000, warmup = 1000, seed = 5,
            backend = "cmdstanr", silent = 2, file = "fits/b08.3.0"
)
            



```


##### 🎛 Parameters

```{r}
#| fig-width: 10
#| fig-height: 2

as_tibble(b8.3) %>% 
  rename("Water" = `b_water_z`,
         "Shade" = `b_shade_z`,,
         ) %>% 
  dplyr::select(c(`Intercept`, `Water`, `Shade`)) %>% 
  pivot_longer(cols = everything(),
               names_to = "Covariate",
               values_to = "Effect") %>% 
  ggplot(aes(x = `Effect`, y = reorder(Covariate, `Effect`))) +
  stat_halfeye(point_interval = median_qi, .width = .95,
               fill = "firebrick4") +
  labs(x = "Effect on Log GDP per capita",
       y = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank())

```


##### Model with Interaction

```{r}
b8.4 <- brm(data = tulips_clean,
            family = gaussian(),
            blooms_s ~ 1 + water_z + shade_z + water_z:shade_z,
            prior = c(
              prior(normal(0.5, 0.12), class = Intercept),
              prior(normal(0, 0.05), class = b, coef = shade_z),
              prior(normal(0, 0.05), class = b, coef = water_z),
              prior(normal(0, 0.05), class = b, coef = water_z:shade_z),
              prior(exponential(2), class = sigma)
              ),
            cores = 4, iter = 2000, warmup = 1000, seed = 5,
            backend = "cmdstanr", silent = 2, file = "fits/b08.4.2"
            )

get_prior(data = tulips_clean,
          family = gaussian(),
          blooms_s ~ 1 + water_z + shade_z + water_z:shade_z)
```


##### 🎛 Parameters

```{r}
#| fig-width: 10
#| fig-height: 2

as_tibble(b8.4) %>% 
  rename("Water" = `b_water_z`,
         "Shade" = `b_shade_z`,
         "Water:Shade" = `b_water_z:shade_z`,
         ) %>% 
  dplyr::select(c(`Intercept`, `Water`, `Shade`, `Water:Shade`)) %>% 
  pivot_longer(cols = everything(),
               names_to = "Covariate",
               values_to = "Effect") %>% 
  ggplot(aes(x = `Effect`, y = reorder(Covariate, `Effect`))) +
  stat_halfeye(point_interval = median_qi, .width = .95,
               fill = "firebrick4") +
  labs(x = "Effect on Log GDP per capita",
       y = NULL) +
  theme_bw() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        panel.grid = element_blank())

```

```{r}
#| fig-width: 7
#| fig-height: 4

figs2 <- list()
figs <- list()


# loop over values of `water_c` and plot predictions
for(s in -1:1) {

  # define the subset of the original data
  tulips_filtered <- tulips_clean[tulips_clean$shade_z == s, ]
  # defining our new data
  nd <- tibble(shade_z = s, water_z = c(-1, 1))
  # use our sampling skills, like before
  f <- 
    fitted(b8.3, 
           newdata = nd,
           summary = F) %>%
    data.frame() %>%
    set_names("-1", "1") %>% 
    mutate(row = 1:n()) %>% 
    pivot_longer(-row,
                 names_to = "water_z",
                 values_to = "blooms_s") %>% 
    mutate(water_z = as.double(water_z))

  fig <- 
    ggplot(data = tulips_filtered,
           aes(x = water_z, y = blooms_s)) +
    geom_line(data = f,
              aes(group = row),
              color = palette_pander(n = 6)[6], alpha = .005) +
    geom_point(color = palette_pander(n = 6)[6]) +
    scale_x_continuous("Water", breaks = c(-1, 0, 1)) +
    labs(title = paste("Shade", s),
         y = "Blooms") +
    coord_cartesian(xlim = c(-1, 1), 
                    ylim = c(0, 1))+
    theme_minimal()
  

  
  f2 <- 
    fitted(b8.4, 
           newdata = nd,
           summary = F) %>%
    data.frame() %>%
    set_names("-1", "1") %>% 
    mutate(row = 1:n()) %>% 
    pivot_longer(-row,
                 names_to = "water_z",
                 values_to = "blooms_s") %>% 
    mutate(water_z = as.double(water_z))

  fig2 <- 
    ggplot(data = tulips_filtered,
           aes(x = water_z, y = blooms_s)) +
    geom_line(data = f2,
              aes(group = row),
              color = palette_pander(n = 6)[6], alpha = .005) +
    geom_point(color = palette_pander(n = 6)[6]) +
    scale_x_continuous("Water", breaks = c(-1, 0, 1)) +
    labs(title = paste("Interaction Shade", s),
         y = "Blooms") +
    coord_cartesian(xlim = c(-1, 1), 
                    ylim = c(0, 1))+
    theme_minimal()
  
  # plot that joint
  
  figs[[s + 2]] <- fig
  figs2[[s + 2]] <- fig2
}


figs[[1]] + figs[[2]] + figs[[3]] + figs2[[1]] + figs2[[2]] + figs2[[3]] 


```




```{r}
# # augment the data
# points <-
#   tulips_filtered %>%
#   expand_grid(fit = c("b8.4", "b8.5")) %>%
#   mutate(x_grid = str_c("shade_cent = ", shade_cent),
#          y_grid = fit)
# 
# # redefine `nd`
# nd <- crossing(shade_cent = -1:1, 
#                water_cent = c(-1, 1))
# 
# # use `fitted()`
# set.seed(8)
# 
# rbind(fitted(b8.4, newdata = nd, summary = F, ndraws = 20),
#       fitted(b8.5, newdata = nd, summary = F, ndraws = 20)) %>%
#   # wrangle
#   data.frame() %>%
#   set_names(mutate(nd, name = str_c(shade_cent, water_cent, sep = "_")) %>% pull()) %>%
#   mutate(row = 1:n(),
#          fit = rep(c("b8.4", "b8.5"), each = n() / 2)) %>%
#   pivot_longer(-c(row:fit), values_to = "blooms_std") %>%
#   separate(name, into = c("shade_cent", "water_cent"), sep = "_") %>%
#   mutate(shade_cent = shade_cent %>% as.double(),
#          water_cent = water_cent %>% as.double()) %>%
#   # these will come in handy for `ggplot2::facet_grid()`
#   mutate(x_grid = str_c("shade_cent = ", shade_cent),
#          y_grid = fit) %>%
#   
#   # plot!
#   ggplot(aes(x = water_cent, y = blooms_std)) +
#   geom_line(aes(group = row),
#             color = palette_pander(n = 6)[6], alpha = 1/5, linewidth = 1/2) +
#   geom_point(data = points,
#              color = palette_pander(n = 6)[6]) +
#   scale_x_continuous("Water (centered)", breaks = c(-1, 0, 1)) +
#   scale_y_continuous("Blooms (standardized)", breaks = c(0, .5, 1)) +
#   ggtitle("Posterior predicted blooms") +
#   coord_cartesian(xlim = c(-1, 1),
#                   ylim = c(0, 1)) +
#   theme(strip.background = element_rect(fill = alpha(palette_pander(n = 2)[2], 1/3))) +
#   facet_grid(y_grid ~ x_grid)
```
```{r}
logit <- function(p){
  return(log(p/(1-p)))
}

invLogit <- function(l){
  return(1/(1+exp(-l)))
}
logit(.5)

invLogit(0)

logit(.07)

invLogit(-2.586689) 

invLogit(-6.23 + 1 + .7 + .6 + .7) 

```