Splines, section 4.6

[sjwild]

Howdy,

One area you asked for help with was for section 4.6, with splines. The code below should reproduce m4.7 using brms. It's not Tidyverse style, but it works.

The trick is to bind d2$doy and B, as each column of B is a predictor. You can see this in McElreath's code, as he uses B %*% w in quap.

Here's the code:

library(brms)
library(splines)
library(tidyverse)

d <- read.csv("https://raw.githubusercontent.com/rmcelreath/rethinking/master/data/cherry_blossoms.csv", 
              sep = ";")
d2 <- subset(d, !is.na(doy))

num_knots <- 15
knot_list <- quantile(d2$year, probs = seq(from = 0, to = 1, length.out = num_knots))

B <- bs(d2$year,
        knots = knot_list[-c(1, num_knots)], 
        degree = 3, 
        intercept = TRUE)

d3 <- data.frame(doy = d2$doy, 
                 B)

m4.7 <- brm(data = d3,
            family = gaussian,
            formula = doy ~ 1 + .,
            prior = c(prior(normal(100, 10), class = Intercept),
                      prior(normal(0, 10), class = b),
                      prior(exponential(1), class = sigma)),
            cores = 4,
            chains = 4,
            seed = 62531,
            iter = 2000,
            warmup = 500)

summary(m4.7)

fits <- fitted(m4.7)
d2 <- cbind(d2, data.frame(fits))

ggplot(d2) + 
  geom_line(mapping = aes(x = year,
                          y = Estimate)) +
  geom_ribbon(mapping = aes(x = year,
                            ymin = Q2.5,
                            ymax = Q97.5),
              alpha = 0.5) +
  labs(x = "year",
       y = "day in year") +
  geom_hline(yintercept = mean(d2$Estimate), linetype = 2) +
  theme_minimal()

I hope that helps.

Steve

[ASKurz]

Hey @sjwild, I won't have time to look at this, for a while. But thanks! This seems helpful.

[ASKurz]

Steve, your solution is excellent. My current implementation lives here.

[sjwild]

Awesome. I am happy to help with the easy ones. I leave the hard ones to @mages!

The new section looks good. I'm always happy to help.

Steve