Interpretation of marginal effects of ordinal predictors

[mhoenicka]

I compute an ordinal logistic regression model with MASS::polr(). One of the predictors is an ordinal parameter AGE with four levels. I understand that the polr-generated model contains three estimates (one less than the number of levels), with AGE.L, AGE.Q, and AGE.C denoting linear, quadratic, and cubic terms, respectively. Now I use avg_comparison() to compute the marginal effects of a +1 unit change, which again returns linear, quadratic, and cubic terms of AGE.
Now I want to quantify the total marginal effect MEt of AGE. Is the math as simple as this:
MEt=ME.L+(ME.Q*u)+(ME.C*u^2), which equals ME.L+ME.Q+ME.C for u=1?
Thanks
Markus

[vincentarelbundock]

The problem is likely that you have hard-coded the polynomials in the dataset before fitting your model. Polynomials should be specified in the fitting formula.

Please read the section on polynomials: https://marginaleffects.com/chapters/interactions.html#polynomial-regression

Also, note that the above is speculation, because you did not supply sufficient information to answer. Please read this page on how to ask good questions about software:

https://marginaleffects.com/bonus/help.html#getting-help

[mhoenicka]

My apologies for assuming that this is a rather generic question. I've created a self-contained example to demonstrate the issue. I have to admit that the last step of the example is merely an assumption.

`library(MASS)
library(marginaleffects)

# create dataset
# predictor: age classes 1 through 4: lo<lomed<himed<hi
AGE_ord<-sample(1:4, 100, replace=TRUE)
AGE_ord<-factor(AGE_ord,ordered=TRUE)

# outcome: acute kidney injury stages 0 through 3 with 0<1<2<3
AKI_ord<-sample(0:3, 100, replace=TRUE)
AKI_ord<-factor(AKI_ord,ordered=TRUE)

my_df<-data.frame(AGE_ord,AKI_ord)

# run ordered logistic regression
olmod<-polr(AKI_ord~AGE_ord,data=my_df,Hess=TRUE)
summary(olmod)

# compute marginal effects
olcomp<-avg_comparisons(olmod)
olcomp$estimate
est<-matrix(olcomp$estimate,ncol=4,dimnames=list(c("L","Q","C"),c("lo","lomed","himed","hi")))
est`

I model the stage of acute kidney injury (0 through 3) as a function of the patients' age class (1 through 4). Both outcome and predictor are declared as ordered factors. The proportional odds logistic regression uses the predictor as is, i.e. I do not specify polynomials or anything. The summary of the model contains the estimates:

Coefficients:
Value Std. Error t value
AGE_ord.L -0.22422 0.3869 -0.57954
AGE_ord.Q 0.21030 0.3717 0.56571
AGE_ord.C 0.02095 0.3559 0.05887
Now this is what I was talking about. polr() automatically uses a polynomial and returns three estimates, one less than the number of levels in AGE_ord. I assume the polynomial to be Lx+Qx^2+Cx^3 (I got that wrong in my initial post, mea culpa). BTW rms::lrm() will do the same, although the names are slightly different: AGE_cat, AGE_cat=3, and AGE_cat=4.

The output of avg_comparison() is a vector of length 12. It is merely a bold assumption of mine that these are four estimates, one per outcome level, for each of the three estimates AGE_ord.L through AGE_ord.C. This is why I converted it to a matrix with the labels showing my assumption:

lo lomed himed hi
L 0.05239002 0.01793816 -0.003243173 -0.06708501
Q 0.07794939 0.02424439 -0.007535115 -0.09465867
C 0.05231553 0.01791755 -0.003232824 -0.06700025

I'll rephrase my questions somewhat more cautiously: How do I compute the total marginal effect of a plus-one unit change of AGE_ord?
Thanks
Markus

[vincentarelbundock]

Thanks for clarifying.

Couple things:

1. I don't think polr() "automatically uses a polynomial". I believe it simply dichotomizes the categorical variable, in the same way other fitting functions like lm() would do.

