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

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)


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?

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.

• Thanks, but I'm still confused how would this work when you don't use the case where coef(m) are both 0 and don't affect the returned value. – thelatemail Sep 22 '14 at 22:53
• Ok, I think I get it now. Something like 1 / (1 + exp(-(sum(coef(m) * c(public=1,pared=1)) - intcpts[2]))) should get there for very likely, where public/pared = 1/1 – thelatemail Sep 22 '14 at 23:09
• I think you need to get over those minus signs. More likely it is 1 / (1 + exp(sum(coef(m) * c(public=1,pared=1) + intcpts[2]))). – Russ Lenth Sep 23 '14 at 0:20
• That seems to not give the result I am looking for, i.e. - it gives 0.00213 instead of the 0.2237 it should give. – thelatemail Sep 23 '14 at 0:32
• Hmmmm. Well, I'm sure about the last plus sign because that is what was needed for the $(0,0)$ case. – Russ Lenth Sep 23 '14 at 0:35