# RMS Ordinal Logistic Regression

**URL:** <https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815>\
**Category:** models\
**Tags:** ordinal\
**Created:** [September 17, 2021, 12:47am UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815 "2021-09-17T00:47:28Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [January 23, 2023, 5:29pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/21 "2023-01-23T17:29:52Z")

</div>

If there is only one non-absorbing state then the previous state is always the same (the non-absorbing one) and is redundant. See if that applies here.

---

<div class="post-metadata">

**Author:** ![Elias\_Eythorsson](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/elias_eythorsson/32/2524_2.png) [@Elias\_Eythorsson](https://discourse.datamethods.org/u/Elias_Eythorsson)\
**Post date:** [May 29, 2023, 5:37pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/22 "2023-05-29T17:37:43Z")

</div>

I have a proportional odds model with four outcome levels: “0-4”, “5-9”, “10-14” and “\>15”. The overall C statistic is 0.766 as given by the `print()` function for the model object. I obtain the same estimate of the C statistic if I use the `validate()` function and transform Summers’ Dxy to the C statistic. If I however use the `val.prob()` function for any of the individual levels of the outcome, I obtain a new estimate of the C-statistic for each level: 0.754 for “5-9”, 0.863 for “10-14” and 0.907 for “\>15”.

Is there anything inappropriate about reporting outcome specific C-statistics if the proportional odds model will predominantly be used to predict one of the four levels in clinical practice? I would do this alongside the global C-statistic.

---

<div class="post-metadata">

**Author:** ![Negin\_Mahmoudi](https://discourse.datamethods.org/letter_avatar_proxy/v4/letter/n/5e9695/32.png) [@Negin\_Mahmoudi](https://discourse.datamethods.org/u/Negin_Mahmoudi)\
**Post date:** [November 30, 2023, 10:48pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/23 "2023-11-30T22:48:51Z")

</div>

I am currently engaged in a research project with an ordinal regression model predicting an ordered response variable within the range of 1-9, I’ve been exploring various methods to estimate the variance of my model’s output (predicted probabilities associated with each class) During this exploration, I came across following methods:

1. Bootstrapping
2. Delta Method
3. Bayesian methods like Markov Chain Monte Carlo (MCMC)

Given that Bayesian methods, such as MCMC, are commonly employed with small sample sizes and Bootstrapping requires generating multiple bootstrap samples, I believe the Delta Method would be a more suitable choice for my specific purpose. However after reading this paper : [Using the Delta Method to Construct Confidence Intervals for Predicted Probabilities, Rates, and Discrete Changes](https://jslsoc.sitehost.iu.edu/stata/ci_computations/spost_deltaci.pdf), I wonder:

1. What factors should guide my decision on the sample size for optimal results?
2. What considerations should I take into account to avoid violating the assumption of asymptotic normality?
3. Are there alternative methods for estimating the uncertainty associated with each predicted class?

P.S. Heartfelt gratitude to Dr. Harrell for providing us with this exceptional platform! @f2harrell

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [December 1, 2023, 1:48pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/24 "2023-12-01T13:48:26Z")

</div>

> [@Elias\_Eythorsson](#):
>
> Is there anything inappropriate about reporting outcome specific C-statistics if the proportional odds model will predominantly be used to predict one of the four levels in clinical practice? I would do this alongside the global C-statistic.

There is nothing wrong with getting c for different levels of Y other than a slight loss of precision. Note that `val.prob` is only for use on a sample that was not used in model development.

> [@Negin\_Mahmoudi](#):
>
> I believe the Delta Method would be a more suitable choice for my specific purpose

I don’t think the delta method has sufficient accuracy in many situations. And always remember that the delta method assumes that the confidence interval is symmetric. That does’t work well for probabilities outside of [0.3, 0.7]. Since Bayesian methods give exact intervals for any transformation of model parameters, it’s harder to get motivated to spend much time on frequentist approximations.

If you stick with a frequentist approach, I’d use the bootstrap instead of the delta method. The problem is which bootstrap from among the many choices. You’ll find differences in confidence interval accuracy among the various bootstraps.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 16, 2024, 7:42pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/25 "2024-02-16T19:42:10Z")

</div>

Hi Professor @f2harrell,

I have a question about getting predicted _means_ and _confidence intervals_ when using the `predict.lrm()` function.

Note: **fit** is my fitted `lrm()` model and **bfit** is the bootstraped model obtained using `bootcov()` (e.g., `bfit <- bootcov(fit, B = 500, coef.reps = T)`).

Using `ggplot(Predict())` seems straightforward:

`ggplot(Predict(bfit, fun = 'mean'))`

Building a nomogram is also working fine:

`nomogram(fit, fun = Mean(fit), funlabel = 'Mean')`

But I’m struggling with the code for getting the predicted _means_ and _confidence intervals_ on new data using the `predict.lrm()` function. I don’t know what the argument `codes` and the ‘design matrix X’ mean. Could you kindly give me an idea of what this code would look like?

Thank you.

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [February 18, 2024, 12:27pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/26 "2024-02-18T12:27:57Z")

</div>

I think that `Predict` or `predict` will give you approximate CIs. To get better CIs that use the bootstrap reps you may need to use `contrast`.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 18, 2024, 8:38pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/27 "2024-02-18T20:38:07Z")

</div>

But if I use:

`Predict(boot, fun = 'mean')`,

then this will give me predictions on the mean scale and the CIs will be nonparametric percentile, right?

If I want predictions on the probability scale (e.g., _Pr{Y\>=j}_, where _j_ is the middle intercept), then I should use:

`fun <- function(x) plogis(x-fit$coef[1] + fit$coef[middle intercept])`,

and then:

`Predict(boot, fun = 'fun')`.

Is that correct?

Thank you!

---

<div class="post-metadata">

**Author:** ![Elias\_Eythorsson](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/elias_eythorsson/32/2524_2.png) [@Elias\_Eythorsson](https://discourse.datamethods.org/u/Elias_Eythorsson)\
**Post date:** [February 18, 2024, 10:31pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/28 "2024-02-18T22:31:15Z")

</div>

Predict (boot, fun = plogis, kint = j)

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 19, 2024, 10:51am UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/29 "2024-02-19T10:51:42Z")

</div>

Hi @Elias_Eythorsson, thank you for your suggestion!

I run both,

`fun <- function(x) plogis(x-fit$coef[1] + fit$coef[middle intercept])`,  
then  
`Predict(boot, fun = 'fun')`,  
and  
`Predict (boot, fun = plogis, kint = j)`

I obtained the same point estimates for the predicted probabilities, but the lower and upper limits of the confidence intervals were different.

I guess the former (`fun = 'fun'`) is using the standard estimated coefficients while the latter (`fun = plogis`) is using the bootstraped coefficients.

Does that make sense?

Thank you!

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [February 20, 2024, 12:33pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/30 "2024-02-20T12:33:18Z")

</div>

Also look at `ExProb`.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 20, 2024, 1:12pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/31 "2024-02-20T13:12:53Z")

</div>

How?  
Like this?

`p <- Predict(boot, x, fun = 'ExProb')`

But my model was fit using `lrm()`, not `orm()`.

What about what I wrote in the post above? Is it correct?

Thank you

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 21, 2024, 8:46pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/32 "2024-02-21T20:46:39Z")

</div>

Calculating contrasts for a `lrm()` model:

`m <- Mean(fit)`

`c <- contrast(fit, list(x = quantile(data$x, 0.75)),`  
`list(x = quantile(data$x, 0.25)))`

`print(c, fun = m)`

The relationship between X and Y is negative, in that going from the first to the third quartile of X results in a decrease in Y (as shown by the partial effect plot - test of association is significant: Chi-square of 13 with P \< 0.001). However, that’s not what the contrast is showing, as the effect is positive. Also, going from first to third quartile or vice-versa results in different effect estimates with `contrast()`, instead of simply changing the signal of the effect (i.e., + and -). Also, there is no SE for the computed contrast (the output from `print(c, fun = m)` shows NA for the SE).

I wonder wheter this has something to do with the `Mean()` function?

@f2harrell?

I would appreciate any help on this.

Thanks

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [February 28, 2024, 3:23pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/33 "2024-02-28T15:23:41Z")

</div>

@f2harrell, please, could you kindly give me some guidance on the above?

Thank you

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [February 29, 2024, 1:17pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/34 "2024-02-29T13:17:13Z")

</div>

I hope others can. I have a heavy teaching load right now.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [March 3, 2024, 12:19pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/35 "2024-03-03T12:19:32Z")

</div>

Hi @Elias_Eythorsson,  
Do you have experience using the `contrast()` function for a `lrm()` model?  
Any comments on the issue above?

Thank you

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [June 12, 2024, 12:48pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/36 "2024-06-12T12:48:11Z")

</div>

Dear all,

I’m using the TBI dataset (complete case analysis) from @Ewout_Steyerberg’s book, _Clinical Prediction Models - 2nd edition_, to fit an _Ordinal (PO) Logistic Regression_ model, using `rms::lrm()`. **My question is at the end of this post**.

Outcome variable is _Glasgow Outcome Scale (GOS, 1-5)_.  
Below is the `Hmisc::describe()` output:  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/e/eee9e801a2628d02278ff67ddeba682b7e552c97.png)

I have collapsed the categories from Motor\_Score and CT\_Class:  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/b/bcb2bfadac31185acd01ecc476dc5f0d4dbd62ee.png)

