In joepowers16/rethinking: Statistical Rethinking book package

8.1. Good King Markov and His island kingdom.

In this version of the code, we've added set.seed(), which helps make the exact results reproducible. If you're new to setting seeds, play around with different values.

set.seed(103)

num_weeks <- 1e5
positions <- rep(0, num_weeks)
current   <- 10
for (i in 1:num_weeks) {
  # record current position
  positions[i] <- current
  # flip coin to generate proposal
  proposal <- current + sample(c(-1, 1), size = 1)
  # now make sure he loops around the archipelago
  if (proposal < 1) proposal <- 10
  if (proposal > 10) proposal <- 1
  # move?
  prob_move <- proposal/current
  current <- ifelse(runif(1) < prob_move, proposal, current)
}

In this document, we'll borrow a theme, theme_ipsum(), from the hrbrthemes package.

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

Figure 8.2.a.

library(tidyverse)

tibble(week = 1:1e5,
       island = positions) %>%
  filter(week < 101) %>%

  ggplot(aes(x = week, y = island)) +
  geom_point(shape = 1) +
  scale_x_continuous(breaks = seq(from = 0, to = 100, by = 20)) +
  scale_y_continuous(breaks = seq(from = 0, to = 10, by = 2)) +
  labs(title = "Behold: The Metropolis algorithm in action!",
       subtitle = "The dots show the king's path over the first 100 weeks.") +
  theme_ipsum()

Figure 8.2.b.

tibble(week = 1:1e5,
       island = positions) %>%
  mutate(island = factor(island)) %>%

  ggplot(aes(x = island)) +
  geom_bar() +
  labs(title = "Old Metropolis shines in the long run.",
       subtitle = "Sure enough, the time the king spent on each island was\nproportional to its population size.") +
  theme_ipsum()

8.3. Easy HMC: ~~map2stan~~ brm()

Here we load the rugged data.

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

Switching from rethinking to brms.

detach(package:rethinking)
library(brms)
rm(rugged)

It takes just a sec to do a little data manipulation.

d <- 
  d %>%
  mutate(log_gdp = log(rgdppc_2000))

dd <-
  d %>%
  filter(complete.cases(rgdppc_2000))

In the context of this chapter, it doesn't make sense to translate McElreath's m8.1 map() code to brm() code. Below, we'll just go directly to the brm() variant of his m8.1stan.

8.3.1. Preparation.

You don't need to do the data processing McElreath does on pages 248 and 249. If you wanted to, however, here's how you might do it within the tidyverse.

dd.trim <-
  dd %>%
  select(log_gdp, rugged, cont_africa)

str(dd.trim)

8.3.2. Estimation.

Finally, we get to work that sweet HMC.

b8.1 <-
  brm(data = dd, family = gaussian,
      log_gdp ~ 1 + rugged + cont_africa + rugged:cont_africa,
      prior = c(set_prior("normal(0, 100)", class = "Intercept"),
                set_prior("normal(0, 10)", class = "b"),
                set_prior("cauchy(0, 2)", class = "sigma")))

w_b8.1 <- waic(b8.1)
l_b8.1 <- loo(b8.1)
l_b8.1_save_psis <- loo(b8.1, save_psis = T)

w_b8.1
l_b8.1

w_b8.1$ pointwise %>% str
l_b8.1$pointwise

w_b8.1$pointwise %>% 
  as_tibble() %>% 
  arrange(p_waic) %>% 
  mutate(i = 1:n()) %>% 

  ggplot(aes(x = i, y = p_waic)) +
  geom_point()

l_b8.1$pointwise %>% 
  as_tibble() %>% 
  arrange(p_loo) %>% 
  mutate(i = 1:n()) %>% 

  ggplot(aes(x = i, y = p_loo)) +
  geom_point()

plot(l_b8.1)

yrep <- posterior_predict(b8.1)

# requires bayesplot version >= 1.5.0
library(bayesplot)

ppc_loo_pit_overlay(
  y = dd$log_gdp,
  yrep = yrep,
  lw = weights(l_b8.1_save_psis$psis_object)
)

The posterior:

print(b8.1)

8.3.3. Sampling again, in parallel.

Here we add chains = 4, cores = 4.

b8.1_4chains_4cores <-
  brm(data = dd, family = gaussian,
      log_gdp ~ 1 + rugged + cont_africa + rugged:cont_africa,
      prior = c(set_prior("normal(0, 100)", class = "Intercept"),
                set_prior("normal(0, 10)", class = "b"),
                set_prior("cauchy(0, 2)", class = "sigma")),
      chains = 4, cores = 4)

This model sampled so fast that it really didn't matter if we sampled in parallel or not. It will for others.

print(b8.1_4chains_4cores)

8.3.4. Visualization.

Unlike the way rethinking's extract.samples() yields a list, brms's posterior_samples() returns a data frame.

post <- posterior_samples(b8.1)
str(post)

