How to perform ordinal LR with rms?

1

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.
This is the beginning of my journey, please help.
I come from medical background.

best regards,
Andrzej

2

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

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
4

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.

5

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 ?

6

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

7

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

8

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

Of course your solution works:

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