Calculating logistic regression coefficients by hand

Viewed 655

I am trying to estimate logistic regression in R calculating everything by hand. I am able to create logit and loglikelihood function however I am not able to solve it using som non-linear solver

I would like to ask for advice

df <- read_csv("http://courses.atlas.illinois.edu/spring2016/STAT/STAT200/RProgramming/data/Default.csv")
df


df$default = ifelse(df$default == "Yes", 1, 0)

logit <- function(x, b0, b1) {
  1/(1 + exp(-b0 - b1*x))
}



Loglikel <- function(y, x, b0, b1) {
  b0 = rep(b0, length(y))
  b1 = rep(b1, length(y))
  p <- logit(x, b0, b1)
  sum(y*log(p)  + (1 - y)*log(1-  p))
}


Loglikel(df$default, df$balance, -10, 0.005)


library(stats4)

mle(Loglikel, 
    start = list(b0 = 0, b1 = 0), 
    fixed = list(y = df$default, x = df$balance))


1 Answers

I took your code and modified it a little bit to pass parameters as a vector:

df <- read_csv("http://courses.atlas.illinois.edu/spring2016/STAT/STAT200/RProgramming/data/Default.csv")
df$default <- ifelse(df$student == "Yes", 1, 0)

logit <- function(x, b0, b1) {
  1/(1 + exp(-b0-b1*x))
}
Loglikel <- function(par, y, x){
  p <- logit(x, par[1], par[2])
  sum(y*log(p)  + (1-y)*log(1-p))
}

We are now ready to use a non-linear solver (e.g., nlm) to obtain estimates for parameters:

nlm_fit <- nlm(Loglikel, p = c(-2,0.001), x=df$balance, y=df$default)

which gives

> nlm_fit
...
$estimate
[1] -2.0002960 -0.2666521
...

nlm uses a Newton-Raphson-type solver to minimize the MLE. At the same time, glm uses the Iteratively Reweighted Least Squares algorithm, which means that output of glm and nlm does not have to agree:

glm_fit <- glm(default ~ balance, family = binomial(link="logit"), data = df)

> glm_fit

Call:  glm(formula = default ~ balance, family = binomial(link = "logit"), 
    data = df)

Coefficients:
(Intercept)      balance  
 -1.7004224    0.0009409 

Check this link out, it gives a good summary of what is going on under the hood of glm.

Related