13 Adventures in Covariance

In this chapter, you’ll see how to… specify varying slopes in combination with the varying intercepts of the previous chapter. This will enable pooling that will improve estimates of how different units respond to or are influenced by predictor variables. It will also improve estimates of intercepts, by borrowing information across parameter types. Essentially, varying slopes models are massive interaction machines. They allow every unit in the data to have its own unique response to any treatment or exposure or event, while also improving estimates via pooling. When the variation in slopes is large, the average slope is of less interest. Sometimes, the pattern of variation in slopes provides hints about omitted variables that explain why some units respond more or less. We’ll see an example in this chapter.

The machinery that makes such complex varying effects possible will be used later in the chapter to extend the varying effects strategy to more subtle model types, including the use of continuous categories, using Gaussian process. (p. 388)

13.1 Varying slopes by construction

How should the robot pool information across intercepts and slopes? By modeling the joint population of intercepts and slopes, which means by modeling their covariance. In conventional multilevel models, the device that makes this possible is a joint multivariate Gaussian distribution for all of the varying effects, both intercepts and slopes. So instead of having two independent Gaussian distributions of intercepts and of slopes, the robot can do better by assigning a two-dimensional Gaussian distribution to both the intercepts (first dimension) and the slopes (second dimension). (p. 389)

In the Rethinking: Why Gaussian? box, McElreath discussed how researchers might use other multivariate distributions to model multiple random effects. The only one he named as an alternative to the Gaussian was the multivariate Student’s \(t\). As it turns out, brms does currently allow users to use multivariate Student’s \(t\) in this way. For details, check out this discussion in the brms GitHub page. Bürkner’s exemplar syntax from his comment on May 13, 2018, was y ~ x + (x | gr(g, dist = "student")). I haven’t experimented with this, but if you do, do consider commenting on how it went.

13.1.1 Simulate the population.

If you follow this section closely, it’s a great template for simulating multilevel code for any of your future projects. You might think of this as an alternative to a frequentist power analysis. Vourre has done some nice work along these lines, too.

a       <-  3.5  # average morning wait time
b       <- -1    # average difference afternoon wait time
sigma_a <-  1    # std dev in intercepts
sigma_b <-  0.5  # std dev in slopes
rho     <- -.7   # correlation between intercepts and slopes

# The next three lines of code simply combine the terms, above
mu <- c(a, b)

cov_ab <- sigma_a * sigma_b * rho
sigma  <- matrix(c(sigma_a^2, cov_ab, 
                   cov_ab, sigma_b^2), ncol = 2)

If you haven’t used matirx() before, you might get a sense of the elements like so.

matrix(c(1, 2, 
         3, 4), nrow = 2, ncol = 2)
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

This next block of code will finally yield our café data.

library(tidyverse)

sigmas <- c(sigma_a, sigma_b)          # standard deviations
rho    <- matrix(c(1, rho,             # correlation matrix
                   rho, 1), nrow = 2)

# now matrix multiply to get covariance matrix
sigma <- diag(sigmas) %*% rho %*% diag(sigmas)

# how many cafes would you like?
n_cafes <- 20

set.seed(5)  # used to replicate example
vary_effects <- 
  MASS::mvrnorm(n_cafes, mu, sigma)

vary_effects <-
  vary_effects %>% 
  as_tibble() %>% 
  rename(a_cafe = V1,
         b_cafe = V2)

head(vary_effects)
## # A tibble: 6 x 2
##   a_cafe b_cafe
##    <dbl>  <dbl>
## 1   4.22 -1.61 
## 2   2.01 -0.752
## 3   4.57 -1.95 
## 4   3.34 -1.19 
## 5   1.70 -0.586
## 6   4.13 -1.14

Let’s make sure we’re keeping this all straight. a_cafe = our café-specific intercepts; b_cafe = our café-specific slopes. These aren’t the actual data, yet. But at this stage, it might make sense to ask What’s the distribution of a_cafe and b_cafe? Our variant of Figure 13.2 contains the answer.

For our plots in this chapter, we’ll use a custom theme. The color palette will come from the “pearl_earring” palette of the dutchmasters package. You can learn more about the original painting, Vermeer’s Girl with a Pearl Earring, here.

# devtools::install_github("EdwinTh/dutchmasters")
library(dutchmasters)

dutchmasters$pearl_earring
##         red(lips)              skin      blue(scarf1)      blue(scarf2)      white(colar) 
##         "#A65141"         "#E7CDC2"         "#80A0C7"         "#394165"         "#FCF9F0" 
##       gold(dress)      gold(dress2) black(background)      grey(scarf3)    yellow(scarf4) 
##         "#B1934A"         "#DCA258"         "#100F14"         "#8B9DAF"         "#EEDA9D" 
##                   
##         "#E8DCCF"

We’ll name our custom theme theme_pearl_earring.

theme_pearl_earring <-
  theme(text       = element_text(color = "#E8DCCF", family = "Courier"),
        strip.text = element_text(color = "#E8DCCF", family = "Courier"),
        axis.text  = element_text(color = "#E8DCCF"),
        axis.ticks = element_line(color = "#E8DCCF"),
        line       = element_line(color = "#E8DCCF"),
        plot.background   = element_rect(fill = "#100F14", color = "transparent"),
        panel.background  = element_rect(fill = "#100F14", color = "#E8DCCF"),
        strip.background  = element_rect(fill = "#100F14", color = "transparent"),
        panel.grid = element_blank(),
        legend.background = element_rect(fill = "#100F14", color = "transparent"),
        legend.key        = element_rect(fill = "#100F14", color = "transparent"),
        axis.line = element_blank())

Now we’re ready to plot Figure 13.2.

vary_effects %>% 
  ggplot(aes(x = a_cafe, y = b_cafe)) +
  geom_point(color = "#80A0C7") +
  geom_rug(color = "#8B9DAF", size = 1/7) +
  theme_pearl_earring

Again, these are not “data.” This is a distribution of parameters.

13.1.2 Simulate observations.

Here we put those simulated parameters to use.

n_visits <- 10
sigma    <-  0.5  # std dev within cafes

set.seed(5)  # used to replicate example
d <-
  vary_effects %>% 
  mutate(cafe      = 1:n_cafes) %>% 
  expand(nesting(cafe, a_cafe, b_cafe), visit = 1:n_visits) %>% 
  mutate(afternoon = rep(0:1, times = n() / 2)) %>% 
  mutate(mu        = a_cafe + b_cafe * afternoon) %>% 
  mutate(wait      = rnorm(n = n(), mean = mu, sd = sigma))

We might peek at the data.

d %>%
  head()
## # A tibble: 6 x 7
##    cafe a_cafe b_cafe visit afternoon    mu  wait
##   <int>  <dbl>  <dbl> <int>     <int> <dbl> <dbl>
## 1     1   4.22  -1.61     1         0  4.22  3.80
## 2     1   4.22  -1.61     2         1  2.61  3.31
## 3     1   4.22  -1.61     3         0  4.22  3.60
## 4     1   4.22  -1.61     4         1  2.61  2.65
## 5     1   4.22  -1.61     5         0  4.22  5.08
## 6     1   4.22  -1.61     6         1  2.61  2.31

Now we’ve finally simulated our data, we are ready to make our version of Figure 13.1, from way back on page 388.

d %>%
  mutate(afternoon = ifelse(afternoon == 0, "M", "A"),
         day       = rep(rep(1:5, each = 2), times = n_cafes)) %>%
  filter(cafe %in% c(3, 5)) %>%
  mutate(cafe = ifelse(cafe == 3, "cafe #3", "cafe #5")) %>%
  
  ggplot(aes(x = visit, y = wait, group = day)) +
  geom_point(aes(color = afternoon), size = 2) +
  geom_line(color = "#8B9DAF") +
  scale_color_manual(values = c("#80A0C7", "#EEDA9D")) +
  scale_x_continuous(breaks = 1:10,
                     labels = rep(c("M", "A"), times = 5)) +
  coord_cartesian(ylim = 0:8) +
  labs(x = NULL, y = "wait time in minutes") +
  theme_pearl_earring +
  theme(legend.position = "none",
        axis.ticks.x    = element_blank()) +
  facet_wrap(~cafe, ncol = 1)

13.1.3 The varying slopes model.

The statistical formula for our varying-slopes model follows the form

\[ \begin{eqnarray} \text{wait}_i & \sim & \text{Normal} (\mu_i, \sigma) \\ \mu_i & = & \alpha_{\text{cafe}_i} + \beta_{\text{cafe}_i} \text{afternoon}_i \\ \begin{bmatrix} \alpha_\text{cafe} \\ \beta_\text{cafe} \end{bmatrix} & \sim & \text{MVNormal} \bigg (\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \mathbf{S} \bigg ) \\ \mathbf S & = & \begin{pmatrix} \sigma_\alpha & 0 \\ 0 & \sigma_\beta \end{pmatrix} \mathbf R \begin{pmatrix} \sigma_\alpha & 0 \\ 0 & \sigma_\beta \end{pmatrix} \\ \alpha & \sim & \text{Normal} (0, 10) \\ \beta & \sim & \text{Normal} (0, 10) \\ \sigma & \sim & \text{HalfCauchy} (0, 1) \\ \sigma_\alpha & \sim & \text{HalfCauchy} (0, 1) \\ \sigma_\beta & \sim & \text{HalfCauchy} (0, 1) \\ \mathbf R & \sim & \text{LKJcorr} (2) \end{eqnarray} \]

Of the notable new parts, \(\mathbf S\) is the covariance matrix and \(\mathbf R\) is the corresponding correlation matrix, which we might more fully express as

\[\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\]

And according to our prior, \(\mathbf R\) is distributed as \(\text{LKJcorr} (2)\). We’ll use rethinking::rlkjcorr() to get a better sense of what that even is.

library(rethinking)

n_sim <- 1e5

set.seed(133)
r_1 <- 
  rlkjcorr(n_sim, K = 2, eta = 1) %>%
  as_tibble()

set.seed(133)
r_2 <- 
  rlkjcorr(n_sim, K = 2, eta = 2) %>%
  as_tibble()

set.seed(133)
r_4 <- 
  rlkjcorr(n_sim, K = 2, eta = 4) %>%
  as_tibble()

Here are the \(\text{LKJcorr}\) distributions of Figure 13.3.

ggplot(data = r_1, aes(x = V2)) +
  geom_density(color = "transparent", fill = "#DCA258", alpha = 2/3) +
  geom_density(data = r_2,
               color = "transparent", fill = "#FCF9F0", alpha = 2/3) +
  geom_density(data = r_4,
               color = "transparent", fill = "#394165", alpha = 2/3) +
  geom_text(data = tibble(x     = c(.83, .62, .46),
                          y     = c(.54, .74, 1),
                          label = c("eta = 1", "eta = 2", "eta = 4")),
            aes(x = x, y = y, label = label),
            color = "#A65141", family = "Courier") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("correlation") +
  theme_pearl_earring