As with McElreath's rethinking, brms allows users to put the post data frame or the brmsfit object directly in pairs().

pairs(b8.1)

Another nice way to customize your pairs plot is with the GGally package.

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

post %>%
  select(b_Intercept:sigma) %>%

  ggpairs()

The above code returned an error: "variables: "x" have non standard format: "b_rugged:cont_africa". Please rename the columns or make a new column." GGally didn't like the way brms named one of our parameters, b_rugged:cont_africa. That's easy enough to fix.

my_pairs_plot <-
  post %>%
  rename(b_rugged_x_cont_africa = `b_rugged:cont_africa`) %>%
  select(b_Intercept:sigma) %>%

  ggpairs()

my_pairs_plot

Since GGally returns a ggplot object, you can customize it as you please.

my_pairs_plot +
  labs(subtitle = "My custom pairs plot") +
  theme_ipsum()

For more ideas on customizing a GGally pairs plot, go here.

8.3.5. Using the samples.

If you want information criteria as a part of your brms model summary, just add loo = T and/or waic = T in the summary() function.

summary(b8.1, loo = T, waic = T)

But beware: the information criteria haven’t taken long to estimate in the models we’ve fit thus far in the text. As your models become more complex and as your data get larger, the time their information criteria will take to compute will increase. So, if you have a complex multilevel model with a large data set, you might be better off computing the information criteria separately from the model summary.

8.3.6. Checking the chain.

Using plot() for a brm() fit returns both density and trace lots for the parameters.

plot(b8.1)

The bayesplot package allows a little more control. Here, we use bayesplot's mcmc_trace() to show only trace plots with our custom theme. Note that mcmc_trace() works with data frames, not brmfit objects. There's a further complication. Recall how we made post (i.e., post <- posterior_samples(b8.1)). Our post data frame carries no information on chains. To retain that information, we'll need to add an add_chain = T argument to our posterior_samples() function.

library(bayesplot)

post <- posterior_samples(b8.1, add_chain = T)

mcmc_trace(post[, c(1:5, 7)], # We need to include column 7 because that contains the chain info 
           facet_args = list(ncol = 3), 
           size = .15) +
  labs(title = "My custom trace plots") +
  scale_color_ipsum() +
  theme_ipsum() +
  theme(legend.position = c(.95, .2))

The bayesplot package offers a variety of diagnostic plots. Here we make autocorrelation plots for all model parameters, one for each HMC chain.

mcmc_acf(post, 
         pars = c("b_Intercept", "b_rugged", "b_cont_africa", "b_rugged:cont_africa", "sigma"),
         lags = 5) +
  scale_color_ipsum() +
  theme_ipsum()

That's just what we like to see--nice L-shaped autocorrelation plots. Those are the kinds of shapes you'd expect when you have reasonably large effective samples. Anyway...

Overthinking: Raw Stan model code.

The stancode() function works in brms much like it does in rethinking.

stancode(b8.1)

You can also get that information with b8.1$model or b8.1$fit@stanmodel.

8.4. Care and feeding of your Markov chain.

8.4.3. Taming a wild chain.

As with rethinking, brms can take data in the form of a list. Recall however, that in order to specify starting values, you need to specify a list of lists with an inits argument, rather than with start, as in rethinking.

b8.2 <-
  brm(data = list(y = c(-1, 1)), 
      family = gaussian,
      y ~ 1,
      prior = c(set_prior("uniform(-1e10, 1e10)", class = "Intercept"),
                set_prior("uniform(0, 1e10)", class = "sigma")),
      inits = list(list(Intercept = 0, sigma = 1),
                   list(Intercept = 0, sigma = 1)),
      chains = 2, iter = 4000, warmup = 1000)

Those were some silly flat priors. Here's the damage.

post <- posterior_samples(b8.2, add_chain = T)

mcmc_trace(post[, c(1:2, 4)],
           size = .25) +
  labs(title = "My version of Figure 8.5.a.",
       subtitle = "These trace plots do not look like the fuzzy caterpillars we usually hope for.") +
  scale_color_ipsum() +
  theme_ipsum() +
  theme(legend.position = c(.85, 1.5),
        legend.direction = "horizontal")

Let's peek at the summary.

print(b8.2)

Holy smokes, those parameters are a mess! Plus we got a nasty warning message, too. Watch our reasonable priors save the day.

b8.3 <-
  brm(data = list(y = c(-1, 1)), 
      family = gaussian,
      y ~ 1,
      prior = c(set_prior("normal(0, 10)", class = "Intercept"),
                set_prior("cauchy(0, 1)", class = "sigma")),
      inits = list(list(Intercept = 0, sigma = 1),
                   list(Intercept = 0, sigma = 1)),
      chains = 2, iter = 4000, warmup = 1000)

print(b8.3)

As in the text, no more warning signs and no more silly estimates. The trace plots look great, too.

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

