GEST    DILATE    EFFACE    CONSIS    CONTR    MEMBRAN    AGE    STRAT    GRAVID    PARIT    DIAB    TRANSF    GEMEL    Y
31           3       100         3        1         2     26         3         1        0       2         2       1     1
28           8         0         3        1         2     25         3         1        0       2         1       2     1
31           3       100         3        2         2     28         3         2        0       2         1       1     1
...

install.packages("optimx")
  library(optimx)

vY = as.matrix(premature['PREMATURE'])
# Recoding the response variable
vY = ifelse(vY == "positif", 1, 0)

mX = as.matrix(premature[c('GEST', 'DILATE', 'EFFACE', 'CONSIS', 'CONTR',
                           'MEMBRAN', 'AGE', 'STRAT', 'GRAVID', 'PARIT',
                           'DIAB', 'TRANSF', 'GEMEL')])

#add an intercept to the predictor variables
mX = cbind(rep(1, nrow(mX)), mX)

#the number of variables and observations
iK = ncol(mX)
iN = nrow(mX)

#define the logistic transformation
logit = function(mX, vBeta) {
  return(exp(mX %*% vBeta)/(1+ exp(mX %*% vBeta)) )
}

# stable parametrisation of the log-likelihood function
logLikelihoodLogitStable = function(vBeta, mX, vY) {
  return(-sum(
    vY*(mX %*% vBeta - log(1+exp(mX %*% vBeta)))
    + (1-vY)*(-log(1 + exp(mX %*% vBeta)))
  )  # sum
  )  # return
}

# score function
likelihoodScore = function(vBeta, mX, vY) {
  return(t(mX) %*% (logit(mX, vBeta) - vY) )
}

# initial set of parameters (arbitrary starting parameters)
vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01, 0.01, 0.001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)

optimLogitLBFGS = optimx(vBeta0, logLikelihoodLogitStable,
                         method = 'L-BFGS-B', gr = likelihoodScore,
                         mX = mX, vY = vY, hessian=TRUE)

optimLogitLBFGS
                p1         p2       p3         p4         p5         p6
L-BFGS-B 9.720242 -0.1652943 0.525449 0.01681583 0.02781123 -0.3921004
                 p7          p8         p9       p10        p11        p12
L-BFGS-B -1.694412 -0.03461208 0.02759248 0.1993573 -0.6718275 0.02537887
                 p13      p14   value fevals gevals niter convcode  kkt1  kkt2
L-BFGS-B -0.8374338 0.625044 187.581    121    121    NA        1 FALSE FALSE
          xtimes
L-BFGS-B  0.044

vBeta0 = c(10, -0.1, -0.3, 0.001, 0.01, 0.01, 0.001, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01)

vBeta0 = rep(0.1, iK)

optimLogitLBFGS
                 p1             p2             p3              p4               p5
L-BFGS-B 0.372672689046 0.206785276091 0.398104550108 0.0175008380158 -0.0460042719084
                 p6             p7               p8            p9            p10
L-BFGS-B 0.139760396213 -1.43192069477 -0.0207666651106 -1.1396642657 0.212186387416
                 p11             p12             p13            p14         value
L-BFGS-B -0.583698421298 0.0576485672766 -0.802789658686 0.993103617257 185.472518798
         fevals gevals niter convcode  kkt1  kkt2 xtimes
L-BFGS-B    121    121    NA        1 FALSE FALSE   0.05