Okay, let’s get ready to model and switch out rethinking for brms.

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

As defined above, our first model has both varying intercepts and afternoon slopes. I should point out that the (1 + afternoon | cafe) syntax specifies that we’d like brm() to fit the random effects for 1 (i.e., the intercept) and the afternoon slope as correlated. Had we wanted to fit a model in which they were orthogonal, we’d have coded (1 + afternoon || cafe).

 b13.1 <- 
  brm(data = d, family = gaussian,
      wait ~ 1 + afternoon + (1 + afternoon | cafe),
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(cauchy(0, 2), class = sd),
                prior(cauchy(0, 2), class = sigma),
                prior(lkj(2), class = cor)),
      iter = 5000, warmup = 2000, chains = 2, cores = 2)

With Figure 13.4, we assess how the posterior for the correlation of the random effects compares to its prior.

post <- posterior_samples(b13.1)

post %>%
  ggplot(aes(x = cor_cafe__Intercept__afternoon)) +
  geom_density(data = r_2, aes(x = V2),
               color = "transparent", fill = "#EEDA9D", alpha = 3/4) +
  geom_density(color = "transparent", fill = "#A65141", alpha = 9/10) +
  annotate("text", label = "posterior", 
           x = -0.2, y = 2.2, 
           color = "#A65141", family = "Courier") +
  annotate("text", label = "prior", 
           x = 0, y = 0.85, 
           color = "#EEDA9D", alpha = 2/3, family = "Courier") +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab("correlation") +
  theme_pearl_earring

McElreath then depicted multidimensional shrinkage by plotting the posterior mean of the varying effects compared to their raw, unpooled estimated. With brms, we can get the cafe-specific intercepts and afternoon slopes with coef(), which returns a three-dimensional list.

# coef(b13.1) %>% glimpse()

coef(b13.1)
## $cafe
## , , Intercept
## 
##    Estimate Est.Error     Q2.5    Q97.5
## 1  4.076657 0.2074001 3.673043 4.485835
## 2  1.930431 0.2083095 1.524455 2.336841
## 3  4.817168 0.2112695 4.400361 5.229888
## 4  3.475913 0.2072789 3.076951 3.879560
## 5  1.775441 0.2182903 1.346382 2.204663
## 6  4.384955 0.2062509 3.980472 4.784476
## 7  3.246504 0.2001336 2.847035 3.646729
## 8  4.008772 0.2107133 3.595458 4.428586
## 9  4.236734 0.2052556 3.836461 4.636787
## 10 3.711203 0.2078821 3.303882 4.109533
## 11 2.168381 0.2143573 1.743727 2.582371
## 12 4.035374 0.2089868 3.628665 4.438421
## 13 4.086637 0.2133844 3.670318 4.497274
## 14 3.591279 0.2054213 3.196173 3.999429
## 15 4.260217 0.2050365 3.857916 4.659703
## 16 3.482390 0.2064012 3.075759 3.880146
## 17 4.040741 0.2113138 3.615844 4.457417
## 18 5.840836 0.2091523 5.438021 6.239810
## 19 3.774387 0.2069485 3.362714 4.179448
## 20 3.856362 0.2105514 3.450898 4.264787
## 
## , , afternoon
## 
##      Estimate Est.Error      Q2.5       Q97.5
## 1  -1.4028304 0.2650106 -1.925897 -0.89745586
## 2  -0.9499469 0.2730138 -1.486945 -0.42107261
## 3  -1.8877371 0.2747891 -2.421923 -1.35021637
## 4  -1.1999865 0.2609638 -1.706151 -0.68144455
## 5  -0.5798662 0.2762329 -1.128942 -0.04157682
## 6  -1.4941387 0.2634284 -2.017196 -0.97103914
## 7  -1.0510192 0.2543420 -1.547794 -0.53539776
## 8  -1.7245287 0.2750980 -2.278482 -1.20547927
## 9  -1.5635515 0.2603850 -2.078395 -1.05338222
## 10 -1.0259160 0.2656137 -1.523594 -0.48770944
## 11 -0.4837134 0.2840094 -1.024349  0.07764580
## 12 -1.2707149 0.2654349 -1.791585 -0.74327996
## 13 -1.8111164 0.2771616 -2.367564 -1.27067278
## 14 -1.6278600 0.2713188 -2.158032 -1.12148827
## 15 -1.6767712 0.2668051 -2.207544 -1.16339863
## 16 -0.9631167 0.2558213 -1.456668 -0.45254115
## 17 -0.6856954 0.2855379 -1.232211 -0.11076829
## 18 -1.5524621 0.2716130 -2.084856 -1.01864804
## 19 -0.9174169 0.2692748 -1.438692 -0.39226069
## 20 -0.9492934 0.2706767 -1.465136 -0.41188546

Here’s the code to extract the relevant elements from the coef() list, convert them to a tibble, and add the cafe index.

partially_pooled_params <-
  # With this line we select each of the 20 cafe's posterior mean (i.e., Estimate)
  # for both `Intercept` and `afternoon`
  coef(b13.1)$cafe[ , 1, 1:2] %>%
  as_tibble() %>%               # convert the two vectors to a tibble
  rename(Slope = afternoon) %>%
  mutate(cafe = 1:nrow(.)) %>%  # add the `cafe` index
  select(cafe, everything())    # simply moving `cafe` to the leftmost position

Like McElreath, we’ll compute the unpooled estimates directly from the data.

# compute unpooled estimates directly from data
un_pooled_params <-
  d %>%
  # With these two lines, we compute the mean value for each cafe's wait time 
  # in the morning and then the afternoon.
  group_by(afternoon, cafe) %>%
  summarise(mean = mean(wait)) %>%
  ungroup() %>%  # Ungrouping allows us to alter afternoon, one of the grouping variables
  mutate(afternoon = ifelse(afternoon == 0, "Intercept", "Slope")) %>%
  spread(key = afternoon, value = mean) %>%  # use `spread()` just as in the previous block
  mutate(Slope = Slope - Intercept)          # Finally, here's our slope!

# Here we combine the partially-pooled and unpooled means into a single data object, 
# which will make plotting easier.
params <-
  # `bind_rows()` will stack the second tibble below the first
  bind_rows(partially_pooled_params, un_pooled_params) %>%
  # index whether the estimates are pooled
  mutate(pooled = rep(c("partially", "not"), each = nrow(.)/2)) 

# Here's a glimpse at what we've been working for
params %>%
  slice(c(1:5, 36:40))
## # A tibble: 10 x 4
##     cafe Intercept  Slope pooled   
##    <int>     <dbl>  <dbl> <chr>    
##  1     1      4.08 -1.40  partially
##  2     2      1.93 -0.950 partially
##  3     3      4.82 -1.89  partially
##  4     4      3.48 -1.20  partially
##  5     5      1.78 -0.580 partially
##  6    16      3.42 -0.836 not      
##  7    17      3.91 -0.348 not      
##  8    18      5.89 -1.50  not      
##  9    19      3.70 -0.733 not      
## 10    20      3.79 -0.774 not

Finally, here’s our code for Figure 13.5.a, showing shrinkage in two dimensions.

ggplot(data = params, aes(x = Intercept, y = Slope)) +
  stat_ellipse(geom = "polygon", type = "norm", level = 1/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 2/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 3/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 4/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 5/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 6/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 7/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 8/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = 9/10, size = 0, alpha = 1/20, fill = "#E7CDC2") +
  stat_ellipse(geom = "polygon", type = "norm", level = .99,  size = 0, alpha = 1/20, fill = "#E7CDC2") +
  geom_point(aes(group = cafe, color = pooled)) +
  geom_line(aes(group = cafe), size = 1/4) +
  scale_color_manual("Pooled?",
                     values = c("#80A0C7", "#A65141")) +
  coord_cartesian(xlim = range(params$Intercept),
                  ylim = range(params$Slope)) +
  theme_pearl_earring

Learn more about stat_ellipse(), here. Let’s prep for Figure 13.5.b.

# retrieve the partially-pooled estimates with `coef()`
partially_pooled_estimates <-
  coef(b13.1)$cafe[ , 1, 1:2] %>%
  as_tibble() %>%                  # convert the two vectors to a tibble
  rename(morning = Intercept) %>%  # the Intercept is the wait time for morning (i.e., `afternoon == 0`)
  mutate(afternoon = morning + afternoon,  # `afternoon` wait time is the `morning` wait time plus the afternoon slope
         cafe      = 1:n()) %>%  # Add the `cafe` index
  select(cafe, everything()) 

# Compute unpooled estimates directly from data
un_pooled_estimates <-
  d %>%
  # As above, with these two lines, we compute each cafe's mean wait value by time of day.
  group_by(afternoon, cafe) %>% 
  summarise(mean = mean(wait)) %>%
  ungroup() %>%  # ungrouping allows us to alter the grouping variable, afternoon
  mutate(afternoon = ifelse(afternoon == 0, "morning", "afternoon")) %>%
  spread(key = afternoon, value = mean)  # this seperates out the values into morning and afternoon columns

estimates <-
  bind_rows(partially_pooled_estimates, un_pooled_estimates) %>%
  mutate(pooled = rep(c("partially", "not"), each = n() / 2))

The code for Figure 13.5.b.

ggplot(data = estimates, aes(x = morning, y = afternoon)) +
  # Nesting `stat_ellipse()` within `mapply()` is a less redundant way to produce the 
  # ten-layered semitransparent ellipses we did with ten lines of `stat_ellipse()` 
  # functions in the previous plot
  mapply(function(level) {
    stat_ellipse(geom  = "polygon", type = "norm",
                 size  = 0, alpha = 1/20, fill = "#E7CDC2",
                 level = level)
    }, 
    # Enter the levels here
    level = c(seq(from = 1/10, to = 9/10, by = 1/10), .99)) +
  geom_point(aes(group = cafe, color = pooled)) +
  geom_line(aes(group = cafe), size = 1/4) +
  scale_color_manual("Pooled?",
                     values = c("#80A0C7", "#A65141")) +
  coord_cartesian(xlim = range(estimates$morning),
                  ylim = range(estimates$afternoon)) +
  labs(x = "morning wait (mins)",
       y = "afternoon wait (mins)") +
  theme_pearl_earring

13.2 Example: Admission decisions and gender

Let’s revisit the infamous UCB admissions data.

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

Here we detach rethinking, reload brms, and augment the data a bit.

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

d <- 
  d %>%
  mutate(male    = ifelse(applicant.gender == "male", 1, 0),
         dept_id = rep(1:6, each = 2))

13.2.1 Varying intercepts.

The statistical formula for our varying-intercepts logistic regression model follows the form

\[ \begin{eqnarray} \text{admit}_i & \sim & \text{Binomial} (n_i, p_i) \\ \text{logit} (p_i) & = & \alpha_{\text{dept_id}_i} + \beta \text{male}_i \\ \alpha_\text{dept_id} & \sim & \text{Normal} (\alpha, \sigma) \\ \alpha & \sim & \text{Normal} (0, 10) \\ \beta & \sim & \text{Normal} (0, 1) \\ \sigma & \sim & \text{HalfCauchy} (0, 2) \\ \end{eqnarray} \]