mcmc_trace(post[, c(1:2, 4)],
           size = .25) +
  labs(title = "My version of Figure 8.5.b",
       subtitle  = "Oh man. This looks so much better.") +
  scale_color_ipsum() +
  theme_ipsum() +
  theme(legend.position = c(.85, 1.5),
        legend.direction = "horizontal")

Our version of Figure 8.6.a.

post %>%
  select(b_Intercept) %>%

  ggplot(aes(x = b_Intercept)) +
  stat_density(geom = "line") +
  geom_line(data = data.frame(x = seq(from = min(post$b_Intercept),
                                      to = max(post$b_Intercept),
                                      length.out = 50)),
            aes(x = x, y = dnorm(x = x, mean = 0, sd = 10)),
            color = ipsum_pal()(1), linetype = 2) +
  theme_ipsum()

And our version of Figure 8.6.b.

post %>%
  select(sigma) %>%

  ggplot(aes(x = sigma)) +
  stat_density(geom = "line") +
  geom_line(data = data.frame(x = seq(from = 0,
                                      to = max(post$sigma),
                                      length.out = 50)),
            aes(x = x, y = dcauchy(x = x, location = 0, scale = 1)*2),
            color = ipsum_pal()(2)[2], linetype = 2) +
  coord_cartesian(xlim = c(0, 10)) +
  theme_ipsum()

8.4.4. Non-identifiable parameters.

It appears that the only way to get a brms version of McElreath's m8.4 and m8.5 is to augment the data. In addition to the Gaussian y vector, we'll add two constants to the data, intercept1 = 1 and intercept2 = 1.

set.seed(8.4)
y <- rnorm(100, mean = 0, sd = 1)

b8.4 <-
  brm(data = list(y = y,
                  intercept1 = 1,
                  intercept2 = 1), 
      family = gaussian,
      y ~ 0 + intercept1 + intercept2,
      prior = c(set_prior("uniform(-1e10, 1e10)", class = "b"),
                set_prior("cauchy(0, 1)", class = "sigma")),
      inits = list(list(intercept1 = 0, intercept2 = 0, sigma = 1),
                   list(intercept1 = 0, intercept2 = 0, sigma = 1)),
      chains = 2, iter = 4000, warmup = 1000,
      seed = 8.4)

Our model results don't perfectly mirror McElreath's, but they're identical in spirit.

print(b8.4)

Note the frightening warning message.

b8.5 <-
  brm(data = list(y = y,
                  intercept1 = 1,
                  intercept2 = 1), 
      family = gaussian,
      y ~ 0 + intercept1 + intercept2,
      prior = c(set_prior("normal(0, 10)", class = "b"),
                set_prior("cauchy(0, 1)", class = "sigma")),
      inits = list(list(intercept1 = 0, intercept2 = 0, sigma = 1),
                   list(intercept1 = 0, intercept2 = 0, sigma = 1)),
      chains = 2, iter = 4000, warmup = 1000)

print(b8.5)

Now we'll do the preparatory work for Figure 8.7. Instead of showing the plots, here, we'll save them as objects, left_column and right_column, in order to combine them below.

post <- posterior_samples(b8.4, add_chain = T)

left_column <-
  mcmc_trace(post[, c(1:3, 5)],
           size = .25,
           facet_args = c(ncol = 1)) +
  scale_color_ipsum() +
  theme_ipsum() +
  theme(legend.position = c(.85, 1.5),
        legend.direction = "horizontal")

post <- posterior_samples(b8.5, add_chain = T)

right_column <-
  mcmc_trace(post[, c(1:3, 5)],
           size = .25,
           facet_args = c(ncol = 1)) +
  scale_color_ipsum() +
  theme_ipsum() +
  theme(legend.position = c(.85, 1.5),
        legend.direction = "horizontal")

With a little help of the multiplot() function we can slap our two plot objects together.

Behold the code for the multiplot() function:

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }

  if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

We're finally ready to use multiplot() to make our version of Figure 8.7.

multiplot(left_column, right_column, cols = 2)

The central message in the text, default to weakly-regularizing priors, holds for brms just as it does in rethinking. For more on the topic, see the recommendations from the Stan team. If you want to dive deeper, check out Dan Simpson's post on Gelman's blog and their corresponding paper with Michael Betancourt.

# rm(i, num_weeks, positions, current, proposal, prob_move, d, dd, b8.1, b8.1_4chains_4cores, post, my_pairs_plot, b8.2, b8.3, b8.4, b8.5, multiplot, left_column, right_column)

Note. The analyses in this document were done with:

• R 3.4.4
• RStudio 1.1.442
• rmarkdown 1.9
• rethinking 1.59
• brms 2.1.9
• rstan 2.17.3
• tidyverse 1.2.1
• hrbrthemes 0.1.0
• GGally 1.3.0
• bayesplot 1.4.0

Reference

McElreath, R. (2016). Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman & Hall/CRC Press.

joepowers16/rethinking documentation built on June 2, 2019, 6:52 p.m.