Assignment 5: Evaluating Bayesian Regression Models
Week 5 — Coffee Ratings
Gerald Hlabiso
1 Introduction
This assignment evaluates Bayesian Normal regression models of coffee
ratings. The response variable is total_cup_points, and the
predictors considered separately are aroma and
aftertaste. The analysis covers Exercises 10.13, 10.14,
10.15, 10.16, and 10.19. Exercise 10.18 is omitted as instructed.
data("coffee_ratings")
coffee_ratings <- coffee_ratings %>%
select(farm_name, total_cup_points, aroma, aftertaste)
head(coffee_ratings)## # A tibble: 6 × 4
## farm_name total_cup_points aroma aftertaste
## <fct> <dbl> <dbl> <dbl>
## 1 "metad plc" 90.6 8.67 8.67
## 2 "metad plc" 89.9 8.75 8.5
## 3 "san marcos barrancas \"san cristobal cuch" 89.8 8.42 8.42
## 4 "yidnekachew dabessa coffee plantation" 89 8.17 8.42
## 5 "metad plc" 88.8 8.25 8.25
## 6 <NA> 88.8 8.58 8.42dim(coffee_ratings)## [1] 1339 42 Exercise 10.13: Getting Started with Coffee Ratings
2.1 10.13.a Independence assumption
The original data contain multiple batches from some farms. Batches from the same farm may share soil, climate, altitude, processing practices, management, and other farm-level characteristics. Their ratings may therefore be more similar than ratings from unrelated farms. Treating all 1,339 batches as independent observations would likely violate the conditional independence assumption of the Bayesian linear regression model because observations are clustered within farms.
The following code checks the number of distinct nonmissing farm names and the number of records without a recorded farm name.
tibble(
distinct_named_farms = n_distinct(coffee_ratings$farm_name, na.rm = TRUE),
missing_farm_names = sum(is.na(coffee_ratings$farm_name))
) %>%
kable(caption = "Farm-name checks")| distinct_named_farms | missing_farm_names |
|---|---|
| 571 | 359 |
There are 571 distinct nonmissing farm names. Missing farm names are
placed into an additional NA group by
group_by(), which explains why the sampling procedure below
returns 572 rows.
2.2 10.13.b Select one batch per farm
set.seed(84735)
new_coffee <- coffee_ratings %>%
group_by(farm_name) %>%
sample_n(1) %>%
ungroup()
dim(new_coffee)## [1] 572 4head(new_coffee)## # A tibble: 6 × 4
## farm_name total_cup_points aroma aftertaste
## <fct> <dbl> <dbl> <dbl>
## 1 - 84 7.67 7.67
## 2 1 76.2 7.33 6.5
## 3 200 farms 67.9 6.75 6.42
## 4 2000 farmers 72.3 6.92 7.08
## 5 2000 farms 80.8 7.42 7.42
## 6 a shu she coffee 阿束社咖啡莊園 80.1 7.25 7.25Selecting one observation from each farm group removes repeated observations within farms from this simplified analysis. This makes the independence assumption more plausible, although it does not guarantee complete independence or remove every possible source of bias.
3 Exercise 10.14: Model Rating with Aroma
3.1 10.14.a Plot and discuss the relationship
ggplot(new_coffee, aes(x = aroma, y = total_cup_points)) +
geom_point(alpha = 0.55, color = "#1f4e79") +
geom_smooth(method = "lm", se = FALSE, color = "#c0392b") +
labs(
title = "Coffee Rating by Aroma Grade",
x = "Aroma grade",
y = "Total cup points"
) +
theme_minimal(base_size = 12)
The scatterplot shows a clear positive, approximately linear association. Coffees with higher aroma grades generally receive higher total cup ratings. The vertical variation around the fitted line appears reasonably consistent across most of the observed aroma range, although the plot alone cannot establish that every model assumption is exactly satisfied.
3.2 10.14.b Fit the Bayesian Normal regression model
Let (Y_i) be the total cup rating and (X_i) be the aroma grade. The model is
\[ Y_i\mid\beta_0,\beta_1,\sigma \stackrel{ind}{\sim} N(\beta_0+\beta_1X_i,\sigma^2). \]
The prior for the centered intercept is (N(75,10^2)). This places approximately 95% of its prior probability between 55 and 95 points. A weakly informative (N(0,10^2)) prior is used for the aroma coefficient.
set.seed(84735)
aroma_model <- stan_glm(
total_cup_points ~ aroma,
data = new_coffee,
family = gaussian(),
prior = normal(0, 10, autoscale = FALSE),
prior_intercept = normal(75, 10, autoscale = FALSE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4,
iter = 5000,
seed = 84735,
refresh = 0
)
print(aroma_model, digits = 3)## stan_glm
## family: gaussian [identity]
## formula: total_cup_points ~ aroma
## observations: 572
## predictors: 2
## ------
## Median MAD_SD
## (Intercept) 35.426 1.986
## aroma 6.158 0.262
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 1.959 0.058
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg3.3 10.14.c Visual and numerical posterior summaries for beta1
aroma_draws <- as.data.frame(aroma_model)
aroma_interval <- posterior_interval(
aroma_model,
pars = "aroma",
prob = 0.95
)
aroma_summary <- tibble(
parameter = "Aroma coefficient (beta1)",
posterior_median = median(aroma_draws$aroma),
posterior_mean = mean(aroma_draws$aroma),
posterior_sd = sd(aroma_draws$aroma),
lower_95 = aroma_interval[1, 1],
upper_95 = aroma_interval[1, 2]
)
kable(aroma_summary, digits = 3,
caption = "Posterior summary of the aroma coefficient")| parameter | posterior_median | posterior_mean | posterior_sd | lower_95 | upper_95 |
|---|---|---|---|---|---|
| Aroma coefficient (beta1) | 6.158 | 6.158 | 0.264 | 5.64 | 6.677 |
ggplot(aroma_draws, aes(x = aroma)) +
geom_density(fill = "#4c78a8", alpha = 0.45) +
geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
geom_vline(
xintercept = median(aroma_draws$aroma),
color = "#c0392b"
) +
labs(
title = "Posterior Distribution of the Aroma Coefficient",
subtitle = "Dashed line: no association; red line: posterior median",
x = expression(beta[1]),
y = "Posterior density"
) +
theme_minimal(base_size = 12)
3.4 10.14.d Interpret the posterior median
The posterior median of the aroma coefficient is 6.16. Thus, a one-point increase in aroma grade is associated with an expected increase of approximately 6.16 points in total cup rating, according to this model.
3.5 10.14.e Evidence of a positive association
The 95% posterior credible interval for the aroma coefficient is (5.64, 6.68). Because this entire interval is above zero, the model provides strong posterior evidence that coffees with better aroma grades tend to receive higher total cup ratings. This is an association and should not be interpreted automatically as a causal effect.
4 Exercise 10.15: Is the Aroma Model Wrong?
4.1 10.15.a Simulate ratings using the first posterior parameter set
first_set <- aroma_draws %>%
slice(1)
first_set %>%
select(`(Intercept)`, aroma, sigma) %>%
kable(digits = 3, caption = "First posterior parameter set")| (Intercept) | aroma | sigma |
|---|---|---|
| 32.882 | 6.502 | 1.96 |
set.seed(84735)
one_simulation <- new_coffee %>%
mutate(
mu = first_set$`(Intercept)` + first_set$aroma * aroma,
simulated_rating = rnorm(
n = n(),
mean = mu,
sd = first_set$sigma
)
)
head(one_simulation)## # A tibble: 6 × 6
## farm_name total_cup_points aroma aftertaste mu simulated_rating
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 - 84 7.67 7.67 82.8 84.1
## 2 1 76.2 7.33 6.5 80.5 80.3
## 3 200 farms 67.9 6.75 6.42 76.8 78.3
## 4 2000 farmers 72.3 6.92 7.08 77.9 78.3
## 5 2000 farms 80.8 7.42 7.42 81.1 82.7
## 6 a shu she coffee 阿束社… 80.1 7.25 7.25 80.0 83.1The code uses the first posterior values of (eta_0), (eta_1), and
(sigma) with all 572 observed aroma grades. It produces one new rating
for every batch in new_coffee.
4.2 10.15.b Compare simulated and observed densities
ggplot(one_simulation) +
geom_density(
aes(x = total_cup_points, color = "Observed"),
linewidth = 1
) +
geom_density(
aes(x = simulated_rating, color = "Simulated"),
linewidth = 1
) +
scale_color_manual(
values = c("Observed" = "#17365d", "Simulated" = "#e67e22")
) +
labs(
title = "Observed and Simulated Coffee Ratings",
x = "Total cup points",
y = "Density",
color = "Sample"
) +
theme_minimal(base_size = 12)
The curves should overlap substantially in their center and overall spread. Any remaining difference in symmetry, peaks, or tails indicates features that a simple Normal regression model does not reproduce perfectly. A single simulated dataset is only one plausible realization, so a broader posterior predictive check is also needed.
4.3 10.15.c Complete posterior predictive check
set.seed(84735)
pp_check(aroma_model, plotfun = "dens_overlay", nreps = 50) +
labs(
title = "Posterior Predictive Check for the Aroma Model",
x = "Total cup points"
)
The observed density should lie largely within the collection of replicated densities. A close match in center and spread indicates that the model can reproduce the broad distribution of ratings. Small differences in shape show that the model remains an approximation.
4.4 10.15.d Assess assumptions 2 and 3
Assumption 2 states that the typical rating is a linear function of aroma. Assumption 3 states that ratings vary Normally around the regression mean with consistent variability across aroma values. The approximately straight pattern in the scatterplot and the broad agreement in the posterior predictive check make these assumptions reasonable approximations. They are not perfectly satisfied because the observed density may contain minor shape features that the simulated Normal densities do not fully reproduce.
5 Exercise 10.16: Posterior Prediction Accuracy
5.1 10.16.a Predict the first batch without posterior_predict()
The first batch in new_coffee has an aroma grade of
7.67. For each posterior parameter set, the code below calculates the
expected rating and then simulates one possible rating.
first_aroma <- new_coffee$aroma[1]
observed_first_rating <- new_coffee$total_cup_points[1]
set.seed(84735)
first_batch_predictions <- aroma_draws %>%
transmute(
mu = `(Intercept)` + aroma * first_aroma,
predicted_rating = rnorm(n(), mean = mu, sd = sigma)
)
ggplot(first_batch_predictions, aes(x = predicted_rating)) +
geom_density(fill = "#4c78a8", alpha = 0.45) +
geom_vline(
xintercept = observed_first_rating,
color = "#c0392b",
linetype = "dashed",
linewidth = 1
) +
labs(
title = "Posterior Predictive Distribution for the First Batch",
subtitle = paste("Observed rating =", observed_first_rating),
x = "Predicted total cup points",
y = "Density"
) +
theme_minimal(base_size = 12)
5.2 10.16.b Raw and standardized prediction errors
Following the chapter, the raw error is the observed rating minus the posterior predictive mean. The standardized error divides the raw error by the posterior predictive standard deviation.
prediction_mean <- mean(first_batch_predictions$predicted_rating)
prediction_sd <- sd(first_batch_predictions$predicted_rating)
raw_error <- observed_first_rating - prediction_mean
standardized_error <- raw_error / prediction_sd
first_batch_error_summary <- tibble(
observed_rating = observed_first_rating,
posterior_predictive_mean = prediction_mean,
posterior_predictive_sd = prediction_sd,
raw_error = raw_error,
standardized_error = standardized_error
)
kable(
first_batch_error_summary,
digits = 3,
caption = "Posterior predictive errors for the first batch"
)| observed_rating | posterior_predictive_mean | posterior_predictive_sd | raw_error | standardized_error |
|---|---|---|---|---|
| 84 | 82.68 | 1.964 | 1.32 | 0.672 |
The raw error is 1.32 rating points. Its sign indicates whether the model underpredicted or overpredicted the observed rating: a positive value means underprediction, whereas a negative value means overprediction. The standardized error is 0.67, meaning the observed rating is 0.67 posterior predictive standard deviations from the predictive mean. A value close to zero indicates an accurate prediction relative to the model’s uncertainty.
5.3 10.16.c Posterior prediction intervals for all batches
set.seed(84735)
all_predictions <- posterior_predict(aroma_model)
ppc_intervals(
y = new_coffee$total_cup_points,
yrep = all_predictions,
x = new_coffee$aroma,
prob = 0.50,
prob_outer = 0.90
) +
labs(
title = "Posterior Prediction Intervals for All Coffee Batches",
x = "Aroma grade",
y = "Total cup points"
) +
theme_minimal(base_size = 12)
The inner segments show 50% posterior prediction intervals, and the outer segments show 90% intervals. Most observed ratings should fall inside the wider 90% intervals. Observations outside an interval are not automatically model failures; even a well-calibrated 90% interval is expected to miss about 10% of outcomes over repeated cases.
5.4 10.16.d Count observations within the 50% intervals
interval_50 <- apply(
all_predictions,
2,
quantile,
probs = c(0.25, 0.75)
)
inside_50 <-
new_coffee$total_cup_points >= interval_50[1, ] &
new_coffee$total_cup_points <= interval_50[2, ]
coverage_50 <- tibble(
batches_inside_50 = sum(inside_50),
total_batches = length(inside_50),
proportion_inside_50 = mean(inside_50)
)
kable(
coverage_50,
digits = 3,
caption = "Observed ratings within their 50% prediction intervals"
)| batches_inside_50 | total_batches | proportion_inside_50 |
|---|---|---|
| 394 | 572 | 0.689 |
The 50% posterior prediction intervals contain 394 of the 572 observed ratings, or 68.9%. Coverage reasonably close to 50% supports calibration of the model’s middle-50% predictive uncertainty.
6 Exercise 10.19: Model Rating with Aftertaste
6.1 10.19.a Fit the aftertaste model
set.seed(84735)
aftertaste_model <- stan_glm(
total_cup_points ~ aftertaste,
data = new_coffee,
family = gaussian(),
prior = normal(0, 10, autoscale = FALSE),
prior_intercept = normal(75, 10, autoscale = FALSE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4,
iter = 5000,
seed = 84735,
refresh = 0
)
print(aftertaste_model, digits = 3)## stan_glm
## family: gaussian [identity]
## formula: total_cup_points ~ aftertaste
## observations: 572
## predictors: 2
## ------
## Median MAD_SD
## (Intercept) 33.160 1.423
## aftertaste 6.619 0.192
##
## Auxiliary parameter(s):
## Median MAD_SD
## sigma 1.566 0.046
##
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanregThe aftertaste model has the same structure and prior distributions as the aroma model. The only change is the predictor.
6.2 10.19.b Quick check for model misspecification
set.seed(84735)
pp_check(aftertaste_model, plotfun = "dens_overlay", nreps = 50) +
labs(
title = "Posterior Predictive Check for the Aftertaste Model",
x = "Total cup points"
)
If the observed density lies comfortably within the replicated densities, with no consistent displacement in center, spread, or tails, the aftertaste model captures the broad rating distribution reasonably well. Minor discrepancies would indicate approximation rather than complete model failure.
6.3 10.19.c Ten-fold cross-validated predictive quality
set.seed(84735)
aroma_cv <- prediction_summary_cv(
data = new_coffee,
model = aroma_model,
k = 10,
prob_inner = 0.50,
prob_outer = 0.95
)
set.seed(84735)
aftertaste_cv <- prediction_summary_cv(
data = new_coffee,
model = aftertaste_model,
k = 10,
prob_inner = 0.50,
prob_outer = 0.95
)
# Different bayesrules versions label the pooled result as either $cv or $pooled.
get_pooled_cv <- function(x) {
if (!is.null(x$cv)) return(x$cv)
if (!is.null(x$pooled)) return(x$pooled)
stop("The pooled cross-validation summary was not found.")
}
aroma_cv_result <- get_pooled_cv(aroma_cv)
aftertaste_cv_result <- get_pooled_cv(aftertaste_cv)
aroma_cv_result## mae mae_scaled within_50 within_95
## 1 0.875705 0.445226 0.676709 0.9580762aftertaste_cv_result## mae mae_scaled within_50 within_95
## 1 0.6894454 0.4390417 0.7166969 0.9684815cv_comparison <- bind_rows(
Aroma = as.data.frame(aroma_cv_result),
Aftertaste = as.data.frame(aftertaste_cv_result),
.id = "Predictor"
)
kable(
cv_comparison,
digits = 3,
caption = "Ten-fold cross-validated posterior predictive quality"
)| Predictor | mae | mae_scaled | within_50 | within_95 |
|---|---|---|---|---|
| Aroma | 0.876 | 0.445 | 0.677 | 0.958 |
| Aftertaste | 0.689 | 0.439 | 0.717 | 0.968 |
The median absolute error (mae) measures the typical
absolute difference between an observed rating and its posterior
predictive center. The scaled MAE expresses that error relative to
predictive uncertainty. Smaller values indicate better predictive
accuracy. The interval coverage statistics should be reasonably close to
their nominal probability levels.
6.4 10.19.d Choose between aroma and aftertaste
mae_column <- intersect(c("mae", "MAE"), names(cv_comparison))[1]
scaled_column <- intersect(
c("mae_scaled", "scaled_mae", "MAE_scaled"),
names(cv_comparison)
)[1]
best_predictor <- cv_comparison %>%
arrange(.data[[mae_column]]) %>%
slice(1) %>%
pull(Predictor)
best_predictor## [1] "Aftertaste"Based on 10-fold cross-validation, Aftertaste is the preferred single predictor because it has the smaller cross-validated MAE. The scaled MAE provides a second comparison after accounting for predictive spread. The size of the difference should also be considered: if the two error measures are nearly equal, the practical predictive advantage is small even when one model ranks first.
7 Conclusion
Sampling one observation per farm reduces the within-farm dependence present in the original data. Both aroma and aftertaste have positive relationships with total cup rating. Posterior predictive checks assess whether the Normal linear models reproduce important features of the observed data, while prediction errors, interval coverage, and 10-fold cross-validation evaluate predictive accuracy. The final choice between aroma and aftertaste is based on out-of-sample cross-validated error rather than the slope estimates alone.
8 Reference
Johnson, A. A., Ott, M. Q., and Dogucu, M. Bayes Rules! An Introduction to Applied Bayesian Modeling. Chapter 10, Exercises 10.13–10.19. https://www.bayesrulesbook.com/chapter-10