Since there’s only one left-hand term in our (1 | dept_id) code, there’s only one random effect.

b13.2 <- 
  brm(data = d, family = binomial,
      admit | trials(applications) ~ 1 + male + (1 | dept_id),
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(cauchy(0, 2), class = sd)),
      iter = 4500, warmup = 500, chains = 3, cores = 3,
      control = list(adapt_delta = 0.99))

Since we don’t have a depth=2 argument in brms::summary(), we’ll have to get creative. One way to look at the parameters is with b13.2$fit:

b13.2$fit
## Inference for Stan model: 986d730d4cab1fb3046dbb2c44ac3aae.
## 3 chains, each with iter=4500; warmup=500; thin=1; 
## post-warmup draws per chain=4000, total post-warmup draws=12000.
## 
##                          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## b_Intercept             -0.57    0.01 0.65  -1.85  -0.94  -0.56  -0.19   0.76  1882    1
## b_male                  -0.10    0.00 0.08  -0.25  -0.15  -0.10  -0.04   0.06  5114    1
## sd_dept_id__Intercept    1.49    0.01 0.58   0.80   1.10   1.37   1.73   3.00  2251    1
## r_dept_id[1,Intercept]   1.24    0.01 0.65  -0.09   0.86   1.23   1.62   2.54  1887    1
## r_dept_id[2,Intercept]   1.20    0.01 0.65  -0.13   0.82   1.19   1.58   2.49  1896    1
## r_dept_id[3,Intercept]  -0.02    0.01 0.65  -1.32  -0.39  -0.02   0.36   1.28  1887    1
## r_dept_id[4,Intercept]  -0.05    0.01 0.65  -1.37  -0.42  -0.06   0.32   1.23  1881    1
## r_dept_id[5,Intercept]  -0.49    0.01 0.65  -1.82  -0.87  -0.50  -0.11   0.78  1900    1
## r_dept_id[6,Intercept]  -2.04    0.01 0.66  -3.38  -2.42  -2.04  -1.65  -0.75  1945    1
## lp__                   -61.77    0.05 2.49 -67.54 -63.22 -61.44 -59.97 -57.89  2375    1
## 
## Samples were drawn using NUTS(diag_e) at Mon Sep 24 19:28:16 2018.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

However, notice that the group-specific parameters don’t match up with those in the text. Though our r_dept_id[1,Intercept] had a posterior mean of 1.27, the number for a_dept[1] in the text is 0.67. This is because the brms package presented the random effects in the non-centered metric. The rethinking package, in contrast, presented the random effects in the centered metric. On page 399, McElreath wrote:

Remember, the values above are the \(\alpha_{DEPT}\) estimates, and so they are deviations from the global mean \(\alpha\), which in this case has posterior mean -0.58. So department A, “[1]” in the table, has the highest average admission rate. Department F, “[6]” in the table, has the lowest.

Here’s another fun fact:

# Numbers taken from the mean column on page 399 in the text
c(0.67, 0.63, -0.59, -0.62, -1.06, -2.61) %>% mean()
## [1] -0.5966667

The average of the rethinking-based centered random effects is within rounding error of the global mean, -0.58. If you want the random effects in the centered metric from brms, you can use the coef() function:

coef(b13.2)
## $dept_id
## , , Intercept
## 
##     Estimate  Est.Error       Q2.5      Q97.5
## 1  0.6748743 0.09869596  0.4799627  0.8678702
## 2  0.6306516 0.11550916  0.4044213  0.8592874
## 3 -0.5832062 0.07462316 -0.7304812 -0.4367849
## 4 -0.6157081 0.08498400 -0.7829412 -0.4495225
## 5 -1.0597297 0.09859278 -1.2584943 -0.8702506
## 6 -2.6065066 0.15474597 -2.9195257 -2.3113986
## 
## , , male
## 
##      Estimate  Est.Error       Q2.5      Q97.5
## 1 -0.09530202 0.07988608 -0.2528348 0.06242025
## 2 -0.09530202 0.07988608 -0.2528348 0.06242025
## 3 -0.09530202 0.07988608 -0.2528348 0.06242025
## 4 -0.09530202 0.07988608 -0.2528348 0.06242025
## 5 -0.09530202 0.07988608 -0.2528348 0.06242025
## 6 -0.09530202 0.07988608 -0.2528348 0.06242025

And just to confirm, the average of the posterior means of the Intercept random effects with brms::coef() is also the global mean within rounding error:

mean(coef(b13.2)$dept_id[ , "Estimate", "Intercept"])
## [1] -0.5932708

Note how coef() returned a three-dimensional list.

coef(b13.2) %>% str()
## List of 1
##  $ dept_id: num [1:6, 1:4, 1:2] 0.675 0.631 -0.583 -0.616 -1.06 ...
##   ..- attr(*, "dimnames")=List of 3
##   .. ..$ : chr [1:6] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:4] "Estimate" "Est.Error" "Q2.5" "Q97.5"
##   .. ..$ : chr [1:2] "Intercept" "male"

If you just want the parameter summaries for the random intercepts, you have to use three-dimensional indexing.

coef(b13.2)$dept_id[ , , "Intercept"]  # this also works: coef(b13.2)$dept_id[ , , 1]
##     Estimate  Est.Error       Q2.5      Q97.5
## 1  0.6748743 0.09869596  0.4799627  0.8678702
## 2  0.6306516 0.11550916  0.4044213  0.8592874
## 3 -0.5832062 0.07462316 -0.7304812 -0.4367849
## 4 -0.6157081 0.08498400 -0.7829412 -0.4495225
## 5 -1.0597297 0.09859278 -1.2584943 -0.8702506
## 6 -2.6065066 0.15474597 -2.9195257 -2.3113986

So to get our brms summaries in a similar format to those in the text, we’ll have to combine coef() with fixef() and VarCorr().

coef(b13.2)$dept_id[, , "Intercept"] %>%
  as_tibble() %>% 
  bind_rows(fixef(b13.2) %>% 
              as_tibble()) %>% 
  bind_rows(VarCorr(b13.2)$dept_id$sd %>% 
              as_tibble())
## # A tibble: 9 x 4
##   Estimate Est.Error   Q2.5   Q97.5
##      <dbl>     <dbl>  <dbl>   <dbl>
## 1   0.675     0.0987  0.480  0.868 
## 2   0.631     0.116   0.404  0.859 
## 3  -0.583     0.0746 -0.730 -0.437 
## 4  -0.616     0.0850 -0.783 -0.450 
## 5  -1.06      0.0986 -1.26  -0.870 
## 6  -2.61      0.155  -2.92  -2.31  
## 7  -0.567     0.649  -1.85   0.756 
## 8  -0.0953    0.0799 -0.253  0.0624
## 9   1.49      0.579   0.795  3.00

And a little more data wrangling will make the summaries easier to read:

coef(b13.2)$dept_id[, , "Intercept"] %>%
  as_tibble() %>% 
  bind_rows(fixef(b13.2) %>% 
              as_tibble()) %>% 
  bind_rows(VarCorr(b13.2)$dept_id$sd %>% 
              as_tibble()) %>% 
  mutate(parameter = c(paste("Intercept [", 1:6, "]", sep = ""), 
                       "Intercept", "male", "sigma")) %>% 
  select(parameter, everything()) %>% 
  mutate_if(is_double, round, digits = 2)
## # A tibble: 9 x 5
##   parameter     Estimate Est.Error  Q2.5 Q97.5
##   <chr>            <dbl>     <dbl> <dbl> <dbl>
## 1 Intercept [1]    0.67      0.1    0.48  0.87
## 2 Intercept [2]    0.63      0.12   0.4   0.86
## 3 Intercept [3]   -0.580     0.07  -0.73 -0.44
## 4 Intercept [4]   -0.62      0.08  -0.78 -0.45
## 5 Intercept [5]   -1.06      0.1   -1.26 -0.87
## 6 Intercept [6]   -2.61      0.15  -2.92 -2.31
## 7 Intercept       -0.570     0.65  -1.85  0.76
## 8 male            -0.1       0.08  -0.25  0.06
## 9 sigma            1.49      0.580  0.8   3

I’m not aware of a slick and easy way to get the n_eff and Rhat summaries into the mix. But if you’re fine with working with the brms-default non-centered parameterization, b13.2$fit gets you those just fine.

One last thing. The broom package offers a very handy way to get those brms random effects. Just throw the model brm() fit into the tidy() function.

library(broom)

tidy(b13.2) %>%
  mutate_if(is.numeric, round, digits = 2)  # This line just rounds the output
##                      term estimate std.error  lower  upper
## 1             b_Intercept    -0.57      0.65  -1.59   0.46
## 2                  b_male    -0.10      0.08  -0.23   0.04
## 3   sd_dept_id__Intercept     1.49      0.58   0.86   2.57
## 4  r_dept_id[1,Intercept]     1.24      0.65   0.22   2.27
## 5  r_dept_id[2,Intercept]     1.20      0.65   0.17   2.23
## 6  r_dept_id[3,Intercept]    -0.02      0.65  -1.04   1.02
## 7  r_dept_id[4,Intercept]    -0.05      0.65  -1.09   0.97
## 8  r_dept_id[5,Intercept]    -0.49      0.65  -1.53   0.54
## 9  r_dept_id[6,Intercept]    -2.04      0.66  -3.09  -1.02
## 10                   lp__   -61.77      2.49 -66.32 -58.36

But note how, just as with b13.2$fit, this approach summarizes the posterior with the non-centered parameterization. Which is a fine parameterization. It’s just a little different from what you’ll get when using precis( m13.2 , depth=2 ), as in the text.

13.2.2 Varying effects of being male.

Now we’re ready to allow our male dummy to varies, too, the statistical model follows the form

\[ \begin{eqnarray} \text{admit}_i & \sim & \text{Binomial} (n_i, p_i) \\ \text{logit} (p_i) & = & \alpha_{\text{dept_id}_i} + \beta_{\text{dept_id}_i} \text{male}_i \\ \begin{bmatrix} \alpha_\text{dept_id} \\ \beta_\text{dept_id} \end{bmatrix} & \sim & \text{MVNormal} \bigg (\begin{bmatrix} \alpha \\ \beta \end{bmatrix}, \mathbf{S} \bigg ) \\ \mathbf S & = & \begin{pmatrix} \sigma_\alpha & 0 \\ 0 & \sigma_\beta \end{pmatrix} \mathbf R \begin{pmatrix} \sigma_\alpha & 0 \\ 0 & \sigma_\beta \end{pmatrix} \\ \alpha & \sim & \text{Normal} (0, 10) \\ \beta & \sim & \text{Normal} (0, 1) \\ (\sigma_\alpha, \sigma_\beta) & \sim & \text{HalfCauchy} (0, 2) \\ \mathbf R & \sim & \text{LKJcorr} (2) \end{eqnarray} \]

