Interpreting odds ratios in ordinal logistic regression

I am trying to predict exam performance (below, average, above) based on whether participants attended a revision class. I am analysing my data in R using a proportional odds logistic regression. I am using the polr function from the MASS package.

Here is an example of my code:

data$$exam_performance <- factor(data$$exam_performance, c("Below", "Average", "Above"))

data$$attended <- factor(data$$attended, c("0", "1"))

model <- polr(formula = exam_performance ~ attended, data = data, Hess = TRUE)

exp(coef(model))


This returns an odds ratio of 18.64. Given my factor orderings, does this tell me that 'the odds of achieving average or above average performance is 18.64 times more if the participant attended than if they did not'?

Also, I notice that MASS does not provide p-values. I've googled this and it seems that LRT is the best way to obtain these. Would this be how I would do it?

intercept_only <- polr(formula = exam_performance ~ 1, data = data, Hess = TRUE)

anova(model, intercept_only)


This results in a Pr(Chi) value of 0. Should this be reported as p < .001?

The table of the factors is:

             Attended
0    1
Exam performance
Below 1482  300
Expected 1094 2822
Above   57 1422

• An OR of 18.64 is most likely a data issue and not a stable result. That would be backed up by the fact the LRT is singular. What's the output of with(data, table(exam_performance, attended))? As a note, it's bad code hygiene to recode a raw variable with factor and overwrite the original output. For instance, start with myval <- 0:1 then run myval <- factor(myval, levels=0:1, labels=c('No', 'Yes') once, you get a nice factor variable, Run it again, and myval is missing. – AdamO Sep 23 '19 at 15:52
• Apologies, these are the same variables. I generated an example here. I will change in the original question now. – Amy Sep 23 '19 at 16:16

That is a huge association. It goes from basically everybody below to everybody above who attend.

Fitting the model:

att <- c(0,1,0,1,0,1)
exam <- factor(c(0,0,1,1,2,2))
w <- c(1482, 300, 1094, 2822, 57, 1422)
f <- polr( exam ~ att, weights=w)


gives

Call:
polr(formula = exam ~ att, weights = w)

Coefficients:
att
2.925251

Intercepts:
0|1       1|2
0.2565983 3.7156750

Residual Deviance: 11686.09
AIC: 11692.09


As noted an OR of 18 ($$\approx(exp(3))$$).

Typing summary for the model gives one way of doing inference:

Call:
polr(formula = exam ~ att, weights = w)

Coefficients:
Value Std. Error t value
att 2.925    0.06634    44.1

Intercepts:
Value   Std. Error t value
0|1  0.2566  0.0390     6.5819
1|2  3.7157  0.0667    55.7015

Residual Deviance: 11686.09
AIC: 11692.09


the two sided Wald $$p$$-value is: 2*pt(44.1, df=3, lower.tail=F) = 0.000025. As you note, MASS does not calculate $$p$$-values this way because the "intercept(s)" terms do not have the same mathematical properties as the intercept in a logistic model, so you don't know what their distribution and standard error might be if the null hypothesis was true. Fitting the reduced model and testing the output with a LRT is the way to overcome this.

If I fit intercept only:

i <- polr( exam ~ 1, weights=w)

then anova(f, i) is:

Likelihood ratio tests of ordinal regression models

Response: exam
Model Resid. df Resid. Dev   Test    Df LR stat. Pr(Chi)
1     1      7175   14382.09
2   att      7174   11686.09 1 vs 2     1 2695.993       0


As you note, huge association and huge test statistic had unsurprising result: a large association. The $$p$$-value of 0 is just a consequence of rounding. $$p$$-values are never exactly 0. Reporting it out to 3 (or even 2) digits using $$p < 0.01$$ suffices, especially since significance testing is more concerned with meeting or exceeding alpha level than the actual precision of the $$p$$-value.

The interpretation of the coefficient is:

The odds of achieving a more desirable exam performance rating for a student that "attended" (clarifying beforehand how attended was defined) was 18 times higher than for a student that did not.

Since the categories are so few, you can also just summarize the predictions:

> round(predict(f, type = 'probs', newdata = data.frame(att=0:1)), 2)
0    1    2
1 0.56 0.41 0.02
2 0.06 0.62 0.31


You can say more than 50% were below average who did not attend, whereas only 6% were below average who did attend. And that only 2% were above expectation who did not attend versus 31% who were above expectation who did attend.

Another implementation of proportional odds that has more "off the shelf" functionality comes from Frank Harrel's rms package, specifically the lrm function. Fitting:

> lrm(exam ~ att, weights = w)
Logistic Regression Model

lrm(formula = exam ~ att, weights = w)

Sum of Weights by Response Category

0    1    2
1782 3916 1479

Model Likelihood     Discrimination    Rank Discrim.
Ratio Test           Indexes           Indexes
Obs             6     LR chi2    2695.99    R2       0.362    C       0.738
0              2     d.f.             1    g        1.755    Dxy     0.476
1              2     Pr(> chi2) <0.0001    gr       5.784    gamma   0.879
2              2                           gp       0.299    tau-a   0.285
Sum of weights7177                          Brier    0.129
max |deriv| 2e-13

Coef    S.E.   Wald Z Pr(>|Z|)
y>=1 -0.2566 0.0390  -6.58 <0.0001
y>=2 -3.7157 0.0667 -55.70 <0.0001
att   2.9253 0.0663  44.10 <0.0001


Gives the same Wald and LRT statistics that I calculated before.