Bayesian Normal Regression Model Evaluation

EXERCISE 10.13

The independence assumption means that each observation should be independent of the other observations. However, the 1,339 coffee batches come from only 571 farms, so some farms have multiple batches in the dataset. Batches from the same farm may share similar growing conditions, soil, climate, and processing methods, which can make their ratings more similar. Therefore, observations from the same farm may be related, violating the independence assumption of the linear regression model.

library(bayesrules)
library(tidyverse)
library(rstanarm)
library(bayesplot)
library(tidybayes)
library(broom.mixed)
set.seed(84735)
new_coffee <- coffee_ratings %>% 
  group_by(farm_name) %>% 
  sample_n(1) %>% 
  ungroup()
dim(new_coffee)
## [1] 572  27

EXERCISE 10.14.a Plot and discuss the observed relationship

ggplot(new_coffee, aes(x = aroma, y = total_cup_points)) +
  geom_point(alpha = 0.55) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Coffee Rating by Aroma Grade",
    x = "Aroma grade",
    y = "Total cup points"
  ) +
  theme_minimal()

The relationship

The scatter plot reveals a positive, approximately linear relationship between aroma and total cup points. While there is natural variation in ratings at any specific aroma grade, coffees with higher aroma scores generally receive higher overall ratings, making a simple Normal regression model a reasonable first approach.

10.14.b Simulate the Normal regression posterior model

set.seed(84735)