Fit the model.

b13.3 <- 
  brm(data = d, family = binomial,
      admit | trials(applications) ~ 1 + male + (1 + male | dept_id),
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(cauchy(0, 2), class = sd),
                prior(lkj(2), class = cor)),
      iter = 5000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = .99,
                     max_treedepth = 12))

McElreath encouraged us to make sure the chains look good. Instead of relying on convenience functions, let’s do it by hand.

post <- posterior_samples(b13.3, add_chain = T)

post %>% 
  select(-lp__) %>% 
  gather(key, value, -chain, -iter) %>% 
  mutate(chain = as.character(chain)) %>% 

  ggplot(aes(x = iter, y = value, group = chain, color = chain)) +
  geom_line(size = 1/15) +
  scale_color_manual(values = c("#80A0C7", "#B1934A", "#A65141", "#EEDA9D")) +
  labs(x = NULL, y = NULL) +
  scale_x_continuous(breaks = c(1001, 5000)) +
  theme_pearl_earring +
  theme(legend.position  = c(.825, .06),
        legend.direction = "horizontal") +
  facet_wrap(~key, ncol = 3, scales = "free_y")

Our chains look great. While we’re at it, let’s examine the \(\hat{R}\) vales in a handmade plot, too.

rhat(b13.3) %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  filter(rowname != "lp__") %>% 
  
  ggplot(aes(x = `.`, y = reorder(rowname, `.`))) + 
  geom_segment(aes(xend = 1, yend = rowname),
               color = "#EEDA9D") +
  geom_point(aes(color = `.` > 1), 
             size = 2) +
  scale_color_manual(values = c("#80A0C7", "#A65141")) +
  labs(x = NULL, y = NULL) +
  theme_pearl_earring +
  theme(legend.position = "none",
        axis.ticks.y    = element_blank(),
        axis.text.y     = element_text(hjust = 0))

Them are some respectable \(\hat{R}\) values. The plot accentuates their differences, but they’re all basically 1 (e.g., see what happens is you set coord_cartesian(xlim = c(0.99, 1.01))). Here are the random effects in the centered metric:

coef(b13.3)
## $dept_id
## , , Intercept
## 
##     Estimate  Est.Error       Q2.5      Q97.5
## 1  1.3060287 0.25584336  0.8242594  1.8202385
## 2  0.7420184 0.32818266  0.1035647  1.3976491
## 3 -0.6470283 0.08457924 -0.8135330 -0.4827617
## 4 -0.6172350 0.10497566 -0.8239028 -0.4132265
## 5 -1.1315001 0.11332819 -1.3580756 -0.9143988
## 6 -2.6007380 0.20242186 -3.0135856 -2.2216006
## 
## , , male
## 
##      Estimate Est.Error       Q2.5      Q97.5
## 1 -0.79382902 0.2692380 -1.3324602 -0.2744163
## 2 -0.21194501 0.3305760 -0.8767613  0.4424014
## 3  0.08194407 0.1391819 -0.1814187  0.3592572
## 4 -0.09240970 0.1394929 -0.3666799  0.1801527
## 5  0.12153994 0.1866106 -0.2342076  0.4933623
## 6 -0.12447792 0.2722629 -0.6711310  0.3969950

We may as well keep our doing-things-by-hand kick going. Instead relying on bayesplog::mcmc_intervals() or tidybayes::pointintervalh() to make our coefficient plot, we’ll combine geom_pointrange() and coord_flip(). But we will need to wrangle a bit to get those brms-based centered random effects into a usefully-formatted tidy tibble.

# As far as I can tell, because `coef()` yields a list, you have to take out the two 
# random effects one at a time, convert them to tibbles, and reassemble them with `bind_rows()`
coef(b13.3)$dept_id[, , 1] %>% 
  as_tibble() %>% 
  bind_rows(
    coef(b13.3)$dept_id[, , 2] %>% 
      as_tibble()
    ) %>% 
  mutate(param   = c(paste("Intercept", 1:6), paste("male", 1:6)),
         reorder = c(6:1, 12:7)) %>% 

  # plot
  ggplot(aes(x = reorder(param, reorder))) +
  geom_hline(yintercept = 0, linetype = 3, color = "#8B9DAF") +
  geom_pointrange(aes(ymin = Q2.5, ymax = Q97.5, y = Estimate, color = reorder < 7),
                  shape = 20, size = 3/4) +
  scale_color_manual(values = c("#394165", "#A65141")) +
  xlab(NULL) +
  coord_flip() +
  theme_pearl_earring +
  theme(legend.position = "none",
        axis.ticks.y    = element_blank(),
        axis.text.y     = element_text(hjust = 0))

Just like in the text, our male slopes are much less dispersed than our intercepts.

13.2.3 Shrinkage.

Figure 13.6.a depicts the correlation between the full UCB model’s varying intercepts and slopes.

library(tidybayes)

post <- posterior_samples(b13.3)

post %>% 
  ggplot(aes(x = cor_dept_id__Intercept__male, y = 0)) +
  geom_halfeyeh(fill = "#394165", color = "#8B9DAF", 
                point_interval = median_qi, .width = .95) +
  scale_x_continuous(breaks = c(-1, median(post$cor_dept_id__Intercept__male), 1),
                     labels = c(-1, "-.35", 1)) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = -1:1) +
  labs(subtitle = "The dot is at the median; the\nhorizontal bar is the 95% CI.",
       x = "correlation") +
  theme_pearl_earring

Much like for Figure 13.5.b, above, it’ll take a little data processing before we’re ready to reproduce Figure 13.6.b.

# Here we put the partially-pooled estimate summaries in a tibble
partially_pooled_params <-
  coef(b13.3)$dept_id[ , 1, ] %>%
  as_tibble() %>%
  rename(intercept = Intercept,
         slope     = male) %>%
  mutate(dept = 1:n()) %>%
  select(dept, everything())

# In order to calculate the unpooled estimates from the data, we'll need a function that 
# can convert probabilities into the logit metric. If you do the algebra, this is just
# a transformation of the `inv_logit_scaled()` function.
prob_to_logit <- function(x){
  -log((1 / x) -1)
  }

# compute unpooled estimates directly from data
un_pooled_params <-
  d %>%
  group_by(male, dept_id) %>%
  summarise(prob_admit = mean(admit / applications)) %>%
  ungroup() %>%
  mutate(male = ifelse(male == 0, "intercept", "slope")) %>%
  spread(key  = male, value = prob_admit) %>%
  rename(dept = dept_id) %>%
  mutate(intercept = prob_to_logit(intercept),  # Here we put our `prob_to_logit()` function to work
         slope     = prob_to_logit(slope)) %>%
  mutate(slope     = slope - intercept)

# Here we combine the partially-pooled and unpooled means into a single data object.
params <-
  bind_rows(partially_pooled_params, un_pooled_params) %>%
  mutate(pooled      = rep(c("partially", "not"), each = n() / 2)) %>%
  mutate(dept_letter = rep(LETTERS[1:6], times = 2))  # This will help with plotting

params
## # A tibble: 12 x 5
##     dept intercept   slope pooled    dept_letter
##    <int>     <dbl>   <dbl> <chr>     <chr>      
##  1     1     1.31  -0.794  partially A          
##  2     2     0.742 -0.212  partially B          
##  3     3    -0.647  0.0819 partially C          
##  4     4    -0.617 -0.0924 partially D          
##  5     5    -1.13   0.122  partially E          
##  6     6    -2.60  -0.124  partially F          
##  7     1     1.54  -1.05   not       A          
##  8     2     0.754 -0.220  not       B          
##  9     3    -0.660  0.125  not       C          
## 10     4    -0.622 -0.0820 not       D          
## 11     5    -1.16   0.200  not       E          
## 12     6    -2.58  -0.189  not       F

Here’s our version of Figure 13.6.b, depicting two-dimensional shrinkage for the partially-pooled multilevel estimates (posterior means) relative to the unpooled coefficients, calculated from the data. The ggrepel::geom_text_repel() function will help us with the in-plot labels.

library(ggrepel)

ggplot(data = params, aes(x = intercept, y = slope)) +
  mapply(function(level){
    stat_ellipse(geom  = "polygon", type = "norm",
                 size  = 0, alpha = 1/20, fill = "#E7CDC2",
                 level = level)
    },  
    level = c(seq(from = 1/10, to = 9/10, by = 1/10), .99)) +
  geom_point(aes(group = dept, color = pooled)) +
  geom_line(aes(group = dept), size = 1/4) +
  scale_color_manual("Pooled?",
                     values = c("#80A0C7", "#A65141")) +
  geom_text_repel(data = params %>% filter(pooled == "partially"),
                  aes(label = dept_letter),
                  color = "#E8DCCF", size = 4, family = "Courier", seed = 13.6) +
  coord_cartesian(xlim = range(params$intercept),
                  ylim = range(params$slope)) +
  labs(x = expression(paste("intercept (", alpha[dept_id], ")")),
       y = expression(paste("slope (", beta[dept_id], ")"))) +
  theme_pearl_earring

13.2.4 Model comparison.

Fit the no-gender model.

b13.4 <- 
  brm(data = d, family = binomial,
      admit | trials(applications) ~ 1 + (1 | dept_id),
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(cauchy(0, 2), class = sd)),
      iter = 5000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = .99,
                     max_treedepth = 12))

Compare the three models by the WAIC.

waic(b13.2, b13.3, b13.4)
##                 WAIC    SE
## b13.2         108.43 16.46
## b13.3          90.88  4.61
## b13.4         105.30 18.00
## b13.2 - b13.3  17.55 13.45
## b13.2 - b13.4   3.12  3.62
## b13.3 - b13.4 -14.42 15.20

The varying slopes model, [b13.3], dominates [the other two]. This is despite the fact that the average slope in [b13.3] is nearly zero. The average isn’t what matters, however. It is the individual slopes, one for each department, that matter. If we wish to generalize to new departments, the variation in slopes suggest that it’ll be worth paying attention to gender, even if the average slope is nearly zero in the population. (pp. 402–403, emphasis in the original)

13.3 Example: Cross-classified chimpanzees with varying slopes

Retrieve the chimpanzees data.

library(rethinking)
data(chimpanzees)
d <- chimpanzees
detach(package:rethinking, unload = T)
library(brms)
rm(chimpanzees)

d <-
  d %>%
  select(-recipient) %>%
  mutate(block_id = block)

My math’s aren’t the best. But if I’m following along correctly, here’s a fuller statistical expression of our cross-classified model.