2. Your question about the "total marginal effect" is under specified. What specific mathematical quantity are you interested in? For example, consider

a) Effect of going from Age 0 to Age 1 on Pr(AKI=1)
b) Effect of going from Age 0 to Age 1 on Pr(AKI=2)
c) Effect of going from Age 0 to Age 1 on Pr(AKI=3)
d) Effect of going from Age 0 to Age 1 on Pr(AKI=4)

Are you interested in the average of a, b, c, d?

Or consider average of:

a) Effect of going from Age 0 to Age 1 on Pr(AKI=1)
b) Effect of going from Age 0 to Age 2 on Pr(AKI=1)
c) Effect of going from Age 0 to Age 3 on Pr(AKI=1)
d) Effect of going from Age 0 to Age 4 on Pr(AKI=1)

These are completely different quantities of interest, and your vocabulary of "total effect" does not distinguish between the two. Unless you can define clearly, mathematically, what you want to estimate, it's hard to provide guidance...

[mhoenicka]

Thanks for bearing with me. I've expanded the example a little to make my question more specific.

library(MASS)
library(marginaleffects)

\# create dataset
\# predictor 1: age classes 1 through 4: lo<lomed<himed<hi
AGE_ord<-sample(1:4, 100, replace=TRUE)
AGE_ord<-factor(AGE_ord,ordered=TRUE)

\# predictor 2: cross-clamp time (numeric)
CCTIME_num<-sample(30:90,100,replace=TRUE)

\# predictor 3: diabetes mellitus (binary)
DM_yn<-sample(0:1,100,replace=TRUE)
DM_yn<-factor(DM_yn)

\# outcome: acute kidney injury stages 0 through 3 with 0<1<2<3
AKI_ord<-sample(0:3, 100, replace=TRUE)
AKI_ord<-factor(AKI_ord,ordered=TRUE)

my_df<-data.frame(AGE_ord,CCTIME_num,DM_yn,AKI_ord)

\# run ordered logistic regression
olmod<-polr(AKI_ord~AGE_ord+CCTIME_num+DM_yn,data=my_df,Hess=TRUE)
summary(olmod)

\# compute marginal effects
olcomp<-avg_comparisons(olmod)
olcomp$estimate
est<-matrix(olcomp$estimate,ncol=4,dimnames=list(c("AGE_ord.L","AGE_ord.Q","AGE_ord.C","CCTIME","DM"),c("0","1","2","3")))
est
`
So now we have a numeric and a binary predictor in addition to the ordinal one. The matrix computed from avg_comparisons() now looks like this:

`                     0             1             2            3`
`AGE_ord.L -0.049222670 -0.0197056298  0.0188674618  0.050060838`
`AGE_ord.Q -0.049629918 -0.0199129159  0.0190020637  0.050540771`
`AGE_ord.C -0.067839350 -0.0300063240  0.0244882261  0.073357447`
`CCTIME     0.001662404  0.0008087751 -0.0005590294 -0.001912149`
`DM         0.044604927  0.0221429384 -0.0148613617 -0.051886504

If I understand the marginaleffects documentation correctly, an increase of the numeric CCTIME parameter by one unit (here: one minute) increases the probability of a patient to end up with AKI 0 by 0.016. Similarly, switching the binary DM predictor from no to yes ("one unit") increases the probability of the outcome AKI 0 by 0.044. Is it possible at all to derive a similar statement from the AGE_ord.L, AGE_ord.Q, and AGE_ord.C marginal effects? Can I calculate what happens if I increase the age class by one unit, e.g. by moving from class 1 to class 2?

BTW my original statement that polr() automagically uses polynomials is derived from posts on stackexchange and stackoverflow:
https://stackoverflow.com/questions/57297771/interpretation-of-l-q-c-4-for-logistic-regression
https://stats.stackexchange.com/questions/233455/interpretation-of-l-q-output-from-a-negative-binomial-glm-with-categorical-d

