Multinomial Logit

mnlogit <- function(data, formula) {
  
  mf <- model.frame(formula, data, na.action=na.exclude)
  y <- model.response(mf)
  X <- model.matrix(formula, mf)
  
  N <- nrow(X)
  K <- ncol(X)
  levels <- unique(y) |> sort()
  M <- length(levels)
  
  Y <- matrix(N*M, nrow=N, ncol=M)
  for (m in 1:M) {
    Y[,m] <- ifelse(y==levels[m], 1, 0)
  }
  b <- numeric(K*(M-1))
  
  mnlogitLL <- function(param) {
    
    b <- matrix(param, nrow=K, ncol=M-1)
    Xb <- cbind(rep(0, N), X %*% b)
    lli <- numeric(N)
    
    for (m in 1:M) {
      lli <- lli + Y[,m]*Xb[,m] - Y[,m]*log(rowSums(exp(Xb)))
    }
    
    return(-sum(lli))
    
  }
  
  out <- optim(b, mnlogitLL, method="BFGS")
  est <- matrix(out$par, nrow=K, ncol=M-1)
  rownames(est) <- colnames(X)
  colnames(est) <- paste0(levels[1], "/", levels[2:M])
  
  return(est)
  
}

Testing the function:

library(mclogit)

## Loading required package: Matrix

housing <- MASS::housing #has an ordinal outcome but we'll ignore that for our purposes
f <- Sat ~ Infl + Freq + Type

t(mblogit(f, data = housing)$coefmat)

## 
## Iteration 1 - deviance = 150.7444 - criterion = 0.5116897
## Iteration 2 - deviance = 150.511 - criterion = 0.001549817
## Iteration 3 - deviance = 150.5107 - criterion = 1.965799e-06
## Iteration 4 - deviance = 150.5107 - criterion = 3.804011e-12
## converged

##                Response categories
## Predictors          Medium         High
##   (Intercept)    1.1352190 -0.740346318
##   InflMedium     0.1291512  0.005033149
##   InflHigh      -0.4256812  0.282489556
##   Freq          -0.0501481  0.026990219
##   TypeApartment  0.8092387 -0.644372087
##   TypeAtrium    -0.4053420  0.264568900
##   TypeTerrace   -0.3689781  0.189799430

mnlogit(housing, f)

##                Low/Medium     Low/High
## (Intercept)    1.13582927 -0.740243313
## InflMedium     0.12920635  0.005028019
## InflHigh      -0.42591879  0.282447109
## Freq          -0.05017388  0.026985033
## TypeApartment  0.80963960 -0.644298089
## TypeAtrium    -0.40569533  0.264391915
## TypeTerrace   -0.36904115  0.189905490