Ordered Logit

inv_logit <- function(x) {
  
  p <- 1 / (1 + exp(-x))
  return(p)
  
}

ologit <- function(data, formula, na.action=na.exclude) {
  
  mf <- model.frame(formula, data, na.action=na.exclude)
  Y <- model.response(mf)
  X <- model.matrix(formula, mf)[,-1]
  N <- nrow(X)
  
  levels <- unique(Y) |> sort()
  M <- length(levels)
  ncut <- M-1
  
  c <- 1:ncut/M
  b <- numeric(ncol(X))
  par <- c(c,b)
  
  ologitLL <- function(param) {
    
    c <- param[1:ncut]
    b <- param[(ncut+1):length(param)]
    
    Xb <-  X %*% b
    lli <- numeric(N)
    
    for (m in 1:M) {
      if (m == 1) {
        lli[Y==levels[m]] <- log(inv_logit(c[m]-Xb[Y==levels[m]]))
      } else if (m < M) {
        lli[Y==levels[m]] <- log(inv_logit(c[m]-Xb[Y==levels[m]]) - inv_logit(c[m-1]-Xb[Y==levels[m]]))
      } else {
        lli[Y==levels[m]] <- log(1 - inv_logit(c[m-1]-Xb[Y==levels[m]]))
      }
    }
    
    return(-sum(lli))
    
  }
  
  out <- optim(par, ologitLL, method="BFGS")
  est <- out$par
  
  names(est)[1:ncut] <- paste0("t", 1:ncut)
  names(est)[(ncut+1):length(est)] <- colnames(X)
  
  return(est)
  
}

Testing the function:

# Example from https://towardsdatascience.com/implementing-and-interpreting-ordinal-logistic-regression-1ee699274cf5

library(carData)
library(MASS)

data(WVS)
f <- poverty~religion+degree+country+age+gender

summary(polr(f, WVS))

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = f, data = WVS)
## 
## Coefficients:
##                  Value Std. Error t value
## religionyes    0.17973   0.077346   2.324
## degreeyes      0.14092   0.066193   2.129
## countryNorway -0.32235   0.073766  -4.370
## countrySweden -0.60330   0.079494  -7.589
## countryUSA     0.61777   0.070665   8.742
## age            0.01114   0.001561   7.139
## gendermale     0.17637   0.052972   3.329
## 
## Intercepts:
##                        Value   Std. Error t value
## Too Little|About Right  0.7298  0.1041     7.0128
## About Right|Too Much    2.5325  0.1103    22.9496
## 
## Residual Deviance: 10402.59 
## AIC: 10420.59

ologit(WVS, f)

## Warning in log(inv_logit(c[m] - Xb[Y == levels[m]]) - inv_logit(c[m - 1] - :
## NaNs produced

## Warning in log(inv_logit(c[m] - Xb[Y == levels[m]]) - inv_logit(c[m - 1] - :
## NaNs produced

## Warning in log(inv_logit(c[m] - Xb[Y == levels[m]]) - inv_logit(c[m - 1] - :
## NaNs produced

##            t1            t2   religionyes     degreeyes countryNorway 
##    0.72937251    2.53210358    0.17978268    0.14087551   -0.32235825 
## countrySweden    countryUSA           age    gendermale 
##   -0.60330084    0.61782030    0.01113128    0.17637937