Why does logistic regression become unstable when classes are well-separated?

Question

Why is it that logistic regression becomes unstable when classes are well-separated? What does well-separated classes mean? I would really appreciate if someone can explain with an example.

Answer 1

It isn't correct that logistic regression in itself becomes unstable when there are separation. Separation means that there are some variables which are very good predictors, which is good, or, separation may be an artifact of too few observations/too many variables. If that is the case, the solution might be to get more data. But separation itself, then, is only a symptom, and not a problem in itself.

So there are really different cases to be treated. First, what is the goal of the analysis? If the final result of the analysis is some classification of cases, separation is no problem at all, it really means that there are very good variables giving very good classification. But if the goal is risk estimation, we need the parameter estimates, and with separation the usual mle (maximum likelihood) estimates do not exist. So we must change estimation method, maybe. There are several proposals in the literature, I will come back to that.

Then there are (as said above) two different possible causes for separation. There might be separation in the full population, or separation might be caused by to few observed cases/too many variables.

What breaks down with separation, is the maximum likelihood estimation procedure. The mle parameter estimates (or at least some of them) becomes infinite. I said in the first version of this answer that that can be solved easily, maybe with bootstrapping, but that does not work, since there will be separation in each bootstrap resample, at least with the usual cases bootstrapping procedure. But logistic regression is still a valid model, but we need some other estimation procedure. Some proposals have been:

1. regularization, like ridge or lasso, maybe combined with bootstrap.
2. exact conditional logistic regression
3. permutation tests, see https://www.ncbi.nlm.nih.gov/pubmed/15515134
4. Firths bias-reduced estimation procedure, see Logistic regression model does not converge
5. surely others ...

If you use R, there is a package on CRAN, SafeBinaryRegression, which help with diagnosing problems with separation, using mathematical optimization methods to check for sure if there is separation or quasiseparation! In the following I will be giving a simulated example using this package, and the elrm package for approximate conditional logistic regression.

First, a simple example with the safeBinaryRegression package. This package just redefines the glm function, overloading it with a test of separation, using linear programming methods. If it detects separation, it exits with an error condition, declaring that the mle does not exist. Otherwise it just runs the ordinary glm function from stats. The example is from its help pages:

library(safeBinaryRegression)   # Some testing of that 
                                # package,
                                # based on its examples
# complete separation:
x  <-  c(-2, -1, 1, 2)
y  <-  c(0, 0, 1, 1)
glm(y ~ x, family=binomial)
glm(y ~ x,  family=binomial,  separation="test")
stats::glm(y ~ x, family=binomial)
# Quasicomplete separation:
x  <-  c(-2, 0, 0, 2)
y  <-  c(0, 0, 1, 1)
glm(y ~ x, family=binomial)
glm(y ~ x,  family=binomial,  separation="test")
stats::glm(y ~ x, family=binomial)

The output from running it:

> # complete separation:
> x  <-  c(-2, -1, 1, 2)
> y  <-  c(0, 0, 1, 1)
> glm(y ~ x, family=binomial)
Error in glm(y ~ x, family = binomial) : 
  The following terms are causing separation among the sample points: (Intercept), x
> glm(y ~ x,  family=binomial,  separation="test")
Error in glm(y ~ x, family = binomial, separation = "test") : 
  Separation exists among the sample points.
    This model cannot be fit by maximum likelihood.
> stats::glm(y ~ x, family=binomial)

Call:  stats::glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
 -9.031e-08    2.314e+01  

Degrees of Freedom: 3 Total (i.e. Null);  2 Residual
Null Deviance:      5.545 
Residual Deviance: 3.567e-10    AIC: 4
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> # Quasicomplete separation:
> x  <-  c(-2, 0, 0, 2)
> y  <-  c(0, 0, 1, 1)
> glm(y ~ x, family=binomial)
Error in glm(y ~ x, family = binomial) : 
  The following terms are causing separation among the sample points: x
> glm(y ~ x,  family=binomial,  separation="test")
Error in glm(y ~ x, family = binomial, separation = "test") : 
  Separation exists among the sample points.
    This model cannot be fit by maximum likelihood.
