Cumulative Ordered Probit / Logit — R (MASS::polr)¶
References:
- Liddell, T.M. & Kruschke, J.K. (2018). Analyzing ordinal data with metric models: What could possibly go wrong? J Exp Psychol Gen 147(12):622–647.
- McCullagh, P. (1980). Regression models for ordinal data. JRSS-B 42(2):109–142.
- Agresti, A. (2002). Categorical Data Analysis (2nd ed.). Wiley.
Model — Proportional Odds (Cumulative Link)¶
For group $g$ and rater $i$:
$$P(y_{ig} \le k) = F(\tau_k - \beta_g), \quad k = 1,\ldots,K$$
Equivalently via a latent continuous variable: $$z_{ig} = \beta_g + \varepsilon_{ig}, \qquad y_{ig} = k \iff \tau_{k-1} < z_{ig} \le \tau_k$$
The link function $F$ determines the error distribution:
polr(method=...) |
$F$ | $\varepsilon_{ig}$ | Latent variance |
|---|---|---|---|
"probit" |
$\Phi$ (std. normal CDF) | $N(0,1)$ | 1 |
"logistic" |
logistic / sigmoid | Logistic$(0,1)$ | $\pi^2/3 \approx 3.29$ |
Key constraint — homoscedastic. $\beta_g$ captures only location; all groups share the
same latent error variance. Per-group scale estimation ($\sigma_g$) requires the heteroscedastic
A&C Gibbs model in ord_probit_python.ipynb.
Notebooks¶
| File | Language | Model | Estimation |
|---|---|---|---|
| ord_probit_r.ipynb (this) | R | Proportional odds — homoscedastic | MASS::polr (MLE, probit + logit) |
| ord_probit_python.ipynb | Python | Heteroscedastic A&C (per-group $\mu_g, \sigma_g$) | Custom Gibbs |
| ord_probit_pymc.ipynb | Python | Heteroscedastic A&C (per-group $\mu_g, \sigma_g$) | PyMC / NUTS |
| ord_probit_gibbs.py | Python | — | A&C Gibbs sampler module |
1. Setup¶
userLib <- file.path(Sys.getenv("USERPROFILE"), "R", "win-library", "4.6")
.libPaths(c(userLib, .libPaths()))
library(MASS)
K <- 5L
SCALE <- pi / sqrt(3) # logit -> probit scale conversion
cat("MASS:", as.character(packageVersion("MASS")), "\n")
cat("K =", K, "ordered categories\n")
cat("Logit->probit scale factor:", round(SCALE, 4), "\n")
MASS: 7.3.65 K = 5 ordered categories Logit->probit scale factor: 1.8138
2. Synthetic Data¶
Four groups crossing low/high latent quality with low/high polarisation (Liddell & Kruschke 2018, Fig. 1). Shared true thresholds: $\boldsymbol{\tau} = (1.5, 2.5, 3.5, 4.5)$. N = 200 per group.
| Group | True $\mu$ | True $\sigma$ | Interpretation |
|---|---|---|---|
| 1 | 3.0 | 0.8 | Mediocre, consistent |
| 2 | 4.5 | 0.8 | High quality, consistent |
| 3 | 3.0 | 2.0 | Same mean as g1, polarising |
| 4 | 4.5 | 2.0 | Same mean as g2, polarising |
polr estimates a single $\hat{\beta}_g$ per group (location shift relative to a reference).
Validation means checking that the ranked order of $\hat{\beta}$ matches the true latent quality
order, and that the probit and logit links agree (after scale conversion).
set.seed(42L)
mu_true <- c(3.0, 4.5, 3.0, 4.5)
sigma_true <- c(0.8, 0.8, 2.0, 2.0)
tau_true <- c(1.5, 2.5, 3.5, 4.5)
G_s <- length(mu_true)
N_per <- 200L
counts_s <- matrix(0L, G_s, K)
for (g in seq_len(G_s)) {
z <- rnorm(N_per, mu_true[g], sigma_true[g])
y <- findInterval(z, tau_true) + 1L
counts_s[g, ] <- tabulate(y, nbins = K)
}
grp_labels <- sprintf("g%d (mu=%.1f sg=%.1f)", seq_len(G_s), mu_true, sigma_true)
emp_mean_s <- as.numeric((counts_s %*% seq_len(K)) / rowSums(counts_s))
rownames(counts_s) <- grp_labels
cat("Counts (rows = groups, cols = categories 1..5):\n")
print(counts_s)
cat("\nEmpirical means:\n")
print(round(setNames(emp_mean_s, grp_labels), 3))
# Long format for polr
y_s <- unlist(lapply(seq_len(G_s), function(g) rep(seq_len(K), counts_s[g, ])))
grp_s <- rep(seq_len(G_s), rowSums(counts_s))
df_s <- data.frame(y = ordered(y_s), g = factor(grp_s))
cat(sprintf("\nLong format: %d obs groups: %d (ref = g1)\n",
nrow(df_s), nlevels(df_s$g)))
Counts (rows = groups, cols = categories 1..5):
[,1] [,2] [,3] [,4] [,5]
g1 (mu=3.0 sg=0.8) 5 44 100 47 4
g2 (mu=4.5 sg=0.8) 0 1 14 90 95
g3 (mu=3.0 sg=2.0) 46 37 47 23 47
g4 (mu=4.5 sg=2.0) 15 23 32 36 94
Empirical means:
g1 (mu=3.0 sg=0.8) g2 (mu=4.5 sg=0.8) g3 (mu=3.0 sg=2.0) g4 (mu=4.5 sg=2.0)
3.005 4.395 2.940 3.855
Long format: 800 obs groups: 4 (ref = g1)
3. MASS::polr — Synthetic Data¶
polr fits the cumulative link model by maximum likelihood (Newton–Raphson).
Output:
coef(fit): group effects $\hat{\beta}_g$ on the latent scale (reference group = 0)fit$zeta: threshold estimates $\hat{\tau}_k$ on the latent scaleconfint(fit): profile-likelihood 95% CI for each coefficient
Fitting both links lets us verify that probit and logit give consistent rankings and that their estimates agree after scaling by $\pi/\sqrt{3} \approx 1.814$.
fit_sp <- polr(y ~ g, data = df_s, method = "probit", Hess = TRUE)
fit_sl <- polr(y ~ g, data = df_s, method = "logistic", Hess = TRUE)
coef_sp <- c(0, coef(fit_sp)) # prepend reference group (g1 = 0)
coef_sl <- c(0, coef(fit_sl))
cat("polr: group effects (probit scale, ref = g1 = 0):\n")
cat(sprintf(" %-28s %8s %8s %8s\n",
"Group", "Probit", "Logit", "L/scale"))
cat(" ", paste(rep("-", 56), collapse = ""), "\n", sep = "")
for (g in seq_len(G_s))
cat(sprintf(" %-28s %8.4f %8.4f %8.4f\n",
grp_labels[g], coef_sp[g], coef_sl[g], coef_sl[g] / SCALE))
cat("\npolr probit thresholds:\n")
print(round(fit_sp$zeta, 4))
cat("\npolr logit thresholds:\n")
print(round(fit_sl$zeta, 4))
cat(sprintf("\nAIC probit=%.1f logit=%.1f\n", AIC(fit_sp), AIC(fit_sl)))
polr: group effects (probit scale, ref = g1 = 0):
Group Probit Logit L/scale
--------------------------------------------------------
g1 (mu=3.0 sg=0.8) 0.0000 0.0000 0.0000
g2 (mu=4.5 sg=0.8) 1.3130 2.1248 1.1715
g3 (mu=3.0 sg=2.0) -0.0331 -0.1287 -0.0709
g4 (mu=4.5 sg=2.0) 0.8189 1.5149 0.8352
polr probit thresholds:
1|2 2|3 3|4 4|5
-1.0721 -0.4191 0.3883 1.1503
polr logit thresholds:
1|2 2|3 3|4 4|5
-1.8850 -0.6975 0.6631 1.9528
AIC probit=2247.8 logit=2242.5
cat("Profile-likelihood 95% CI for group effects (probit):\n")
ci_sp <- suppressMessages(confint(fit_sp)) # profile CI, may print trace
cat(sprintf(" %-28s %8s %8s %8s\n",
"Group", "Coef", "CI_lo", "CI_hi"))
cat(" ", paste(rep("-", 56), collapse = ""), "\n", sep = "")
# g1 = 0 (reference, no CI)
cat(sprintf(" %-28s %8.4f %8s %8s (reference)\n",
grp_labels[1L], 0, "-", "-"))
for (g in seq(2L, G_s)) {
cf <- coef(fit_sp)[g - 1L]
cat(sprintf(" %-28s %8.4f %8.4f %8.4f\n",
grp_labels[g], cf, ci_sp[g - 1L, 1L], ci_sp[g - 1L, 2L]))
}
cat("\nGroup ranking (probit):",
paste(grp_labels[order(-coef_sp)], collapse = " > "), "\n")
cat("True quality (mu) order: g2=g4 > g1=g3 (tied within each mu pair)\n")
cat("polr conflates location AND scale: consistent films rank above polarising\n")
cat("ones with the same true mu, because polarisation creates more 1-star ratings\n")
cat("which shift cumulative probabilities downward.\n")
Profile-likelihood 95% CI for group effects (probit): Group Coef CI_lo CI_hi -------------------------------------------------------- g1 (mu=3.0 sg=0.8) 0.0000 - - (reference) g2 (mu=4.5 sg=0.8) 1.3130 1.0921 1.5350 g3 (mu=3.0 sg=2.0) -0.0331 -0.2395 0.1733 g4 (mu=4.5 sg=2.0) 0.8189 0.6069 1.0314 Group ranking (probit): g2 (mu=4.5 sg=0.8) > g4 (mu=4.5 sg=2.0) > g1 (mu=3.0 sg=0.8) > g3 (mu=3.0 sg=2.0) True quality (mu) order: g2=g4 > g1=g3 (tied within each mu pair) polr conflates location AND scale: consistent films rank above polarising ones with the same true mu, because polarisation creates more 1-star ratings which shift cumulative probabilities downward.
Synthetic results (validated)¶
polr: group effects (probit scale, ref = g1 = 0):
Group Probit Logit L/scale
--------------------------------------------------------
g1 (mu=3.0 sg=0.8) 0.0000 0.0000 0.0000
g2 (mu=4.5 sg=0.8) 1.3130 2.1248 1.1715
g3 (mu=3.0 sg=2.0) -0.0331 -0.1287 -0.0709
g4 (mu=4.5 sg=2.0) 0.8189 1.5149 0.8352
polr probit thresholds:
1|2 2|3 3|4 4|5
-1.0721 -0.4191 0.3883 1.1503
AIC probit=2247.8 logit=2242.5
Profile-likelihood 95% CI for group effects (probit):
Group Coef CI_lo CI_hi
g1 (mu=3.0 sg=0.8) 0.0000 - - (reference)
g2 (mu=4.5 sg=0.8) 1.3130 1.0921 1.5350
g3 (mu=3.0 sg=2.0) -0.0331 -0.2395 0.1733
g4 (mu=4.5 sg=2.0) 0.8189 0.6069 1.0314
Ranking: g2 > g4 > g1 > g3
Analysis — synthetic¶
What polr gets right:
- Both high-quality groups (g2, g4; true $\mu=4.5$) rank above both mediocre groups (g1, g3; true $\mu=3.0$) with non-overlapping 95% CIs. The two quality tiers are cleanly separated.
- g1 and g3 (true $\mu=3.0$) are statistically indistinguishable: $\hat{\beta}_{g3}=-0.033$, CI $[-0.24, 0.17]$ includes zero. Correct.
What polr gets wrong — the scale conflation:
- g2 ($\mu=4.5$, $\sigma=0.8$) outranks g4 ($\mu=4.5$, $\sigma=2.0$) by 0.49 probit units, even though both have the same true quality. The consistent film (g2) appears substantially better because its ratings concentrate in categories 4–5; the polarising film (g4) scatters 15–23% of ratings into categories 1–3, which the homoscedastic model interprets as lower location.
- Symmetrically, g3 ($\mu=3.0$, $\sigma=2.0$) is ranked slightly below g1 ($\mu=3.0$, $\sigma=0.8$). g3's 23% 1-star share (vs g1's 2.5%) pulls the estimated location down; the compensating 23% 5-star share cannot restore it under the proportional odds constraint.
Logit vs probit: logit AIC (2242.5) beats probit (2247.8) by 5.3 — marginal but consistent. Logit coefficients are 1.62–1.63× probit (slightly less than the $\pi/\sqrt{3}=1.81$ constant-variance prediction), reflecting the different tail shapes interacting with the data.
Key takeaway: polr correctly identifies which quality tier each movie belongs to, but systematically misranks movies within a tier when their audience dispersion differs. For movies that are genuinely controversial (high $\sigma$), polr underestimates their quality. The heteroscedastic model in ord_probit_python.ipynb resolves this by estimating $\sigma_g$ explicitly.
4. Movies Data¶
Dataset description¶
Source: Liddell, T.M. & Kruschke, J.K. (2018), supplementary data. Amazon Prime video ratings collected circa 2017–2018 for a cross-section of 36 titles available in the US catalogue.
File: MoviesData.csv — one row per movie, columns ID, Descrip, n1–n5
(integer counts of 1-star through 5-star ratings).
Content: 36 titles spanning theatrical films and Amazon Original series. The mix is intentionally varied — prestige dramas, genre thrillers, limited series, classic re-releases — covering the full range from niche arthouse (N < 600) to Amazon Original flagship shows (N > 100K).
| Characteristic | Value |
|---|---|
| Movies | 36 |
| Total ratings | 284,671 |
| N per movie: min / median / max | 501 / 1,271 / 111,232 |
| Star scale | 1–5 (ordinal) |
Note on N: The dataset is deliberately unbalanced. Amazon Originals (Sneaky Pete, Man in the
High Castle) have orders-of-magnitude more ratings than limited theatrical releases (Megan Leavey,
Bosch Season 3). This creates heterogeneous precision across movies — a feature the heteroscedastic
model in ord_probit_python.ipynb handles naturally through the posterior width of $\mu_g$.
L&K headline pair¶
Liddell & Kruschke's central illustration of the metric model's failure:
| Movie | N | Empirical mean | % 1–2 star | Counts (1→5) |
|---|---|---|---|---|
| Priceless | 745 | 4.412 | 11.9% | 67, 22, 22, 60, 574 |
| Doctor Thorne S1 | 21,079 | 4.420 | 5.0% | 422, 632, 1897, 4848, 13280 |
Means differ by just 0.008 stars, yet Priceless has 2.4× the low-star share. Priceless is bimodal: 77% five-star fans vs 12% one-two-star detractors — a genuinely divided audience. Doctor Thorne is right-skewed and consistent: 63% five-star, 5% low-star. A metric analysis would call them equivalent; the ordinal model separates them.
movies <- read.csv("MoviesData.csv")
counts_m <- as.matrix(movies[, c("n1","n2","n3","n4","n5")])
mnames <- movies$Descrip
G_m <- nrow(counts_m)
n_m <- as.integer(rowSums(counts_m))
emp_mean_m <- as.numeric((counts_m %*% seq_len(K)) / n_m)
cat(sprintf("Movies: %d Total: %s ratings\n",
G_m, format(sum(counts_m), big.mark = ",")))
cat(sprintf("N/movie: min=%d median=%d max=%s\n",
min(n_m), as.integer(median(n_m)), format(max(n_m), big.mark = ",")))
pair_idx <- which(mnames %in% c("Priceless",
"Julian Fellowes Presents Doctor Thorne Season 1"))
cat("\nL&K headline pair:\n")
for (g in pair_idx) {
pct <- 100 * sum(counts_m[g, 1:2]) / n_m[g]
cat(sprintf(" %-50s N=%6d mean=%.3f %%low=%.1f\n",
mnames[g], n_m[g], emp_mean_m[g], pct))
cat(" Counts:", counts_m[g, ], "\n")
}
Movies: 36 Total: 284,671 ratings N/movie: min=501 median=1271 max=111,232 L&K headline pair: Priceless N= 745 mean=4.412 %low=11.9 Counts: 67 22 22 60 574 Julian Fellowes Presents Doctor Thorne Season 1 N= 21079 mean=4.420 %low=5.0 Counts: 422 632 1897 4848 13280
# Full movie summary table
emp_mean_m <- as.numeric((counts_m %*% seq_len(K)) / n_m)
pct_low <- 100 * (counts_m[, 1L] + counts_m[, 2L]) / n_m
pct_top <- 100 * counts_m[, 5L] / n_m
cat(sprintf("%-3s %-50s %7s %5s %6s %6s\n",
"ID", "Movie", "N", "Mean", "%Low", "%5star"))
cat(paste(rep("-", 82), collapse=""), "\n")
for (g in seq_len(G_m)) {
cat(sprintf("%3d %-50s %7d %5.3f %6.1f %6.1f\n",
g, mnames[g], n_m[g], emp_mean_m[g], pct_low[g], pct_top[g]))
}
ID Movie N Mean %Low %5star ---------------------------------------------------------------------------------- 1 The Whole Truth 700 3.770 17.0 35.0 2 Priceless 745 4.412 11.9 77.0 3 Allied 846 3.970 16.0 48.0 4 The Infiltrator 4319 4.140 9.0 48.0 5 Miss Sloane 1201 4.071 20.0 66.0 6 Room 13897 4.490 4.0 66.0 7 The Sea of Trees 5418 3.590 23.0 34.0 8 Megan Leavey 1127 4.730 3.0 84.0 9 Silence 523 3.428 33.1 40.0 10 The Only Living Boy in New York 587 4.262 12.9 65.1 11 Brawl in Cell Block 99 713 3.628 28.1 47.0 12 Freedom Writers 1160 4.658 4.1 79.9 13 The Curious Case of Benjamin Button 1296 4.379 10.0 70.0 14 Creed 11993 4.330 6.0 57.0 15 Triple 9 543 3.630 21.0 35.0 16 Taking Chance 5793 4.810 1.0 86.0 17 Rise 1012 3.820 17.0 41.0 18 Revolutionary Road 501 3.643 25.0 41.1 19 Sneaky Pete Season 1 57281 4.650 3.0 76.0 20 Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 5.0 63.0 21 Lifted 1610 4.461 6.0 69.0 22 Bosch Season 3 1531 4.758 4.0 89.0 23 The Man in the High Castle Season 1 111232 4.480 6.0 69.0 24 Ladies in Lavender 774 4.421 7.0 67.1 25 Me Before You 4966 4.310 12.0 69.0 26 Hunted Season 1 3468 4.280 7.0 56.0 27 The Guardian Season 1 2081 4.599 2.0 71.0 28 Poldark Season 2 3298 4.740 4.0 87.0 29 The Inheritance 711 4.236 11.1 61.9 30 The Marvelous Mrs Maisel Season 1 15594 4.840 3.0 93.0 31 Babel 1247 3.200 36.0 29.0 32 Tin Star Season 1 1626 3.301 39.0 41.0 33 In Secret 509 3.472 23.0 28.1 34 Philip K. Dick's Electric Dreams - Season 1 822 3.779 27.0 56.0 35 The Choice 3546 4.290 12.0 67.0 36 Fortitude Season 2 922 3.858 24.1 56.0
y_m <- unlist(lapply(seq_len(G_m), function(g) rep(seq_len(K), counts_m[g, ])))
grp_m <- rep(seq_len(G_m), n_m)
df_m <- data.frame(y = ordered(y_m), g = factor(grp_m))
cat("Fitting polr probit...\n")
t0 <- proc.time()["elapsed"]
fit_mp <- polr(y ~ g, data = df_m, method = "probit", Hess = FALSE)
cat(sprintf(" %.1fs\n", proc.time()["elapsed"] - t0))
cat("Fitting polr logit...\n")
t0 <- proc.time()["elapsed"]
fit_ml <- polr(y ~ g, data = df_m, method = "logistic", Hess = FALSE)
cat(sprintf(" %.1fs\n", proc.time()["elapsed"] - t0))
cat("\nPolr probit thresholds:\n")
print(round(fit_mp$zeta, 4))
cat(sprintf("\nAIC probit=%d logit=%d\n",
round(AIC(fit_mp)), round(AIC(fit_ml))))
coef_mp <- c(0, coef(fit_mp))
coef_ml <- c(0, coef(fit_ml))
Fitting polr probit...
23.6s
Fitting polr logit...
20.5s
Polr probit thresholds:
1|2 2|3 3|4 4|5
-1.1592 -0.8207 -0.4193 0.2649
AIC probit=524182 logit=523799
5. Rankings¶
rank_metric <- rank(-emp_mean_m, ties.method = "first")
rank_probit <- rank(-coef_mp, ties.method = "first")
rank_logit <- rank(-coef_ml, ties.method = "first")
comp <- data.frame(
movie = mnames,
N = n_m,
emp_mean = round(emp_mean_m, 3),
beta_probit = round(coef_mp, 4),
r_metric = as.integer(rank_metric),
r_probit = as.integer(rank_probit),
r_logit = as.integer(rank_logit),
shift = as.integer(rank_metric) - as.integer(rank_probit)
)
cat("Full ranking by polr probit:\n")
print(comp[order(comp$r_probit),
c("movie","N","emp_mean","beta_probit","r_metric","r_probit","shift")],
row.names = FALSE)
cat("\nLargest disagreements metric vs probit (|shift| >= 3):\n")
big <- comp[abs(comp$shift) >= 3L, ]
print(big[order(-abs(big$shift)),
c("movie","N","emp_mean","beta_probit","r_metric","r_probit","shift")],
row.names = FALSE)
cat("\nHeadline pair:\n")
print(comp[pair_idx,
c("movie","N","emp_mean","beta_probit","r_metric","r_probit","shift")],
row.names = FALSE)
Full ranking by polr probit:
movie N emp_mean beta_probit
The Marvelous Mrs Maisel Season 1 15594 4.840 1.6319
Bosch Season 3 1531 4.758 1.3777
Taking Chance 5793 4.810 1.3567
Poldark Season 2 3298 4.740 1.3045
Megan Leavey 1127 4.730 1.2206
Freedom Writers 1160 4.658 1.0684
Sneaky Pete Season 1 57281 4.650 0.9935
The Guardian Season 1 2081 4.599 0.8727
Priceless 745 4.412 0.8112
The Man in the High Castle Season 1 111232 4.480 0.7559
Lifted 1610 4.461 0.7397
Room 13897 4.490 0.7249
The Curious Case of Benjamin Button 1296 4.379 0.6919
Ladies in Lavender 774 4.421 0.6876
Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 0.6395
Me Before You 4966 4.310 0.6285
The Choice 3546 4.290 0.5928
The Only Living Boy in New York 587 4.262 0.5545
Creed 11993 4.330 0.5163
The Inheritance 711 4.236 0.4999
Hunted Season 1 3468 4.280 0.4716
Miss Sloane 1201 4.071 0.4259
The Infiltrator 4319 4.140 0.3135
Allied 846 3.970 0.2108
Fortitude Season 2 922 3.858 0.2022
Philip K. Dick's Electric Dreams - Season 1 822 3.779 0.1466
Rise 1012 3.820 0.0653
The Whole Truth 700 3.770 0.0000
Brawl in Cell Block 99 713 3.628 -0.0191
Revolutionary Road 501 3.643 -0.0414
Triple 9 543 3.630 -0.0896
The Sea of Trees 5418 3.590 -0.1142
Silence 523 3.428 -0.1802
In Secret 509 3.472 -0.2182
Tin Star Season 1 1626 3.301 -0.2504
Babel 1247 3.200 -0.3761
r_metric r_probit shift
1 1 0
3 2 1
2 3 -1
4 4 0
5 5 0
6 6 0
7 7 0
8 8 0
14 9 5
10 10 0
11 11 0
9 12 -3
15 13 2
12 14 -2
13 15 -2
17 16 1
18 17 1
20 18 2
16 19 -3
21 20 1
19 21 -2
23 22 1
22 23 -1
24 24 0
25 25 0
27 26 1
26 27 -1
28 28 0
31 29 2
29 30 -1
30 31 -1
32 32 0
34 33 1
33 34 -1
35 35 0
36 36 0
Largest disagreements metric vs probit (|shift| >= 3):
movie N emp_mean beta_probit r_metric r_probit shift
Priceless 745 4.412 0.8112 14 9 5
Room 13897 4.490 0.7249 9 12 -3
Creed 11993 4.330 0.5163 16 19 -3
Headline pair:
movie N emp_mean beta_probit
Priceless 745 4.412 0.8112
Julian Fellowes Presents Doctor Thorne Season 1 21079 4.420 0.6395
r_metric r_probit shift
14 9 5
13 15 -2
Movies results (validated)¶
Full ranking by polr probit (top 10 and bottom 5 shown):
movie N mean beta r_met r_prob shift
The Marvelous Mrs Maisel S1 15594 4.840 1.632 1 1 0
Bosch S3 1531 4.758 1.378 3 2 +1
Taking Chance 5793 4.810 1.357 2 3 -1
Poldark S2 3298 4.740 1.305 4 4 0
Megan Leavey 1127 4.730 1.221 5 5 0
Freedom Writers 1160 4.658 1.068 6 6 0
Sneaky Pete S1 57281 4.650 0.994 7 7 0
The Guardian S1 2081 4.599 0.873 8 8 0
Priceless 745 4.412 0.811 14 9 +5 ← bimodal
Man in the High Castle S1 111232 4.480 0.756 10 10 0
...
Babel 1247 3.200 -0.376 36 36 0
Largest disagreements (|shift| >= 3):
Priceless 745 4.412 0.811 r_metric=14 r_probit= 9 shift=+5
Room 13897 4.490 0.725 r_metric= 9 r_probit=12 shift=-3
Creed 11993 4.330 0.516 r_metric=16 r_probit=19 shift=-3
L&K headline pair:
Priceless 745 mean=4.412 beta=0.811 r_metric=14 r_probit= 9 shift=+5
Doctor Thorne 21079 mean=4.420 beta=0.640 r_metric=13 r_probit=15 shift=-2
→ polr ranks Priceless ABOVE Doctor Thorne (rank 9 vs 15)
Analysis — movies¶
Overall agreement. Most movies are ranked consistently by the metric mean and polr probit. Twenty-six of 36 have shift = 0 or ±1; the top four (Mrs Maisel, Taking Chance, Bosch, Poldark) and the bottom five are identical across both rankings. The ordinal model adds the most information for movies with unusual rating distributions.
The headline pair — a counterintuitive result.
| Movie | N | Mean | % 1–2★ | % 5★ | polr β | r_metric | r_probit | shift |
|---|---|---|---|---|---|---|---|---|
| Priceless | 745 | 4.412 | 11.9% | 77.0% | 0.811 | 14 | 9 | +5 |
| Doctor Thorne S1 | 21,079 | 4.420 | 5.0% | 63.0% | 0.640 | 13 | 15 | −2 |
polr ranks Priceless above Doctor Thorne (rank 9 vs 15). Why?
The proportional odds model fits cumulative probabilities $P(y \le k)$ with a single location shift $\hat\beta_g$. For Priceless vs Doctor Thorne:
| Cumulative threshold | Priceless | Doctor Thorne | Verdict |
|---|---|---|---|
| $P(y \le 1)$ | 9.0% | 2.0% | Doctor Thorne better |
| $P(y \le 2)$ | 11.9% | 5.0% | Doctor Thorne better |
| $P(y \le 3)$ | 14.9% | 14.0% | Nearly identical |
| $P(y \le 4)$ | 23.0% | 37.0% | Priceless better |
The top threshold ($P(y \le 4)$) carries the most information. Priceless's 77% five-star share means only 23% of its ratings fall at or below 4 stars, compared to 37% for Doctor Thorne. This 14-percentage-point gap at the top outweighs the 7-point gap in low-star rates. The single-location model treats the 5-star concentration as the dominant signal — and concludes Priceless has higher "quality."
The homoscedastic model's fundamental error. It cannot distinguish:
- Priceless: fans love it (77% ★★★★★) but detractors hate it (12% ★–★★) — high $\mu$, high $\sigma$
- Doctor Thorne: almost everyone likes it (63% ★★★★★, only 5% low-star) — moderate $\mu$, low $\sigma$
To polr, both look like "high-quality films," but Priceless's more extreme 5-star share makes
it appear better. The heteroscedastic Gibbs sampler in ord_probit_python.ipynb estimates a
separate $\sigma_g$ per movie, revealing that Priceless ($\mu=7.67$, $\sigma=4.45$) has genuinely
high latent quality for fans but extreme polarisation, while Doctor Thorne ($\mu=5.04$, $\sigma=1.56$)
is reliably liked by nearly everyone. The choice between them depends on risk preference, not quality.
Room (shift = −3) and Creed (shift = −3). Both are downranked by polr relative to the metric mean. Room (N = 13,897, mean = 4.490) and Creed (N = 11,993, mean = 4.330) have rating distributions that concentrate around 4–5 stars without a dominant 5-star spike. A high empirical mean driven partly by 4-star ratings converts to a lower $\hat\beta$ than a film with the same mean but more 5-star mass.
Summary. MASS::polr does not solve the L&K problem — it still conflates quality and consistency,
just via cumulative probabilities rather than raw means. Both the metric model and polr rank Priceless
above Doctor Thorne. The heteroscedastic model in ord_probit_python.ipynb gives the full picture:
two separate parameters per movie, $\mu_g$ and $\sigma_g$, letting the viewer decide.
6. Save¶
write.csv(comp, "movies_comparison_polr.csv", row.names = FALSE)
cat("Saved: movies_comparison_polr.csv\n")
cat(sprintf("Rows: %d Cols: %s\n",
nrow(comp), paste(names(comp), collapse = ", ")))
Saved: movies_comparison_polr.csv Rows: 36 Cols: movie, N, emp_mean, beta_probit, r_metric, r_probit, r_logit, shift