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.
score
forbloodValue = 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. $\endgroup$