> stats::glm(y ~ x, family=binomial)

Call:  stats::glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
  5.009e-17    9.783e+00  

Degrees of Freedom: 3 Total (i.e. Null);  2 Residual
Null Deviance:      5.545 
Residual Deviance: 2.773    AIC: 6.773

Now we simulate from a model which can be closely approximated by a logistic model, except that above a certain cutoff the event probability is exactly 1.0. Think about a bioassay problem, but above the cutoff the poison always kills:

pl  <-  function(a, b, x) 1/(1+exp(-a-b*x))
a  <-  0
b  <-  1.5
x_cutoff  <-  uniroot(function(x) pl(0,1.5,x) - 
                       0.98,lower=1,upper=3.5)$root
### circa 2.6
pltrue  <-  function(a, b, x) ifelse(x < x_cutoff, pl(a, b, 
                 x), 1.0)

x <- -3:3

### Let us simulate many times from this model,  and try to 
    # estimate it
### with safeBinaryRegression::glm  That way we can estimate 
    #the probability
### of separation from this model

set.seed(31415926)  ### May I have a large container of 
                       # coffee 
replications  <-  1000
p  <-  pltrue(a, b, x)
err  <-  0
good  <- 0

for (i in 1:replications) {
    y  <- rbinom(length(x), 1, p)
    res  <-  try(glm(y~x, family=binomial), silent=TRUE)
    if (inherits(res,"try-error")) err <-  err+1 else 
                      good <- good+1
}
P_separation  <-  err/replications
P_separation

When running this code, we estimate the probability of separation as 0.759. Run the code yourself, it is fast!

Then we extend this code to try different estimations procedures, mle and approximate conditional logistic regression from elrm. Running this simulation take around 40 minutes on my computer.

library(elrm)  # from CRAN
set.seed(31415926)  ### May I have a large container of 
                        # coffee
replications  <-  1000
GOOD  <-  numeric(length=replications) ### will be set to one 
           # when MLE exists!
COEFS <- matrix(as.numeric(NA), replications, 2)
COEFS.elrm <- matrix(as.numeric(NA), replications, 2) # But 
                 # we'll only use second col for x
p  <-  pltrue(a, b, x)
err   <-  0
good  <- 0

for (i in 1:replications) {
    y  <- rbinom(length(x), 1, p)
    res  <-  try(glm(y~x, family=binomial), silent=TRUE)
    if (inherits(res,"try-error")) err <-  err+1 else{ 
          good <- good+1
          GOOD[i] <- 1 }
    # Using stats::glm
    mod  <-  stats::glm(y~x, family=binomial)
    COEFS[i, ]  <-  coef(mod)
    # Using elrm:
    DATASET  <-  data.frame(x=x, y=y, n=1)
    mod.elrm  <-  elrm(y/n ~ x,  interest= ~ x -1, r=4, 
                       iter=10000, burnIn=1000,
                       dataset=DATASET)
    COEFS.elrm[i, 2 ]  <-  mod.erlm$coeffs       
}
### Now we can compare coefficient estimates of x,
###  when there are separation,  and when not:

non  <-  which(GOOD==1)
cof.mle.non  <-  COEFS[non, 2, drop=TRUE]
cof.mle.sep  <-  COEFS[-non, 2, drop=TRUE]
cof.elrm.non  <-  COEFS.elrm[non, 2, drop=TRUE]
cof.elrm.sep  <-  COEFS.elrm[-non, 2, drop=TRUE]

Now we want to plot the results, but before that, note that ALL the conditional estimates are equal! That is really strange and should need an explanation ... The common value is 0.9523975. But at least we obtained finite estimates, with confidence intervals which contains the true value (not shown here). So I will only show a histogram of the mle estimates in the cases without separation:

hist(cof.mle.non, prob=TRUE)

