> 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!