Binomial vs. proportional odds logistic regression

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:

but they don't seem to answer my question.

• Ordinal regression has a higher likelihood of fitting the data, and provides a bigger variety of types of predictions. Don't interpret the proportional odds model in terms of probabilities of transitions but rather in terms of covariate-specific cumulative probabilities. It's not so hard to interpret. If you have a specific question about interpretation please state it. Sep 23 '21 at 11:42
• @FrankHarrell Can you maybe elaborate (in an answer?) why POLR should be more likely to fit the data? For my concrete dataset, binomial regression produced lower errors. Regarding interpretation: Say I want to predict the score for bloodValue = 1. POLR gives me: 0.18 0.14 0.14 0.30 0.17 0.04 for the levels 0-5. Now, I know that score 3 is the most probable, but it is confusing that 0 is more probable than 1, 2, and 4. With binomial regression I can easily see the expected value and the standard error of the score. Sep 23 '21 at 13:47

1 Answer

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.).

• Thanks, I got that. I believe it boils down to the question of the generating process. For school grades, calculated from the number of correct answers on an exam, it is certainly possible that the questions are independent. I may answer any three questions out of five correctly and get a "C". In problems I encounter in practice, we often don't know the true process, we just observe certain phenomena and try to summarise them in a score. Sep 24 '21 at 7:35
• Besides the generating process, plot the cumulative distributions from (1) the data, (2) that dictated by the binomial distribution, and (3) that dictated by the proportional odds model. See how (2) and (3) match (1) or not. Sep 24 '21 at 11:35
• Yes, in the exam example assuming multiple binomials for each question could be a reasonable model (especially if the final score is number of correctly answered questions). For a clinical score from Dead -> in ICU -> in hospital -> is at home with symptoms -> has no symptoms that seems very dubious, and something that treats these as general ordered categories makes much more sense to me. Sep 24 '21 at 11:40
• Thanks, @Björn, it makes sense. So, for qSOFA (mdcalc.com/qsofa-quick-sofa-score-sepsis), binomial may be justified, but not for GCS (mdcalc.com/glasgow-coma-scale-score-gcs): You may have low blood pressure independently of you mental state, but you cannot be able to spontaneously open your eyes if you cannot open them in response to pain. Sep 29 '21 at 14:00
• I've got no idea whether qSOFA will be well described by a binomial model. I'm slightly skeptical, because my default assumption is that ordinal models tend to be a good idea, but don't base a decision on my gut feeling. E.g. before trying a new method on a study, I like to try it out on some previous data (if I can get any) to see whether it behaves sensibly / manages to find parameters that can even describe the data well. Sep 30 '21 at 6:54