Calculate coefficients in a ordinal logistic regression with R

Question

Following the question about manually fitting logistic regression, can someone provide the same 'manual' way to fit a ordinal logistic regression with ordered categorical response?

Answer 1

Let's assume you have an ordinal response with 4 categories, so you want to estimate 3 parameters which will be the thresholds of the latent continuous response.

First think the log-likelihood function for a given set of $\mathbf{\theta}$ parameters you want to estimate.

ll.ord.log <- function(theta, x, y){
    t1    <- theta[1]
    t2    <- theta[2]
    t3    <- theta[3]
    b     <- matrix(theta[4:5])
    one   <- cbind(log( 1/(1+exp(x%*%b-t1)) - 0 ), y)
    two   <- cbind(log( 1/(1+exp(x%*%b-t2)) - 1/(1+exp(x%*%b-t1)) ), y)
    three <- cbind(log( 1/(1+exp(x%*%b-t3)) - 1/(1+exp(x%*%b-t2)) ), y)
    four  <- cbind(log( 1                   - 1/(1+exp(x%*%b-t3)) ), y)

    ll <- sum(one[,1][one[,2]==1]) + 
          sum(two[,1][two[,2]==2]) + 
          sum(three[,1][three[,2]==3]) + 
          sum(four[,1][four[,2]==4])
    -ll
  }

Then set some initial values for the optimization and some constrains. Thus, we use the constrOptim function which enables us to perform optimization given some constrains of the estimated parameters, since the thresholds need to be ordered.

theta0 <- c(-2,1,4,1,-1)
ui     =  rbind(c(-1,1,0,0,0), c(0,-1,1,0,0), c(-1,0,1,0,0))
ci     =  c(.1,.1,.1)

Then, you can use the function from R to estimate the maximum likelihood and the estimates at this value.

constrOptim(theta0, ll.ord.log, NULL, x=x, y=y, ui=ui, ci=ci)