Statistical Rethinking 2026, A08

R
Bayes
item response model
MCMC
Using ulam and MCMC diagnostics on wine-tasting data.
Author

Paulina Sell

Published

Mar 2026

TipTL;DR

In this eighth exercise from the beginner section, we’re using ulam and MCMC diagnostics for the first time. We modify the item response model from the lecture to find out wether French and American judges differ in how they rate wine by estimating (total) causal contrast. Even though I’ve never heard of item response models before, this was not too difficult, because we really just modify what McElreath showed in his lecture.

Link to full course on GitHub
Link to online lecture

Exercise

“Modify the item response model from the lecture to ask if there are any average differences between French and American judges. Use ulam() and be sure to check chain convergence. Do the French and American judges differ in either their harshness or their discrimination?””

Load data

First, we load the package and data.

library(rethinking)
data(Wines2012)
data <- Wines2012

Let’s have a look:

'data.frame':   180 obs. of  6 variables:
 $ judge     : Factor w/ 9 levels "Daniele Meulder",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ flight    : Factor w/ 2 levels "red","white": 2 2 2 2 2 2 2 2 2 2 ...
 $ wine      : Factor w/ 20 levels "A1","A2","B1",..: 1 3 5 7 9 11 13 15 17 19 ...
 $ score     : num  10 13 14 15 8 13 15 11 9 12 ...
 $ wine.amer : int  1 1 0 0 1 1 1 0 1 0 ...
 $ judge.amer: int  0 0 0 0 0 0 0 0 0 0 ...

The data is from a wine-tasting event, where 20 wines (10 from France, 10 from New Jersey) have been tasted by 9 French and American judges.
Now, we want to find out whether there are any differences in harshness or discrimination between French and American judges.

DAG

Let’s use a DAG to think about this.

Q = Wine quality
S = Score
X = Wine origin
J = Judge preferences
Z = Judge origin

Statistical Model

The score \(S\) assigned to a wine is modeled as normally distributed around a mean \(\mu\):

\[ \text{S} \sim \mathcal{N}(\mu, \sigma) \] \[ \mu = (\text{Q}_W + \text{O}_X - \text{H}_Z) * \text{D}_Z \]

Inside the brackets, three additive components determine the underlying score before judge effects are applied: \(\text{Q}_W\) is the intrinsic quality of wine \(\text{W}\); \(\text{O}_X\)​ is a bonus or penalty depending on whether the wine’s origin \(\text{X}\) matches the judge’s prior expectations; and \(\text{H}_Z\) is the harshness of judge \(\text{Z}\), which shifts scores downward. This sum is then multiplied by \(\text{D}_Z\)​, the discrimination of judge \(\text{Z}\).

Priors

I don’t know what to expect from the model yet, so I’ll first set some default priors but I know that discrimination (\(\text{D}\)) and \(\sigma\) must be positive.

\[ \text{Q} \sim \mathcal{N}(0, 1) \] \[ \text{O} \sim \mathcal{N}(0,1) \] \[ \text{H} \sim \mathcal{N}(0,1) \] \[ \text{D} \sim \text{exponential}(1) \] \[ \sigma \sim \text{exponential}(1) \]

Prepare data

Now we need to prepare the data. This includes standardizing the outcome variable, \(\text{S}\).

data <- list(
  S = standardize(data$score),
  W = data$wine,
  X = ifelse(data$wine.amer == 1, 1, 2), # 1 = american wine, 2 = french wine
  Z = ifelse(data$judge.amer == 1, 1, 2) # 1 = american judge, 2 = french judge
)

Fit model

Let’s try and fit the model:

m.wines <- ulam(
  alist(
    S ~ dnorm(mu, sigma),
    mu <- (Q[W] + O[X] - H[Z]) * D[Z],
    Q[W] ~ dnorm(0, 1),
    O[X] ~ dnorm(0, 1),
    H[Z] ~ dnorm(0, 1),
    D[Z] ~ dexp(1),
    sigma ~ dexp(1)
  ),
  data = data,
  chains = 4
)
show(m.wines)
Hamiltonian Monte Carlo approximation
2000 samples from 4 chains

