Ordered Logit
<- function(x) {
inv_logit
<- 1 / (1 + exp(-x))
p return(p)
}
<- function(data, formula, na.action=na.exclude) {
ologit
<- model.frame(formula, data, na.action=na.exclude)
mf <- model.response(mf)
Y <- model.matrix(formula, mf)[,-1]
X <- nrow(X)
N
<- unique(Y) |> sort()
levels <- length(levels)
M <- M-1
ncut
<- 1:ncut/M
c <- numeric(ncol(X))
b <- c(c,b)
par
<- function(param) {
ologitLL
<- param[1:ncut]
c <- param[(ncut+1):length(param)]
b
<- X %*% b
Xb <- numeric(N)
lli
for (m in 1:M) {
if (m == 1) {
==levels[m]] <- log(inv_logit(c[m]-Xb[Y==levels[m]]))
lli[Yelse if (m < M) {
} ==levels[m]] <- log(inv_logit(c[m]-Xb[Y==levels[m]]) - inv_logit(c[m-1]-Xb[Y==levels[m]]))
lli[Yelse {
} ==levels[m]] <- log(1 - inv_logit(c[m-1]-Xb[Y==levels[m]]))
lli[Y
}
}
return(-sum(lli))
}
<- optim(par, ologitLL, method="BFGS")
out <- out$par
est
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)
<- poverty~religion+degree+country+age+gender
f
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