\[\begin{eqnarray} \text{pulled_left}_i & \sim & \text{Binomial} (n = 1, p_i) \\ \text{logit} (p_i) & = & \alpha_i + (\beta_{1i} + \beta_{2i} \text{condition}_i) \text{prosoc_left}_i \\ \alpha_i & = & \alpha + \alpha_{\text{actor}_i} + \alpha_{\text{block_id}_i} \\ \beta_{1i} & = & \beta_1 + \beta_{1, \text{actor}_i} + \beta_{1, \text{block_id}_i} \\ \beta_{2i} & = & \beta_2 + \beta_{2, \text{actor}_i} + \beta_{2, \text{block_id}_i} \\ \begin{bmatrix} \alpha_\text{actor} \\ \beta_{1, \text{actor}} \\ \beta_{2, \text{actor}} \end{bmatrix} & \sim & \text{MVNormal} \begin{pmatrix} \begin{bmatrix}0 \\ 0 \\ 0 \end{bmatrix} , \mathbf{S}_\text{actor} \end{pmatrix} \\ \begin{bmatrix} \alpha_\text{block_id} \\ \beta_{1, \text{block_id}} \\ \beta_{2, \text{block_id}} \end{bmatrix} & \sim & \text{MVNormal} \begin{pmatrix} \begin{bmatrix}0 \\ 0 \\ 0 \end{bmatrix} , \mathbf{S}_\text{block_id} \end{pmatrix} \\ \mathbf S_\text{actor} & = & \begin{pmatrix} \sigma_{\alpha_\text{actor}} & 0 & 0 \\ 0 & \sigma_{\beta_{1_\text{actor}}} & 0 \\ 0 & 0 & \sigma_{\beta_{2_\text{actor}}} \end{pmatrix} \mathbf R_\text{actor} \begin{pmatrix} \sigma_{\alpha_\text{actor}} & 0 & 0 \\ 0 & \sigma_{\beta_{1_\text{actor}}} & 0 \\ 0 & 0 & \sigma_{\beta_{2_\text{actor}}} \end{pmatrix} \\ \mathbf S_\text{block_id} & = & \begin{pmatrix} \sigma_{\alpha_\text{block_id}} & 0 & 0 \\ 0 & \sigma_{\beta_{1_\text{block_id}}} & 0 \\ 0 & 0 & \sigma_{\beta_{2_\text{block_id}}} \end{pmatrix} \mathbf R_\text{block_id} \begin{pmatrix} \sigma_{\alpha_\text{block_id}} & 0 & 0 \\ 0 & \sigma_{\beta_{1_\text{block_id}}} & 0 \\ 0 & 0 & \sigma_{\beta_{2_\text{block_id}}} \end{pmatrix} \\ (\sigma_{\alpha_\text{actor}}, \sigma_{\beta_{1_\text{actor}}}, \sigma_{\beta_{2_\text{actor}}}) & \sim & \text{HalfCauchy} (0, 2) \\ (\sigma_{\alpha_\text{block_id}}, \sigma_{\beta_{1_\text{block_id}}}, \sigma_{\beta_{2_\text{block_id}}}) & \sim & \text{HalfCauchy} (0, 2) \\ \mathbf R_\text{actor} & \sim & \text{LKJcorr} (4) \\ \mathbf R_\text{block_id} & \sim & \text{LKJcorr} (4) \end{eqnarray}\]

And now each \(\mathbf R\) is a \(3 \times 3\) correlation matrix.

Let’s fit this beast.

b13.6 <- 
  brm(data = d, family = binomial,
      pulled_left ~ 1 + prosoc_left + condition:prosoc_left +
        (1 + prosoc_left + condition:prosoc_left | actor) +
        (1 + prosoc_left + condition:prosoc_left | block_id),
      prior = c(prior(normal(0, 1), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(cauchy(0, 2), class = sd),
                prior(lkj(4), class = cor)),
      iter = 5000, warmup = 1000, chains = 3, cores = 3)

Even though it’s not apparent in the syntax, our model b13.6 was already fit using the non-centered parameterization. Behind the scenes, Bürkner has brms do this automatically. It’s been that way all along.

If you recall from last chapter, we can compute the number of effective samples for our parameters like so.

ratios_cp <- neff_ratio(b13.6)

neff <-
  ratios_cp %>% 
  as_tibble %>% 
  rename(neff_ratio = value) %>% 
  mutate(neff       = neff_ratio * 12000)

head(neff)
## # A tibble: 6 x 2
##   neff_ratio  neff
##        <dbl> <dbl>
## 1      0.174 2088.
## 2      0.518 6214.
## 3      0.627 7519.
## 4      0.320 3836.
## 5      0.465 5577.
## 6      0.550 6595.

Now we’re ready for our variant of Figure 13.7. The handy ggbeeswarm package and its geom_quasirandom() function will give a better sense of the distribution.

library(ggbeeswarm)

neff %>%
  ggplot(aes(x = factor(0), y = neff)) +
  geom_boxplot(fill = "#394165", color = "#8B9DAF") +
  geom_quasirandom(method = "tukeyDense",
                   size = 2/3, color = "#EEDA9D", alpha = 2/3) +
  scale_x_discrete(NULL, breaks = NULL,
                   expand = c(.75, .75)) +
  scale_y_continuous(breaks = c(0, 6000, 12000)) +
  coord_cartesian(ylim = 0:12000) +
  labs(y = "effective samples",
       subtitle = "The non-centered\nparameterization is the\nbrms default. No fancy\ncoding required.") +
  theme_pearl_earring

As in the last chapter, we’ll use the bayesplot::mcmc_neff() function to examine the ratio of n.eff and the fill number of post-warm-up iterations, N. Ideally, that ratio is closer to 1 than not.

library(bayesplot)

color_scheme_set(c("#DCA258", "#EEDA9D", "#394165", "#8B9DAF", "#A65141", "#A65141"))

mcmc_neff(ratios_cp, size = 2) +
  theme_pearl_earring

Here are our standard deviation parameters.

tidy(b13.6) %>%
  filter(str_detect(term , "sd_")) %>%
  mutate_if(is.numeric, round, digits = 2)
##                                 term estimate std.error lower upper
## 1                sd_actor__Intercept     2.37      0.91  1.31  4.08
## 2              sd_actor__prosoc_left     0.45      0.37  0.04  1.12
## 3    sd_actor__prosoc_left:condition     0.52      0.48  0.04  1.42
## 4             sd_block_id__Intercept     0.22      0.20  0.02  0.59
## 5           sd_block_id__prosoc_left     0.57      0.40  0.07  1.31
## 6 sd_block_id__prosoc_left:condition     0.51      0.44  0.04  1.29

McElreath discussed rethinking::link() in the middle of page 407. He showed how his link(m13.6NC) code returned a list of four matrices, of which the p matrix was of primary interest. The brms::fitted() function doesn’t work quite the same way, here.

fitted(b13.6,
       summary = F,
       nsamples = 1000) %>% 
  str()
## Using the maximum response value as the number of trials.
## Warning: Using 'binomial' families without specifying 'trials' on the left-hand side of the model
## formula is deprecated.
##  num [1:1000, 1:504] 0.365 0.303 0.321 0.248 0.27 ...

First off, recall that fitted() returns summary values, by default. If we want individual values, set summary = FALSE. It’s also the fitted() default to use all posterior iterations, which is 12,000 in this case. To match the text, we need to set nsamples = 1000. But those are just details. The main point is that fitted() only returns one matrix, which is the analogue to the p matrix in the text.

Moving forward, before we can follow along with McElreath’s R code 13.27, we need to refit the simpler model from way back in Chapter 12.

b12.5 <- 
  brm(data = d, family = binomial,
      pulled_left ~ 1 + prosoc_left + condition:prosoc_left + 
        (1 | actor) + (1 | block_id),
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(cauchy(0, 1), class = sd)),
      iter = 5000, warmup = 1000, chains = 3, cores = 3)

Now we can compare them with waic().

waic(b13.6, b12.5)
##                 WAIC    SE
## b13.6         534.77 19.89
## b12.5         532.76 19.69
## b13.6 - b12.5   2.00  4.10
model_weights(b13.6, b12.5, weights = "waic")
##     b13.6     b12.5 
## 0.2685155 0.7314845

In this example, no matter which varying effect structure you use, you’ll find that actors vary a lot in their baseline preference for the left-hand lever. Everything else is much less important. But using the most complex model, [b13.6], tells the correct story. Because the varying slopes are adaptively regularized, the model hasn’t overfit much, relative to the simpler model that contains only the important intercept variation. (p. 408)

13.4 Continuous categories and the Gaussian process

There is a way to apply the varying effects approach to continuous categories… The general approach is known as Gaussian process regression. This name is unfortunately wholly uninformative about what it is for and how it works.

We’ll proceed to work through a basic example that demonstrates both what it is for and how it works. The general purpose is to define some dimension along which cases differ. This might be individual differences in age. Or it could be differences in location. Then we measure the distance between each pair of cases. What the model then does is estimate a function for the covariance between pairs of cases at different distances. This covariance function provides one continuous category generalization of the varying effects approach. (p. 410)

13.4.1 Example: Spatial autocorrelation in Oceanic tools.

# load the distance matrix
library(rethinking)
data(islandsDistMatrix)

# display short column names, so fits on screen
d_mat <- islandsDistMatrix
colnames(d_mat) <- c("Ml", "Ti", "SC", "Ya", "Fi", 
                     "Tr", "Ch", "Mn", "To", "Ha")
round(d_mat, 1)
##             Ml  Ti  SC  Ya  Fi  Tr  Ch  Mn  To  Ha
## Malekula   0.0 0.5 0.6 4.4 1.2 2.0 3.2 2.8 1.9 5.7
## Tikopia    0.5 0.0 0.3 4.2 1.2 2.0 2.9 2.7 2.0 5.3
## Santa Cruz 0.6 0.3 0.0 3.9 1.6 1.7 2.6 2.4 2.3 5.4
## Yap        4.4 4.2 3.9 0.0 5.4 2.5 1.6 1.6 6.1 7.2
## Lau Fiji   1.2 1.2 1.6 5.4 0.0 3.2 4.0 3.9 0.8 4.9
## Trobriand  2.0 2.0 1.7 2.5 3.2 0.0 1.8 0.8 3.9 6.7
## Chuuk      3.2 2.9 2.6 1.6 4.0 1.8 0.0 1.2 4.8 5.8
## Manus      2.8 2.7 2.4 1.6 3.9 0.8 1.2 0.0 4.6 6.7
## Tonga      1.9 2.0 2.3 6.1 0.8 3.9 4.8 4.6 0.0 5.0
## Hawaii     5.7 5.3 5.4 7.2 4.9 6.7 5.8 6.7 5.0 0.0

If you wanted to use color to more effectively visualize the values in the matirx, you might do something like this.

