7 Interactions

Every model so far in [McElreath’s text] 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?… To model deeper conditionality—where the importance of one predictor depends upon another predictor—we need interaction. Interaction is a kind of conditioning, a way of allowing parameters (really their posterior distributions) to be conditional on further aspects of the data. (p. 210)

7.1 Building an interaction.

“Africa is special” (p. 211). Let’s load the rugged data to see one of the reasons why.

library(rethinking)
data(rugged)
d <- rugged

And here we switch out rethinking for brms.

detach(package:rethinking, unload = T)
library(brms)
rm(rugged)

We’ll continue to use tidyverse-style syntax to wrangle the data.

library(tidyverse)

# make the log version of criterion
d <- 
  d %>%
  mutate(log_gdp = log(rgdppc_2000))

# extract countries with GDP data
dd <-
  d %>%
  filter(complete.cases(rgdppc_2000))

# split countries into Africa and not-Africa
d.A1 <-
  dd %>%
  filter(cont_africa == 1)

d.A0 <-
  dd %>%
  filter(cont_africa == 0)

The first two models predicting log_gdp are univariable.

b7.1 <-
  brm(data = d.A1, family = gaussian,
      log_gdp ~ 1 + rugged,
      prior = c(prior(normal(8, 100), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(uniform(0, 10), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4)

b7.2 <-
  update(b7.1, 
         newdata = d.A0)

In the text, McElreath more or less dared us to figure out how to make Figure 7.2. Here’s the brms-relevant data wrangling.

nd <- 
  tibble(rugged = seq(from = 0, to = 6.3, length.out = 30))

fitd_b7.1 <-
  fitted(b7.1, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

fitd_b7.2 <-
  fitted(b7.2, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

# Here we'll put both in a single data object, with `fitd_b7.1` stacked atop `fitd_b7.2`
fitd_both <-
  full_join(fitd_b7.1, fitd_b7.2) %>%
  mutate(cont_africa = rep(c("Africa", "not Africa"), each = 30))

For this chapter, we’ll take our plot theme from the ggthemes package, which you can learn more about here.

# install.packages("ggthemes", dependencies = T)
library(ggthemes)

Here’s the plot code for our version of Figure 7.2.

dd %>%
  mutate(cont_africa = ifelse(cont_africa == 1, "Africa", "not Africa")) %>%
  
  ggplot(aes(x = rugged)) +
  geom_ribbon(data = fitd_both,
              aes(ymin = Q2.5, 
                  ymax = Q97.5,
                  fill = cont_africa),
              alpha = 1/4) +
  geom_line(data = fitd_both,
              aes(y = Estimate, 
                  color = cont_africa)) +
  geom_point(aes(y = log_gdp, color = cont_africa),
             size = 2/3) +
  scale_colour_pander() +
  scale_fill_pander() +
  scale_x_continuous("Terrain Ruggedness Index", expand = c(0, 0)) +
  ylab("log GDP from year 2000") +
  theme_pander() + 
  theme(text = element_text(family = "Times"),
        legend.position = "none") +
  facet_wrap(~cont_africa)

7.1.1 Adding a dummy variable doesn’t work.

Here’s our model with all the countries, but without the cont_africa dummy.

b7.3 <-
  update(b7.1, 
         newdata = dd)

Now we’ll add the dummy.

b7.4 <-
  update(b7.3, 
         newdata = dd,
         formula = log_gdp ~ 1 + rugged + cont_africa)

And here we can compare the two models with information criteria.

waic(b7.3, b7.4)
##               WAIC    SE
## b7.3        539.48 12.99
## b7.4        475.94 14.83
## b7.3 - b7.4  63.54 14.61
loo(b7.3, b7.4)
##              LOOIC    SE
## b7.3        539.49 12.99
## b7.4        475.97 14.84
## b7.3 - b7.4  63.52 14.61

Happily, the WAIC and the LOO are in agreement. The model with the dummy, b7.4, fit the data much better. Here are the WAIC model weights.

model_weights(b7.3, b7.4,
              weights = "waic") %>% 
  round(digits = 3)
## b7.3 b7.4 
##    0    1

As in the text, almost all the weight goes to the multivariable model, b7.4. Before we can plot that model, we need to wrangle a bit.

nd <- 
  tibble(rugged      = seq(from = 0, to = 6.3, length.out = 30) %>% 
           rep(., times = 2),
         cont_africa = rep(0:1, each = 30))

fitd_b7.4 <-
  fitted(b7.4, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd) %>%
  mutate(cont_africa = ifelse(cont_africa == 1, "Africa", "not Africa"))

Behold our Figure 7.3.

dd %>%
  mutate(cont_africa = ifelse(cont_africa == 1, "Africa", "not Africa")) %>%
  
ggplot(aes(x = rugged)) +
  geom_ribbon(data = fitd_b7.4,
              aes(ymin  = Q2.5, 
                  ymax  = Q97.5,
                  fill  = cont_africa,
                  group = cont_africa),
              alpha = 1/4) +
  geom_line(data = fitd_b7.4,
              aes(y     = Estimate, 
                  color = cont_africa,
                  group = cont_africa)) +
  geom_point(aes(y = log_gdp, color = cont_africa),
             size = 2/3) +
  scale_colour_pander() +
  scale_fill_pander() +
  scale_x_continuous("Terrain Ruggedness Index", expand = c(0, 0)) +
  ylab("log GDP from year 2000") +
  theme_pander() + 
  theme(text = element_text(family = "Times"),
        legend.position  = c(.69, .94),
        legend.title     = element_blank(),
        legend.direction = "horizontal")

7.1.2 Adding a linear interaction does work.

Yes, it sure does. But before we fit, here’s the equation:

\[ \begin{eqnarray} \text{log_gdp}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha + \gamma_i \text{rugged}_i + \beta_2 \text{cont_africa}_i \\ \gamma_i & = & \beta_1 + \beta_3 \text{cont_africa}_i \end{eqnarray} \]

Fit the model.

b7.5 <-
  update(b7.4,
         formula = log_gdp ~ 1 + rugged*cont_africa)

For kicks, we’ll just use the LOO to compare the last three models.

loo(b7.3, b7.4, b7.5,
    cores = 4)  # We can speed things up using the `cores` argument
##              LOOIC    SE
## b7.3        539.49 12.99
## b7.4        475.97 14.84
## b7.5        469.30 14.57
## b7.3 - b7.4  63.52 14.61
## b7.3 - b7.5  70.20 14.68
## b7.4 - b7.5   6.68  5.86

And we can weigh the models based on the LOO rather than the WAIC, too.

model_weights(b7.3, b7.4, b7.5,
              weights = "loo") %>% 
  round(digits = 3)
##  b7.3  b7.4  b7.5 
## 0.000 0.034 0.966

7.1.2.1 Overthinking: Conventional form of interaction.

The conventional equation for the interaction model might look like:

\[ \begin{eqnarray} \text{log_gdp}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha + \beta_1 \text{rugged}_i + \beta_2 \text{cont_africa}_i + \beta_3 \text{rugged}_i \times \text{cont_africa}_i \end{eqnarray} \]

Instead of the y ~ 1 + x1*x2 approach, which will work fine with brm(), you can use this more explicit syntax.

b7.5b <-
  update(b7.5,
         formula = log_gdp ~ 1 + rugged + cont_africa + rugged:cont_africa)

From here on, I will default to this style of syntax for interactions.

Since this is the same model, it yields the same information criteria estimates within simulation error.

waic(b7.5, b7.5b)
##                WAIC    SE
## b7.5         469.15 14.55
## b7.5b        469.28 14.53
## b7.5 - b7.5b  -0.13  0.14
loo(b7.5, b7.5b)
##               LOOIC    SE
## b7.5         469.30 14.57
## b7.5b        469.42 14.55
## b7.5 - b7.5b  -0.12  0.14

7.1.3 Plotting the interaction.

Here’s our prep work for the figure.

fitd_b7.5 <-
  fitted(b7.5, newdata = nd) %>%  # we can use the same `nd` data from last time
  as_tibble() %>%
  bind_cols(nd) %>%
  mutate(cont_africa = ifelse(cont_africa == 1, "Africa", "not Africa"))

And here’s the code for our version of Figure 7.4.

dd %>%
  mutate(cont_africa = ifelse(cont_africa == 1, "Africa", "not Africa")) %>%
  
ggplot(aes(x = rugged)) +
  geom_ribbon(data = fitd_b7.5,
              aes(ymin  = Q2.5, 
                  ymax  = Q97.5,
                  fill  = cont_africa,
                  group = cont_africa),
              alpha = 1/4) +
  geom_line(data = fitd_b7.5,
              aes(y     = Estimate, 
                  color = cont_africa,
                  group = cont_africa)) +
  geom_point(aes(y = log_gdp, color = cont_africa),
             size = 2/3) +
  scale_colour_pander() +
  scale_fill_pander() +
  scale_x_continuous("Terrain Ruggedness Index", expand = c(0, 0)) +
  ylab("log GDP from year 2000") +
  theme_pander() + 
  theme(text = element_text(family = "Times"),
        legend.position = "none") +
  facet_wrap(~cont_africa)

7.1.4 Interpreting an interaction estimate.

Interpreting interaction estimates is tricky. It’s trickier than interpreting ordinary estimates. And for this reason, I usually advise against trying to understand an interaction from tables of numbers along. Plotting implied predictions does far more for both our own understanding and for our audience’s.

7.1.4.1 Parameters change meaning.

In a simple linear regression with no interactions, each coefficient says how much the average outcome, \(\mu\), changes when the predictor changes by one unit. And since all of the parameters have independent influences on the outcome, there’s no trouble in interpreting each parameter separately. Each slope parameter gives us a direct measure of each predictor variable’s influence.

Interaction models ruin this paradise. (p. 220)

Return the parameter estimates.

posterior_summary(b7.5)
##                          Estimate  Est.Error          Q2.5         Q97.5
## b_Intercept             9.1823382 0.13894772    8.90659529    9.45228293
## b_rugged               -0.1832468 0.07735517   -0.33161340   -0.03199895
## b_cont_africa          -1.8407667 0.22221419   -2.25894413   -1.40785674
## b_rugged:cont_africa    0.3433508 0.13148168    0.08390145    0.59336056
## sigma                   0.9514394 0.05071302    0.85997921    1.05643893
## lp__                 -244.3853233 1.53536617 -248.18162453 -242.34387796

“Since \(\gamma\) (gamma) doesn’t appear in this table—it wasn’t estimated—we have to compute it ourselves” (p. 221). Like in the text, we’ll do so first by working with the point estimates.

# within Africa
fixef(b7.5)[2, 1] + fixef(b7.5)[4, 1] * 1
## [1] 0.160104
# outside Africa
fixef(b7.5)[2, 1] + fixef(b7.5)[4, 1] * 0
## [1] -0.1832468

7.1.4.2 Incorporating uncertainty.

To get some idea of the uncertainty around those \(\gamma\) values, we’ll need to use the whole posterior. Since \(\gamma\) depends upon parameters, and those parameters have a posterior distribution, \(\gamma\) must also have a posterior distribution. Read the previous sentence again a few times. It’s one of the most important concepts in processing Bayesian model fits. Anything calculated using parameters has a distribution. (p. 212)

Like McElreath, we’ll avoid integral calcus in favor of working with the posterior_samples().

post <- posterior_samples(b7.5) 

post %>%
  transmute(gamma_Africa    = b_rugged + `b_rugged:cont_africa`,
            gamma_notAfrica = b_rugged) %>%
  gather(key, value) %>%
  group_by(key) %>%
  summarise(mean = mean(value))
## # A tibble: 2 x 2
##   key               mean
##   <chr>            <dbl>
## 1 gamma_Africa     0.160
## 2 gamma_notAfrica -0.183

And here is our version of Figure 7.5.

post %>%
  transmute(gamma_Africa    = b_rugged + `b_rugged:cont_africa`,
            gamma_notAfrica = b_rugged) %>%
  gather(key, value) %>%
  
  ggplot(aes(x = value, group = key, color = key, fill = key)) +
  geom_density(alpha = 1/4) +
  scale_colour_pander() +
  scale_fill_pander() +
  scale_x_continuous(expression(gamma), expand = c(0, 0)) +
  scale_y_continuous(NULL, breaks = NULL) +
  ggtitle("Terraine Ruggedness slopes",
          subtitle = "Blue = African nations, Green = others") +
  theme_pander() + 
  theme(text = element_text(family = "Times"),
        legend.position = "none")

What proportion of these differences is below zero?

post %>%
  mutate(gamma_Africa    = b_rugged + `b_rugged:cont_africa`,
         gamma_notAfrica = b_rugged) %>% 
  mutate(diff            = gamma_Africa -gamma_notAfrica) %>%
  summarise(Proportion_of_the_difference_below_0 = sum(diff < 0) / length(diff))
##   Proportion_of_the_difference_below_0
## 1                              0.00675

7.2 Symmetry of the linear interaction.

Consider for example the GDP and terrain ruggedness problem. The interaction there has two equally valid phrasings.

  1. How much does the influence of ruggedness (on GDP) depend upon whether the nation is in Africa?
  2. How much does the influence of being in Africa (on GDP) depend upon ruggedness?

While these two possibilities sound different to most humans, your golem thinks they are identical. (p. 223)

7.2.1 Buridan’s interaction.

Recall the original equation.

\[ \begin{eqnarray} \text{log_gdp}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha + \gamma_i \text{rugged}_i + \beta_2 \text{cont_africa}_i \\ \gamma_i & = & \beta_1 + \beta_3 \text{cont_africa}_i \end{eqnarray} \]

Next McElreath replaced \(\gamma_i\) with the expression for \(\mu_i\).

\[ \begin{eqnarray} \mu_i & = & \alpha + (\beta_1 + \beta_3 \text{cont_africa}_i) \times \text{rugged}_i + \beta_2 \text{cont_africa}_i \\ & = & \alpha + \beta_1 \text{rugged}_i + \beta_3 \text{rugged}_i \times \text{cont_africa}_i + \beta_2 \text{cont_africa}_i \end{eqnarray} \]

And new we’ll factor together the terms containing \(\text{cont_africa}_i\).

\[\mu_i = \alpha + \beta_1 \text{rugged}_i + \underbrace{(\beta_2 + \beta_3 \text{rugged}_i)}_G \times \text{cont_africa}_i\]

And just as in the text, our \(G\) term looks a lot like the original \(\gamma_i\) term.

7.2.2 Africa depends upon ruggedness.

Here is our version of McElreath’s Figure 7.6.

# new predictor data for `fitted()`
nd <- 
  tibble(rugged      = rep(range(dd$rugged), times = 2),
         cont_africa = rep(0:1,              each = 2))

# `fitted()`
fitd_b7.5 <-
  fitted(b7.5, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd) %>% 
  mutate(ox = rep(c(-0.05, 0.05), times = 2))

# augment the `dd` data a bit
dd %>% 
  mutate(ox = ifelse(rugged > median(rugged), 0.05, -0.05),
         cont_africa = cont_africa + ox) %>%
  select(cont_africa, everything()) %>%
  
  # plot
  ggplot(aes(x = cont_africa)) +
  geom_ribbon(data = fitd_b7.5,
              aes(ymin  = Q2.5, 
                  ymax  = Q97.5,
                  fill  = factor(ox),
                  group = factor(ox)),
              alpha = 1/4) +
  geom_line(data = fitd_b7.5,
              aes(y        = Estimate,
                  color    = factor(ox),
                  group    = factor(ox),
                  linetype = factor(ox))) +
  geom_point(aes(y = log_gdp, color = factor(ox)),
             alpha = 1/2, shape = 1) +
  scale_colour_pander() +
  scale_fill_pander() +
  scale_x_continuous("Continent", breaks = 0:1, 
                     labels = c("other", "Africa")) +
  coord_cartesian(xlim = c(-.2, 1.2)) +
  ylab("log GDP from year 2000") +
  theme_pander() +
  theme(text = element_text(family = "Times"),
        legend.position = "none")

7.3 Continuous interactions

Though continuous interactions can be more challenging to interpret, they’re just as easy to fit as interactions including dummies.

7.3.1 The data.

Look at the tulips.

library(rethinking)
data(tulips)
d <- tulips
str(d)
## 'data.frame':    27 obs. of  4 variables:
##  $ bed   : Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 2 ...
##  $ water : int  1 1 1 2 2 2 3 3 3 1 ...
##  $ shade : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ blooms: num  0 0 111 183.5 59.2 ...

7.3.2 The un-centered models.

The likelihoods for the next two models are

\[ \begin{eqnarray} \text{blooms}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha + \beta_1 \text{water}_i + \beta_2 \text{shade}_i \\ \alpha & \sim & \text{Normal} (0, 100) \\ \beta_1 & \sim & \text{Normal} (0, 100) \\ \beta_2 & \sim & \text{Normal} (0, 100) \\ \sigma & \sim & \text{Uniform} (0, 100) \\ \end{eqnarray} \] and

\[ \begin{eqnarray} \text{blooms}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha + \beta_1 \text{water} + \beta_2 \text{shade}_i + \beta_3 \text{water}_i \times \text{shade}_i\\ \alpha & \sim & \text{Normal} (0, 100) \\ \beta_1 & \sim & \text{Normal} (0, 100) \\ \beta_2 & \sim & \text{Normal} (0, 100) \\ \beta_3 & \sim & \text{Normal} (0, 100) \\ \sigma & \sim & \text{Uniform} (0, 100) \\ \end{eqnarray} \]

Load brms.

detach(package:rethinking, unload = T)
library(brms)
rm(tulips)

Here we continue with McElreath’s very-flat priors for the multivariable and interaction models.

b7.6 <-
  brm(data = d, family = gaussian,
      blooms ~ 1 + water + shade,
      prior = c(prior(normal(0, 100), class = Intercept),
                prior(normal(0, 100), class = b),
                prior(uniform(0, 100), class = sigma)),
      iter = 2000, warmup = 1000, cores = 4, chains = 4)
## Warning: There were 44 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
b7.7 <- 
  update(b7.6, 
         formula = blooms ~ 1 + water + shade + water:shade)
## Warning: There were 6 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: Examine the pairs() plot to diagnose sampling problems

Much like in the text, these models yielded divergent transitions. Here, we’ll try to combat them by following Stan’s advice and “[increase] adapt_delta above 0.8.” While we’re at it, we’ll put better priors on \(\sigma\).

b7.6 <-
  update(b7.6,
         prior = c(prior(normal(0, 100), class = Intercept),
                   prior(normal(0, 100), class = b),
                   prior(cauchy(0, 10), class = sigma)),
         control = list(adapt_delta = 0.9))

b7.7 <- 
  update(b7.6, 
         formula = blooms ~ 1 + water + shade + water:shade)

Increasing adapt_delta did the trick. Instead of coeftab(), we can also use posterior_summary(), which gets us most of the way there.

posterior_summary(b7.6) %>% round(digits = 2)
##             Estimate Est.Error    Q2.5   Q97.5
## b_Intercept    60.98     41.53  -22.32  140.61
## b_water        74.12     14.22   46.06  101.96
## b_shade       -40.98     14.70  -69.54  -11.86
## sigma          61.58      9.07   46.80   81.61
## lp__         -169.76      1.49 -173.40 -167.86
posterior_summary(b7.7) %>% round(digits = 2)
##               Estimate Est.Error    Q2.5   Q97.5
## b_Intercept    -106.32     64.25 -224.44   26.34
## b_water         159.54     29.81   96.06  214.18
## b_shade          43.50     29.84  -17.30   99.20
## b_water:shade   -42.98     13.71  -68.31  -14.15
## sigma            49.85      7.57   37.48   66.87
## lp__           -170.69      1.70 -174.88 -168.39

This is an example where HMC yielded point estimates notably different from MAP. However, look at the size of those posterior standard deviations (i.e., ‘Est.Error’ column)! The MAP estimates are well within a fraction of those \(SD\)s.

Anyway, let’s look at WAIC.

waic(b7.6, b7.7)
##               WAIC   SE
## b7.6        304.02 7.63
## b7.7        293.84 8.06
## b7.6 - b7.7  10.18 5.31

Why not compute the WAIC weights?

model_weights(b7.6, b7.7, weights = "waic")
##        b7.6        b7.7 
## 0.006122731 0.993877269

As in the text, almost all the weight went to the interaction model, b7.7.

7.3.3 Center and re-estimate.

Here’s a tidyverse way to center the predictors.

d <-
  d %>%
  mutate(shade_c = shade - mean(shade),
         water_c = water - mean(water))

Now refit the models with our shiny new centered predictors.

b7.8 <-
  brm(data = d, family = gaussian,
      blooms ~ 1 + water_c + shade_c,
      prior = c(prior(normal(130, 100), class = Intercept),
                prior(normal(0, 100), class = b),
                prior(cauchy(0, 10), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9))

b7.9 <- 
  update(b7.8, 
         formula = blooms ~ 1 + water_c + shade_c + water_c:shade_c)
posterior_summary(b7.8) %>% round(digits = 2)
##             Estimate Est.Error    Q2.5   Q97.5
## b_Intercept   128.97     11.65  106.14  152.10
## b_water_c      74.14     14.00   46.46  101.42
## b_shade_c     -40.86     14.46  -69.78  -12.13
## sigma          61.27      8.65   47.27   81.31
## lp__         -168.85      1.45 -172.60 -167.03
posterior_summary(b7.9) %>% round(digits = 2)
##                   Estimate Est.Error    Q2.5   Q97.5
## b_Intercept         128.87      9.34  110.17  146.90
## b_water_c            74.71     11.55   51.74   97.84
## b_shade_c           -41.01     11.88  -63.99  -17.91
## b_water_c:shade_c   -51.77     14.45  -80.74  -23.92
## sigma                49.65      7.48   37.57   66.62
## lp__               -168.56      1.70 -172.72 -166.30

And okay fine, if you really want a coeftab()-like summary, here’s a grotesque way to do it.

# First, we reformat `b7.8`
posterior_summary(b7.8) %>% 
  data.frame() %>% 
  rownames_to_column() %>% 
  # Now we bind to it the rows of `b7.9`
  bind_rows(
    posterior_summary(b7.9) %>% 
      data.frame() %>%
      rownames_to_column()
  ) %>% 
  # Miscellaneous wrangling
  mutate(model = c(rep("b7.8", times = 5),
                   rep("b7.9", times = 6))) %>% 
  select(rowname, Estimate, model) %>% 
  filter(rowname != "lp__") %>% 
  spread(key = model, value = Estimate) %>% 
  rename(parameter = rowname) %>% 
  mutate_if(is.double, round, digits = 2)
##           parameter   b7.8   b7.9
## 1       b_Intercept 128.97 128.87
## 2         b_shade_c -40.86 -41.01
## 3         b_water_c  74.14  74.71
## 4 b_water_c:shade_c     NA -51.77
## 5             sigma  61.27  49.65

Anyway, centering helped a lot. Now, not only do the results in the text match up better than those from Stan, but the ‘Est.Error’ values are uniformly smaller.

7.3.3.1 Estimation worked better.

Nothing to add, here.

7.3.3.2 Estimates changed less across models.

k <- fixef(b7.7)
k[1] + k[2]*2 + k[3]*2 + k[4]*2*2
## [1] 127.8318
k <- fixef(b7.9)
k[1] + k[2]*0 + k[3]*0 + k[4]*0*0
## [1] 128.8718
print(b7.9)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: blooms ~ water_c + shade_c + water_c:shade_c 
##    Data: d (Number of observations: 27) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept         128.87      9.34   110.17   146.90       4000 1.00
## water_c            74.71     11.55    51.74    97.84       4000 1.00
## shade_c           -41.01     11.88   -63.99   -17.91       4000 1.00
## water_c:shade_c   -51.77     14.45   -80.74   -23.92       4000 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma    49.65      7.48    37.57    66.62       4000 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

7.3.4 Plotting implied predictions.

Now we’re ready for the bottom row of Figure 7.7. Here’s our variation on McElreath’s tryptych loop code, adjusted for brms and ggplot2.

# loop over values of water.c and plot predictions
shade_seq <- -1:1

for(w in -1:1){
  # defining the subset of the original data
  dt <- d[d$water_c == w, ]
  # defining our new data
  nd <- tibble(water_c = w, shade_c = shade_seq)
  # using our sampling skills, like before
  fitd_7.9 <- 
    fitted(b7.9, newdata = nd) %>%
    as_tibble() %>%
    bind_cols(nd)
  
  # specifying our custom plot
  fig <- 
    ggplot() +
    geom_ribbon(data = fitd_7.9, 
                aes(x    = shade_c,
                    ymin = Q2.5,
                    ymax = Q97.5), 
                fill = "#CC79A7", alpha = 1/5) +
    geom_line(data = fitd_7.9, 
              aes(x = shade_c, y = Estimate), 
              color = "#CC79A7") +
    geom_point(data = dt, aes(x = shade_c, y = blooms),
               shape = 1, color = "#CC79A7") +
    coord_cartesian(xlim = range(d$shade_c), 
                    ylim = range(d$blooms)) +
    scale_x_continuous("Shade (centered)", breaks = c(-1, 0, 1)) +
    labs("Blooms", 
         title = paste("Water (centered) =", w)) +
    theme_pander() + 
    theme(text = element_text(family = "Times"))
  
  # plotting that joint
  plot(fig)
}

But we don’t necessarily need a loop. We can achieve all of McElreath’s Figure 7.7 with fitted(), some data wrangling, and a little help from ggplot2::facet_grid().

# fitted() for model b7.8
fitted(b7.8) %>%
  as_tibble() %>%
  # adding fitted() for model b7.9
  bind_rows(
    fitted(b7.9) %>% 
      as_tibble()
  ) %>% 
  # We'll want to index the models
  mutate(fit  = rep(c("b7.8", "b7.9"), each = 27)) %>% 
  # Here we add the data, `d`
  bind_cols(bind_rows(d, d)) %>% 
  # These will come in handy for `ggplot2::facet_grid()`
  mutate(x_grid = paste("water_c =", water_c),
         y_grid = paste("model: ", fit)) %>% 
  
  ggplot(aes(x = shade_c)) +
  geom_ribbon(aes(ymin = Q2.5,
                  ymax = Q97.5), 
              fill = "#CC79A7", alpha = 1/5) +
  geom_line(aes(y = Estimate), 
            color = "#CC79A7") +
  geom_point(aes(y = blooms, group = x_grid), 
             shape = 1, color = "#CC79A7") +
  coord_cartesian(xlim = range(d$shade_c), 
                  ylim = range(d$blooms)) +
  scale_x_continuous("Shade (centered)", breaks = c(-1, 0, 1)) +
  ylab("Blooms") +
  theme_pander() + 
  theme(text = element_text(family = "Times"),
        panel.background = element_rect(color = "black")) +
  facet_grid(y_grid ~ x_grid)

7.4 Interactions in design formulas

The brms syntax generally follows the design formulas typical of lm(). Hopefully this is all old hat.

7.5 Summary Bonus: marginal_effects()

The brms package includes the marginal_effects() function as a convenient way to look at simple effects and two-way interactions. Recall the simple univariable model, b7.3:

b7.3$formula
## log_gdp ~ 1 + rugged

We can look at the regression line and its percentile-based intervals like so:

marginal_effects(b7.3)

If we nest marginal_effects() within plot() with a points = T argument, we can add the original data to the figure.

plot(marginal_effects(b7.3), points = T)

We can further customize the plot. For example, we can replace the intervals with a spaghetti plot. While we’re at it, we can use point_args to adjust the geom_jitter() parameters.

plot(marginal_effects(b7.3,
                      spaghetti = T, nsamples = 200),
     points = T,
     point_args = c(alpha = 1/2, size = 1))

With multiple predictors, things get more complicated. Consider our multivariable, non-interaction model, b7.4.

b7.4$formula
## log_gdp ~ rugged + cont_africa
marginal_effects(b7.4)

We got one plot for each predictor, controlling the other predictor at zero. Note how the plot for cont_africa treated it as a continuous variable. This is because the variable was saved as an integer in the original data set:

b7.4$data %>% 
  glimpse()
## Observations: 170
## Variables: 3
## $ log_gdp     <dbl> 7.492609, 8.216929, 9.933263, 9.407032, 7.792343, 9.212541, 10.143191, 10.2...
## $ rugged      <dbl> 0.858, 3.427, 0.769, 0.775, 2.688, 0.006, 0.143, 3.513, 1.672, 1.780, 0.388...
## $ cont_africa <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, ...

One way to fix that is to adjust the data set and refit the model.

d_factor <-
  b7.4$data %>% 
  mutate(cont_africa = factor(cont_africa))

b7.4_factor <- update(b7.4, newdata = d_factor)

Using the update() syntax often speeds up the re-fitting process.

marginal_effects(b7.4_factor)

Now our second marginal plot more clearly expresses the cont_africa predictor as categorical.

Things get more complicated with the interaction model, b7.5.

b7.5$formula
## log_gdp ~ rugged + cont_africa + rugged:cont_africa
marginal_effects(b7.5)

The marginal_effects() function defaults to expressing interactions such that the first variable in the term–in this case, rugged–is on the x axis and the second variable in the term–cont_africa, treated as an integer–is depicted in three lines corresponding its mean and its mean +/- one standard deviation. This is great for continuous variables, but incoherent for categorical ones. The fix is, you guessed it, to refit the model after adjusting the data.

d_factor <-
  b7.5$data %>% 
  mutate(cont_africa = factor(cont_africa))

b7.5_factor <- update(b7.5, newdata = d_factor)

Just for kicks, we’ll use probs = c(.25, .75) to return 50% intervals, rather than the conventional 95%.

marginal_effects(b7.5_factor,
                 probs = c(.25, .75))

With the effects argument, we can just return the interaction effect, which is where all the action’s at. While we’re at it, we’ll use plot() to change some of the settings.

plot(marginal_effects(b7.5_factor,
                      effects = "rugged:cont_africa", 
                      spaghetti = T, nsamples = 150),
     points = T,
     point_args = c(alpha = 2/3, size = 1), mean = F)

Note, the ordering of the variables matters for the interaction term. Consider our interaction model for the tulips data.

b7.9$formula
## blooms ~ water_c + shade_c + water_c:shade_c

The plot tells a slightly different story, depending on whether you specify effects = "shade_c:water_c" or effects = "water_c:shade_c".

plot(marginal_effects(b7.9, 
                      effects = "shade_c:water_c"),
     points = T)

plot(marginal_effects(b7.9, 
                      effects = "water_c:shade_c"),
     points = T)

One might want to evaluate the effects of the second term in the interaction–water_c, in this case–at values other than the mean and the mean +/- one standard deviation. When we reproduced the bottom row of Figure 7.7, we expressed the interaction based on values -1, 0, and 1 for water_c. We can do that, here, by using the int_conditions argument. It expects a list, so we’ll put our desired water_c values in just that.

ic <- 
  list(water.c = c(-1, 0, 1))

plot(marginal_effects(b7.9, 
                      effects = "shade_c:water_c",
                      int_conditions = ic),
     points = T)

Session info

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggthemes_3.5.0     forcats_0.3.0      stringr_1.3.1      dplyr_0.7.6        purrr_0.2.5       
##  [6] readr_1.1.1        tidyr_0.8.1        tibble_1.4.2       tidyverse_1.2.1    brms_2.5.0        
## [11] Rcpp_0.12.18       rstan_2.17.3       StanHeaders_2.17.2 ggplot2_3.0.0     
## 
## loaded via a namespace (and not attached):
##   [1] pacman_0.4.6              utf8_1.1.4                ggstance_0.3             
##   [4] tidyselect_0.2.4          htmlwidgets_1.2           munsell_0.5.0            
##   [7] codetools_0.2-15          nleqslv_3.3.2             DT_0.4                   
##  [10] miniUI_0.1.1.1            withr_2.1.2               Brobdingnag_1.2-5        
##  [13] colorspace_1.3-2          highr_0.7                 knitr_1.20               
##  [16] rstudioapi_0.7            stats4_3.5.1              Rttf2pt1_1.3.7           
##  [19] bayesplot_1.6.0           labeling_0.3              mnormt_1.5-5             
##  [22] bridgesampling_0.4-0      rprojroot_1.3-2           coda_0.19-1              
##  [25] xfun_0.3                  R6_2.2.2                  markdown_0.8             
##  [28] HDInterval_0.2.0          reshape_0.8.7             assertthat_0.2.0         
##  [31] promises_1.0.1            scales_0.5.0              beeswarm_0.2.3           
##  [34] gtable_0.2.0              rlang_0.2.1               extrafontdb_1.0          
##  [37] lazyeval_0.2.1            broom_0.4.5               inline_0.3.15            
##  [40] yaml_2.1.19               reshape2_1.4.3            abind_1.4-5              
##  [43] modelr_0.1.2              threejs_0.3.1             crosstalk_1.0.0          
##  [46] backports_1.1.2           httpuv_1.4.4.2            rsconnect_0.8.8          
##  [49] extrafont_0.17            tools_3.5.1               bookdown_0.7             
##  [52] psych_1.8.4               RColorBrewer_1.1-2        ggridges_0.5.0           
##  [55] plyr_1.8.4                base64enc_0.1-3           progress_1.2.0           
##  [58] prettyunits_1.0.2         zoo_1.8-2                 LaplacesDemon_16.1.1     
##  [61] haven_1.1.2               magrittr_1.5              colourpicker_1.0         
##  [64] mvtnorm_1.0-8             matrixStats_0.54.0        hms_0.4.2                
##  [67] shinyjs_1.0               mime_0.5                  evaluate_0.10.1          
##  [70] arrayhelpers_1.0-20160527 xtable_1.8-2              shinystan_2.5.0          
##  [73] readxl_1.1.0              gridExtra_2.3             rstantools_1.5.0         
##  [76] compiler_3.5.1            maps_3.3.0                crayon_1.3.4             
##  [79] htmltools_0.3.6           later_0.7.3               lubridate_1.7.4          
##  [82] MASS_7.3-50               Matrix_1.2-14             cli_1.0.0                
##  [85] bindr_0.1.1               igraph_1.2.1              pkgconfig_2.0.1          
##  [88] foreign_0.8-70            xml2_1.2.0                svUnit_0.7-12            
##  [91] dygraphs_1.1.1.5          vipor_0.4.5               rvest_0.3.2              
##  [94] digest_0.6.15             rmarkdown_1.10            cellranger_1.1.0         
##  [97] shiny_1.1.0               gtools_3.8.1              nlme_3.1-137             
## [100] jsonlite_1.5              bindrcpp_0.2.2            mapproj_1.2.6            
## [103] viridisLite_0.3.0         pillar_1.2.3              lattice_0.20-35          
## [106] loo_2.0.0                 httr_1.3.1                glue_1.2.0               
## [109] xts_0.10-2                shinythemes_1.1.1         pander_0.6.2             
## [112] stringi_1.2.3