Binomial vs. proportional odds logistic regression

Question

I have some clinical data, blood lab values (continuous) and clinical scores (ordinal: 0, 1, ..., 5), and I'd like to model the score as a function of the blood values. I know I could perform proportional odds logistic regression, but I find it hard to interpret. I am not interested in the probabilities of transitions between the levels, I can't find a nice way of visualising it, and computing the p-values requires extra effort.

I thought I could treat the score as a binomial variable, the number of successes in five trials, where the probability of a success would depend on the blood lab value. I've tried it out and I've been pretty satisfied with the results:

Call:
glm(formula = cbind(score, score_5) ~ bloodValue, family = "binomial", 
    data = tb1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8426  -1.0013   0.0366   0.9358   3.5244  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.4733     0.3224  -7.670 1.71e-14 ***
bloodValue    2.2781     0.3059   7.447 9.56e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 224.30  on 69  degrees of freedom
Residual deviance: 143.22  on 68  degrees of freedom
AIC: 231.45

Number of Fisher Scoring iterations: 4

However, the predictions differ from those obtained using POLR:

Obviously, they cannot be both right and, since POLR exists and I've never seen binomial regression applied this way before, I suspect my approach is wrong. On the other hand, binomial model has lower errors, both MAE (0.94 vs. 1.71) and MSE (1.8 vs. 4.6), so it cannot be that wrong. I'd appreciate it if someone could point out the mistake in my thinking. Or are both approaches legitimate, only under different assumptions? What could they be?

Below is my code for reproducing the results:

library(tidyverse)

tb1 = structure(list(
  bloodValue = c(0.81, 0.43, 0.6, 1.74, 1.33, 0.43, 
    0.64, 1.18, 1.06, 0, 1.52, 1.03, 0.86, 1.23, 1.45, 0.76, 0.93, 
    0.84, 2.6, 0.77, 1.23, 0, 1.23, 0.83, 1.1, 0.93, 0.69, 0, 1.74, 
    1.68, 0.69, 1, 0.89, 2.24, 0.59, 1.16, 0, 1.06, 0, 1.21, 1.62, 
    1.31, 1.26, 0, 1.09, 0.76, 0, 0.81, 1.27, 0.63, 1.19, 1.39, 0.44, 
    1.4, 0.85, 1.68, 0.44, 1.51, 0, 1.12, 0.44, 0.76, 0.69, 0.53, 
    1.12, 0.35, 0.73, 1.12, 0.81, 0.88),
  score = c(2, 1, 0, 5, 3, 
    0, 2, 2, 4, 0, 4, 3, 3, 3, 3, 0, 1, 3, 4, 2, 3, 0, 5, 5, 4, 0, 
    0, 0, 4, 2, 5, 1, 2, 4, 0, 3, 0, 2, 0, 4, 0, 4, 4, 0, 4, 3, 0, 
    0, 5, 3, 3, 4, 0, 4, 3, 4, 1, 1, 0, 0, 1, 1, 0, 0, 3, 0, 3, 5, 
    3, 1)),
  row.names = c(NA, -70L),
  class = c("tbl_df", "tbl", "data.frame"))

tb1 = tb1 %>% mutate(score_5 = 5-score)
bm1 = glm(cbind(score, score_5) ~ bloodValue, data=tb1, family = "binomial")
summary(bm1)
tb1 %>%
  ggplot(aes(bloodValue, score/5, succ=score, fail=score_5)) +
  geom_jitter(width=0, height=.01, alpha=.5) +
  geom_smooth(
    method = "glm",
    method.args = list(family="binomial"),
    formula=cbind(succ, fail) ~ x
  )

po1 = MASS::polr(factor(score, levels=seq(0, 5)) ~ bloodValue, data=tb1)
pred.polr = as.numeric(predict(po1))-1
pred.bm   = round(5*predict(bm1, type="response"))
tibble(truth = tb1$score, polr=pred.polr, binom=pred.bm) %>%
  gather("method", "predicted", 2:3) %>%
  ggplot(aes(truth, predicted, colour=method)) +
    geom_jitter(width=.1, height=.1, alpha=.5, size=2) +
  geom_smooth(method="lm")

# MAE
sum(abs(as.numeric(pred.polr)-1 - tb1$score)) / nrow(tb1)
sum(abs(pred.bm - tb1$score)) / nrow(tb1)

# MSE
sum((as.numeric(pred.polr)-1 - tb1$score)^2)  / nrow(tb1)
sum((pred.bm - tb1$score)^2) / nrow(tb1)

P.S. I've read the following questions:

• Which regression model to use for a probability as dependent variable?
• Regression for an outcome (ratio or fraction) between 0 and 1
• How to do logistic regression in R when outcome is fractional (a ratio of two counts)?
• Which glm family to use for ordinal DV?

but they don't seem to answer my question.

Answer 1

With the binomial regression, you are assuming that when you see level 3, you had 5 independent binary yes/no outcomes that each could have been yes or no independent (conditional on the observed covariate values for that subject) from each other. You are assuming that the probabilities of each binary outcome are the same (conditional on the observed covariates values, which are identical for a subject), which may mismatch your observed data rather badly. This model also mismatches the data generating process: you can only get to level 3, if you are already above levels 0, 1 and 2, but if you did not get from level 3 to 4, then whether you hypothetically would have mde it from 4 to 5 does not matter. In a Bayesian setting something like posterior predictive checks might be able to illustrate just how much this does or does not mismatch the actually data distribution.

So, on a whole, that model produces a smaller standard error by making lots and lots of assumptions and these assumptions seem to be rather questionable to me. A small standard error certainly does not mean anything is good about a model (counter example: "my model" always estimates 42 with a standard error of 0).

Proportional odds logistic regression, of course, also makes assumptions. These are less strong than the ones spelled out above and can, if truly necessary, be relaxed further. It would be a rather obvious starting point though.

Perhaps one additional perspective would be to do cross-validation and to see which of the two models predicts better on the out-of-fold data (admittedly not quite the same task as inference), perhaps simply in terms of mean absolute error (0=correct category predicted, 1=missed by one category, 2=missed by two categories etc.).