Sampling durations (seconds):
  chain_id warmup sampling total
1        1   0.29     0.15  0.43
2        2   0.22     0.17  0.39
3        3   0.23     0.16  0.39
4        4   0.23     0.16  0.39

Formula:
S ~ dnorm(mu, sigma)
mu <- (Q[W] + O[X] - H[Z]) * D[Z]
Q[W] ~ dnorm(0, 1)
O[X] ~ dnorm(0, 1)
H[Z] ~ dnorm(0, 1)
D[Z] ~ dexp(1)
sigma ~ dexp(1)
precis(m.wines, depth = 2)
              mean        sd        5.5%     94.5%     rhat  ess_bulk
Q[1]   0.126309575 0.9813506 -1.48226008 1.6981076 1.002720 2329.6184
Q[2]  -0.001095050 0.9385775 -1.51257624 1.4836040 1.000985 2581.4028
Q[3]   0.324823919 0.9295976 -1.17433415 1.7428587 1.006169 2913.5069
Q[4]   0.477692484 0.9289642 -1.04378094 1.9363799 1.000278 2897.4570
Q[5]  -0.209456757 0.9525154 -1.75571999 1.3328744 1.006877 3565.8587
Q[6]  -0.373732296 0.9236920 -1.85636866 1.1090363 1.003826 2836.3192
Q[7]   0.201993557 0.9221609 -1.30439300 1.6124867 1.004106 2467.9735
Q[8]   0.342879260 0.9437186 -1.17680077 1.8220423 1.004328 3093.0414
Q[9]   0.137084205 0.9446527 -1.34153194 1.6540625 1.000235 3530.6184
Q[10]  0.181911543 0.9116711 -1.28089111 1.6455186 1.003975 4138.7504
Q[11]  0.036758881 0.9171739 -1.43957058 1.5535230 1.002499 3033.9597
Q[12]  0.005378341 0.8694459 -1.39511658 1.4011498 1.002552 2832.7288
Q[13] -0.032693025 0.9604719 -1.61198397 1.5118294 1.002274 3160.3759
Q[14] -0.059601173 0.9456081 -1.57312773 1.5072017 1.000432 3128.2079
Q[15] -0.261984245 0.9423224 -1.79330216 1.2890327 1.002022 2907.0542
Q[16] -0.173822430 0.9697528 -1.70007329 1.3574755 1.004209 3585.5744
Q[17] -0.075508392 0.9230980 -1.52468738 1.3564489 1.000282 2965.1625
Q[18] -0.820807311 0.9159865 -2.30992737 0.6573951 0.998802 2707.5380
Q[19] -0.239046355 0.9321102 -1.70309866 1.2662142 1.000224 3055.0710
Q[20]  0.338630163 0.9320098 -1.16353990 1.8218597 1.000476 3323.7440
O[1]  -0.349069178 0.7619026 -1.56222536 0.8415831 1.000332 2208.1820
O[2]   0.319352696 0.8018144 -0.97331430 1.6193383 1.001018 2667.5292
H[1]  -0.462736442 0.8049304 -1.72348470 0.8848093 1.000234 2419.2475
H[2]   0.428964554 0.7999668 -0.85274317 1.7342868 1.003506 2353.2180
D[1]   0.118090698 0.0894043  0.01142418 0.2817639 1.001971  912.3974
D[2]   0.138417930 0.1077334  0.01295865 0.3354079 1.001895  841.5341
sigma  0.992958732 0.0543589  0.90805449 1.0802554 1.003431 2689.2884
trankplot(m.wines)

Waiting to draw page 2 of 2

Causal contrast

Total causal contrast in mean: wine discrimination

Code 
post <- extract.samples(m.wines)

mu_contrast <- post$D[, 2] - post$D[, 1] # french - american discrimination

dens(
  mu_contrast,
  lwd = 2,
  xlab = "posterior mean discrimination contrast"
)

Total causal contrast in predicted score: wine discrimination

