Model 7.6

[ASKurz]

In Section 7.1.1, McElreath fit a series of higher-order polynomial models. Most work fine with brms. However, the sixth-order polynomial model 7.6 seems elusive. The catch is it's based on seven data points. McElreath got it to work with rethinking::quap() like this:

# load
library(tidyverse)
library(rethinking)

# define the data
d <- 
  tibble(species = c("afarensis", "africanus", "habilis", "boisei", "rudolfensis", "ergaster", "sapiens"), 
         brain   = c(438, 452, 612, 521, 752, 871, 1350), 
         mass    = c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)) %>% 
  mutate(mass_std  = (mass - mean(mass)) / sd(mass),
         brain_std = brain / max(brain))

# fit the the model
m7.6 <- quap( 
  alist(
    brain_std ~ dnorm(mu, 0.001), 
    mu <- a + b[1] * mass_std + b[2] * mass_std^2 + b[3] * mass_std^3 + b[4] * mass_std^4 + b[5] * mass_std^5 + b[6] * mass_std^6,
    a ~ dnorm(0.5, 1),
    b ~ dnorm(0, 10)
  ), 
  data = d, start = list(b = rep(0, 6)))

Here's what McElreath wrote about that model on page 198: "That last model, m7.6, has one trick in it. The standard deviation is replaced with a constant value 0.001. The model will not work otherwise, for a very important reason that will become clear as we plot these monsters."

If you have the brms solution, I'd love to see it.

[sjwild]

For this one, we need to use a custom family. We can use the Stan's normal_lpdf and normal_rng if we replace sigma with 0.001.

This code should work. It's based off the brms vignette on custom families.

d  <- 
  tibble(species = c("afarensis", "africanus", "habilis", "boisei", "rudolfensis", "ergaster", "sapiens"), 
         brain   = c(438, 452, 612, 521, 752, 871, 1350), 
         mass    = c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)) %>% 
  mutate(mass_std  = (mass - mean(mass)) / sd(mass),
         brain_std = brain / max(brain))

custom_normal <- custom_family(
  "custom_normal", dpars = "mu",
  links = "identity",
  type = "real"
)

stan_funs  <- "real custom_normal_lpdf(real y, real mu) {
  return normal_lpdf(y | mu, 0.001);
}
real custom_normal_rng(real mu) {
  return normal_rng(mu, 0.001);
}
" 

stanvars <- stanvar(scode = stan_funs, block = "functions")

m7.6 <- brm(data = d,
            family = custom_normal,
            formula = brain_std ~ 1 + poly(mass_std, 6),
            prior = c(prior(normal(0.5, 1), class = Intercept),
                      prior(normal(0, 10), class = b)),
            stanvars = stanvars,
            chains = 4,
            cores = 4, 
            iter = 2000,
            warmup = 1000)

expose_functions(m7.6, vectorize = TRUE)

posterior_epred_custom_normal <- function(prep) {
  mu <- prep$dpars$mu
  mu 
}

mass_std <- seq(-1.5, 1.5, length.out = 1000)
fits <- fitted(m7.6, newdata = list(mass_std = mass_std))
fits <- data.frame(fits, 
                   mass_std = mass_std)

ggplot() + geom_line(data = fits, 
                     mapping = aes(x = mass_std, 
                                   y = Estimate),
                     colour = "blue") +
  geom_ribbon(data = fits, mapping = aes(x = mass_std,
                                         ymin = Q2.5,
                                         ymax = Q97.5),
              fill = "blue",
              alpha = 0.5) +
  geom_point(d, mapping = aes(x = mass_std,
                              y = brain_std),
             colour = "blue") +
  theme_minimal()

I hope that helps!

Steve

[ASKurz]

Dude, you're on fire. So I can acknowledge you in the thank you section of the next update, what's your author name?

[sjwild]

My name is Stephen Wild. I'm happy to contribute a little. I found your first edition helpful when I started learning bayes and Stan, so I figured a little help is the least I can do. If all goes well, I may (or may not, most likely) succeed with the ODE models.

Steve

[ASKurz]

Right on. If you can get those ODE's, it'd be a major help.

[ASKurz]

The current implementation lives here.

[sjwild]

That's awesome! The ODE problem is interesting, because we're actually dealing with two issues: ODEs and a state-space model (at least for lynx-hare). This one's going to require some thinking.

[ASKurz]

I'd also be interested in an example or two of a much simpler ODE, something easier for us novices to grapple with.

[sjwild]

I agree. My knowledge of ODEs is limited to 2 things: 1) they exist, and 2) I couldn't solve one if my life depended on it.