C05 · Likert scales are not numbers — Bayesian Modelling in Practice

The brief

A survey-research team ran a moral-judgement study. Each participant read 30 short stories in which someone's action leads to a death, and rated how morally permissible the action was from 1 (absolutely not) to 7 (absolutely). The stories vary three features:

action - the harm comes from doing something rather than from failing to act;
intention - the death is a means to the agent's goal, not a side effect;
contact - the agent physically touches the victim (in this coding contact stories are not also flagged as action).

The team analyses the ratings the way most survey teams do: ordinary linear regression on the numbers 1-7. A reviewer has told them that this is indefensible for Likert-type data. The team lead is unconvinced: "Everyone does it and the effects are huge. Show me one conclusion that would change."

Your job is to find out who is right - and to be honest about it in both directions.

How this notebook works

Each task states what to deliver, not how. Write your code in the YOUR CODE HERE cells.
Stuck? h.hint("task2") reveals hints one level at a time: nudge → approach → code skeleton. Try to get by on nudges.
h.check("task2", cut_first=...) compares your numbers with the reference solution.
A full worked solution lives in notebooks/solutions/. Open it only when you are done (or truly stuck).

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
from scipy import stats
from scipy.special import expit, logit

from pymc_challenges import Hints, data

RANDOM_SEED = 2006
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")

h = Hints("C05")
h.tasks()

C05 - Trolley problems: is a 1-7 rating a number?

task0 - Look at the outcome (2 hints)
task1 - The team's linear model and a check that exposes it (3 hints, 2 check(s))
task2 - Ordered logistic, cutpoints only (3 hints, 2 check(s))
task3 - Ordinal regression with predictors (3 hints, 4 check(s))
task4 - Predictions on the response scale (3 hints, 3 check(s))
task5 - Participant and story effects (3 hints, 2 check(s))
task6 - Does the metric shortcut change a conclusion? (3 hints, 4 check(s))

The data

One row per rating. response is the 1-7 rating, id the participant, story one of 12 story settings (a boat, a burning building, a switch...), and action, intention, contact the 0/1 features of the version of the story that was shown. age, male and edu describe the participant.

data.describe("trolley")
trolley = data.load("trolley")
trolley.head()

trolley: 9930 ordinal (1-7) permissibility ratings from 331 participants.
  source: Cushman, Young & Hauser (2006) moral-judgement experiment, via `rethinking`.
  url:    https://raw.githubusercontent.com/rmcelreath/rethinking/master/data/Trolley.csv

	case	response	order	id	age	edu	intention	contact	story	action2
0	cfaqu	4	2	96;434	14	Middle School	0	1	aqu	1
1	cfbur	3	31	96;434	14	Middle School	0	1	bur	1
2	cfrub	4	16	96;434	14	Middle School	0	1	rub	1
3	cibox	3	32	96;434	14	Middle School	1	1	box	1
4	cibur	3	4	96;434	14	Middle School	1	1	bur	1

Only six combinations of the three features occur in the experiment. We will call them scenarios and refer to them by these names throughout:

scenarios = pd.DataFrame(
    {"action": [0, 1, 0, 0, 1, 0], "intention": [0, 0, 0, 1, 1, 1], "contact": [0, 0, 1, 0, 0, 1]},
    index=pd.Index(
        ["none", "action", "contact", "intention", "action+intention", "contact+intention"],
        name="scenario",
    ),
)
key = ["action", "intention", "contact"]
trolley["scenario"] = pd.Categorical(
    trolley[key].merge(scenarios.reset_index(), on=key, how="left")["scenario"],
    categories=scenarios.index,
)
scenarios.assign(n=trolley.scenario.value_counts().reindex(scenarios.index))

	action	intention	contact	n
scenario
none	0	0	0	1655
action	1	0	0	2648
contact	0	0	1	993
intention	0	1	0	1986
action+intention	1	1	0	1655
contact+intention	0	1	1	993

Task 0 · Look at the outcome

Deliver