Code 
# sample from posterior
D_A <- rnorm(1000, post$D[, 1], post$sigma)
D_F <- rnorm(1000, post$D[, 2], post$sigma)

# calculate contrast
S_contrast <- D_F - D_A

dens(
  S_contrast,
  lwd = 2,
  xlab = "posterior contrast in discrimination (French - American)"
)

Here we can see that 49.9% of the times we randomly sample a new judge, the French judge is more discriminating than the American judge. And 50.1% of the time, the American judge is more discriminating than the French judge. So it seems like there is no difference between American and French wine judges in terms of wine discrimination.
Let’s now look at harshness.

Total causal contrast in mean: harshness

Code 
mu_contrast_H <- post$H[, 2] - post$H[, 1] # french - american discrimination

dens(
  mu_contrast_H,
  lwd = 2,
  xlab = "posterior mean harshness contrast"
)

Total causal contrast in predicted score: harshness

Code 
# sample from posterior
H_A <- rnorm(1000, post$H[, 1], post$sigma)
H_F <- rnorm(1000, post$H[, 2], post$sigma)

# calculate contrast
S_contrast_H <- H_F - H_A

dens(
  S_contrast_H,
  lwd = 2,
  xlab = "posterior contrast in harshness (French - American)"
)

Here we can see that 69.9% of the times we randomly sample a new judge, the French judge is more harsh in their judgement than the American judge. And 30.1% of the time, the American judge is more harsh in their judgement than the French judge.

Lineplot

Now I want to try to recreate the lineplot that McElreath showed in his slides.

First, we need to create a sequence of wine quality (Q).

library(tidyverse)
Warning: package 'dplyr' was built under R version 4.4.3
# Simulation settings
n_lines <- 40
Q_seq <- seq(-2, 2, len = 60)
# 1 = American, 2 = French
Z_american <- 1
Z_french <- 2

Now, we need a function that gets the right values out of the posteriors. We first sample indices from the posterior, so that we “decide” a priori what row from the posterior we’ll use. In the for loop, we then fill the empty list we initiated with values from the full posterior in order to get the trajectories for a given Z.

make_lines <- function(Z_val, n_lines, Q_seq, post) {
  idx <- sample(nrow(post$sigma), n_lines, replace = FALSE)

  results <- vector("list", length(idx))

  for (i in seq_along(idx)) {
    s <- idx[i] # Full posterior draw: O, H & D come from the same sample row
    X_draw <- sample(1:2, 1) # random wine origin
    O <- post$O[s, X_draw]
    H <- post$H[s, Z_val]
    D <- post$D[s, Z_val]

    results[[i]] <- data.frame(
      Q = Q_seq,
      mu = (Q_seq + O - H) * D,
      line = i
    )
  }

  # Combine all the individual data frames into one
  do.call(rbind, results)
}

Now, we just call the function twice (for French and American judges) and then combine both data frames the data by row-binding it.

df_plot <- bind_rows(
  make_lines(Z_american, n_lines, Q_seq, post) |>
    mutate(judge_origin = "American"),
  make_lines(Z_french, n_lines, Q_seq, post) |>
    mutate(judge_origin = "French")
) |>
  mutate(judge_origin = factor(judge_origin, levels = c("American", "French")))

And now we’re ready to plot!

ggplot(df_plot, aes(x = Q, y = mu, group = line)) +
  geom_vline(xintercept = 0, col = "lightcoral") +
  geom_hline(yintercept = 0, col = "lightcoral") +
  geom_line(alpha = 0.35, linewidth = 0.5, colour = "#2c7bb6") +
  facet_wrap(~judge_origin, ncol = 2) +
  labs(
    x = "Wine quality (Q)",
    y = "Predicted score (μ)",
    title = "Predicted scores by judge origin",
    subtitle = paste0(n_lines, " posterior draws per panel")
  ) +
  theme_classic()

So: do the French and American judges actually differ in either their harshness or discrimination? Yes, they do differ in harshness. French judges tend to give lower scores in general compared to American judges. But their discrimination between wines does not differ.