[

What is remarkable is that all the estimates is lesser than the true value 1.5. That can have to do with the fact that we simulated from a modified model, needs investigation.

Answer 2

There are good answers here from @sean501 and @kjetilbhalvorsen. You asked for an example. Consider the figure below. You might come across some situation in which the data generating process is like that depicted in panel A. If so, it is quite possible that the data you actually gather look like those in panel B. Now, when you use data to build a statistical model, the idea is to recover the true data generating process or at least come up with an approximation that is reasonably close. Thus, the question is, will fitting a logistic regression to the data in B yield a model that approximates the blue line in A? If you look at panel C, you can see that the gray line better approximates the data than the true function does, so in seeking the best fit, the logistic regression will 'prefer' to return the gray line rather than the blue one. It doesn't stop there, however. Looking at panel D, the black line approximates the data better than the gray one—in fact, it is the best fit that could possibly occur. So that is the line the logistic regression model is pursuing. It corresponds to an intercept of negative infinity and a slope of infinity. That is, of course, very far from the truth that you are hoping to recover. Complete separation can also cause problems with the calculation of the p-values for your variables that come standard with logistic regression output (the explanation there is slightly different and more complicated). Moreover, trying to combine the fit here with other attempts, for example with a meta-analysis, will just make the other findings less accurate.

Answer 3

It means that there is a hyperplane such that on one side there are all the positive points and on the other side all the negative. The maximum likelihood solution is then flat 1 on one side and flat 0 on other side, which is 'achieved' with the logistic function by having the coefficients at infinity.

Answer 4

What you are calling "separation" (not 'seperation') covers two different situations that end up causing the same issue – which I would not call, however, an issue of "instability" as you do.

An illustration: Survival on the Titanic

• Let $DV \in (0, 1)$ be a binary dependent variable, and $SV$ a separating, independent variable.

  Let's suppose that $SV$ is the class of the passengers on the Titanic, and that $DV$ indicates whether they survived the wreckage, with $0$ indicating death and $1$ indicating survival.

• Complete separation is the situation where $SV$ predicts all values of $DV$.

  That would be the case if all first-class passengers on the Titanic had survived the wreckage, and none of the second-class passengers had survived.

• Quasi-complete separation is the situation where $SV$ predicts either all cases where $DV = 0$, or all cases where $DV = 1$, but not both.

  That would be the case if some first-class passengers on the Titanic had survived the wreckage, and none of the second-class passengers had survived. In that case, passenger class $SV$ predicts all cases where $DV = 1$, but not all cases where $DV = 0$.

  Reversely, if only some second-class passengers on the Titanic had died in the wreckage, then passenger class $SV$ predicts all cases where $DV = 0$, but not all cases where $DV = 1$, which includes both first-class and second-class passengers.

What you are calling "well-separated classes" is the situation where a binary outcome variable $DV$ (e.g. survival on the Titanic) can be completely or quasi-completely mapped to a predictor $SV$ (e.g. passenger class membership; $SV$ need not be binary as it is in my example).

Why is logistic regression "unstable" in these cases?

This is well explained in Rainey 2016 and Zorn 2005.

• Under complete separation, your logistic model is going to look for a logistic curve that assigns, for example, all probabilities of $DV$ to $1$ when $SV = 1$, and all probabilities to $DV$ to $0$ when $SV = 0$.

  This corresponds to the aforementioned situation where only and all first-class passengers of the Titanic survive, with $SV = 1$ indicating first-class passenger membership.

  This is problematic because the logistic curve lies strictly between $0$ and $1$, which means that, to model the observed data, the maximisation is going to push some of its terms towards infinity, in order, if you like, to make $SV$ "infinitely" predictive of $DV$.

• The same problem arises under quasi-complete separation, as the logistic curve will still need to assign only values of either $0$ or $1$ to $DV$ in one of two cases, $SV = 0$ or $SV = 1$.

In both cases, the likelihood function of your model will be unable to find a maximum likelihood estimate: it will only find an approximation of that value by approaching it asymptotically.

What you are calling "instability" is the fact that, in cases of complete or quasi-complete separation, there is no finite likelihood for the logistic model to reach. I would not use that term, however: the likelihood function is, in fact, being pretty "stable" (monotonic) in its assignment of coefficient values towards infinity.

---

Note: my example is fictional. Survival on the Titanic did not boil down just to passenger class membership. See Hall (1986).