The distribution of response over its seven values, overall and for each scenario.
The cumulative proportions $P(\text{response} \le k)$ for $k = 1..6$ and the same numbers on the log-odds scale. Keep them - you will meet them again in Task 2.
Two sentences: which features of this outcome will a Normal likelihood be unable to reproduce?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A Normal distribution has one hump, symmetric thin tails and support on the whole real line. Go through those three properties one at a time with the bar chart of response in front of you - overall, and then separately for the most extreme scenario.

2 Approach

trolley.response.value_counts(normalize=True).sort_index() gives the seven shares, .cumsum() the cumulative proportions (drop the last one - it is 1 by construction) and scipy.special.logit puts them on the log-odds scale. pd.crosstab(trolley.scenario, trolley.response, normalize="index") gives the same seven shares per scenario.

Look at the spacing of the cumulative log-odds. If the seven categories were equally spaced slices of a smooth latent distribution, what would that spacing look like?

Task 1 · The team's model

This is the model in the team's paper:

$$\text{response}_i \sim \text{Normal}(\mu_i, \sigma), \qquad \mu_i = a + b_A A_i + b_I I_i + b_C C_i + b_{AI} A_i I_i + b_{CI} C_i I_i$$

Deliver

The fitted model with clean diagnostics.
A posterior predictive check that you design to show the team what, if anything, is wrong with it. The responses are seven integers and the scenarios are very different from one another - a single smoothed density overlay will not do. Make the check reusable: you will want to run the same check on a better model later.
One number for the abstract: what share of the model's predicted ratings are impossible (below 1 or above 7)?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Fitting is the easy part. For the check, ask: what can the data do that the model cannot, and what can the model do that the data cannot? The data can only be 1, 2, ..., 7. What does a draw from the model's posterior predictive look like? And is the misfit the same in every scenario, or worst where the mean is closest to an end of the scale?

2 Approach

Compare category shares: for each replicated dataset and each of the six scenarios, the share of ratings equal to 1, 2, ..., 7, shown as intervals on top of the observed bars. The linear model predicts real numbers, so give it its best shot: round to the nearest integer and clip to 1-7 (np.clip(np.rint(y_rep), 1, 7)), and report separately how many raw predictions fell below 1 or above 7.

You do not need pm.sample_posterior_predictive on 9,930 rows x 4,000 draws. All rows of one scenario share the same mu, so for a few hundred posterior draws (az.extract(idata, num_samples=400)) simulate n_scenario Normal draws directly with NumPy. Write the plotting part as a function of an array of replicated shares with shape (replicates, scenarios, categories) - Task 3 can then reuse it.

3 Skeleton

X = np.column_stack([A, I, C, A * I, C * I]).astype(float)
with pm.Model(coords={"coef": COEFS}) as linear_model:
    a = pm.Normal("a", ___, ___)
    b = pm.Normal("b", 0, ___, dims="coef")
    sigma = pm.HalfNormal("sigma", ___)
    pm.Normal("response", a + pm.math.dot(X, b), sigma, observed=trolley.response.values)
    linear_idata = pm.sample(random_seed=RANDOM_SEED)

post = az.extract(linear_idata, num_samples=400, random_seed=RANDOM_SEED)
mu = post["a"].values[:, None] + post["b"].values.T @ X_scen.T        # (400, 6)
raw = [rng.normal(mu[:, [s]], post["sigma"].values[:, None], size=(400, n_s[s]))
       for s in range(6)]
share_impossible = np.mean([___ for r in raw for value in ...])        # below 1 or above 7
freq = ...   # (400, 6, 7): np.bincount of the rounded, clipped ratings / n

Self-checks in the notebook (h.check("task1", …)):

sigma — posterior mean of the residual standard deviation of the linear model, about 1.8 rating points
share_impossible — share of posterior predictive draws below 1 or above 7 (a proportion, not a percentage), about 0.12

Task 2 · An ordinal likelihood

Leave the predictors aside for a moment. Model the seven categories as an ordered logistic outcome: a latent continuous "permissibility" with a logistic distribution, cut into seven bins by six cutpoints.

Deliver

