Logit

inv_logit <- function(x) {
  
  p <- 1 / (1 + exp(-x))
  return(p)
  
}

logit <- function(data, formula) {
  
  mf <- model.frame(formula, data, na.action=na.exclude)
  Y <- model.response(mf, type="double")
  X <- model.matrix(formula, mf)
  
  b <- numeric(ncol(X))
  
  logitll <- function(b) {
    
    p <- inv_logit((X %*% b))
    lli <- Y*log(p) + (1-Y)*log(1-p)
    
    return(-sum(lli))
    
  }
  
  est <- optim(b, logitll, method="BFGS")
  b <- est$par
  names(b) <- colnames(X)
  
  return(b)
  
}

Testing the function:

set.seed(42)

N <- 1000
X1 <- rnorm(N)
X2 <- rnorm(N)
ystar <- 1.5 + 2*X1 + 3*X2 + rnorm(N)
y <- rbinom(N, 1, inv_logit(ystar))

df <- data.frame(
  y, X1, X2
)

f <- y ~ X1 + X2

coef(glm(f,binomial,df))

## (Intercept)          X1          X2 
##    1.413335    1.839204    2.472301

logit(df, f)

## (Intercept)          X1          X2 
##    1.413368    1.839245    2.472373