Manual calculation of ordinal logistic regression probabilities - starting at different comparison

Question

I am having a little issue with manually calculating some predicted probabilities from a ordinal logistic regression (I am doing this as a learning exercise - I am aware I don't have to do this manually).
Take the following example from the UCLA statistics training website:

Import the data:

library(foreign)
library(MASS)
dat <- read.dta("http://www.ats.ucla.edu/stat/data/ologit.dta")

Run an ordinal logistic regression:

m <- polr(apply ~ public + pared, data = dat, Hess=TRUE)

Get the predicted probabilities:

newdat <- expand.grid(public=unique(dat$public),pared=unique(dat$pared))
cbind(newdat,predict(m,newdat,type="probs"))

#  public pared  unlikely somewhat likely very likely
#1      0     0 0.5960409       0.3255546  0.07840448 # reproduce this row of results
#2      1     0 0.5720066       0.3421339  0.08585947
#3      0     1 0.3246644       0.4682876  0.20704799
#4      1     1 0.3033554       0.4728797  0.22376497

Now, manually reproduce these values:

  # store the intercept values
intcpts <- coef(summary(m))[3:4]

  # unlikely only
unl <- (1 / (1 + exp(-(intcpts[1] + sum(coef(m) * c(public=0,pared=0) )  ))))
# 0.5960409
  # somewhat likely ((unlikely + somewhat likely) - unlikely)
swl <- (1 / (1 + exp(-(intcpts[2] + sum(coef(m) * c(public=0,pared=0) )  )))) - unl
# 0.3255546
  # very likely (1 - somewhat likely  - unlikely)
vrl <- 1 - swl - unl
# 0.07840448

So this all works fine.
However, in the process of figuring this out, I found another page on the UCLA site (http://www.ats.ucla.edu/stat/spss/dae/ologit.htm) which seemed to suggest that one should calculate the values the other way around (very likely, then somewhat, then unlikely - reference category last). As in:

$ P(Y = 2) = \left(\frac{1}{1 + e^{-(a_{2}+b_{1}x_{1} + b_{2}x_{2} + b_{3}x_{3})}}\right) $

$ P(Y = 1) = \left(\frac{1}{1 + e^{-(a_{1}+b_{1}x_{1} + b_{2}x_{2} + b_{3}x_{3})}}\right) - P(Y = 2) $

$P(Y = 0) = 1 - P(Y = 1) - P(Y = 2) $

Now I realise this should give the exact same results as the probabilities always add to 1, but I can't seem to get their logic to work in terms of actual R code or manual calculations. I am stuck always starting with unlikely and then moving on up to very likely.

Am I missing something in the logic here? Can anyone suggest a simple calculation method?

Answer 1

Well, I get the correct results by removing the minus signs from the exponential...

> (P2 <- 1 / (1 + exp(intcpts[2])))
[1] 0.07840448

> (P1 <- 1 / (1 + exp(intcpts[1])) - P2)
[1] 0.3255546

> (P0 <- 1 - P1 - P2)
[1] 0.5960409

I streamlined these calculations because in the first row of your table, the values of public and pared are both zero.

Perhaps SPSS doesn't parameterize the model in the same way? Note in the help page for polr, the model is stated as: $$ \mbox{logit}\; P(Y \le k \;|\; x) = \zeta_k - \beta'x $$ where $\beta'x = \eta$ is the linear predictor. Note the minus sign there.