Before fitting: what does your prior on the cutpoints imply for the seven category probabilities? Show it.
The fitted intercept-only model with clean diagnostics. If your first attempt dies before a single draw is taken, do not reach for a different sampler - find out which term of the log-probability is to blame at the starting point, and fix the model.
Compare the posterior cutpoints with a number you already have from Task 0. What are the cutpoints, in words?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

The category probabilities of an ordered-logistic model are differences of neighbouring cumulative probabilities. What has to be true of the cutpoints for all seven differences to be positive? Now look at where the sampler starts: model.initial_point() and model.point_logps() tell you the starting values and which term of the log-probability is infinite there.

2 Approach

Two things are needed and neither is enough alone:

transform=pm.distributions.transforms.ordered on the cutpoint prior, so the sampler moves in an unconstrained space where every point maps to increasing cutpoints;
a strictly increasing starting value. The default start for a vector of Normal(0, s) variables is all zeros, and the ordered transform takes log of the gaps between neighbours. Pass e.g. initvals={"cut": np.linspace(-2, 2, 6)} to pm.sample. (Setting initval= on the distribution also works for sampling, but then pm.compute_log_likelihood refuses the model later - you will want it in Task 5.)

pm.OrderedLogistic expects the outcome coded 0..6 and uses P(y <= k) = logistic(cut_k - eta). With eta=0 the cutpoints are simply the log-odds of the cumulative proportions - compare with Task 0.

For the prior: pm.sample_prior_predictive ignores transforms, so it will give you unsorted cutpoints. The prior you are actually using is the distribution of six sorted Normal draws, which is a one-liner in NumPy.

3 Skeleton

ordered = pm.distributions.transforms.ordered
y = trolley.response.values - 1
with pm.Model(coords={"cutpoint": CUTS}) as cut_model:
    cut = pm.Normal("cut", 0, ___, dims="cutpoint", transform=ordered)
    pm.OrderedLogistic("response", eta=0.0, cutpoints=cut, observed=y)
    cut_idata = pm.sample(initvals={"cut": ___}, random_seed=RANDOM_SEED)

def ordinal_probs(eta, cut):            # eta (...,), cut (..., 6) -> (..., 7)
    cdf = expit(cut - np.asarray(eta)[..., None])
    cdf = np.concatenate([np.zeros_like(cdf[..., :1]), cdf, np.ones_like(cdf[..., :1])], -1)
    return np.diff(cdf, axis=-1)

prior_cut = np.sort(rng.normal(0, ___, size=(500, 6)), axis=1)
prior_p = ordinal_probs(np.zeros(500), prior_cut)

Self-checks in the notebook (h.check("task2", …)):

cut_first — posterior mean of the lowest cutpoint (between responses 1 and 2), the log-odds of P(response <= 1), about -1.9
cut_last — posterior mean of the highest cutpoint (between 6 and 7), the log-odds of P(response <= 6), about 1.8

Task 3 · Ordinal regression

Add the team's predictors - action, intention, contact and the two interactions with intention - to the ordinal model, with the convention that a positive coefficient means "more permissible".

Deliver

Priors for the coefficients that you can defend on the log-odds scale.
The fitted model with clean diagnostics. (Memory warning: by default pm.OrderedLogistic stores the seven category probabilities of every row for every draw. With 9,930 rows that is gigabytes. Find the switch.)
Your Task 1 check, run on this model and shown next to the linear one.

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Where does the linear predictor go, and with which sign? In pm.OrderedLogistic, P(y <= k) = logistic(cut_k - eta): a larger eta makes every "at most k" less likely, i.e. pushes ratings up. And think before you add an intercept to eta: what would happen to the posterior if you added 1 to the intercept and to all six cutpoints?

2 Approach

eta = pm.math.dot(X, b) with no intercept - the cutpoints are the intercepts. On the log-odds scale a coefficient of 1 is a big effect, so something like Normal(0, 0.5) is a reasonable weakly-informative prior.

Pass compute_p=False to pm.OrderedLogistic, otherwise a (chain, draw, 9930, 7) array of probabilities is stored in the posterior.

