Binary logit modelling with R - Issue finding same results

Question

I need your help to figure out something about the estimation of simple binary logit model in R.

As nicely explained on the following website (https://stats.idre.ucla.edu/r/dae/logit-regression/) the GLM command can be used for the estimation of such model (when specifying the binomial(link='logit') link function). I did manage to find same results as those from the example, but then I tried to code my own log-likelihood function and to optimize it with both "optim" and "maxLik" optimization packages, but I could not find same results! This should not be the case (and the results also differ across these two optimization packages!). Am I missing something obvious here?

1. R code for the "reference" model (from the example)

data = NULL
data = read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

data$rank2 = ifelse(data$rank==2,1,0)
data$rank3 = ifelse(data$rank==3,1,0)
data$rank4 = ifelse(data$rank==4,1,0)
data$cst = 1

logit = glm(admit ~ -1 + cst + gre + gpa + rank2 + rank3 + rank4,
              data=data, family=binomial(link='logit'))

2. R code for the log-likelihood (LL) function

BL_LL = function(param, Data, X, Y){
  num = as.matrix(Data[,X]) %*% as.vector(param[1:length(X)])
  prb = exp(num) / (1+exp(num))
  llik = Data[,Y]*log(prb) + (1-Data[,Y])*log(1-prb)
  return(-sum(llik))}

3. R code for the estimation of the binary logit model with the user-specified LL function

3.1. Model specification and data declaration

D = data
X = c('cst','gre','gpa','rank2','rank3','rank4')
Y ='admit'
sv = matrix(0, ncol=length(X))

3.2. Model estimation with "maxLik" and 2 different optimization routines

m1 = sapply(lx = c('NM','BFGS'), function(x){
  k = NULL
  k = maxLik(BL_LL, start=sv, Data=D, X=X, Y=Y, method=x)
  round(100*coef(k)/as.numeric(coef(logit)),2)})
m1

3.3. Model estimation with "optim" and 2 different optimization routines

m2 = sapply(c("Nelder-Mead","BFGS"), function(x){
  k = NULL
  k = optim(par=sv, fn=BL_LL, Data=D, X=X, Y=Y, method=x)
  round(100*k$par/as.numeric(coef(logit)),2)})
m2

The two output tables (i.e., m1 and m2) measure the differences between the "reference" estimates (from the example) and those obtained with the user-specified LL function. In case of perfect matching, all the values should be close to 100%, but if you run the code you will notice some very large differences. I don't understand what might be causing this issue. My understanding is that the binary logit model has a closed-form solution and then both NM and BFGS algorithms should return same results disregarding the specification of the starting values. I don't think providing a user-specified gradient function would make a difference (I haven't tried though). Any thoughts? (Thanks in advance for your help!)

Answer 1

Alright - It seems that the last sentence "I don't think providing a user-specified gradient function would make a difference" was the issue. But I don't really understand why (it must have somethign to do with the coding of the optimization routines as GLM is not based on BFGS but on IWLS). I've tried to estimate the binary logit model with user-specified log-likelihood, gradient and hessian functions. I still obtained different results when using the "maxLik" optimization package, but at least the results match when using "optimx" (Newer version of "optim" built-in package) [Remark: I have used "optimx" instead of "optim" because the last one does not allow for user-specified hessian function].

Here is the R code if you want to replicate/double check the results:

### BL modelling functions

BL_LL = function(param, Data, X, Y){
  num = as.matrix(Data[,X]) %*% as.vector(param[1:length(X)])
  prb = exp(num) / (1+exp(num))
  llik = Data[,Y]*log(prb) + (1-Data[,Y])*log(1-prb)
  return(-sum(llik))}

BL_GRAD = function(param, Data, X, Y){
  num = as.matrix(Data[,X]) %*% as.vector(param[1:length(X)])
  prb = exp(num) / (1+exp(num))
  v = Data[,Y] - prb
  vx = as.matrix(Data[,X])*as.vector(v)
  grad = colSums(vx)
  return(-grad)}

BL_HESS = function(param, Data, X, Y){
  mX = as.matrix(Data[,X])
  num = mX %*% as.vector(param[1:length(X)])
  prb = exp(num) / (1+exp(num))
  v = diag(as.vector(prb*(1-prb)))
  hess = -t(mX) %*% v %*% mX
  return(-hess)}

### Get data from web example
# https://stats.idre.ucla.edu/r/dae/logit-regression/ 

data = NULL
data = read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")

data$rank2 = ifelse(data$rank==2,1,0)
data$rank3 = ifelse(data$rank==3,1,0)
data$rank4 = ifelse(data$rank==4,1,0)
data$cst = 1

### Model estimation with GLM

logit = glm(admit ~ -1 + cst + gre + gpa + rank2 + rank3 + rank4,
              data=data, family=binomial(link='logit'))

### Model specification

D = data
X = c('cst','gre','gpa','rank2','rank3','rank4')
Y ='admit'
sv = matrix(0, ncol=length(X))

### Model estimation with "maxlik"

l.maxlik = as.list(NULL)
l.maxlik[[1]] = coef(maxLik(start=sv, Data=D, X=X, Y=Y, method="BFGS",
                            logLik=BL_LL,
                            grad=NULL,
                            hess=NULL))
l.maxlik[[2]] = coef(maxLik(start=sv, Data=D, X=X, Y=Y, method="BFGS",
                            logLik=BL_LL,
                            grad=BL_GRAD,
                            hess=NULL))
l.maxlik[[3]] = coef(maxLik(start=sv, Data=D, X=X, Y=Y, method="BFGS",
                            logLik=BL_LL,
                            grad=BL_GRAD,
                            hess=BL_HESS))
m1 = matrix(unlist(l.maxlik), ncol=length(l.maxlik))
round(100*m1/coef(logit),2)

### Model estimation with "optimx"

l.optimx = as.list(NULL)
l.optimx[[1]] = coef(optimx(par=as.vector(sv), Data=D, X=X, Y=Y, method="BFGS",
                            fn=BL_LL,
                            gr=NULL,
                            hess=NULL))
l.optimx[[2]] = coef(optimx(par=as.vector(sv), Data=D, X=X, Y=Y, method="BFGS",
                            fn=BL_LL,
                            gr=BL_GRAD,
                            hess=NULL))
l.optimx[[3]] = coef(optimx(par=as.vector(sv), Data=D, X=X, Y=Y, method="BFGS",
                            fn=BL_LL,
                            gr=BL_GRAD,
                            hess=BL_HESS))
m3 = matrix(unlist(l.optimx), ncol=length(l.optimx))
round(100*m3/coef(logit),2)