Prediction Interval Cherry Blossoms

cherryblossoms.R

library(ggplot2)
library(rethinking)

data("cherry_blossoms")

d <- cherry_blossoms[complete.cases(cherry_blossoms$temp), ]
d <- d[complete.cases(d$doy), ]

# posterior model 
flist <- alist(
  doy ~ dnorm(mu, sigma),
  mu <- a + b*temp,
  # intercept around day 100
  a ~ dnorm(100, 20),
  b ~ dnorm(0, 10),
  sigma ~ dunif(0, 10)
)

model <- quap(flist, data = d)

# simulate posterior

temp.seq <- seq(from = 2, 10, length.out = 1e4)
mu <- link(model, data = list(temp = temp.seq))

mu.mean <- apply(mu, 2, mean)
mu.PI <- apply(mu, 2, PI, 0.97)

doy.sim <- sim(model, data = list(temp = temp.seq))
doy.PI <- apply(doy.sim, 2, PI, 0.97)

# plot MAP, confidence interval and prediction interval

ggplot() +
  xlim(2, 10) +
  geom_point(aes(x = d$temp, y = d$doy), color = "pink") +
  geom_function(fun = function(x) {precis(model)$mean[1] + precis(model)$mean[2]*x}, color = "red") +
  geom_ribbon(aes(x = temp.seq, ymin = mu.PI[1, ], ymax = mu.PI[2, ]), fill = "grey", alpha = 0.5) +
  geom_ribbon(aes(x = temp.seq, ymin = doy.PI[1, ], ymax = doy.PI[2, ]), fill = "lightblue", alpha = 0.2) +
  theme_minimal() + 
  labs(x = "March temperature",
       y = "Day of first blossom")