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.22     0.17  0.39
2        2   0.27     0.11  0.38
3        3   0.26     0.14  0.40
4        4   0.25     0.09  0.34

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.123465749 0.97575650 -1.440763325 1.6495920 1.0003460 2570.3033
Q[2]   0.030662599 0.88132125 -1.356920072 1.4844007 1.0010745 2854.2619
Q[3]   0.327323118 0.95541807 -1.181727540 1.8013415 1.0014439 2983.2897
Q[4]   0.483026143 0.93642324 -1.001230703 1.9416481 1.0035814 2571.5348
Q[5]  -0.179842103 0.93373831 -1.723549261 1.3074550 1.0015178 2402.6979
Q[6]  -0.358023365 0.99750423 -1.982210917 1.1942719 1.0016167 2617.2790
Q[7]   0.205577491 0.94350133 -1.258083744 1.7145641 1.0015935 2597.1744
Q[8]   0.284033674 0.94930158 -1.278533821 1.7724713 1.0039858 3019.5553
Q[9]   0.110260261 0.92143379 -1.363597686 1.6065035 0.9996029 3231.1725
Q[10]  0.167190122 0.93393834 -1.310138719 1.5893373 1.0030825 3103.6843
Q[11]  0.041208426 0.94666437 -1.403878414 1.5749833 1.0014884 2702.7889
Q[12]  0.009942758 0.91949389 -1.445394309 1.4586559 1.0074807 3297.8660
Q[13] -0.024210488 0.93301530 -1.505807118 1.4799362 1.0051855 2409.2746
Q[14] -0.062376830 0.93070628 -1.480602800 1.4836175 0.9993438 2848.0866
Q[15] -0.262846971 0.94407006 -1.746395629 1.2486929 1.0014280 2383.0151
Q[16] -0.161475764 0.94117503 -1.634898402 1.3060047 1.0052496 2756.9533
Q[17] -0.078791764 0.92193321 -1.533535371 1.4449256 1.0005394 2730.6148
Q[18] -0.815084367 0.95226301 -2.298834491 0.7079721 1.0016923 2231.5821
Q[19] -0.226559537 0.89111522 -1.626016512 1.2478942 1.0034786 2588.9268
Q[20]  0.319105704 0.88721589 -1.109397795 1.7469635 1.0003520 2411.1255
O[1]  -0.340617529 0.78585758 -1.591898832 0.8974470 1.0032088 2106.2612
O[2]   0.297991781 0.80233085 -0.972382835 1.5586870 1.0035348 2073.8805
H[1]  -0.442185490 0.81060866 -1.663883983 0.9034959 1.0011919 2125.1326
H[2]   0.433571013 0.80061706 -0.865779704 1.6676488 0.9991557 1875.4711
D[1]   0.114872077 0.08842752  0.009598625 0.2791226 1.0066272  910.5951
D[2]   0.136709755 0.10046497  0.013884204 0.3235602 1.0038967  820.9098
sigma  0.993843641 0.05314539  0.912985762 1.0802499 1.0017855 2348.1853
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 51.6% of the times we randomly sample a new judge, the French judge is more discriminating than the American judge. And 48.4% 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.4% of the times we randomly sample a new judge, the French judge is more harsh in their judgement than the American judge. And 30.6% 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.