Different outcome groups with MASS:polr and rms::lrm

[mhoenicka]

Short version:
when using marginaleffects::avg_slopes() or avg_comparisons() on MASS::polr() and rms::lrm() fits of identical ordinal logistic models, I seem to get one estimate per predictor and outcome variable level (4 in our case) from the polr() fit but only three estimates per predictor from the lrm() fit. Am I doing something wrong, or are my expectations wrong?

Full description:
We model the severity of postoperative acute kidney injury as a function of several preoperative and intraoperative predictors. The outcome parameter AKI_Kreatinin_KDIGO is an ordered factor. 0 (zero) denotes absence of acute kidney injury. 1 through 3 indicate the three stages of acute kidney injury with increasing severity.

str(rhe_amelia_complete$AKI_Kreatinin_KDIGO)
 Ord.factor w/ 4 levels "0"<"1"<"2"<"3": 1 2 2 1 1 2 2 1 1 1 ...

The model formula looks like this, with a total of 26 binary, categorical, and numeric factors:

f<-"AKI_Kreatinin_KDIGO~AGE+ARTFLOW_1.5_incl_to_2.75_excl_cont+..."

I fit ordinal logistic regression models with MASS::polr() and rms::lrm() using the same formula.

olmodel<-polr(as.formula(f),data=rhe_amelia_complete,Hess=TRUE) # method="logistic" (default)

dd <- datadist(rhe_amelia_complete)
options(datadist='dd')
olmodel_lrm<-lrm(as.formula(f),data=rhe_amelia_complete,x=TRUE,y=TRUE)

I calculate the slopes of both models and display the estimates:

olslopes<-avg_slopes(olmodel)
olslopes$estimate

olslopes_lrm<-avg_slopes(olmodel_lrm)
olslopes_lrm$estimate
 
olslopes$estimate is a vector of length 104. I figured these are the slopes of our 26 predictors for each of the four levels of the outcome parameter, for example:


                 0             1             2             3
AGE  -0.0035496610  0.0026187686  5.994254e-04  3.314670e-04


olslopes_lrm$estimate is a vector of length 78. This I interpreted as three columns for 26 predictors:

                0             1             2
AGE  0.0035490765  0.0009307322  3.314111e-04

The values of columns 0 (polr) and 0 (lrm) as well as of columns 3 (polr) and 2 (lrm) of these outputs are identical within a reasonable precision. However, lrm() estimates in column 1 are somewhere inbetween the polr() estimates in columns 1 and 2. Why is that, and how do I interpret the estimates of the lrm()-derived slopes in column 1?

Thanks
Markus

[vincentarelbundock]

Thanks for the report.

marginaleffects does not determine the groupings itself, and completely defers to the upstream modeling package to do so. In this case, you'll note that if we call the base R predict() function on models from these two packages, we'll get matrices of different shapes, and with different labels.

avg_slopes() is completely agnostic with respect to the meaning of the quantities generated by predict(). It gives you the average slope for each column, with respect to the predictor of interest. So the interpretation depends on what the upstream packages have decided to put in those columns, and you should refer to the respective packages documentation for the predict() methods.

Here's a minimal reproducible and self-contained example:

library(rms)
library(MASS)
mod1 <- lrm(factor(carb) ~ hp + mpg, mtcars)
mod2 <- polr(factor(carb) ~ hp + mpg, mtcars)

predict(mod1, type = "fitted") |> head()
#>                        y>=2      y>=3      y>=4        y>=6         y>=8
#> Mazda RX4         0.7367057 0.2061653 0.1168384 0.004276733 0.0013229833
#> Mazda RX4 Wag     0.7367057 0.2061653 0.1168384 0.004276733 0.0013229833
#> Datsun 710        0.6219511 0.1324721 0.0721720 0.002519033 0.0007782996
#> Hornet 4 Drive    0.7357330 0.2053468 0.1163226 0.004255457 0.0013163820
#> Hornet Sportabout 0.9526702 0.6513585 0.4876268 0.029971840 0.0094398634
#> Valiant           0.7141569 0.1882451 0.1056493 0.003820538 0.0011814879

predict(mod2, type = "probs") |> head()
#>                            1         2          3          4           6
#> Mazda RX4         0.26306260 0.5310115 0.08923438 0.11242212 0.002949111
#> Mazda RX4 Wag     0.26306260 0.5310115 0.08923438 0.11242212 0.002949111
#> Datsun 710        0.37783157 0.4898954 0.06021444 0.06954455 0.001737523
#> Hornet 4 Drive    0.26403733 0.5308567 0.08893103 0.11192677 0.002934419
#> Hornet Sportabout 0.04723621 0.3015372 0.16373520 0.45754371 0.020518345
#> Valiant           0.28559232 0.5263809 0.08250878 0.10170390 0.002634991
#>                              8
#> Mazda RX4         0.0013203269
#> Mazda RX4 Wag     0.0013203269
#> Datsun 710        0.0007765262
#> Hornet 4 Drive    0.0013137215
#> Hornet Sportabout 0.0094292987
#> Valiant           0.0011791551

[mhoenicka]

Thanks Vincent for the explanation. The output makes sense now.
Markus