d_mat %>%
  as_tibble() %>%
  gather() %>%
  rename(column = key,
         distance = value) %>%
  mutate(row          = rep(rownames(d_mat), times = 10),
         row_order    = rep(9:0,            times = 10),
         column_order = rep(0:9,            each  = 10)) %>%
  
  ggplot(aes(x = reorder(column, column_order), 
             y = reorder(row,    row_order))) + 
  geom_raster(aes(fill = distance)) + 
  geom_text(aes(label = round(distance, digits = 1)),
            size = 3, family = "Courier", color = "#100F14") +
  scale_fill_gradient(low = "#FCF9F0", high = "#A65141") +
  scale_x_discrete(position = "top", expand = c(0, 0)) +
  scale_y_discrete(expand = c(0, 0)) +
  labs(x = NULL, y = NULL) +
  theme_pearl_earring +
  theme(axis.ticks  = element_blank(),
        axis.text.y = element_text(hjust = 0))

Figure 13.8 shows the “shape of the function relating distance to the covariance \(\mathbf{K}_{ij}\).”

tibble(
  x       = seq(from = 0, to = 4, by = .01),
  linear  = exp(-1 * x),
  squared = exp(-1 * x^2)) %>%
  
  ggplot(aes(x = x)) +
  geom_line(aes(y = linear),
            color = "#B1934A", linetype = 2) +
  geom_line(aes(y = squared),
            color = "#DCA258") +
  scale_x_continuous("distance", expand = c(0, 0)) +
  scale_y_continuous("correlation", 
                     breaks = c(0, .5, 1),
                     labels = c(0, ".5", 1)) +
  theme_pearl_earring

data(Kline2) # load the ordinary data, now with coordinates

d <- 
  Kline2 %>%
  mutate(society = 1:10)

rm(Kline2)

d %>% glimpse()
## Observations: 10
## Variables: 10
## $ culture     <fct> Malekula, Tikopia, Santa Cruz, Yap, Lau Fiji, Trobriand, Chuuk, Manus, Tong...
## $ population  <int> 1100, 1500, 3600, 4791, 7400, 8000, 9200, 13000, 17500, 275000
## $ contact     <fct> low, low, low, high, high, high, high, low, high, low
## $ total_tools <int> 13, 22, 24, 43, 33, 19, 40, 28, 55, 71
## $ mean_TU     <dbl> 3.2, 4.7, 4.0, 5.0, 5.0, 4.0, 3.8, 6.6, 5.4, 6.6
## $ lat         <dbl> -16.3, -12.3, -10.7, 9.5, -17.7, -8.7, 7.4, -2.1, -21.2, 19.9
## $ lon         <dbl> 167.5, 168.8, 166.0, 138.1, 178.1, 150.9, 151.6, 146.9, -175.2, -155.6
## $ lon2        <dbl> -12.5, -11.2, -14.0, -41.9, -1.9, -29.1, -28.4, -33.1, 4.8, 24.4
## $ logpop      <dbl> 7.003065, 7.313220, 8.188689, 8.474494, 8.909235, 8.987197, 9.126959, 9.472...
## $ society     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

Switch out rethinking for brms.

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

Okay, it appears this is going to be a bit of a ride. It’s not entirely clear to me if we can fit a Gaussian process model in brms that’s a direct equivalent to what McElreath did with rethinking. But we can try. First, note our use of the gp() syntax in the brm() function, below. We’re attempting to tell brms that we would like to include latitude and longitude (i.e., lat and long2, respectively) in a Gaussian process. Also note how our priors are a little different than those in the text. I’ll explain, below. Let’s just move ahead and fit the model.

b13.7 <- 
  brm(data = d, family = poisson,
      total_tools ~ 1 + gp(lat, lon2) + logpop,
      prior = c(prior(normal(0, 10), class = Intercept),
                prior(normal(0, 1), class = b),
                prior(inv_gamma(2.874624, 0.393695), class = lscale),
                prior(cauchy(0, 1), class = sdgp)),
      iter = 1e4, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.999,
                     max_treedepth = 12))

Here’s the model summary.

posterior_summary(b13.7) %>%
  round(digits = 2)
##                   Estimate Est.Error   Q2.5  Q97.5
## b_Intercept           1.45      1.11  -0.72   3.78
## b_logpop              0.23      0.10   0.02   0.44
## sdgp_gplatlon2        0.52      0.36   0.16   1.46
## lscale_gplatlon2      0.23      0.13   0.07   0.58
## zgp_gplatlon2[1]     -0.60      0.79  -2.16   0.92
## zgp_gplatlon2[2]      0.44      0.85  -1.26   2.08
## zgp_gplatlon2[3]     -0.63      0.70  -1.98   0.88
## zgp_gplatlon2[4]      0.88      0.71  -0.48   2.32
## zgp_gplatlon2[5]      0.25      0.77  -1.25   1.78
## zgp_gplatlon2[6]     -1.00      0.79  -2.56   0.59
## zgp_gplatlon2[7]      0.13      0.72  -1.39   1.52
## zgp_gplatlon2[8]     -0.18      0.88  -1.91   1.61
## zgp_gplatlon2[9]      0.40      0.92  -1.51   2.11
## zgp_gplatlon2[10]    -0.32      0.83  -1.95   1.33
## lp__                -51.61      3.16 -58.80 -46.46

Our Gaussian process parameters are different than McElreath’s. From the brms reference manual, here’s the brms parameterization:

\[k(x_{i},x_{j}) = sdgp^2 \text{exp} (-||x_{i} - x_{j}||/2lscale^2)\]

What McElreath called \(\eta\), Bürkner called \(sdgp\). While McElreath estimated \(\eta^2\), brms simply estimated \(sdgp\). So we’ll have to square our sdgp_gplatlon2 before it’s on the same scale as etasq in the text. Here it is.

posterior_samples(b13.7) %>% 
  transmute(sdgp_squared = sdgp_gplatlon2^2) %>% 
  mean_hdi(sdgp_squared, .width = .89) %>% 
  mutate_if(is.double, round, digits = 3)
##   sdgp_squared .lower .upper .width .point .interval
## 1        0.399  0.002  0.773   0.89   mean       hdi

Now we’re in the ballpark. In our model brm() code, above, we just went with the flow and kept the cauchy(0, 1) prior on sdgp.

Now look at the denominator of the inner part of Bürkner equation, \(2lscale^2\). This appears to be the brms equivalent to McElreath’s \(\rho^2\). Or at least it’s what we’ve got. Anyway, also note that McElreath estimated \(\rho^2\) directly as rhosq. If I’m doing the algebra correctly–and that may well be a big if–, we might expect:

\[\rho^2 = 1/(2 \times (lscale^2))\]

But that doesn’t appear to be the case. Sigh.

posterior_samples(b13.7) %>% 
  transmute(rho_squared = 1/(2*(lscale_gplatlon2^2))) %>% 
  mean_hdi(rho_squared, .width = .89) %>% 
  mutate_if(is.double, round, digits = 3)
##   rho_squared .lower .upper .width .point .interval
## 1      21.995  0.513 47.817   0.89   mean       hdi

Oh man, that isn’t even close to the 2.67 McElreath reported in the text. The plot deepens. If you look back, you’ll see we used a very different prior for \(lscale\). Here is it: inv_gamma(2.874624, 0.393695). Use get_prior() to discover where that came from.

get_prior(data = d, family = poisson,
      total_tools ~ 1 + gp(lat, lon2) + logpop)
##                           prior     class         coef group resp dpar nlpar bound
## 1                                       b                                         
## 2                                       b       logpop                            
## 3           student_t(3, 3, 10) Intercept                                         
## 4                normal(0, 0.5)    lscale                                         
## 5 inv_gamma(2.874624, 0.393695)    lscale gp(lat,lon2)                            
## 6           student_t(3, 0, 10)      sdgp                                         
## 7                                    sdgp gp(lat,lon2)

That is, we used the brms default prior for \(lscale\). In a GitHub exchange, Bürkner pointed out that brms uses special priors for \(lscale\) parameters based on Michael Betancourt [of the Stan team]’s vignette on the topic. Though it isn’t included in this document, I also ran the model with the cauchy(0, 1) prior and the results were quite similar. So the big discrepancy between our model and the one in the text isn’t based on that prior.

Now that we’ve started, we may as well keep going down the comparison train. Let’s reproduce McElreath’s model with rethinking.

Switch out brms for rethinking.

detach(package:brms, unload = T)
library(rethinking)
m13.7 <- map2stan(
    alist(
        total_tools ~ dpois(lambda),
        log(lambda) <- a + g[society] + bp*logpop,
        g[society] ~ GPL2( Dmat , etasq , rhosq , 0.01 ),
        a ~ dnorm(0,10),
        bp ~ dnorm(0,1),
        etasq ~ dcauchy(0,1),
        rhosq ~ dcauchy(0,1)
    ),
    data=list(
        total_tools=d$total_tools,
        logpop=d$logpop,
        society=d$society,
        Dmat=islandsDistMatrix),
    warmup=2000 , iter=1e4 , chains=4)

Alright, now we’ll work directly with the posteriors to make some visual comparisons.

# rethinking-based posterior
post_m13.7 <- rethinking::extract.samples(m13.7)[2:5] %>% as_tibble()

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

# brms-based posterior
post_b13.7 <- posterior_samples(b13.7)

Here’s the model intercept posterior, by package:

post_m13.7[, "a"] %>% 
  bind_rows(post_b13.7[, "b_Intercept"] %>% 
              as_tibble() %>% 
              rename(a = value)) %>% 
  mutate(package = rep(c("rethinking", "brms"), each = nrow(post_m13.7))) %>% 
  
  ggplot(aes(x = a, fill = package)) +
  geom_density(size = 0, alpha = 1/2) +
  scale_fill_manual(NULL, values = c("#80A0C7", "#A65141")) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "Not identical, but pretty close",
       x     = "intercept") +
  theme_pearl_earring

The slope:

post_m13.7[, "bp"] %>% 
  bind_rows(post_b13.7[, "b_logpop"] %>% 
              as_tibble() %>% 
              rename(bp = value)) %>% 
  mutate(package = rep(c("rethinking", "brms"), each = nrow(post_m13.7))) %>% 
  
  ggplot(aes(x = bp, fill = package)) +
  geom_density(size = 0, alpha = 1/2) +
  scale_fill_manual(NULL, values = c("#80A0C7", "#A65141")) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "Again, pretty close",
       x     = "slope") +
  theme_pearl_earring

This one, \(\eta^2\), required a little transformation:

post_m13.7[, "etasq"] %>% 
  bind_rows(post_b13.7[, "sdgp_gplatlon2"] %>% 
              as_tibble() %>%
              mutate(value = value^2) %>% 
              rename(etasq = value)) %>% 
  mutate(package = rep(c("rethinking", "brms"), each = nrow(post_m13.7))) %>% 
  
  ggplot(aes(x = etasq, fill = package)) +
  geom_density(size = 0, alpha = 1/2) +
  scale_fill_manual(NULL, values = c("#80A0C7", "#A65141")) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "Still in the same ballpark",
       x     = expression(eta^2)) +
  coord_cartesian(xlim = 0:3) +
  theme_pearl_earring

\(\rho^2\) required more extensive transformation of the brms posterior:

post_m13.7[, "rhosq"] %>%
  bind_rows(post_b13.7[, "lscale_gplatlon2"] %>% 
              as_tibble() %>%
              transmute(value = 1/(2*(value^2))) %>%
              # transmute(value = value^2) %>% 
              rename(rhosq = value)) %>% 
  mutate(package = rep(c("rethinking", "brms"), each = nrow(post_m13.7))) %>% 
  
  ggplot(aes(x = rhosq, fill = package)) +
  geom_density(size = 0) +
  scale_fill_manual(NULL, values = c("#80A0C7", "#A65141")) +
  labs(title    = "Holy smokes are those not the same!",
       subtitle = "Notice how differently the y axes got scaled. Also, that brms density is\nright skewed for days.",
       x        = expression(rho^2)) +
  coord_cartesian(xlim = 0:50) +
  theme_pearl_earring +
  theme(legend.position = "none") +
  facet_wrap(~package, scales = "free_y")

I’m in clinical psychology. Folks in my field don’t tend to use Gaussian processes, so getting to the bottom of this is low on my to-do list. Perhaps one of y’all are more experienced with Gaussian processes and see a flaw somewhere in my code. Please hit me up if you do.

Anyways, here’s our brms + ggplot2 version of Figure 13.9.

ggplot(data = tibble(x = c(0, 50.2)), aes(x = x)) +
  mapply(function(etasq, rhosq) {
    stat_function(fun = function(x, etasq, rhosq) etasq*exp(-rhosq*x^2), 
                  args = list(etasq = etasq, rhosq = rhosq), 
                  size = 1/4,
                  alpha = 1/4,
                  color = "#EEDA9D")
  }, 
  etasq = post_b13.7[1:100, "sdgp_gplatlon2"]^2,
  rhosq = post_b13.7[1:100, "lscale_gplatlon2"]^2*.5
  ) +
  stat_function(fun = function(x) median(post_b13.7$sdgp_gplatlon2)^2 *exp(-median(post_b13.7[1:100, "lscale_gplatlon2"] )^2*.5*x^2),
                color = "#EEDA9D", size = 1.1) +
  coord_cartesian(ylim = 0:1) +
  scale_x_continuous(breaks = seq(from = 0, to = 50, by = 10),
                     expand = c(0, 0)) +
  labs(x = "distance (thousand km)", 
       y = "covariance") +
  theme_pearl_earring

Do note the scale on which we placed our x axis. Our brms parameterization resulted in a gentler decline in spatial covariance.

Let’s finish this up and “push the parameters back through the function for \(\mathbf{K}\), the covariance matrix” (p. 415).

# compute posterior median covariance among societies
k <- matrix(0, nrow = 10, ncol = 10)
for (i in 1:10)
    for (j in 1:10)
        k[i, j] <- median(post_b13.7$sdgp_gplatlon2^2) * exp(-median(post_b13.7$lscale_gplatlon2^2) * islandsDistMatrix[i, j]^2)

diag(k) <- median(post_b13.7$sdgp_gplatlon2^2) + 0.01

k %>% round(2)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,] 0.19 0.18 0.18 0.08 0.17 0.16 0.12 0.13 0.16  0.05
##  [2,] 0.18 0.19 0.18 0.09 0.17 0.16 0.13 0.14 0.16  0.06
##  [3,] 0.18 0.18 0.19 0.10 0.17 0.16 0.14 0.15 0.15  0.06
##  [4,] 0.08 0.09 0.10 0.19 0.06 0.14 0.17 0.17 0.04  0.02
##  [5,] 0.17 0.17 0.17 0.06 0.19 0.12 0.10 0.10 0.18  0.07
##  [6,] 0.16 0.16 0.16 0.14 0.12 0.19 0.16 0.18 0.10  0.03
##  [7,] 0.12 0.13 0.14 0.17 0.10 0.16 0.19 0.17 0.07  0.05
##  [8,] 0.13 0.14 0.15 0.17 0.10 0.18 0.17 0.19 0.08  0.03
##  [9,] 0.16 0.16 0.15 0.04 0.18 0.10 0.07 0.08 0.19  0.07
## [10,] 0.05 0.06 0.06 0.02 0.07 0.03 0.05 0.03 0.07  0.19

And we’ll continue to follow suit and change these to a correlation matrix.

# convert to correlation matrix
rho <- round(cov2cor(k), 2)
# add row/col names for convenience
colnames(rho) <- c("Ml","Ti","SC","Ya","Fi","Tr","Ch","Mn","To","Ha")
rownames(rho) <- colnames(rho)

rho %>% round(2)
##      Ml   Ti   SC   Ya   Fi   Tr   Ch   Mn   To   Ha
## Ml 1.00 0.94 0.93 0.44 0.89 0.80 0.63 0.69 0.82 0.26
## Ti 0.94 1.00 0.94 0.47 0.89 0.81 0.68 0.71 0.81 0.30
## SC 0.93 0.94 1.00 0.52 0.86 0.84 0.72 0.76 0.77 0.29
## Ya 0.44 0.47 0.52 1.00 0.29 0.74 0.86 0.85 0.20 0.12
## Fi 0.89 0.89 0.86 0.29 1.00 0.62 0.49 0.51 0.93 0.36
## Tr 0.80 0.81 0.84 0.74 0.62 1.00 0.83 0.92 0.51 0.16
## Ch 0.63 0.68 0.72 0.86 0.49 0.83 1.00 0.89 0.37 0.24
## Mn 0.69 0.71 0.76 0.85 0.51 0.92 0.89 1.00 0.40 0.15
## To 0.82 0.81 0.77 0.20 0.93 0.51 0.37 0.40 1.00 0.34
## Ha 0.26 0.30 0.29 0.12 0.36 0.16 0.24 0.15 0.34 1.00

The correlations in our rho matrix look a little higher than those in the text. Before we get see them in a plot, let’s consider psize. If you really want to scale the points in Figure 13.10.a like McElreath did, you can make the psize variable in a tidyverse sort of way as follows. However, if you compare the psize method and the default ggplot2 method using just logpop, you’ll see the difference is negligible. In that light, I’m going to be lazy and just use logpop in my plots.

d %>% 
  transmute(psize = logpop / max(logpop)) %>% 
  transmute(psize = exp(psize * 1.5) - 2)
##        psize
## 1  0.3134090
## 2  0.4009582
## 3  0.6663711
## 4  0.7592196
## 5  0.9066890
## 6  0.9339560
## 7  0.9834797
## 8  1.1096138
## 9  1.2223112
## 10 2.4816891

As far as I can figure, you still have to get rho into a tidy data frame before feeding it into ggplot2. Here’s my attempt at doing so.

tidy_rho <-
  rho %>%
  as.data.frame() %>% 
  rownames_to_column() %>% 
  bind_cols(d %>% select(culture, logpop, total_tools, lon2, lat)) %>% 
  gather(colname, correlation, -rowname, -culture, -logpop, -total_tools, -lon2, -lat) %>% 
  mutate(group = paste(pmin(rowname, colname), pmax(rowname, colname))) %>% 
  select(rowname, colname, group, culture, everything())

head(tidy_rho)
##   rowname colname group    culture   logpop total_tools  lon2   lat correlation
## 1      Ml      Ml Ml Ml   Malekula 7.003065          13 -12.5 -16.3        1.00
## 2      Ti      Ml Ml Ti    Tikopia 7.313220          22 -11.2 -12.3        0.94
## 3      SC      Ml Ml SC Santa Cruz 8.188689          24 -14.0 -10.7        0.93
## 4      Ya      Ml Ml Ya        Yap 8.474494          43 -41.9   9.5        0.44
## 5      Fi      Ml Fi Ml   Lau Fiji 8.909235          33  -1.9 -17.7        0.89
## 6      Tr      Ml Ml Tr  Trobriand 8.987197          19 -29.1  -8.7        0.80

Okay, here’s our version of Figure 13.10.a.

tidy_rho %>%       
  ggplot(aes(x = lon2, y = lat)) +
  geom_line(aes(group = group, alpha = correlation^2),
            color = "#80A0C7") +
  geom_point(data = d, aes(size = logpop), color = "#DCA258") +
  geom_text_repel(data = d, aes(label = culture), 
                  seed = 0, point.padding = .3, size = 3, color = "#FCF9F0") +
  scale_alpha_continuous(range = c(0, 1)) +
  labs(x = "longitude",
       y = "latitude") +
  coord_cartesian(xlim = range(d$lon2),
                  ylim = range(d$lat)) +
  theme(legend.position = "none") +
  theme_pearl_earring 

Yep, as expressed by the intensity of the colors of the connecting lines, those correlations are more pronounced than those in the text. Here’s our version of Figure 13.10.b.

# new data for fitted()
nd <- 
  tibble(logpop = seq(from = 6, to = 14, length.out = 30),
         lat    = median(d$lat),
         lon2   = median(d$lon2))

# fitted()
ftd <-
  fitted(b13.7, newdata = nd) %>% 
  as_tibble() %>% 
  bind_cols(nd)
  
# plot
tidy_rho %>% 
  ggplot(aes(x = logpop)) +
  geom_ribbon(data = ftd,
              aes(ymin = Q2.5, ymax = Q97.5),
              fill = "#394165", alpha = .5) +
  geom_line(data = ftd,
            aes(y = Estimate), color = "#100F14", linetype = 1, size = 1.1) + #  80A0C7 100F14
  geom_line(aes(y = total_tools, group = group, alpha = correlation^2),
            color = "#80A0C7") +
  geom_point(data = d, 
             aes(y = total_tools, size = logpop), color = "#DCA258") +
  geom_text_repel(data = d, 
                  aes(y = total_tools, label = culture), 
                  seed = 0, point.padding = .3, size = 3, color = "#FCF9F0") +
  scale_alpha_continuous(range = c(0, 1)) +
  labs(x = "log population",
       y = "total tools") +
  coord_cartesian(xlim = range(d$logpop),
                  ylim = range(d$total_tools)) +
  theme(legend.position = "none") +
  theme_pearl_earring

Same deal. Our higher correlations make for a more intensely-webbed plot. To learn more on Bürkner’s thoughts on this model in brms, check out the thread on this issue.

13.5 Summary Bonus: Another Berkley-admissions-data-like example.

McElreath uploaded recordings of him teaching out of his text for a graduate course during the 2017/2018 fall semester. In the beginning of lecture 13 from week 7, he discussed a paper from van der Lee and Ellemers (2015) published an article in PNAS. Their paper suggested male researchers were more likely than female researchers to get research funding in the Netherlands. In their initial analysis (p. 12350) they provided a simple \(\chi^2\) test to test the null hypothesis there was no difference in success for male versus female researchers, for which they reported \(\chi_{df = 1}^2 = 4.01, p = .045\). Happily, van der Lee and Ellemers provided their data values in their supplemental material (i.e., Table S1.), which McElreath also displayed in his video.