For the check you need replicated category counts per scenario. Compute the seven probabilities per scenario and draw (six distinct eta values only - use the six-row design matrix of the scenarios), then rng.multinomial(n_scenario, p).

3 Skeleton

with pm.Model(coords={"cutpoint": CUTS, "coef": COEFS}) as ordinal_model:
    cut = pm.Normal("cut", 0, 1.5, dims="cutpoint", transform=ordered)
    b = pm.Normal("b", 0, ___, dims="coef")
    pm.OrderedLogistic("response", eta=___, cutpoints=cut, observed=y, compute_p=False)
    ordinal_idata = pm.sample(initvals={"cut": CUT_INIT}, random_seed=RANDOM_SEED)

post = az.extract(ordinal_idata, var_names=["b", "cut"], num_samples=400)
eta = post["b"].values.T @ X_scen.T                                # (400, 6)
p = ordinal_probs(eta, post["cut"].values.T[:, None, :])           # (400, 6, 7)
freq = np.stack([[rng.multinomial(n_s[s], p_d[s]) / n_s[s] for s in range(6)] for p_d in p])

Self-checks in the notebook (h.check("task3", …)):

b_action — posterior mean of the action coefficient (positive = more permissible), about -0.47
b_intention — posterior mean of the intention coefficient, about -0.29
b_contact — posterior mean of the contact coefficient, about -0.34
b_contact_intention — posterior mean of the contact x intention interaction, about -1.2 (the largest effect)

Task 4 · Back to the scale people answered on

"The contact:intention coefficient is -1.2 on the latent log-odds scale" is not a sentence you can put in front of the survey team.

Deliver, from the ordinal model, with 94% intervals:

A figure of the predicted distribution over the seven categories for each scenario.
A table with the expected rating and $P(\text{response} \le 3)$ ("leans impermissible") for each scenario.
The headline contrast: going from none to contact+intention, how much does the expected rating fall, and by how many percentage points does $P(\text{response} \le 3)$ rise?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

A coefficient on the latent scale does not translate into a fixed change in the rating: the same shift in eta moves a lot of probability when the cutpoints around it are close together, and little when you are already near 1 or 7. So do not interpret coefficients - compute the seven probabilities for each scenario, for every posterior draw, and summarise whatever you care about from those.

2 Approach

For each posterior draw and scenario: eta = b @ x_scenario, then P(y <= k) = expit(cut_k - eta), pad with 0 and 1 and take np.diff to get the seven probabilities. The expected rating is (p * [1..7]).sum(-1), P(response <= 3) is p[..., :3].sum(-1). Contrasts between scenarios are differences per draw, then summarised - not differences of summaries.

3 Skeleton

post = az.extract(ordinal_idata, var_names=["b", "cut"])
eta = post["b"].values.T @ X_scen.T                                # (draws, 6)
P = ordinal_probs(eta, post["cut"].values.T[:, None, :])           # (draws, 6, 7)
mean_rating = (P * np.arange(1, 8)).sum(-1)                        # (draws, 6)
p_le3 = P[..., :3].sum(-1)
shift = mean_rating[:, ___] - mean_rating[:, ___]                  # per draw
shift.mean(), np.quantile(shift, [0.03, 0.97])

Self-checks in the notebook (h.check("task4", …)):

mean_none — posterior mean of the expected rating in the 'none' scenario, about 4.8
p_le3_contact_intention — posterior mean of P(response <= 3) in 'contact+intention', about 0.63 (a probability)
shift_mean — change in expected rating from 'none' to 'contact+intention', about -1.9 (negative - less permissible)

Task 5 · People differ

Every participant rated 30 stories, and so far the models have treated their 30 answers as independent. Some people are strict, some lenient; some stories are simply more upsetting than others whatever their features.

To stay inside a laptop runtime budget this task uses a random subsample of 150 participants (4,500 ratings), defined below. Use it for every model in this task so that they can be compared.

Deliver

