Why does the glm function converge and not give an error when all y's are equal to the same value?

Question

I need to fit a univariate logistic model with few observations (between 10 and 20). In some cases, y is equal to the same value (example 1) for all observations. Theoretically, the model should not converge. But, when I use the glm function in R it doesn't show me an error or a warning! Here is an example of code that I tested:

x=sample(c(0,1),20,replace = TRUE)
y=rep(1,20)
summary(glm(y~x, family = binomial))

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

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
3.971e-06  3.971e-06  3.971e-06  3.971e-06  3.971e-06  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.557e+01  7.200e+04       0        1
x           -4.549e-10  9.708e+04       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 0.000e+00  on 19  degrees of freedom
Residual deviance: 3.154e-10  on 18  degrees of freedom
AIC: 4

Number of Fisher Scoring iterations: 24

What do I need to change so that the glm function gives me an error instead of a result? I don't understand how the model can converge! Contrary to this question I getting "algorithm did not converge" I don't get a warning with my model.

I find it strange that R does not show me an error while with another software (SAS) I get an error and the program stops!

Note also that even with 3 observations the glm function converge and gives results without displaying a warning !

y=rep(1,3)
x=sample(c(0,1),3,replace = TRUE)
summary(glm(y~x, family = binomial))

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

Deviance Residuals: 
        1          2          3  
1.079e-05  1.079e-05  1.079e-05  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.357e+01  7.946e+04       0        1
x           -2.545e-08  9.732e+04       0        1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 0.0000e+00  on 2  degrees of freedom
Residual deviance: 3.4957e-10  on 1  degrees of freedom
AIC: 4

Number of Fisher Scoring iterations: 22

I also want to specify that when we increase the number of observations y>= 101, I get this warning :

x=sample(c(0,1),101,replace = TRUE)
y=rep(1,101)

Warning message:
glm.fit: the algorithm did not converge

Whereas with less observations I don't get this warning

I'd be grateful for your help.

Answer 1

In some cases, y is equal to the same value (example 1) for all observations. Theoretically, the model should not converge.

Nonsense. This is a very simple dataset for which the maximum likelihood results are known in closed form. Convergence (in terms of the fitted values) is quite simple.

But, when I use the glm function in R it doesn't show me an error or a warning!

Of course not. Instead it has converged and given you the correct results. If all the y=1, then the maximum likelihood solution occurs when:

• all fitted values are equal to 1
• residual deviance is 0
• intercept is positive infinity
• slope is zero

and that is exactly what R's glm function has given you. To working precision, the intercept value is 25.57 is large enough to be effectively infinite, because larger values would only change the fitted values in the 10th decimal place.

Note that, although it is hard to measure convergence of the intercept to infinity, it is very simple to measure convergence of the fitted values to 1 and convergence of the deviance to zero. It is exactly to handle cases like this that the convergence criterion for the glm iterative algorithm is defined in terms of the fitted values rather than in terms of the coefficient estimates.

The definition of the convergence criterion in terms of fitted values goes back to the earliest papers on generalized linear models nearly 50 years ago. Indeed, the whole glm iterative algorithm is defined in terms of the fitted values rather than the coefficient estimates. The algorithm is described in full in the original 1972 paper by Nelder and Wedderburn and has been repeated in many books and references since then. The algorithm was implemented in the original glm software GLIM back in the 1970s and GLIM did not return errors for this type of dataset either.

What do I need to change so that the glm function gives me an error instead of a result?

Why would you want to stop the function from working properly and giving correct results?

another software (SAS) I get an error and the program stops!

Why don't you write to SAS and ask them why such an expensive piece of software as SAS can't cope with such a small simple dataset?

Answer to similar question on R-devel mailing list

A very similar question was asked on the R-devel mailing list nearly 20 years ago, when Brian Ripley and I gave the same answer as above:

• https://stat.ethz.ch/pipermail/r-devel/2003-June/026812.html
• https://stat.ethz.ch/pipermail/r-devel/2003-June/026816.html

Example with ``manual'' calculations

Here is a simple example where I implement the iteratively reweighted least square (IRWLS) algorithm proposed by Nelder and Wedderburn (1972) directly, instead of using the implementation in the glm function.

Let's suppose that there are four observations and two groups:

n <- 4
y <- rep(1,n)
g <- gl(2,n/2)
X <- model.matrix(~g)

We need a starting value for the success probabilities (which for binary regression are the same as the fitted values). I've started the success probabilities at 0.5 for every observation. In practice, it would obviously be possible to do better than that.

p <- rep(0.5,n)
eta <- qlogis(p)

Now run the IRWLS algorithm for 10 iterations, outputing the fitted probabilities after each iteration.

for (i in 1:10) {
  w <- p*(1-p)
  linkderiv <- p*(1-p)
  z <- (y-p)/linkderiv + eta
  fit <- lm.wfit(X,z,w)
  eta <- fit$fitted.value
  p <- plogis(eta)
  cat(p,"\n")
}

The output shows that the fitted values converge quickly to 1:

0.8807971 0.8807971 0.8807971 0.8807971 
0.958327 0.958327 0.958327 0.958327 
0.9849145 0.9849145 0.9849145 0.9849145 
0.9944816 0.9944816 0.9944816 0.9944816 
0.997974 0.997974 0.997974 0.997974 
0.9992552 0.9992552 0.9992552 0.9992552 
0.9997261 0.9997261 0.9997261 0.9997261 
0.9998992 0.9998992 0.9998992 0.9998992 
0.9999629 0.9999629 0.9999629 0.9999629 
0.9999864 0.9999864 0.9999864 0.9999864 

No drama!

I could easily refine the code to prevent subtractive cancelation when p is close to 1 by defining

w <- linkderiv <- plogis(eta) * plogis(eta, lower.tail=FALSE)

but the more readable code above already handles the example well enough.

Answer 2

In a regression framework, if the dependent variable shows no variation (say it is 5 for each observation), the regression produces results, it sets the constant to match the dependent variable (in this case 5), and the coefficients on other independent variables (if any) to 0. This happens in the glm case too. The high value of the constant is there so that the resulting predicted value is 1 (to match the dependent).

It is not an error that the dependent variable shows no variation. Estimation routines work. The issue is not with the estimation techniques and with the algorithms. The real question is what do you want to explain if there is no variation to explain?

I hope it helps