In this article we’ll analyze data from a blind test. Blind tests are fun experiments that can determine, for instance, if family members or friends can distinguish between:

• Pepsi and Cola
• Red and white wines
• Brand-name and generic foods

Analysis of such data is interesting, especially when you’ve designed the experiment and collected the data yourself. You might learn some statistics along the way, and that knowledge could prove useful in industry settings as well as in everyday situations.

The beer tasting experiment

A friend wanted to determine her true beer preferences, so we decided to run a blind test:

• We bought 6 different brands of beer.
• She tasted each brand on 3 different occasions.
• The ordering of the 18 tastings was randomized.
• Each of the 18 tastings was conducted on a separate day.

In every tasting, the unknown beer brand was given a rating on a scale from 1 to 10. Knowledge of which brands were given which rating was not disclosed until the experiment was finished. The results of the full experiment are shown in the figure below, and the full dataset is available at the end of this article.

The Oslo-based beer brand “Schous” is the winner, with an average rating of

Can my friend tell beers apart?

Let us examine the figure above a bit more critically. The ratings are somewhat consistent, in the sense that they do not look like random integers between 1 and 10. However, it is surprising that the very same beer is awarded both the lowest observed rating (3) and the highest observed rating (8).

This raises the question: can my friend consistently tell beers apart? Let’s test this with a permutation hypothesis test. We’ll test if the explained variance \(r^2\) in the dataset is higher than what could expect if knowledge of beer brands contained no information.

Let \(y_i\) be the rating at tasting \(i\). Defint the variable \(x_{ij}\) to be a one-hot encoding of the brands: it is equal to \(1\) if beer brand \(j\) was tasted in tasting \(i\), otherwise it is equal to \(0\). We’ll assume a linear model with structure

where \(\epsilon_i\) is normally distributed noise. The coefficients \(\beta_j\) have a natural interpretation: they are the average ratings of each brand \(j\).

Fitting this model to the data gives an explained variance \(r^2 = 0.375\). The beer brand explains around one third of the total variane in the ratings. To put this number into perspective, consider that:

• A perfectly consistent beer-rater would achieve \(r^2 = 1\) on a training set (and on a test set), since the beer brand would explain their ratings perfectly and account for all variation.
• A completely random beer-rater would on average achieve \(r^2 \approx 0\) on a test set, since the beer brand would not explain the ratings at all. However, on a training set they would achieve \(r^2 > 0\) since there are six degrees of freedom in the brand coefficients \(\beta_j\), and these will overfit and capture some of the variance.

Let the null hypothesis \(H_0\) be that the beer brands have no effect on the ratings at all. Permuting the rows of the data set randomly (so there is no relationship between brands and ratings) produces the following distribution of \(r^2\) under the null hypothesis.

• If \(H_0\) is true (brands have no effect), we will on average observe an \(r^2\) of around \(0.292\). The variation in \(r^2\)-values in the permuted data sets is large, and anything between \(0.08\) and \(0.57\) is relatively common (coverage of 90%).
• The observed \(r^2\) is \(0.375\), which is certainly higher than the average, but by no means extreme.
• The \(p\)-value, which is the probability of observing a value of \(0.375\) or higher, given that \(H_0\) is true, is around \(0.289\).

Scientists typically reject \(H_0\) if the \(p\)-value is lower than \(0.05\). Here the \(p\)-value is much larger; this level of consistency in the ratings is well within the realm of what could happen even if my friend cannot tell beers apart. It might be that any apparent consistency is attributable to chance rather than genuine tasting ability.

Summary

The brand “Schous” is the winner, but there is a concerning amount of variation in the ratings within each brand. It’s not clear that my friend can consistently tell beers apart in this experiment. The \(p\)-value is so high that any scientific paper would report a “no effect” conclusion, failing to discard \(H_0\). That being said, I am also reluctant to decisively claim that my friend cannot tell beers apart. While \(H_0\) is mathematically convenient, it imposes a somewhat one-sided burden of proof.

People’s tasting abilities may be less consistent than they believe. The paper An Examination of Judge Reliability at a major U.S. Wine Competition investigates this to some extent. A larger, more involved experiment might be needed to figure out if my friend can tell beers apart. For now, the data suggests Schous remains the preferred choice.

Dataset

Here is the dataset used in this article in .csv format:

order,beer_brand,rating
1,Frydenlund,3
2,CB,7
3,Pokal,5
4,Tuborg,3
5,Frydenlund,5
6,Tuborg,8
7,Frydenlund,5
8,Hansa,6
9,CB,7
10,CB,5
11,Pokal,3
12,Schous,7
13,Hansa,3
14,Schous,8
15,Hansa,8
16,Tuborg,4
17,Schous,7
18,Pokal,5

Python code

Here’s a sketch of Python code to perform the statistical test. Python 3.12.8 was used with Pandas 2.2.3, NumPy 1.26.4 and scikit-learn 1.5.2.

import pandas as pd
from sklearn.preprocessing import OneHotEncoder
from sklearn.linear_model import LinearRegression
import numpy as np

# Read dataset
df = pd.read_csv("dataset.csv")

# Create X and y
encoder = OneHotEncoder(drop=None)
X = encoder.fit_transform(df[["beer_brand"]])
y = df["rating"].values

# Fit linear model, compute r2
linreg = LinearRegression(fit_intercept=False)
linreg.fit(X, y)
r2 = linreg.score(X, y)
print(f"r^2={r2:.3f}")  # r^2=0.375

# Randomly permute data and fit models
rng = np.random.default_rng(42)
r2_scores = []
for simulation in range(999):
    y_permuted = rng.permutation(y)
    linreg.fit(X, y_permuted)
    r2_scores.append(linreg.score(X, y_permuted))

p_value = np.mean(np.array(r2_scores) >= r2)
print(f"p_value={p_value:.3f}")  # p_value=0.311