aroma_model <- stan_glm(
  total_cup_points ~ aroma,
  data = new_coffee,
  family = gaussian,
  prior_intercept = normal(75, 10),
  prior = normal(0, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4,
  iter = 10000,
  seed = 84735
)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 6.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.191 seconds (Warm-up)
## Chain 1:                0.397 seconds (Sampling)
## Chain 1:                0.588 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.2 seconds (Warm-up)
## Chain 2:                0.389 seconds (Sampling)
## Chain 2:                0.589 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.187 seconds (Warm-up)
## Chain 3:                0.525 seconds (Sampling)
## Chain 3:                0.712 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.8e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.196 seconds (Warm-up)
## Chain 4:                0.438 seconds (Sampling)
## Chain 4:                0.634 seconds (Total)
## Chain 4:
# Prior information
prior_summary(aroma_model)
## Priors for model 'aroma_model' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 75, scale = 10)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 22)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.36)
## ------
## See help('prior_summary.stanreg') for more details
# Posterior numerical results
print(aroma_model)
## stan_glm
##  family:       gaussian [identity]
##  formula:      total_cup_points ~ aroma
##  observations: 572
##  predictors:   2
## ------
##             Median MAD_SD
## (Intercept) 35.4    2.0  
## aroma        6.2    0.3  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 2.0    0.1   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg
# Detailed results including Rhat and effective sample size
summary(aroma_model)
## 
## Model Info:
##  function:     stan_glm
##  family:       gaussian [identity]
##  formula:      total_cup_points ~ aroma
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 572
##  predictors:   2
## 
## Estimates:
##               mean   sd   10%   50%   90%
## (Intercept) 35.4    2.0 32.8  35.4  38.0 
## aroma        6.2    0.3  5.8   6.2   6.5 
## sigma        2.0    0.1  1.9   2.0   2.0 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 82.1    0.1 82.0  82.1  82.2 
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##               mcse Rhat n_eff
## (Intercept)   0.0  1.0  19018
## aroma         0.0  1.0  19162
## sigma         0.0  1.0  18509
## mean_PPD      0.0  1.0  20023
## log-posterior 0.0  1.0   8803
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
# Trace plots showing all 4 MCMC chains
plot(
  aroma_model,
  plotfun = "trace",
  pars = c("(Intercept)", "aroma", "sigma")
)

# Posterior regression lines
new_coffee %>%
  add_fitted_draws(aroma_model, n = 100) %>%
  ggplot(aes(x = aroma, y = total_cup_points)) +
  geom_point(alpha = 0.4) +
  geom_line(
    aes(y = .value, group = .draw),
    alpha = 0.15
  ) +
  labs(
    title = "Posterior Regression Model: Coffee Rating by Aroma",
    x = "Aroma grade",
    y = "Total cup points"
  ) +
  theme_minimal()

Interpretation of results

The model shows a clear positive relationship between aroma and coffee rating. The posterior median for the aroma coefficient is 6.2, with an 80% posterior interval from 5.8 to 6.5. This means that a one-point increase in aroma is associated with about a 6.2-point increase in total cup points, on average. Since the interval is entirely above zero, there is strong posterior evidence of a positive association. The model also converged well because all Rhat values are 1.00 and the effective sample sizes are high. The posterior regression lines show the same positive pattern, supporting the numerical results.

10.14.c Visual and numerical posterior summaries for \(\beta_1\)

aroma_draws <- as.data.frame(aroma_model)

ggplot(aroma_draws, aes(x = aroma)) +
  geom_density(alpha = 0.45) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(
    title = "Posterior Distribution of the Aroma Coefficient",
    x = expression(beta[1]),
    y = "Density"
  ) +
  theme_minimal()

aroma_summary <- tidy(
  aroma_model,
  effects = "fixed",
  conf.int = TRUE,
  conf.level = 0.95
) %>%
  filter(term == "aroma")

# Define the missing variables so the inline text works
aroma_median <- aroma_summary$estimate
aroma_lower <- aroma_summary$conf.low
aroma_upper <- aroma_summary$conf.high

aroma_summary
## # A tibble: 1 × 5
##   term  estimate std.error conf.low conf.high
##   <chr>    <dbl>     <dbl>    <dbl>     <dbl>
## 1 aroma     6.16     0.261     5.65      6.68
aroma_prob_positive <- mean(aroma_draws$aroma > 0)
aroma_prob_positive
## [1] 1

EXPLANATION:

The posterior median of the aroma coefficient is 6.164. This means that, according to the fitted model, a one-point increase in aroma grade is associated with an estimated increase of about 6.164 points in the typical total coffee rating.

10.14.e Posterior evidence of a positive association

The 95% posterior credible interval for the aroma coefficient is approximately (5.647, 6.677). In addition, the posterior probability that \(\beta_1>0\) is 1. If the credible interval is entirely above 0, this provides strong posterior evidence that coffees with better aroma grades tend to have higher total ratings. This describes an association and does not by itself establish causation.

EXERCISE 10.15: Is the Aroma Model Wrong?

10.15.a Simulate 572 ratings using the first posterior parameter set

The first posterior draw gives one plausible parameter set \((\beta_0,\beta_1,\sigma)\). I use that set together with the observed aroma grades to simulate one complete dataset.

first_parameters <- aroma_draws[1, ]

b0_first <- first_parameters[["(Intercept)"]]
b1_first <- first_parameters[["aroma"]]
sigma_first <- first_parameters[["sigma"]]

set.seed(84735)

simulated_first <- new_coffee %>%
  mutate(
    mu_first = b0_first + b1_first * aroma,
    simulated_rating = rnorm(
      n = n(),
      mean = mu_first,
      sd = sigma_first
    )
  )

head(simulated_first)
## # A tibble: 6 × 29
##   owner farm_name mill  in_country_partner country_of_origin altitude_low_meters
##   <fct> <fct>     <fct> <fct>              <fct>                           <dbl>
## 1 tayl… -         <NA>  Kenya Coffee Trad… Kenya                            1650
## 2 will… 1         1     Specialty Coffee … Indonesia                        3500
## 3 myri… 200 farms coeb… Specialty Coffee … Haiti                             350
## 4 myri… 2000 far… coop… Specialty Coffee … Haiti                             640
## 5 myri… 2000 far… utca… Specialty Coffee … Haiti                            1000
## 6 cqi … a shu sh… a sh… Blossom Valley In… Taiwan                            150
## # ℹ 23 more variables: altitude_high_meters <dbl>, altitude_mean_meters <dbl>,
## #   number_of_bags <dbl>, bag_weight <fct>, species <fct>, variety <fct>,
## #   processing_method <fct>, aroma <dbl>, flavor <dbl>, aftertaste <dbl>,
## #   acidity <dbl>, body <dbl>, balance <dbl>, uniformity <dbl>,
## #   clean_cup <dbl>, sweetness <dbl>, moisture <dbl>,
## #   category_one_defects <dbl>, category_two_defects <dbl>, color <fct>,
## #   total_cup_points <dbl>, mu_first <dbl>, simulated_rating <dbl>
nrow(simulated_first)
## [1] 572

10.15.b Compare the simulated and observed densities

density_data <- bind_rows(
  tibble(
    rating = new_coffee$total_cup_points,
    source = "Observed ratings"
  ),
  tibble(
    rating = simulated_first$simulated_rating,
    source = "Simulated ratings"
  )
)

ggplot(density_data, aes(x = rating, linetype = source)) +
  geom_density(linewidth = 1) +
  labs(
    title = "Observed and Simulated Coffee Ratings",
    x = "Total cup points",
    y = "Density",
    linetype = "Sample"
  ) +
  theme_minimal()

density_summary <- density_data %>%
  group_by(source) %>%
  summarise(
    mean = mean(rating),
    median = median(rating),
    sd = sd(rating),
    .groups = "drop"
  )

density_summary
## # A tibble: 2 × 4
##   source             mean median    sd
##   <chr>             <dbl>  <dbl> <dbl>
## 1 Observed ratings   82.1   82.3  2.75
## 2 Simulated ratings  82.2   82.1  2.78

Discusion

The simulated ratings should resemble the observed ratings if the model is capturing the important features of the data. The two density curves show how closely one posterior-plausible dataset matches the real data. Their central locations should be similar, but differences in the shape or tails show that the Normal regression model is only an approximation. A single simulated dataset is limited, so the next posterior predictive check uses many posterior draws.

10.15.c More complete posterior predictive check

set.seed(84735)

pp_check(aroma_model) +
  labs(
    title = "Posterior Predictive Check: Aroma Model",
    x = "Total cup points"
  )

In this plot, the observed outcome distribution is compared with distributions simulated from the fitted Bayesian model. Substantial overlap indicates that the model reproduces the main behavior of the observed ratings. Any repeated differences in the tails or overall shape indicate places where the model is imperfect.

10.15.d Evaluate assumptions 2 and 3

Assumption 2 states that the typical value of \(Y\) changes linearly with \(X\). The scatterplot from Exercise 10.14 shows a clear positive and roughly linear trend between aroma and total cup points, so this assumption appears reasonable as a first approximation.

Assumption 3 states that ratings at a given predictor value vary approximately Normally around the regression mean with a common residual standard deviation. The posterior predictive check shows whether Normal datasets generated by the fitted model resemble the observed ratings. If the simulated curves generally overlap the observed curve, this assumption is reasonably supported for the main part of the data. Any mismatch in the extreme ratings shows that the assumption is not perfect. Overall, I would treat assumptions 2 and 3 as useful approximations rather than exact descriptions of the coffee ratings.

EXERCISE 10.16: Posterior Predictive Accuracy of the Aroma Model

10.16.a Posterior predictive model for a batch with aroma = 7.67

This part is simulated manually and does not use posterior_predict().

x_new <- 7.67

set.seed(84735)

first_batch_predictive <- aroma_draws %>%
  mutate(
    mu_new = `(Intercept)` + aroma * x_new,
    y_new = rnorm(
      n = n(),
      mean = mu_new,
      sd = sigma
    )
  )

# Posterior predictive distribution
ggplot(first_batch_predictive, aes(x = y_new)) +
  geom_density(alpha = 0.45) +
  labs(
    title = "Posterior Predictive Distribution",
    subtitle = "For a coffee batch with aroma grade = 7.67",
    x = "Predicted total cup points",
    y = "Density"
  ) +
  theme_minimal()

# Show where aroma = 7.67 is in the observed data
ggplot(new_coffee, aes(x = aroma, y = total_cup_points)) +
  geom_point(alpha = 0.5) +
  geom_vline(
    xintercept = 7.67,
    linetype = "dashed",
    linewidth = 1
  ) +
  annotate(
    "text",
    x = 7.67,
    y = max(new_coffee$total_cup_points),
    label = "Aroma = 7.67",
    hjust = -0.1
  ) +
  labs(
    title = "Observed Coffee Ratings and Aroma = 7.67",
    x = "Aroma grade",
    y = "Total cup points"
  ) +
  theme_minimal()

# Numerical posterior predictive summary
first_batch_predictive %>%
  summarise(
    predictive_mean = mean(y_new),
    predictive_sd = sd(y_new),
    predictive_median = median(y_new),
    lower_50 = quantile(y_new, 0.25),
    upper_50 = quantile(y_new, 0.75),
    lower_95 = quantile(y_new, 0.025),
    upper_95 = quantile(y_new, 0.975)
  )
##   predictive_mean predictive_sd predictive_median lower_50 upper_50 lower_95
## 1        82.66198      1.954686          82.64533 81.33092 83.97605 78.89919
##   upper_95
## 1 86.51351

The density represents plausible ratings for a new coffee batch with an aroma grade of 7.67 while accounting for both posterior uncertainty in the regression parameters and the natural residual variation among coffee ratings.

10.16.b Raw and standardized posterior predictive error

actual_rating <- 84

predictive_mean <- mean(first_batch_predictive$y_new)
predictive_sd <- sd(first_batch_predictive$y_new)

raw_error <- actual_rating - predictive_mean
standardized_error <- raw_error / predictive_sd

tibble(
  actual_rating = actual_rating,
  predictive_mean = predictive_mean,
  predictive_sd = predictive_sd,
  raw_error = raw_error,
  standardized_error = standardized_error
)
## # A tibble: 1 × 5
##   actual_rating predictive_mean predictive_sd raw_error standardized_error
##           <dbl>           <dbl>         <dbl>     <dbl>              <dbl>
## 1            84            82.7          1.95      1.34              0.685
raw_direction <- ifelse(raw_error >= 0, "above", "below")
raw_direction
## [1] "above"

EXPLANATION

The raw error is 1.338 points. Therefore, the observed rating of 84 is about 1.338 points above the model’s posterior predictive mean. The standardized error is 0.685, meaning the observed rating is about 0.685 posterior predictive standard deviations above the mean prediction. A standardized error with a large absolute value would indicate an unusually poor prediction for this case.

10.16.c ppc_intervals() for all batches

For readability, I order the batches by aroma grade before plotting. The model and observed outcomes remain paired correctly.

new_coffee_ordered <- new_coffee %>%
  arrange(aroma)

set.seed(84735)

all_predictions_aroma <- posterior_predict(
  aroma_model,
  newdata = new_coffee_ordered
)

ppc_intervals(
  y = new_coffee_ordered$total_cup_points,
  yrep = all_predictions_aroma,
  prob = 0.50,
  prob_outer = 0.95
) +
  labs(
    title = "Posterior Predictive Intervals: Aroma Model",
    x = "Coffee batches ordered by aroma",
    y = "Total cup points"
  )

DISCUSSION

The plot compares each observed rating with its posterior predictive intervals. Observed points that fall inside the intervals are consistent with values anticipated by the model, while points outside the 95% intervals represent cases that the model predicts poorly. The 50% intervals are narrower, so fewer observed values are expected to fall inside them. For this specific dataset, ordering the batches by aroma visually confirms the positive linear trend between aroma scores and overall ratings. Overall, the model successfully captures the primary behavior of the data, but several severe outliers fall far below the lower prediction bounds. This indicates the model underperforms in the lower tail, reinforcing that the Normal model assumptions are useful approximations rather than a perfect fit.

10.16.d Number of ratings within the 50% posterior prediction interval

interval_50 <- apply(
  all_predictions_aroma,
  2,
  quantile,
  probs = c(0.25, 0.75)
)

inside_50 <- (
  new_coffee_ordered$total_cup_points >= interval_50[1, ]
) & (
  new_coffee_ordered$total_cup_points <= interval_50[2, ]
)

inside_50_count <- sum(inside_50)
inside_50_proportion <- mean(inside_50)

tibble(
  number_inside_50 = inside_50_count,
  proportion_inside_50 = inside_50_proportion
)
## # A tibble: 1 × 2
##   number_inside_50 proportion_inside_50
##              <int>                <dbl>
## 1              394                0.689

There are 394 batches whose observed ratings fall within their 50% posterior prediction intervals. This is 68.9% of the batches.

EXERCISE 10.19: Model Coffee Ratings by Aftertaste

For this model, I keep the same prior structure used for the aroma model so that the two predictors can be compared under the same prior assumptions.

10.19.a Fit the aftertaste posterior model

set.seed(84735)

aftertaste_model <- stan_glm(
  total_cup_points ~ aftertaste,
  data = new_coffee,
  family = gaussian,
  prior_intercept = normal(75, 10),
  prior = normal(0, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4,
  iter = 10000,
  seed = 84735
)
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 3.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 1: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 1: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 1: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 1: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 1: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 1: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 1: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 1: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 1: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 1: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 1: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.43 seconds (Warm-up)
## Chain 1:                0.799 seconds (Sampling)
## Chain 1:                1.229 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.8e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 2: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 2: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 2: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 2: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 2: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 2: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 2: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 2: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 2: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 2: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 2: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.311 seconds (Warm-up)
## Chain 2:                0.628 seconds (Sampling)
## Chain 2:                0.939 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 2.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 3: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 3: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 3: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 3: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 3: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 3: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 3: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 3: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 3: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 3: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 3: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.3 seconds (Warm-up)
## Chain 3:                0.495 seconds (Sampling)
## Chain 3:                0.795 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.19 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 10000 [  0%]  (Warmup)
## Chain 4: Iteration: 1000 / 10000 [ 10%]  (Warmup)
## Chain 4: Iteration: 2000 / 10000 [ 20%]  (Warmup)
## Chain 4: Iteration: 3000 / 10000 [ 30%]  (Warmup)
## Chain 4: Iteration: 4000 / 10000 [ 40%]  (Warmup)
## Chain 4: Iteration: 5000 / 10000 [ 50%]  (Warmup)
## Chain 4: Iteration: 5001 / 10000 [ 50%]  (Sampling)
## Chain 4: Iteration: 6000 / 10000 [ 60%]  (Sampling)
## Chain 4: Iteration: 7000 / 10000 [ 70%]  (Sampling)
## Chain 4: Iteration: 8000 / 10000 [ 80%]  (Sampling)
## Chain 4: Iteration: 9000 / 10000 [ 90%]  (Sampling)
## Chain 4: Iteration: 10000 / 10000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.258 seconds (Warm-up)
## Chain 4:                0.547 seconds (Sampling)
## Chain 4:                0.805 seconds (Total)
## Chain 4:
prior_summary(aftertaste_model)
## Priors for model 'aftertaste_model' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 75, scale = 10)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 20)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.36)
## ------
## See help('prior_summary.stanreg') for more details
print(aftertaste_model)
## stan_glm
##  family:       gaussian [identity]
##  formula:      total_cup_points ~ aftertaste
##  observations: 572
##  predictors:   2
## ------
##             Median MAD_SD
## (Intercept) 33.1    1.4  
## aftertaste   6.6    0.2  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 1.6    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

10.19.b Quick posterior predictive check

set.seed(84735)

pp_check(aftertaste_model) +
  labs(
    title = "Posterior Predictive Check: Aftertaste Model",
    x = "Total cup points"
  )

This plot compares the observed distribution of coffee ratings with distributions generated from the aftertaste model. If the observed density follows the same general shape as the simulated densities, the model captures the main features of the outcome distribution. Consistent differences indicate where the model is wrong or incomplete. In this case, the model is slightly wrong: the observed data has a much sharper, higher peak than the model’s simulations, indicating the model expects a wider variance in middle-tier scores than actually exists. The model also smooths over several distinct outliers in the lower tail.

10.19.c 10-fold cross-validated posterior predictive quality

set.seed(84735)

aftertaste_cv <- prediction_summary_cv(
  model = aftertaste_model,
  data = new_coffee,
  k = 10
)

# Results for each of the 10 folds
aftertaste_cv$folds
##    fold       mae mae_scaled within_50 within_95
## 1     1 0.5335567  0.3308355 0.7586207 0.9827586
## 2     2 0.7035429  0.4500871 0.6842105 0.9473684
## 3     3 0.7806366  0.4803939 0.7719298 1.0000000
## 4     4 0.4117100  0.2557490 0.8596491 0.9649123
## 5     5 0.5425611  0.3899349 0.7192982 0.9473684
## 6     6 0.6962022  0.4365138 0.5964912 0.9649123
## 7     7 0.7929982  0.5064621 0.6842105 0.9298246
## 8     8 0.6723484  0.4289112 0.7368421 0.9649123
## 9     9 0.9493724  0.5928659 0.6315789 0.9824561
## 10   10 0.8169294  0.5164621 0.7068966 0.9655172
# Overall cross-validation results
aftertaste_cv$cv
##         mae mae_scaled within_50 within_95
## 1 0.6899858  0.4388215 0.7149728  0.965003

mae: the mean absolute error, or the average absolute difference between the observed ratings and the model’s predicted ratings. Smaller values mean better predictive accuracy.

mae_scaled: the mean absolute error on a standardized scale, relative to the model’s predictive uncertainty. Smaller values indicate better predictive performance.

within_50: the proportion of observed ratings that fall inside the 50% posterior prediction intervals. A well-calibrated model should be reasonably close to 0.50.

within_95: the proportion of observed ratings that fall inside the 95% posterior prediction intervals. A well-calibrated model should be reasonably close to 0.95.

The overall cross-validation results show an MAE of 0.690, meaning that the model’s predictions differ from the actual coffee ratings by about 0.69 points on average. The within_50 value of 0.715 indicates that about 71.5% of observed ratings were within the 50% prediction intervals, while the within_95 value of 0.965 indicates that about 96.5% were within the 95% prediction intervals.

10.19.d Compare aroma and aftertaste

Because Exercise 10.19 asks for a final comparison, I also calculate the same 10-fold cross-validated summaries for the previously fitted aroma model. This allows the two predictors to be compared using the same type of out-of-sample assessment.

set.seed(84735)

aroma_cv <- prediction_summary_cv(
  model = aroma_model,
  data = new_coffee,
  k = 10
)

cv_comparison <- bind_rows(
  Aroma = aroma_cv$cv,
  Aftertaste = aftertaste_cv$cv,
  .id = "predictor"
)

cv_comparison
##    predictor       mae mae_scaled within_50 within_95
## 1      Aroma 0.8817667  0.4503440 0.6854507 0.9598004
## 2 Aftertaste 0.6899858  0.4388215 0.7149728 0.9650030
preferred_predictor <- cv_comparison %>%
  slice_min(order_by = mae, n = 1, with_ties = FALSE) %>%
  pull(predictor)

preferred_mae <- cv_comparison %>%
  filter(predictor == preferred_predictor) %>%
  pull(mae)

other_predictor <- cv_comparison %>%
  filter(predictor != preferred_predictor) %>%
  pull(predictor)

other_mae <- cv_comparison %>%
  filter(predictor != preferred_predictor) %>%
  pull(mae)

If I could only pick one predictor of coffee bean ratings, I would choose aftertaste. The cross-validation results show that the aftertaste model produces a lower mean absolute error (0.690) compared to the aroma model (0.882), indicating that aftertaste yields more consistently accurate predictions of the actual ratings.

AI Usage

I used Gemini AI tool as a support resource while completing this assignment. I used it to help me understand some of the R code, organize and format parts of my work, and check whether my answers and interpretations were on the right track. I also used it to help me understand errors I came across while running my code and to suggest possible corrections. After receiving the feedback, I went back to my own work, made the necessary changes, reran the code, and reviewed the results. I did not simply copy the output; I used the explanations to improve my understanding and then modified my work before writing the final explanations.