# How to perform ordinal LR with rms?

Hi,
I am beginner in LR subjects. I would like to perform ordinal logistic regression with this data:

``````ordinal_data <- structure(list(ODIPain = structure(c(3L, 1L, 2L, 3L, 4L, 3L,
4L, 2L, 3L, 5L, 2L, 5L, 5L, 6L, 2L, 3L, 1L, 2L, 3L, 3L, 1L, 3L,
3L, 2L, 2L, 5L, 5L, 2L, 5L, 3L, 5L, 1L, 3L, 3L, 3L, 1L, 5L, 3L,
5L, 1L, 1L, 2L, 1L, 2L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 4L, 6L, 4L,
3L, 3L, 3L, 3L, 1L, 4L, 5L, 4L, 3L, 3L, 1L, 3L, 1L, 4L, 3L, 3L,
2L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 2L, 2L, 1L, 3L, 4L, 4L, 3L,
2L, 2L, 2L, 2L, 2L, 1L, 1L, 3L, 1L, 3L, 1L, 3L, 4L, 4L, 3L, 3L,
1L, 2L, 3L, 3L, 3L, 3L, 5L, 2L, 2L), levels = c("[0] No pain",
"[1] Very mild pain", "[2] Moderate pain", "[3] Fairly severe pain",
"[4] Very severe pain", "[5] Worst imaginable pain"), class = c("ordered",
"factor")), Arm = structure(c(2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L,
1L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L,
2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L,
1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L,
1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 2L,
1L, 1L, 2L, 1L, 2L, 2L, 2L), levels = c("A", "B"), class = "factor"),
Age_centered = c(-6.15315315315316, 12.8468468468468, -9.15315315315316,
14.8468468468468, 12.8468468468468, 2.84684684684684, -10.1531531531532,
-18.1531531531532, -1.15315315315316, 8.84684684684684, -17.1531531531532,
13.8468468468468, 9.84684684684684, 17.8468468468468, -19.1531531531532,
-7.15315315315316, -10.1531531531532, -19.1531531531532,
-7.15315315315316, 0.846846846846844, -17.1531531531532,
5.84684684684684, -25.1531531531532, -1.15315315315316, -15.1531531531532,
4.84684684684684, 1.84684684684684, 12.8468468468468, -11.1531531531532,
5.84684684684684, -6.15315315315316, -0.153153153153156,
20.8468468468468, 5.84684684684684, -0.153153153153156, 12.8468468468468,
-19.1531531531532, -11.1531531531532, 1.84684684684684, 0.846846846846844,
-21.1531531531532, 9.84684684684684, 15.8468468468468, 14.8468468468468,
-12.1531531531532, -11.1531531531532, -9.15315315315316,
5.84684684684684, -4.15315315315316, 12.8468468468468, 1.84684684684684,
-7.15315315315316, -3.15315315315316, 7.84684684684684, 0.846846846846844,
-4.15315315315316, 5.84684684684684, -0.153153153153156,
1.84684684684684, -7.15315315315316, 1.84684684684684, -9.15315315315316,
6.84684684684684, 9.84684684684684, 17.8468468468468, 5.84684684684684,
9.84684684684684, -10.1531531531532, -5.15315315315316, 18.8468468468468,
21.8468468468468, -0.153153153153156, 2.84684684684684, -8.15315315315316,
-5.15315315315316, 5.84684684684684, 2.84684684684684, -15.1531531531532,
2.84684684684684, 25.8468468468468, -11.1531531531532, 27.8468468468468,
2.84684684684684, 20.8468468468468, -0.153153153153156, -2.15315315315316,
12.8468468468468, -0.153153153153156, 0.846846846846844,
11.8468468468468, -8.15315315315316, 3.84684684684684, 22.8468468468468,
5.84684684684684, 12.8468468468468, 4.84684684684684, 11.8468468468468,
-5.15315315315316, -17.1531531531532, -7.15315315315316,
-16.1531531531532, 0.846846846846844, -13.1531531531532,
-13.1531531531532, -19.1531531531532, -15.1531531531532,
-11.1531531531532, -4.15315315315316, -0.153153153153156,
-4.15315315315316, -7.15315315315316)), row.names = c(NA,
-111L), class = "data.frame")
``````

