Ordinal Regression
The following demonstrates a standard cumulative link ordinal regression model via maximum likelihood. Default is with probit link function. Alternatively you can compare it with a logit link, which will result in values roughly 1.7*parameters estimates from the probit.
Data
This data generation is from the probit perspective, where the underlying continuous latent variable is normally distributed.
library(tidyverse)
set.seed(808)
= 1000 # Sample size
N = cbind(x1 = rnorm(N), x2 = rnorm(N)) # predictor variables
x = c(1,-1) # coefficients
beta = rnorm(N, mean = x %*% beta) # the underlying latent variable
y_star = y_star > -1.5 # -1.50 first cutpoint
y_1 = y_star > .75 # 0.75 second cutpoint
y_2 = y_star > 1.75 # 1.75 third cutpoint
y_3 = y_1 + y_2 + y_3 + 1 # target
y
table(y)
y
1 2 3 4
175 495 182 148
= data.frame(x, y = factor(y)) d
Function
<- function(par, X, y, probit = TRUE) {
ordinal_ll = length(unique(y)) # number of classes K
K = K-1 # number of cutpoints/thresholds
ncuts = par[(1:ncuts)] # cutpoints
cuts = par[-(1:ncuts)] # regression coefficients
beta = X %*% beta # linear predictor
lp = rep(0, length(y)) # log likelihood
ll = ifelse(probit, pnorm, plogis) # which link to use
pfun
for(k in 1:K){
if (k==1) {
==k] = pfun((cuts[k] - lp[y==k]), log = TRUE)
ll[y
}else if (k < K) {
==k] = log(pfun(cuts[k] - lp[y==k]) - pfun(cuts[k-1] - lp[y==k]))
ll[y
}else {
==k] = log(1 - pfun(cuts[k-1] - lp[y==k]))
ll[y
}
}
-sum(ll)
}
Estimation
= c(-1, 1, 2, 0, 0) # initial values init
= optim(
fit_probit
init,
ordinal_ll,y = y,
X = x,
probit = TRUE,
control = list(reltol = 1e-10)
)
= optim(
fit_logit
init,
ordinal_ll,y = y,
X = x,
probit = FALSE,
control = list(reltol = 1e-10)
)
Comparison
We can compare our results with the ordinal package.
library(ordinal)
= clm(y ~ x1 + x2, data = d, link = 'probit')
fit_ordpack_probit = clm(y ~ x1 + x2, data = d, link = 'logit') fit_ordpack_logit
method | cut_1 | cut_2 | cut_3 | beta1 | beta2 | |
---|---|---|---|---|---|---|
probit | ordinal_ll | -1.609 | 0.731 | 1.798 | 1.019 | -1.051 |
probit | ordpack | -1.609 | 0.731 | 1.798 | 1.019 | -1.051 |
logit | ordinal_ll | -2.872 | 1.284 | 3.168 | 1.806 | -1.877 |
logit | ordpack | -2.872 | 1.284 | 3.168 | 1.806 | -1.877 |
Source
Original code available at https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/ordinal_regression.R