How to get residuals from an ordinal logit/probit, and which ones to get

Question

I stumbled upon this question on Stack, where someone asked how to get the residuals from a polr regression, to which Ben Bolker answers as follows:

My question now is:

I would like to use a control function (CF) / two stage residual inclusion (2SRI) with a polr in the first stage. I was wondering if any of these different suggestions on how to calculate the residuals from the polr would be suitable for this purpose (if someone could show me how to do that in R, that would be amazing).

The normal 2SLS approach would be to include the fitted values of polr (which exist as the probabilities of each level of the ordinal var), but the result would be inconsistent when mimicking the 2SLS (by using the fitted values), and I don't know any software that does an ordinal-2SLS directly.

Answer 1

There are actually plenty of ways to get residuals from an ordinal probit/logit. Although polr does not provide any residuals, vglm provides several. See ?residualsvglm from the VGAMpackage (see also below).

NOTE: However, for the specific situation at hand Wooldridge (2014) suggests using the generalised residuals as described in Vella (1993). These are as far as I know currently not available in R, although I am working on that, but they are in Stata (using predict gr, score)

Residuals in VLGM

Surrogate residuals for polr

You can use the package sure (link), to calculate surrogate residuals with resids. The package is based on this paper, in the Journal of the American Statistical Association.

library(sure) # for residual function and sample data sets
library(MASS) # for polr function

df1 <- df1
df1$x1 <- df1$x
df1$x <- NULL
df1$y <- df2$y
df1$x2 <- df2$x
df1$x3 <- df3$x

options(contrasts = c("contr.treatment", "contr.poly"))
mod1 <- polr(as.ordered(y) ~ x1 + x2 + x3, data=df1, method='probit')

fit <- mod1$fitted.values
res <- resids(mod1)

EDIT: One big issue is the following (from ?resids):

"Note: Surrogate residuals require sampling from a continuous distribution; consequently, the result will be different with every call to resids. The internal functions used for sampling from truncated distributions when method = "latent" are based on modified versions of rtrunc and qtrunc."

Even when running resids(mod1, nsim=1000, method="latent"), there is no convergence of the outcome.