ModelPredictions35.R

A <- function(samples){
  as.matrix(samples[,c("b_Intercept" ,"b_temp")])
}
x <- c(1, 35)
prob <- 0.975 
lin.samples <- posterior_samples(lin.mod)
n <- nrow(lin.samples)
mu <- A(lin.samples) %*% x
sigma <- lin.samples[,"sigma_units"]
(lin.q <- quantile(rnorm(n, mu, sigma), prob))
#    97.5% 
# 1025.077
log.lin.samples <- posterior_samples(log.lin.mod)
mu <- A(log.lin.samples) %*% x
sigma <- log.lin.samples[,"sigma_log_units"]
(log.lin.q <- quantile(exp(rnorm(n, mu +  0.5*sigma^2, sigma)), prob))
#    97.5% 
# 2460.108
pois.samples <- posterior_samples(pois.mod)
mu <- exp(A(pois.samples) %*% x)
(pois.q <- quantile(rpois(n, mu) , prob))
#  97.5% 
#   1515 
bin.samples <- posterior_samples(bin.mod)
inv.logit <- function(u) exp(u)/(1+exp(u))
mu <- inv.logit( A(bin.samples) %*% x)
(bin.q <- quantile(rbinom(n, size = 800, mu), prob))
#   97.5% 
# 761.025
percentiles <- c(lin.q, log.lin.q, pois.q, bin.q)
b <- barplot(percentiles, 
             names.arg = c("Linear", "Log-transformed",
                           "Poisson", "Binomial"),
             ylab="Predicted ice cream units",
             main="Predicted 97.5%ile at 35ºC")
text(b, percentiles-75, round(percentiles))