Both deal with logistic, not ordinal logistic regression, but the main point is the representation of the effect of an ordinal predictor in a statistical model. The answer on stackoverflow seems to explain the issue, but the author loses me midway. However, I gather from his explanations that I cannot simply plug the original AGE_cat values into the polynomial defined by the estimates (boldfaced in the original reply), but does this hold true for marginal effects too? The answer also states that one should focus on the sign of the estimates rather than on their magnitude to describe the effect.

I'd appreciate if you had some additional clues.
Markus

[vincentarelbundock]

I recommend you forget about this whole idea of polynomial, stop reshaping the data, and try to interpret the results directly as printed by marginaleffects.

I also think it might be beneficial for you to read I have published online for free. Especially chapter 9 on categorical outcomes, and chapters 3, 4, 5, 6 on the conceptual framework: https://marginaleffects.com/

Scroll all the way down for two sentences of interpretation.

library(MASS)
library(marginaleffects)
set.seed(48103)

AGE_ord <- sample(1:4, 100, replace = TRUE)
AGE_ord <- factor(AGE_ord, ordered = TRUE)
CCTIME_num <- sample(30:90, 100, replace = TRUE)
DM_yn <- sample(0:1, 100, replace = TRUE)
DM_yn <- factor(DM_yn)
AKI_ord <- sample(0:3, 100, replace = TRUE)
AKI_ord <- factor(AKI_ord, ordered = TRUE)
my_df <- data.frame(AGE_ord, CCTIME_num, DM_yn, AKI_ord)
olmod <- polr(AKI_ord ~ AGE_ord + CCTIME_num + DM_yn, data = my_df, Hess = TRUE)

cmp <- avg_comparisons(olmod, variables = list("AGE_ord" = "sequential"))
cmp
#> 
#>  Group Contrast Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
#>      0    2 - 1  0.13036     0.0950  1.372    0.170 2.6 -0.0558 0.3166
#>      0    3 - 2 -0.05162     0.1018 -0.507    0.612 0.7 -0.2511 0.1479
#>      0    4 - 3  0.04175     0.1089  0.383    0.701 0.5 -0.1717 0.2552
#>      1    2 - 1  0.03777     0.0326  1.160    0.246 2.0 -0.0261 0.1016
#>      1    3 - 2 -0.00918     0.0192 -0.478    0.633 0.7 -0.0468 0.0285
#>      1    4 - 3  0.00793     0.0207  0.383    0.701 0.5 -0.0326 0.0485
#>      2    2 - 1 -0.05819     0.0446 -1.305    0.192 2.4 -0.1456 0.0292
#>      2    3 - 2  0.02638     0.0521  0.506    0.613 0.7 -0.0758 0.1285
#>      2    4 - 3 -0.02119     0.0558 -0.380    0.704 0.5 -0.1306 0.0882
#>      3    2 - 1 -0.10994     0.0845 -1.302    0.193 2.4 -0.2755 0.0556
#>      3    3 - 2  0.03443     0.0684  0.503    0.615 0.7 -0.0997 0.1685
#>      3    4 - 3 -0.02850     0.0735 -0.388    0.698 0.5 -0.1725 0.1155
#> 
#> Term: AGE_ord
#> Type:  probs

First row of output:

On average, moving from AGE_ord 1 to AGE_ord 2 changes the predicted probability of AKI_ord=0 by 0.13036 percentage point.

Last row of output:

On average, moving from AGE_ord 3 to AGE_ord 4 changes the predicted probability of AKI_ord=3 by -0.02850 percentage point.

[mhoenicka]

Jeez. I have always analyzed the estimates as I had to pool them from several imputations anyway. The marginaleffects::print() output explains it all. Thanks for rubbing this under my nose!
Markus

[vincentarelbundock]

:)

Glad we could figure this out. These things are surprisingly tricky, and I'm really hoping for this package to simplify things for people.