Unexpected behaviour of logit regression with glm in R I recently was puzzled by the behaviour of R's glm when trying to compute a logistic regression
data <- data.frame(
    response = rbinom(600, 1, prob=rep(c(1,0.5,0),each=200)),
    predictor = rep(c("A","B","C"), each = 200)
)

llm1 <- glm(response ~ predictor, data = data, family=binomial(link='logit'))
summary(llm1)

The output tells me that there is apparently no difference in my predictors
Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    20.57    1253.73   0.016    0.987
predictorB    -20.81    1253.73  -0.017    0.987
predictorC    -41.13    1773.04  -0.023    0.981


When I compute glm with repsonse values where not all values are ones or zeros in predictors "A" and "C" respectively, I receive results as I had expected them.
data2 <- data.frame(
    response = rbinom(600, 1, prob=rep(c(0.95,0.5,0.05),each=200)),
    predictor = rep(c("A","B","C"), each = 200)
)

llm2 <- glm(response ~ predictor, data = data2, family=binomial(link='logit'))
summary(llm2)

# output:
...
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.9444     0.3244   9.075  < 2e-16 ***
predictorB   -2.8243     0.3540  -7.978 1.49e-15 ***
predictorC   -6.2612     0.5033 -12.441  < 2e-16 ***
...

I suspect it has something to do with the link function, but I can't really pinpoint it. Has anyone experienced this and can offer an explanation for the observed behaviour?
 A: This is due to perfect separation of the response variable, where the covariates perfectly predict the outcome.
This is the case in your example when predictor = "A"(always predict 1), and since A goes into the intercept you get an intercept of 20, pushing the probability toward 1. It's also the case when prediction = "C" but now the opposite way, since you now get $P(y = 1) = \exp(20 - 40) / (1 + \exp(20 - 40))$ giving you a prediction of 0.
The maximum likelihood estimate is not valid when you have perfect separation and the estimates should, in theory, be $\infty$, since the bigger the coefficient, the closer the probability is to 1 (or 0) and the better your prediction. Why is the estimate not infinity? Probably due to some stoping rule inside glm. The likelihood becomes very flat (see figure 2 here) and there is not a big enough change in the likelihood when doing an additional iteration in the optimization, even though you are not at a global maximum.
The flatness of the likelihood is also the reason behind the big standard errors, since these are based on the inverse Hessian (matrix of the second derivatives of the likelihood wrt the parameters), so a small curvature (flat likelihood) leads to large elements in the inverse Hessian.
As you can see when you don't have perfect separation, the parameter estimates are more sensible.
EDIT
One way to solve the problem of perfect separation is to add a prior to your parameters/penalize big estimates. This can be done using Ridge or Lasso, where you penalize the likelihood when $\beta$ becomes big. For your example:
data <- data.frame(
  response = rbinom(600, 1, prob=rep(c(1,0.5,0),each=200)),
  predictor = rep(c("A","B","C"), each = 200)
)

mod3 = ridge::logisticRidge(response ~ predictor, data = data)

summary(mod)

Call:
ridge::logisticRidge(formula = response ~ predictor, data = data)


Coefficients:
            Estimate Scaled estimate Std. Error (scaled) t value (scaled) Pr(>|t|)    
(Intercept)    2.438              NA                  NA               NA       NA    
predictorB    -2.471         -28.531               3.349            -8.52   <2e-16 ***
predictorC    -4.446         -51.339               3.835           -13.39   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ridge paramter: 0.0005737279, chosen automatically, computed using 2 PCs

Degrees of freedom: model 2.97 , variance 2.94 

where now the parameter estimates are a lot closer to zero, but you still get problems with the standard errors.
