In joepowers16/rethinking: Statistical Rethinking book package
options(width = 100)
For our plots in this document, 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
We'll name our custom theme theme_pearl_earring. Here we make it.
library(tidyverse)
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())
13.1. Varying slopes by construction
13.1.1. Simulate the population.
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 <- (-0.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 with this
matrix(c(1, 2,
3, 4), nrow = 2, ncol = 2)
sigmas <- c(sigma_a, sigma_b) # standard deviations
Rho <- matrix(c(1, rho,
rho, 1), nrow = 2) # correlation matrix
# now matrix multiply to get covariance matrix
Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas)
N_cafes <- 20
library(MASS)
set.seed(5) # used to replicate example
vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
What's the distribution of a_cafe and b_cafe?, you might ask.
tibble(a_cafe = vary_effects[ , 1],
b_cafe = vary_effects[ , 2]) %>%
ggplot(aes(x = a_cafe, y = b_cafe)) +
geom_point(color = "#80A0C7") +
geom_rug(color = "#8B9DAF", size = 1/7) +
scale_x_continuous(expand = c(.5, .5)) +
theme_pearl_earring
13.1.2. Simulate observations.
Before we simulate our observations, we'll need to detach the MASS package and reload the tidyverse in order to use the tidyverse::select() function.
detach(package:MASS, unload = T)
library(tidyverse)
N_visits <- 10
sigma <- 0.5 # std dev within cafes
set.seed(5) # used to replicate example
d <-
tibble(cafe = rep(1:N_cafes, each = N_visits),
afternoon = rep(0:1, N_visits*N_cafes/2),
mu = rep(vary_effects[ , 1], each = N_visits) + rep(vary_effects[ , 2], each = N_visits)*afternoon,
wait = rnorm(N_visits*N_cafes, mu, sigma)) %>%
select(-mu)
Here's a look at the data.
d %>%
glimpse()
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),
x_order = rep(1:10, times = N_cafes)) %>%
filter(cafe %in% c(3, 5)) %>%
mutate(cafe = ifelse(cafe == 3, "cafe #3", "cafe #5")) %>%
ggplot(aes(x = x_order, y = wait, group = day)) +
geom_point(aes(color = afternoon), size = 2) +
scale_color_manual(values = c("#80A0C7", "#EEDA9D")) +
geom_line(color = "#8B9DAF") +
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") +
facet_wrap(~cafe, ncol = 1) +
theme_pearl_earring +
theme(legend.position = "none")
Here are the simulations for Figure 13.3.
library(rethinking)
set.seed(133)
R_1 <-
rlkjcorr(1e5, K = 2, eta = 1) %>%
as_tibble()
set.seed(133)
R_2 <-
rlkjcorr(1e5, K = 2, eta = 2) %>%
as_tibble()
set.seed(133)
R_4 <-
rlkjcorr(1e5, K = 2, eta = 4) %>%
as_tibble()
The code for 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) +
labs(x = "correlation") +
theme_pearl_earring
Our first model has both varying intercepts and afternoon slopes.
detach(package:rethinking, unload = T)
library(brms)
b13.1 <-
brm(data = d, family = gaussian,
wait ~ 1 + afternoon + (1 + afternoon | cafe),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 2)", class = "sd"),
set_prior("cauchy(0, 2)", class = "sigma"),
set_prior("lkj(2)", class = "cor")),
iter = 5000, warmup = 2000, chains = 2, cores = 2)
Figure 13.4.
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 = 1/2) +
geom_density(color = "transparent", fill = "#A65141", alpha = 3/4) +
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) +
labs(x = "correlation") +
theme_pearl_earring
McElreath then depicts 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)
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() %>% # Converting the two vectors to a tibble
rename(Slope = afternoon) %>%
mutate(cafe = 1:nrow(.)) %>% # Adding the `cafe` index
select(cafe, everything()) # simply moving `cafe` to the left-most column of the tibble
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) %>% # using 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) %>%
mutate(pooled = rep(c("partially", "not"), each = nrow(.)/2)) # indexing whether the estimates are pooled
# Here's a glimpse at what we've been working for
params %>%
slice(c(1:5, 36:40))
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
Here we 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() %>% # Converting 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:nrow(.)) %>% # Adding 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 = nrow(.)/2))
The code for Figure 13.5.b.
ggplot(data = estimates, aes(x = morning, y = afternoon)) +
# Nesting stat_ellipse() with 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.
We start by only letting the intercepts vary in this one.
b13.2 <-
brm(data = d, family = binomial,
admit | trials(applications) ~ 1 + male + (1 | dept_id),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b"),
set_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
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, presents 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()
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)
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"])
Note how coef() returned a three-dimensional list.
coef(b13.2) %>% str()
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]
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())
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)
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
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 our male dummy varies, too.
b13.3 <-
brm(data = d, family = binomial,
admit | trials(applications) ~ 1 + male + (1 + male | dept_id),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b"),
set_prior("cauchy(0, 2)", class = "sd"),
set_prior("lkj(2)", class = "cor")),
iter = 5000, warmup = 1000, chains = 4, cores = 4,
control = list(adapt_delta = .99,
max_treedepth = 12))
The random effects in the centered metric:
coef(b13.3)
It just takes a little data wrangling to put the brms-based centered random effects into a tidy tibble with which we might make a coefficient plot, like McElreath did on page 401.
# 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 into tibbles, and then 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)) %>%
# The plot
ggplot(aes(x = reorder(param, reorder))) +
geom_hline(yintercept = 0, linetype = 3, color = "#E8DCCF") +
geom_linerange(aes(ymin = `2.5%ile`, ymax = `97.5%ile`),
color = "#E7CDC2") +
geom_point(aes(y = Estimate),
color = "#A65141") +
xlab(NULL) +
coord_flip() +
theme_pearl_earring +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0))
13.2.3. Shrinkage.
Figure 13.6.a., the correlation between the full UCB model's varying intercepts and slopes.
post <- posterior_samples(b13.3)
post %>%
ggplot(aes(x = cor_dept_id__Intercept__male)) +
geom_density(color = "transparent", fill = "#8B9DAF") +
geom_vline(xintercept = median(post$cor_dept_id__Intercept__male), color = "#394165") +
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 line is at the median.",
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:nrow(.)) %>%
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 Gelman and Hill's invlogit() 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 custom 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 = nrow(.)/2)) %>%
mutate(dept_letter = rep(LETTERS[1:6], times = 2)) # This will help with plotting
params
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 package and its geom_text_repel() function will help us with the in-plot labels.
library(ggrepel)
set.seed(6457240) # for ggrepel::geom_text_repel()
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") +
coord_cartesian(xlim = range(params$intercept),
ylim = range(params$slope)) +
theme_pearl_earring
13.2.4. Model comparison.
Our no-gender model:
b13.4 <-
brm(data = d, family = binomial,
admit | trials(applications) ~ 1 + (1 | dept_id),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("cauchy(0, 2)", class = "sd")),
iter = 5000, warmup = 1000, chains = 4, cores = 4,
control = list(adapt_delta = .99,
max_treedepth = 12))
waic(b13.2, b13.3, b13.4)
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)
Here's our cross-classified model. For you SEM lovers, this reminds me of a factor model with a method effect (e.g., a bifactor model).
b13.6 <-
brm(data = d, family = binomial,
pulled_left ~ 1 + prosoc_left + condition:prosoc_left +
(1 + prosoc_left + condition:prosoc_left | block_id) +
(1 + prosoc_left + condition:prosoc_left | actor),
prior = c(set_prior("normal(0, 1)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b"),
set_prior("cauchy(0, 2)", class = "sd"),
set_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.
ratios_cp <- neff_ratio(b13.6)
neff <-
ratios_cp %>%
as_tibble %>%
rename(neff_ratio = value) %>%
mutate(neff = neff_ratio*12000)
head(neff)
Our variant of Figure 13.7. The handy ggbeeswarm package and it's 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
The bayesplot package contains a sweet of handy diagnostic features. mcmc_neff(), for example, makes it easy 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)
Here we 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 | block_id) +
(1 | actor),
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 10)", class = "b"),
set_prior("cauchy(0, 1)", class = "sd")),
iter = 5000, warmup = 1000, chains = 3, cores = 3)
The WAIC comparison:
waic(b13.6, b12.5)
13.4. Continuous categories and the Gaussian process
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
Dmat <- islandsDistMatrix
colnames(Dmat) <- c("Ml", "Ti", "SC", "Ya", "Fi",
"Tr", "Ch", "Mn", "To", "Ha")
round(Dmat, 1)
If you wanted to use color to more effectively visualize the values in the matirx, you might do something like this.
Dmat %>%
as_tibble() %>%
gather() %>%
rename(Column = key,
distance = value) %>%
mutate(Row = rep(rownames(Dmat), 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., 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(expand = c(0, 0)) +
scale_y_continuous(breaks = c(0, .5, 1),
labels = c(0, ".5", 1)) +
labs(x = "distance",
y = "correlation") +
theme_pearl_earring
data(Kline2) # load the ordinary data, now with coordinates
d <-
Kline2 %>%
mutate(society = 1:10)
d %>% glimpse()
Okay, it appears this is going to be a bit of a ride. It's not entire 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. Note how our priors are a little different. I'll explain, below. Let's just move ahead and fit the model.
detach(package:rethinking, unload = T)
library(brms)
b13.7 <-
brm(data = d, family = poisson,
total_tools ~ 1 + gp(lat, lon2) + logpop,
prior = c(set_prior("normal(0, 10)", class = "Intercept"),
set_prior("normal(0, 1)", class = "b"),
set_prior("inv_gamma(2.874624, 0.393695)", class = "lscale"),
set_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.
posterior_summary(b13.7) %>%
round(digits = 2)
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^2exp(-||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) %>%
summarise_at(vars(sdgp_squared), funs(mean, sd)) %>%
round(digits = 2)
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 whereas 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*(lscale^2))$$
But that doesn't appear to be the case. Sigh.
posterior_samples(b13.7) %>%
transmute(rho_squared = 1/(2*(lscale_gplatlon2^2))) %>%
summarise_at(vars(rho_squared), funs(mean, sd)) %>%
round(digits = 2)
Oh man, that isn't even close to the 2.67 (51.60) 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). Here's where that came from:
get_prior(data = d, family = poisson,
total_tools ~ 1 + gp(lat, lon2) + logpop)
That is, we used the brms defualt 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 the prior.
Now that we’ve started, we may as well keep going down the comparison train. Let’s reproduce McElreath’s model with 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.
post_m13.7 <- rethinking::extract.samples(m13.7)[2:5] %>% as_tibble()
detach(package:rethinking, unload = T)
library(brms)
post_b13.7 <- posterior_samples(b13.7)
Here's the model intercept, by package:
post_m13.7[, "a"] %>%
bind_rows(post_b13.7[, "b_Intercept"] %>%
as_tibble() %>%
rename(a = value)) %>%
mutate(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>%
ggplot(aes(x = a, fill = model)) +
geom_density(size = 0, alpha = 1/2) +
scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>%
ggplot(aes(x = bp, fill = model)) +
geom_density(size = 0, alpha = 1/2) +
scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>%
ggplot(aes(x = etasq, fill = model)) +
geom_density(size = 0, alpha = 1/2) +
scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>%
ggplot(aes(x = rhosq, fill = model)) +
geom_density(size = 0) +
scale_fill_manual(values = c("#80A0C7", "#A65141")) +
labs(title = "Holy smokes are those not the same!",
subtitle = "Notice how differently the y axes got scaled. Also, the 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(~model, 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. The 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)
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)
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)
Okay, here's our Figure 13.10.a.
library(ggrepel)
# 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 hackish attempt at doing so. I’m sure there are more elegant ways to do this.
# If you’ve got ideas, share your code.
Rho %>%
as_tibble() %>%
gather() %>%
mutate(culture = rep(unique(d$culture), each = 10) %>% as.character(),
culture_2 = rep(unique(d$culture), times = 10) %>% as.character()) %>%
mutate(group = paste(pmin(culture, culture_2), pmax(culture, culture_2)),
culture = culture %>% as.factor(),
rev_value = max(value) - value) %>%
left_join(d, by = "culture") %>%
# The plot
ggplot(aes(x = lon2, y = lat)) +
geom_line(aes(group = group, alpha = value^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),
y = 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. Here's our 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)
# making Rho tidy
Rho %>%
as_tibble() %>%
gather() %>%
mutate(culture = rep(unique(d$culture), each = 10) %>% as.character(),
culture_2 = rep(unique(d$culture), times = 10) %>% as.character()) %>%
mutate(group = paste(pmin(culture, culture_2), pmax(culture, culture_2)),
culture = culture %>% as.factor(),
rev_value = max(value) - value) %>%
left_join(d, by = "culture") %>%
# The plot
ggplot(aes(x = logpop)) +
geom_ribbon(data = ftd,
aes(ymin = `2.5%ile`, ymax = `97.5%ile`),
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 = value^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),
y = 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.
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^2$ (1) = 4.01, P = 0.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
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(set_prior("normal(0, 4)", class = "Intercept"),
set_prior("normal(0, 4)", class = "b")),
iter = 5000, warmup = 1000, chains = 4, cores = 4)
The chains look great. Here are the posterior summaries:
tidy(b13.bonus_0) %>%
filter(term != "lp__") %>%
mutate_if(is.numeric, round, digits = 2)
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(.))
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(set_prior("normal(0, 4)", class = "Intercept"),
set_prior("normal(0, 4)", class = "b"),
set_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(set_prior("normal(0, 4)", class = "Intercept"),
set_prior("normal(0, 4)", class = "b"),
set_prior("cauchy(0, 1)", class = "sd"),
set_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)
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 = discipline, y = Estimate,
ymin = `2.5%ile`,
ymax = `97.5%ile`)) +
geom_hline(yintercept = 0, color = "#E8DCCF", size = 1/10) +
geom_hline(yintercept = fixef(b13.bonus_2)[2], linetype = 3, color = "#A65141") +
geom_pointrange(shape = 20, color = "#A65141") +
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)
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(.))
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.
Note. The analyses in this document were done with:
- R 3.4.4
- RStudio 1.1.442
- rmarkdown 1.9
- dutchmasters 0.1.0
- tidyverse 1.2.1
- MASS 7.3-47
- rethinking 1.59
- brms 2.2.0
- rstan 2.17.3
- broom 0.4.2
- ggrepel 0.7.0
- ggbeeswarm 0.5.3
- bayesplot 1.5.0
Reference
McElreath, R. (2016). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman & Hall/CRC Press.
rm(d, Rho, Sigma, theme_pearl_earring, vary_effects, a, b, cov_ab, islandsDistMatrix, Mu, N_cafes, N_visits, neff, ratios_cp, rho, R_1, R_2, R_4, sigma, sigma_a, sigma_b, sigmas, post, b13.1, b13.2, b13.3, b13.4, b13.6, b12.5, partially_pooled_params, un_pooled_params, params, partially_pooled_estimates, un_pooled_estimates, estimates, Dmat, K, nd, ftd, b13.7, m13.7, post_b13.7, post_m13.7, prob_to_logit, funding, b13.bonus_0, b13.bonus_1, b13.bonus_2)
joepowers16/rethinking documentation built on June 2, 2019, 6:52 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
options(width = 100)
For our plots in this document, 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
We'll name our custom theme theme_pearl_earring. Here we make it.
library(tidyverse) 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())
13.1. Varying slopes by construction
13.1.1. Simulate the population.
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 <- (-0.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 with this matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2)
sigmas <- c(sigma_a, sigma_b) # standard deviations Rho <- matrix(c(1, rho, rho, 1), nrow = 2) # correlation matrix # now matrix multiply to get covariance matrix Sigma <- diag(sigmas) %*% Rho %*% diag(sigmas) N_cafes <- 20 library(MASS) set.seed(5) # used to replicate example vary_effects <- mvrnorm(N_cafes, Mu, Sigma)
What's the distribution of a_cafe and b_cafe?, you might ask.
tibble(a_cafe = vary_effects[ , 1], b_cafe = vary_effects[ , 2]) %>% ggplot(aes(x = a_cafe, y = b_cafe)) + geom_point(color = "#80A0C7") + geom_rug(color = "#8B9DAF", size = 1/7) + scale_x_continuous(expand = c(.5, .5)) + theme_pearl_earring
13.1.2. Simulate observations.
Before we simulate our observations, we'll need to detach the MASS package and reload the tidyverse in order to use the tidyverse::select() function.
detach(package:MASS, unload = T) library(tidyverse) N_visits <- 10 sigma <- 0.5 # std dev within cafes set.seed(5) # used to replicate example d <- tibble(cafe = rep(1:N_cafes, each = N_visits), afternoon = rep(0:1, N_visits*N_cafes/2), mu = rep(vary_effects[ , 1], each = N_visits) + rep(vary_effects[ , 2], each = N_visits)*afternoon, wait = rnorm(N_visits*N_cafes, mu, sigma)) %>% select(-mu)
Here's a look at the data.
d %>% glimpse()
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), x_order = rep(1:10, times = N_cafes)) %>% filter(cafe %in% c(3, 5)) %>% mutate(cafe = ifelse(cafe == 3, "cafe #3", "cafe #5")) %>% ggplot(aes(x = x_order, y = wait, group = day)) + geom_point(aes(color = afternoon), size = 2) + scale_color_manual(values = c("#80A0C7", "#EEDA9D")) + geom_line(color = "#8B9DAF") + 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") + facet_wrap(~cafe, ncol = 1) + theme_pearl_earring + theme(legend.position = "none")
Here are the simulations for Figure 13.3.
library(rethinking) set.seed(133) R_1 <- rlkjcorr(1e5, K = 2, eta = 1) %>% as_tibble() set.seed(133) R_2 <- rlkjcorr(1e5, K = 2, eta = 2) %>% as_tibble() set.seed(133) R_4 <- rlkjcorr(1e5, K = 2, eta = 4) %>% as_tibble()
The code for 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) + labs(x = "correlation") + theme_pearl_earring
Our first model has both varying intercepts and afternoon slopes.
detach(package:rethinking, unload = T) library(brms) b13.1 <- brm(data = d, family = gaussian, wait ~ 1 + afternoon + (1 + afternoon | cafe), prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("normal(0, 10)", class = "b"), set_prior("cauchy(0, 2)", class = "sd"), set_prior("cauchy(0, 2)", class = "sigma"), set_prior("lkj(2)", class = "cor")), iter = 5000, warmup = 2000, chains = 2, cores = 2)
Figure 13.4.
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 = 1/2) + geom_density(color = "transparent", fill = "#A65141", alpha = 3/4) + 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) + labs(x = "correlation") + theme_pearl_earring
McElreath then depicts 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)
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() %>% # Converting the two vectors to a tibble rename(Slope = afternoon) %>% mutate(cafe = 1:nrow(.)) %>% # Adding the `cafe` index select(cafe, everything()) # simply moving `cafe` to the left-most column of the tibble
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) %>% # using 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) %>% mutate(pooled = rep(c("partially", "not"), each = nrow(.)/2)) # indexing whether the estimates are pooled # Here's a glimpse at what we've been working for params %>% slice(c(1:5, 36:40))
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
Here we 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() %>% # Converting 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:nrow(.)) %>% # Adding 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 = nrow(.)/2))
The code for Figure 13.5.b.
ggplot(data = estimates, aes(x = morning, y = afternoon)) + # Nesting stat_ellipse() with 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.
We start by only letting the intercepts vary in this one.
b13.2 <- brm(data = d, family = binomial, admit | trials(applications) ~ 1 + male + (1 | dept_id), prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("normal(0, 1)", class = "b"), set_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
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, presents 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()
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)
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"])
Note how coef() returned a three-dimensional list.
coef(b13.2) %>% str()
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]
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())
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)
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
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 our male dummy varies, too.
b13.3 <- brm(data = d, family = binomial, admit | trials(applications) ~ 1 + male + (1 + male | dept_id), prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("normal(0, 1)", class = "b"), set_prior("cauchy(0, 2)", class = "sd"), set_prior("lkj(2)", class = "cor")), iter = 5000, warmup = 1000, chains = 4, cores = 4, control = list(adapt_delta = .99, max_treedepth = 12))
The random effects in the centered metric:
coef(b13.3)
It just takes a little data wrangling to put the brms-based centered random effects into a tidy tibble with which we might make a coefficient plot, like McElreath did on page 401.
# 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 into tibbles, and then 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)) %>% # The plot ggplot(aes(x = reorder(param, reorder))) + geom_hline(yintercept = 0, linetype = 3, color = "#E8DCCF") + geom_linerange(aes(ymin = `2.5%ile`, ymax = `97.5%ile`), color = "#E7CDC2") + geom_point(aes(y = Estimate), color = "#A65141") + xlab(NULL) + coord_flip() + theme_pearl_earring + theme(axis.ticks.y = element_blank(), axis.text.y = element_text(hjust = 0))
13.2.3. Shrinkage.
Figure 13.6.a., the correlation between the full UCB model's varying intercepts and slopes.
post <- posterior_samples(b13.3) post %>% ggplot(aes(x = cor_dept_id__Intercept__male)) + geom_density(color = "transparent", fill = "#8B9DAF") + geom_vline(xintercept = median(post$cor_dept_id__Intercept__male), color = "#394165") + 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 line is at the median.", 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:nrow(.)) %>% 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 Gelman and Hill's invlogit() 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 custom 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 = nrow(.)/2)) %>% mutate(dept_letter = rep(LETTERS[1:6], times = 2)) # This will help with plotting params
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 package and its geom_text_repel() function will help us with the in-plot labels.
library(ggrepel) set.seed(6457240) # for ggrepel::geom_text_repel() 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") + coord_cartesian(xlim = range(params$intercept), ylim = range(params$slope)) + theme_pearl_earring
13.2.4. Model comparison.
Our no-gender model:
b13.4 <- brm(data = d, family = binomial, admit | trials(applications) ~ 1 + (1 | dept_id), prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("cauchy(0, 2)", class = "sd")), iter = 5000, warmup = 1000, chains = 4, cores = 4, control = list(adapt_delta = .99, max_treedepth = 12))
waic(b13.2, b13.3, b13.4)
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)
Here's our cross-classified model. For you SEM lovers, this reminds me of a factor model with a method effect (e.g., a bifactor model).
b13.6 <- brm(data = d, family = binomial, pulled_left ~ 1 + prosoc_left + condition:prosoc_left + (1 + prosoc_left + condition:prosoc_left | block_id) + (1 + prosoc_left + condition:prosoc_left | actor), prior = c(set_prior("normal(0, 1)", class = "Intercept"), set_prior("normal(0, 1)", class = "b"), set_prior("cauchy(0, 2)", class = "sd"), set_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.
ratios_cp <- neff_ratio(b13.6) neff <- ratios_cp %>% as_tibble %>% rename(neff_ratio = value) %>% mutate(neff = neff_ratio*12000) head(neff)
Our variant of Figure 13.7. The handy ggbeeswarm package and it's 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
The bayesplot package contains a sweet of handy diagnostic features. mcmc_neff(), for example, makes it easy 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)
Here we 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 | block_id) + (1 | actor), prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("normal(0, 10)", class = "b"), set_prior("cauchy(0, 1)", class = "sd")), iter = 5000, warmup = 1000, chains = 3, cores = 3)
The WAIC comparison:
waic(b13.6, b12.5)
13.4. Continuous categories and the Gaussian process
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 Dmat <- islandsDistMatrix colnames(Dmat) <- c("Ml", "Ti", "SC", "Ya", "Fi", "Tr", "Ch", "Mn", "To", "Ha") round(Dmat, 1)
If you wanted to use color to more effectively visualize the values in the matirx, you might do something like this.
Dmat %>% as_tibble() %>% gather() %>% rename(Column = key, distance = value) %>% mutate(Row = rep(rownames(Dmat), 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., 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(expand = c(0, 0)) + scale_y_continuous(breaks = c(0, .5, 1), labels = c(0, ".5", 1)) + labs(x = "distance", y = "correlation") + theme_pearl_earring
data(Kline2) # load the ordinary data, now with coordinates d <- Kline2 %>% mutate(society = 1:10) d %>% glimpse()
Okay, it appears this is going to be a bit of a ride. It's not entire 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. Note how our priors are a little different. I'll explain, below. Let's just move ahead and fit the model.
detach(package:rethinking, unload = T) library(brms) b13.7 <- brm(data = d, family = poisson, total_tools ~ 1 + gp(lat, lon2) + logpop, prior = c(set_prior("normal(0, 10)", class = "Intercept"), set_prior("normal(0, 1)", class = "b"), set_prior("inv_gamma(2.874624, 0.393695)", class = "lscale"), set_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.
posterior_summary(b13.7) %>% round(digits = 2)
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^2exp(-||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) %>% summarise_at(vars(sdgp_squared), funs(mean, sd)) %>% round(digits = 2)
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 whereas 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*(lscale^2))$$
But that doesn't appear to be the case. Sigh.
posterior_samples(b13.7) %>% transmute(rho_squared = 1/(2*(lscale_gplatlon2^2))) %>% summarise_at(vars(rho_squared), funs(mean, sd)) %>% round(digits = 2)
Oh man, that isn't even close to the 2.67 (51.60) 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). Here's where that came from:
get_prior(data = d, family = poisson, total_tools ~ 1 + gp(lat, lon2) + logpop)
That is, we used the brms defualt 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 the prior.
Now that we’ve started, we may as well keep going down the comparison train. Let’s reproduce McElreath’s model with 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.
post_m13.7 <- rethinking::extract.samples(m13.7)[2:5] %>% as_tibble() detach(package:rethinking, unload = T) library(brms) post_b13.7 <- posterior_samples(b13.7)
Here's the model intercept, by package:
post_m13.7[, "a"] %>% bind_rows(post_b13.7[, "b_Intercept"] %>% as_tibble() %>% rename(a = value)) %>% mutate(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>% ggplot(aes(x = a, fill = model)) + geom_density(size = 0, alpha = 1/2) + scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>% ggplot(aes(x = bp, fill = model)) + geom_density(size = 0, alpha = 1/2) + scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>% ggplot(aes(x = etasq, fill = model)) + geom_density(size = 0, alpha = 1/2) + scale_fill_manual(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(model = rep(c("m13.7", "b13.7"), each = nrow(post_m13.7))) %>% ggplot(aes(x = rhosq, fill = model)) + geom_density(size = 0) + scale_fill_manual(values = c("#80A0C7", "#A65141")) + labs(title = "Holy smokes are those not the same!", subtitle = "Notice how differently the y axes got scaled. Also, the 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(~model, 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. The 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)
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)
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)
Okay, here's our Figure 13.10.a.
library(ggrepel) # 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 hackish attempt at doing so. I’m sure there are more elegant ways to do this. # If you’ve got ideas, share your code. Rho %>% as_tibble() %>% gather() %>% mutate(culture = rep(unique(d$culture), each = 10) %>% as.character(), culture_2 = rep(unique(d$culture), times = 10) %>% as.character()) %>% mutate(group = paste(pmin(culture, culture_2), pmax(culture, culture_2)), culture = culture %>% as.factor(), rev_value = max(value) - value) %>% left_join(d, by = "culture") %>% # The plot ggplot(aes(x = lon2, y = lat)) + geom_line(aes(group = group, alpha = value^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), y = 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. Here's our 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) # making Rho tidy Rho %>% as_tibble() %>% gather() %>% mutate(culture = rep(unique(d$culture), each = 10) %>% as.character(), culture_2 = rep(unique(d$culture), times = 10) %>% as.character()) %>% mutate(group = paste(pmin(culture, culture_2), pmax(culture, culture_2)), culture = culture %>% as.factor(), rev_value = max(value) - value) %>% left_join(d, by = "culture") %>% # The plot ggplot(aes(x = logpop)) + geom_ribbon(data = ftd, aes(ymin = `2.5%ile`, ymax = `97.5%ile`), 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 = value^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), y = 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.
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^2$ (1) = 4.01, P = 0.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
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(set_prior("normal(0, 4)", class = "Intercept"), set_prior("normal(0, 4)", class = "b")), iter = 5000, warmup = 1000, chains = 4, cores = 4)
The chains look great. Here are the posterior summaries:
tidy(b13.bonus_0) %>% filter(term != "lp__") %>% mutate_if(is.numeric, round, digits = 2)
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(.))
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(set_prior("normal(0, 4)", class = "Intercept"), set_prior("normal(0, 4)", class = "b"), set_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(set_prior("normal(0, 4)", class = "Intercept"), set_prior("normal(0, 4)", class = "b"), set_prior("cauchy(0, 1)", class = "sd"), set_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)
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 = discipline, y = Estimate, ymin = `2.5%ile`, ymax = `97.5%ile`)) + geom_hline(yintercept = 0, color = "#E8DCCF", size = 1/10) + geom_hline(yintercept = fixef(b13.bonus_2)[2], linetype = 3, color = "#A65141") + geom_pointrange(shape = 20, color = "#A65141") + 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)
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(.))
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.
Note. The analyses in this document were done with:
- R 3.4.4
- RStudio 1.1.442
- rmarkdown 1.9
- dutchmasters 0.1.0
- tidyverse 1.2.1
- MASS 7.3-47
- rethinking 1.59
- brms 2.2.0
- rstan 2.17.3
- broom 0.4.2
- ggrepel 0.7.0
- ggbeeswarm 0.5.3
- bayesplot 1.5.0
Reference
McElreath, R. (2016). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman & Hall/CRC Press.
rm(d, Rho, Sigma, theme_pearl_earring, vary_effects, a, b, cov_ab, islandsDistMatrix, Mu, N_cafes, N_visits, neff, ratios_cp, rho, R_1, R_2, R_4, sigma, sigma_a, sigma_b, sigmas, post, b13.1, b13.2, b13.3, b13.4, b13.6, b12.5, partially_pooled_params, un_pooled_params, params, partially_pooled_estimates, un_pooled_estimates, estimates, Dmat, K, nd, ftd, b13.7, m13.7, post_b13.7, post_m13.7, prob_to_logit, funding, b13.bonus_0, b13.bonus_1, b13.bonus_2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.