Their data follows the same structure as the Berkley admissions data. In his lecture, McElreath suggested their \(\chi^2\) test is an example of Simpson’s paradox, just as with the Berkley data. He isn’t the first person to raise this criticism (see Volker and SteenBeek’s critique, which McElreath also pointed to in the lecture).

Here are the data:

funding <- 
  tibble(
    discipline   = rep(c("Chemical sciences", "Physical sciences", "Physics", "Humanities", "Technical sciences",  "Interdisciplinary", "Earth/life sciences", "Social sciences", "Medical sciences"),
                     each = 2),
    gender       = rep(c("m", "f"), times = 9),
    applications = c(83, 39, 135, 39, 67, 9, 230, 166, 189, 62, 105, 78, 156, 126, 425, 409, 245, 260) %>% as.integer(),
    awards       = c(22, 10, 26, 9, 18, 2, 33, 32, 30, 13, 12, 17, 38, 18, 65, 47, 46, 29) %>% as.integer(),
    rejects      = c(61, 29, 109, 30, 49, 7, 197, 134, 159, 49, 93, 61, 118, 108, 360, 362, 199, 231) %>% as.integer(),
    male         = ifelse(gender == "f", 0, 1) %>% as.integer()
  )

funding
## # A tibble: 18 x 6
##    discipline          gender applications awards rejects  male
##    <chr>               <chr>         <int>  <int>   <int> <int>
##  1 Chemical sciences   m                83     22      61     1
##  2 Chemical sciences   f                39     10      29     0
##  3 Physical sciences   m               135     26     109     1
##  4 Physical sciences   f                39      9      30     0
##  5 Physics             m                67     18      49     1
##  6 Physics             f                 9      2       7     0
##  7 Humanities          m               230     33     197     1
##  8 Humanities          f               166     32     134     0
##  9 Technical sciences  m               189     30     159     1
## 10 Technical sciences  f                62     13      49     0
## 11 Interdisciplinary   m               105     12      93     1
## 12 Interdisciplinary   f                78     17      61     0
## 13 Earth/life sciences m               156     38     118     1
## 14 Earth/life sciences f               126     18     108     0
## 15 Social sciences     m               425     65     360     1
## 16 Social sciences     f               409     47     362     0
## 17 Medical sciences    m               245     46     199     1
## 18 Medical sciences    f               260     29     231     0

Let’s fit a few models.

First, we’ll fit an analogue to the initial van der Lee and Ellemers \(\chi^2\) test. Since we’re Bayesian modelers, we’ll use a simple logistic regression, using male (dummy coded 0 = female, 1 = male) to predict admission (i.e., awards).

b13.bonus_0 <- 
  brm(data = funding, family = binomial,
      awards | trials(applications) ~ 1 + male,
      # Note our continued use of weakly-regularizing priors
      prior = c(prior(normal(0, 4), class = Intercept),
                prior(normal(0, 4), class = b)),
      iter = 5000, warmup = 1000, chains = 4, cores = 4)

If you inspect them, the chains look great. Here are the posterior summaries:

tidy(b13.bonus_0) %>%
  filter(term != "lp__") %>%
  mutate_if(is.numeric, round, digits = 2)
##          term estimate std.error lower upper
## 1 b_Intercept    -1.75      0.08 -1.88 -1.61
## 2      b_male     0.21      0.10  0.04  0.38

Yep, the 95% intervals for male dummy exclude zero. If you wanted a one-sided Bayesian \(p\)-value, you might do something like:

posterior_samples(b13.bonus_0) %>%
  summarise(One_sided_Bayesian_p_value = filter(., b_male <= 0) %>% nrow()/nrow(.))
##   One_sided_Bayesian_p_value
## 1                  0.0209375

Pretty small. But recall how Simpson’s paradox helped us understand the Berkley data. Different departments in Berkley had different acceptance rates AND different ratios of male and female applicants. Similarly, different academic disciplines in the Netherlands might have different award rates for funding AND different ratios of male and female applications.

Just like in section 13.2, let’s fit two more models. The first model will allow intercepts to vary by discipline. The second model will allow intercepts and the male dummy slopes to vary by discipline.

b13.bonus_1 <- 
  brm(data = funding, family = binomial,
      awards | trials(applications) ~ 1 + male + (1 | discipline),
      prior = c(prior(normal(0, 4), class = Intercept),
                prior(normal(0, 4), class = b),
                prior(cauchy(0, 1), class = sd)),
      iter = 5000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = .99))

b13.bonus_2 <- 
  brm(data = funding, family = binomial,
      awards | trials(applications) ~ 1 + male + (1 + male | discipline),
      prior = c(prior(normal(0, 4), class = Intercept),
                prior(normal(0, 4), class = b),
                prior(cauchy(0, 1), class = sd),
                prior(lkj(4), class = cor)),
      iter = 5000, warmup = 1000, chains = 4, cores = 4,
      control = list(adapt_delta = .99))

We’ll compare the models with information criteria.

waic(b13.bonus_0, b13.bonus_1, b13.bonus_2)
##                             WAIC   SE
## b13.bonus_0               129.69 8.90
## b13.bonus_1               125.74 7.35
## b13.bonus_2               116.61 5.60
## b13.bonus_0 - b13.bonus_1   3.95 6.29
## b13.bonus_0 - b13.bonus_2  13.09 5.60
## b13.bonus_1 - b13.bonus_2   9.14 2.76

The WAIC suggests the varying intercepts/varying slopes model made the best sense of the data. Here’s what the random intercepts look like in a coefficient plot.

coef(b13.bonus_2)$discipline[, , 2] %>% 
  as_tibble() %>% 
  mutate(discipline = c("Chemical sciences", "Physical sciences",
                        "Physics", "Humanities", "Technical sciences",
                        "Interdisciplinary", "Earth/life sciences",
                        "Social sciences", "Medical sciences")) %>%
  
  ggplot(aes(x = reorder(discipline, Estimate), y = Estimate,
             ymin = Q2.5,
             ymax = Q97.5)) +
  geom_hline(yintercept = 0, color = "#E8DCCF", alpha = 1/4) +
  geom_hline(yintercept = fixef(b13.bonus_2)[2], linetype = 3, color = "#80A0C7") +
  geom_pointrange(shape = 20, size = 3/4, color = "#8B9DAF") +
  labs(title    = "Random slopes for the male dummy",
       subtitle = "The vertical dotted line is the posterior mean of the fixed effect for the\nmale dummy. The dots and horizontal lines are the posterior means and\npercentile-based 95% intervals, respectively. The values are on the log scale.",
       x        = NULL, y = NULL) +
  coord_flip(ylim = -1:1) +
  theme_pearl_earring +
  theme(axis.ticks.y = element_blank(),
        axis.text.y  = element_text(hjust = 0))

Note how the 95% intervals for all the random male slopes contain zero within their bounds. Here are the fixed effects:

tidy(b13.bonus_2) %>%
  filter(str_detect(term , "b_")) %>%
  mutate_if(is.numeric, round, digits = 2)
##          term estimate std.error lower upper
## 1 b_Intercept    -1.63      0.14 -1.85 -1.38
## 2      b_male     0.15      0.17 -0.13  0.42

And if you wanted a one-sided Bayesian \(p\)-value for the male dummy for the full model:

posterior_samples(b13.bonus_2) %>%
  summarise(One_sided_Bayesian_p_value = filter(., b_male <= 0) %>% nrow()/nrow(.))
##   One_sided_Bayesian_p_value
## 1                    0.17975

So, the estimate of the gender bias is small and consistent with the null hypothesis. Which is good! We want gender equality for things like funding success.

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] brms_2.5.0         bayesplot_1.6.0    ggbeeswarm_0.6.0   ggrepel_0.8.0      tidybayes_1.0.1   
##  [6] broom_0.4.5        Rcpp_0.12.18       rstan_2.17.3       StanHeaders_2.17.2 dutchmasters_0.1.0
## [11] forcats_0.3.0      stringr_1.3.1      dplyr_0.7.6        purrr_0.2.5        readr_1.1.1       
## [16] tidyr_0.8.1        tibble_1.4.2       ggplot2_3.0.0      tidyverse_1.2.1   
## 
## 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] labeling_0.3              mnormt_1.5-5              bridgesampling_0.4-0     
##  [22] rprojroot_1.3-2           coda_0.19-1               xfun_0.3                 
##  [25] R6_2.2.2                  markdown_0.8              HDInterval_0.2.0         
##  [28] reshape_0.8.7             assertthat_0.2.0          promises_1.0.1           
##  [31] scales_0.5.0              beeswarm_0.2.3            gtable_0.2.0             
##  [34] rlang_0.2.1               extrafontdb_1.0           lazyeval_0.2.1           
##  [37] inline_0.3.15             yaml_2.1.19               reshape2_1.4.3           
##  [40] abind_1.4-5               modelr_0.1.2              threejs_0.3.1            
##  [43] crosstalk_1.0.0           backports_1.1.2           httpuv_1.4.4.2           
##  [46] rsconnect_0.8.8           extrafont_0.17            tools_3.5.1              
##  [49] bookdown_0.7              psych_1.8.4               RColorBrewer_1.1-2       
##  [52] ggridges_0.5.0            plyr_1.8.4                base64enc_0.1-3          
##  [55] progress_1.2.0            prettyunits_1.0.2         zoo_1.8-2                
##  [58] LaplacesDemon_16.1.1      haven_1.1.2               magrittr_1.5             
##  [61] colourpicker_1.0          mvtnorm_1.0-8             matrixStats_0.54.0       
##  [64] hms_0.4.2                 shinyjs_1.0               mime_0.5                 
##  [67] evaluate_0.10.1           arrayhelpers_1.0-20160527 xtable_1.8-2             
##  [70] shinystan_2.5.0           readxl_1.1.0              gridExtra_2.3            
##  [73] rstantools_1.5.0          compiler_3.5.1            maps_3.3.0               
##  [76] crayon_1.3.4              htmltools_0.3.6           later_0.7.3              
##  [79] lubridate_1.7.4           MASS_7.3-50               Matrix_1.2-14            
##  [82] cli_1.0.0                 bindr_0.1.1               igraph_1.2.1             
##  [85] pkgconfig_2.0.1           foreign_0.8-70            xml2_1.2.0               
##  [88] svUnit_0.7-12             dygraphs_1.1.1.5          vipor_0.4.5              
##  [91] rvest_0.3.2               digest_0.6.15             rmarkdown_1.10           
##  [94] cellranger_1.1.0          shiny_1.1.0               gtools_3.8.1             
##  [97] nlme_3.1-137              jsonlite_1.5              bindrcpp_0.2.2           
## [100] mapproj_1.2.6             viridisLite_0.3.0         pillar_1.2.3             
## [103] lattice_0.20-35           loo_2.0.0                 httr_1.3.1               
## [106] glue_1.2.0                xts_0.10-2                shinythemes_1.1.1        
## [109] pander_0.6.2              stringi_1.2.3