I have been trying to use GLM.jl to fit some logistic regression models and the fit failed. I don’t remember the exact error message but it was something about singularity.
That annoyed me for a couple of years but I never figured why because using RCall.jl was straight-forward enough and I just use it to fit the logistic regression.
Recently, I decided to look into the scikit-learn implementation and it occurred to me why that is the case!
The GLM.jl implementation using Julia’s powerful numerical abilities implemented the EM algorithms where as scikit just invoked a numerical solver that doesn’t use EM.
The EM must have some technical condition where the fit will fail while scikit-learn and R’s implementation will just treat it like any other numerical problem and try to find the best coefficients!
Admittedly the R’s coefficients for the problem I tried to solve contained some NA coefficients which I needed to discard but I actually prefer it to GLM.jl’s approach.
In Julia, is there an implementation of GLM that uses a solver like they do in R? Can’t seem to find one but I guess it’s easy to structure the problem as an optimisation problem and invoke a solver.
Not an expert and don’t have an answer but I wonder whether it has anything to do with precision/scale. I had a similar problem before between R and Matlab. Turns out the default scale of Matlab was higher than in R resulting in singularity in R and a result in Matlab.
Edit: Agree that Julia should aim to obtain the same results as R when it comes to statistical analysis. Currently this doesn’t seem to be the case for other models also
If you can give a reproducible example and the exact error message, peopel can give you more helpful advice. Without that, it’s hard for this thread to go anywhere productive. Please read: make it easier to help you
Both the R and the GLM.jl implementation of Generalized Linear Models use Iteratively Reweighted Least Squares (IRLS), not EM. The least squares part will fail if the coefficients are undefined due to a singular model matrix (i.e. the X matrix). A general optimizer is less likely to detect the singularity and, depending on the convergence criteria, may well declare convergence near the subspace of possible solutions.
Most analysts, including me, prefer to know when the model is computationally singular. GLM.jl uses a Cholesky factorization of X’WX to solve for the new coefficient vector at each iteration. I have an unregistered package at GitHub - dmbates/GLMMng.jl: Experimental version of GLM and GLMM fitting in Julia that uses a QR decomposition of √W * X, which will be slightly slower but better able to handle near singularity, If you want to try that I will add some documentation (right now it is a test-bed with perfunctory documentation). Or, if you could provide a MWE then we can check the condition number of the weighted model matrix.