So far I usually work with SPSS and trying to use R as well.
I hope my question is not trivial in this respected forum.
I would like to use rms package and perform ordinal logistic regression.
I would like to get an output in the form of data frame not a lists if possible.
I also work a lot with Likert data and want to learn the proper way of analysing it.
I would be very grateful for any help here.

What I have tried so far:

``````ormys <- orm(ODIPain ~ Age_centered + Arm, data = ordinal_data)
summary_orm <- summary(ormys)
``````

I want to learn how to use: rms::lrtest() function and how to interpret the results.
I come from medical background.

best regards,
Andrzej

I assume with LR you mean likelihood-ratio test? You would simply fit two nested models and apply `lrtest`. For example: For a LR test for the variable `Arm`, you could do the following:

``````ormys <- orm(ODIPain ~ Age_centered + Arm, data = ordinal_data)
ormys0 <- orm(ODIPain ~ Age_centered, data = ordinal_data)
lrtest(ormys0, ormys)
``````

Thanks, by LR I meant Logistic Regression, but likelihood-ratio test you mentioned is also very useful.
I wonder why this is not working:

``````model <- lrm(ODIPain ~ Age_centered + Arm, data = ordinal_data)

model_no_proportional <- lrm(ODIPain ~ Age_centered + Arm + rcs(Age_centered, 2), data = ordinal_data)

lrt_result <- rms::lrtest(model_no_proportional, model)

``````

and this is not working as well:

``````library(MASS)
model <- polr(ODIPain ~ Age_centered + Arm, data = ordinal_data, Hess = TRUE)
model_no_proportional <- polr(ODIPain ~ Age_centered + Arm + I(Age_centered^2),
data = ordinal_data, Hess = TRUE)
lrt_result <-rms::lrtest(model_no_proportional, model)

Error in command 'if (df1 == df2) stop("models are not nested")':
value missing where TRUE/FALSE is required
``````

`lrtest` is from the `rms` package and doesnâ€™t support `polr` models, afaik. To perform a likelihood ratio test with `polr`, you would use `anova(model_no_proportional, model, test = "Chisq")`, I presume.

Also, please note that in the `lrm` model, you fit restricted cubic splines but get an error, because the number of knots needs to be 3 or higher. In the `polr` model, you fit a quadratic polynomial, which is not the same as a spline. To get equal results, you would do the following:

``````# Using rms
model <- lrm(ODIPain ~ Age_centered + Arm, data = ordinal_data)
model_no_proportional <- lrm(ODIPain ~ Arm + poly(Age_centered, 2), data = ordinal_data)
rms::lrtest(model_no_proportional, model)

# With polr
model <- polr(ODIPain ~ Age_centered + Arm, data = ordinal_data, Hess = TRUE)
model_no_proportional <- polr(ODIPain ~ Age_centered + Arm + I(Age_centered^2),
data = ordinal_data, Hess = TRUE)
anova(model_no_proportional, model, test = "Chisq")
``````

The p-values are identical.

Thank you very much indeed, with every answer I am learning a lot.
So what would be the proper way of analysing my ordinal_data ?

Is it better to use orm() or lrm() , and what does this test say:

``````lrtest(ormys0, ormys)
``````

Is it also a way to do some plotting when using ordinal logistic regression ?

This discussion is on the wrong page. It should have been on https://discourse.datamethods.org/t/rms-ordinal-logistic-regression

So what do we do ? Can you please answer my questions or transfer this topic to as proper place and answer my question ?

Thank you , I want to fully understand why this attempts are failing:

``````model_no_proportional <- lrm(ODIPain ~ Age_centered + Arm + rcs(Age_centered, 3), data = ordinal_data)

knots <- quantile(ordinal_data\$Age_centered, c(0.25, 0.5, 0.75))
model_no_proportional <- lrm(ODIPain ~ Age_centered + Arm + rcs(Age_centered, knots = knots), data = ordinal_data)

number of knots in rcs defaulting to 5
Error in 'if (asc[i] == 8) next' command:
value missing where TRUE/FALSE is required
``````

``````model_no_proportional <- lrm(ODIPain ~ Arm + poly(Age_centered, 2), data = ordinal_data)