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)


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:

                                    0    1
Exam performance 
                           Below 1482  300 
                        Expected 1094 2822 
                           Above   57 1422
  • $\begingroup$ 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. $\endgroup$ – AdamO Sep 23 '19 at 15:52
  • $\begingroup$ Apologies, these are the same variables. I generated an example here. I will change in the original question now. $\endgroup$ – 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)


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


      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:

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

    Value Std. Error t value
att 2.925    0.06634    44.1

    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.