The model fitted for predicting GOS\>=4 at 6 months was the following:

`lrm(GOS ~ Age + Motor_Score + Pupillary_Reac + Hypoxia + Hypotension + CT_Class + TSAH, data = data, x = T, y = T)`

Here’s the output from `rms::validate()`:  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/c/ca881316eb0009cffb35bc3bd27477fed186923d.png)

And from `rms::calibrate()`:  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/e/e842283dbc304a7d1b7f2a7bcddb225666e562c6.png)

For comparison purposes, I’ve also fitted a _Binary Logistic Regression_ model (_GOS 1,2,3_ vs _GOS 4,5_), which surprisingly (for me at least) seemed to outperform the _Ordinal_ model:  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/9/93dfcbeb721e372a9206b6d9510bdb716ed52b3a.png)

![image](https://discourse.datamethods.org/uploads/default/original/2X/d/d3043abc0ec0ebf3322b58760a72565ea4e017bf.png)

**Could the reason for this be the fact that there’s evidence for non-PO in my _Ordinal_ model?**  
 ![image](https://discourse.datamethods.org/uploads/default/original/2X/b/b517361bedfec97cd3a3faa00be6b5df51a59f74.png)

Thank you!

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 12, 2024, 3:12pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/37 "2024-06-12T15:12:16Z")

</div>

There is evidence for non-PO but the degree of non-PO is only “semi-consequential”. The larger issue is that when you compare quantities such as c and D\_{xy}, predicting an ordinal response is a much more difficult task than predicting a binary outcome.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [June 12, 2024, 6:07pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/38 "2024-06-12T18:07:02Z")

</div>

Thank you for the reply Professor!  
Would it be ok then for me to proceed with the ordinal model?  
I ask not only because of the evidence for non-PO but also because of the _mean absolute prediction error_, which is 5 times smaller for the binary logistic model.

---

<div class="post-metadata">

**Author:** ![f2harrell](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/f2harrell/32/165_2.png) [@f2harrell](https://discourse.datamethods.org/u/f2harrell)\
**Post date:** [June 13, 2024, 12:51pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/39 "2024-06-13T12:51:31Z")

</div>

This is an extremely difficult decision. Both calibration curves are good. But AIC favors the multinomial model for prediction accuracy. I would slightly lean towards the multinomial model but not by much.

---

<div class="post-metadata">

**Author:** ![FurlanLeo](https://discourse.datamethods.org/user_avatar/discourse.datamethods.org/furlanleo/32/3812_2.png) [@FurlanLeo](https://discourse.datamethods.org/u/FurlanLeo)\
**Post date:** [June 13, 2024, 1:35pm UTC](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815/40 "2024-06-13T13:35:23Z")

</div>

I was considering only the Ordinal and Binary models, and completely forgot about the Multinomial one (despite comparing it against the Ordinal). The calibration curves you see above are from the Ordinal and Binary models. From your last comment, it seems that you are not favouring the Binary model. Any particular reason? Many thanks!

[Previous page](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815.md?page=1)

[Next page](https://discourse.datamethods.org/t/rms-ordinal-logistic-regression/4815.md?page=3)