An ordinal model with a varying intercept per participant and an effect per story, with clean diagnostics. 150 correlated intercepts next to six cutpoints is a location-identification problem waiting to happen - deal with it in the model.
How large is the between-participant standard deviation compared with the biggest treatment effect? Put it in words.
A LOO comparison with the Task 3 model refitted to the same subsample. How much do the participant and story effects matter for prediction?
Put the coefficients of the two models side by side. They all moved away from zero. Explain why that is expected and not a contradiction.

# a dedicated generator, so that everybody gets the same 150 people whatever was run before
sub_ids = np.sort(np.random.default_rng(RANDOM_SEED).choice(trolley.id.unique(), size=150, replace=False))
sub = trolley[trolley.id.isin(sub_ids)].reset_index(drop=True)
print(f"{sub.id.nunique()} participants, {len(sub)} ratings")

150 participants, 4500 ratings

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

The six cutpoints already pin down where the latent scale sits. If all 150 participant intercepts could drift up together while the cutpoints drift up with them, nothing in the likelihood would notice. How can you give each participant an offset while making that joint drift impossible? You met the tool in E02 and C02.

2 Approach

u = pm.ZeroSumNormal("u", sigma=sigma_participant, dims="participant") and the same for story, added to eta via integer indices from pd.factorize(..., sort=True). With 30 ratings per participant the centred form samples well. Keep compute_p=False, keep the ordered transform and initvals.

For LOO both models need a log_likelihood group: pm.compute_log_likelihood(idata) inside the model context (this is the step that fails if you used initval= on the cutpoints instead of initvals= in pm.sample). Then az.compare({...}). Divide the elpd difference by the number of ratings to get a per-rating figure you can explain.

For part 4, think of what a population-level logit model is: an average over people of curves that are steep but shifted relative to each other.

3 Skeleton

pid, participants = pd.factorize(sub.id, sort=True)
sid, stories = pd.factorize(sub.story, sort=True)
coords = {"cutpoint": CUTS, "coef": COEFS, "participant": participants, "story": stories}
with pm.Model(coords=coords) as hier_model:
    cut = pm.Normal("cut", 0, 1.5, dims="cutpoint", transform=ordered)
    b = pm.Normal("b", 0, 0.5, dims="coef")
    sigma_participant = pm.HalfNormal("sigma_participant", ___)
    sigma_story = pm.HalfNormal("sigma_story", ___)
    u = pm.ZeroSumNormal("u", ___, dims="participant")
    v = pm.ZeroSumNormal("v", ___, dims="story")
    eta = pm.math.dot(X_sub, b) + u[___] + v[___]
    pm.OrderedLogistic("response", eta=eta, cutpoints=cut, observed=y_sub, compute_p=False)
    hier_idata = pm.sample(initvals={"cut": CUT_INIT}, nuts_sampler="nutpie",
                           random_seed=RANDOM_SEED)
    pm.compute_log_likelihood(hier_idata)
az.compare({"population": pop_sub_idata, "hierarchical": hier_idata})

Self-checks in the notebook (h.check("task5", …)):

sigma_participant — posterior mean of the between-participant sd on the log-odds scale, about 2.05
elpd_diff — elpd_loo of the hierarchical model minus the population model on the 150-participant subsample (positive, about 1500)

Task 6 · So was the reviewer right?

Time to answer the team lead. Use the two population models fitted to the full data (Tasks 1 and 3) and, for the last part, the hierarchical model.

Deliver

For every pair of scenarios (15 pairs) the effect of moving from one to the other according to each model and according to the raw data, on three quantities: the expected rating, $P(\text{response} \le 3)$ and $P(\text{response} = 1)$ ("absolutely impermissible"). For the linear model, take its predictive distribution rounded to the nearest category, with everything below 1.5 counted as 1 and everything above 6.5 as 7. For which pair, and which quantity, do the two models disagree most - and which of them is closer to the data?
The linear model says the none → contact+intention drop is the same number of rating points for everybody. What does the hierarchical ordinal model say for a typical participant, a strict one (2 sd below) and a lenient one (2 sd above)? Is there any sign of this in the raw data?
A paragraph for the team lead: which of their conclusions survive, which do not, and what should they change?

Your turn Work this task in the notebook.

Hints for this task — open one at a time

1 Nudge

Count parameters before you compute anything: six scenario means, and a linear model with an intercept and five coefficients. What does that imply for how well the linear model reproduces the six mean ratings - whatever the likelihood? Then ask which quantities depend on more than the mean.

2 Approach

Linear model: per draw and scenario, mu is the expected rating; the rounded predictive distribution is np.diff(norm.cdf((edges - mu) / sigma)) with edges = [-inf, 1.5, 2.5, ..., 6.5, inf]. Ordinal model: the (draws, 6, 7) probability array from Task 4. Reduce both to posterior means per scenario, build all 15 pairwise differences, and take abs(linear - ordinal).

Participant part: in the hierarchical model a participant z standard deviations from the average has eta = b @ x + z * sigma_participant (average story: story effect 0). Convert to an expected rating for 'none' and 'contact+intention' and difference per draw. For the raw-data comparison, group participants by their mean rating on the other four scenarios so the grouping does not use the two cells you are comparing.

3 Skeleton

edges = np.array([-np.inf, ___, np.inf])        # where does a Normal 'round' to each of 1..7?
P_lin = np.diff(stats.norm.cdf((edges - mu_lin[..., None]) / sigma[:, None, None]), axis=-1)
p1_lin, p1_ord = P_lin[..., 0].mean(0), P_ord[..., 0].mean(0)       # per scenario
pairs = [(i, j) for i in range(6) for j in range(i + 1, 6)]
gap_p1 = [abs((p1_lin[j] - p1_lin[i]) - (p1_ord[j] - p1_ord[i])) for i, j in pairs]

eta_z = b_h @ X_scen[[NONE, CI]].T + ___                                # (draws, 2): shift by a participant z sd's from average
m = (ordinal_probs(eta_z, cut_h[:, None, :]) * np.arange(1, 8)).sum(-1)
drop = m[:, 1] - m[:, 0]

Self-checks in the notebook (h.check("task6", …)):

max_gap_mean — largest absolute disagreement between the two models about a shift in EXPECTED RATING over the 15 pairs, in rating points - tiny, below 0.06
max_gap_p1 — largest absolute disagreement about a shift in P(response = 1), as a proportion, about 0.066
drop_typical — drop in expected rating none -> contact+intention for a typical participant under the hierarchical ordinal model, about -2.2
drop_lenient — the same drop for a participant 2 sd more lenient than average, about -1.0 (much smaller - ceiling)

Going further

The ordinal model over-predicts the midpoint 4 in the two harshest scenarios. Relax proportional odds - for example let intention shift the cutpoints rather than the latent mean, or add a "fence-sitter" mixture component at 4 - and see whether the Task 1 check comes clean.
Give every participant their own contact:intention slope (correlated with their intercept via pm.LKJCholeskyCov). Does it reproduce the gradient in the observed drops across the leniency fifths that the intercept-only model missed?
edu is itself ordinal. Model it as a monotonic predictor: a total effect times the cumulative sum of a pm.Dirichlet simplex over the education steps. Add age and male.
Fit the hierarchical model to all 331 participants (a few minutes with nutpie). Do the subsample's conclusions hold up?

h.progress()

C05: 0/20 hints used, 17/17 checks passed.

Reference solution

A complete, executed solution with commentary on every task — including the ways the obvious model fails. Try the hint ladder first.

Open the worked solution →

#The brief

#How this notebook works

#The data

#Task 0 · Look at the outcome

#Task 1 · The team's model

#Task 2 · An ordinal likelihood

#Task 3 · Ordinal regression

#Task 4 · Back to the scale people answered on

#Task 5 · People differ

#Task 6 · So was the reviewer right?

#Going further

Reference solution

The brief

How this notebook works

The data

Task 0 · Look at the outcome

Task 1 · The team's model

Task 2 · An ordinal likelihood

Task 3 · Ordinal regression

Task 4 · Back to the scale people answered on

Task 5 · People differ

Task 6 · So was the reviewer right?

Going further