# Singular fit with simplest random structure in glmer (lme4)?

I am trying to run mixed models (logistic regression) on a dataframe with the glmer function from lme4 but I always receive this message: "boundary (singular) fit: see ?isSingular"

Even if I create a model with just an intercept and and the simplest random part (random intercept for one factor), the variance for this random factor is 0.

 Family: binomial  ( logit )
Formula: PointGagneparleServeur ~ 1 + (1 | Tour)
Data: DataModel_Logit_allRF_AusOpen
AIC       BIC    logLik  deviance  df.resid
480.7822  488.5765 -238.3911  476.7822       362
Random effects:
Groups Name        Std.Dev.
Tour   (Intercept) 0
Number of obs: 364, groups:  Tour, 6
Fixed Effects:
(Intercept)
0.5639
convergence code 0; 1 optimizer warnings; 0 lme4 warnings


Though I have observations for all the values of the factor :

table(DataModel_Logit_allRF_AusOpen$$PointGagneparleServeur,DataModel_Logit_allRF_AusOpen$$Tour)

1erTour 2emeTour 3emeTour 8eme Quart Demi
0      26       24       12   35    20   15
1      40       36       37   59    32   28



and the dependent variable PointGagneparleServeur is actually numeric.

(FYI, i recently "upgraded" my os to Catalina 10.15. Experiencing several bugs with other(non programming) softwares since. So, I am mentionning it just in case it could play a role...)

Does anyone have an idea on why I have this issue ?

• Please can you post the output from glm() (ie without random effects) Jan 12 '20 at 19:56
• please see below the comment of JTH Jan 13 '20 at 14:58

Indeed, your data do not support that hypothesis that there is significant variation in the outcome between the levels of the grouping factor.

library(lme4)
library(tidyverse)

dat <- data.frame(f = rep(letters[1:6], c(26 + 40, 24 + 36, 12 + 37, 35 + 59, 20 + 32, 15 + 28)),
y = rep(c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1), c(26, 40, 24, 36, 12, 37, 35, 59, 20, 32, 15, 28)))

fit <- glmer(y ~ 1 + (1 | f), family = binomial, data = dat)

summary(fit)

> summary(fit)
Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: y ~ 1 + (1 | f)
Data: dat

AIC      BIC   logLik deviance df.resid
480.8    488.6   -238.4    476.8      362

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.3257 -1.3257  0.7543  0.7543  0.7543

Random effects:
Groups Name        Variance Std.Dev.
f      (Intercept) 0        0
Number of obs: 364, groups:  f, 6

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.5639     0.1090   5.173 2.31e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
convergence code: 0
boundary (singular) fit: see ?isSingular


This is because you don't have enough data to conclude that the groups significantly differ with respect to the outcome measure:

> dat %>% group_by(f) %>% summarize(m = mean(y))
# A tibble: 6 x 2
f         m
<fct> <dbl>
1 a     0.606
2 b     0.6
3 c     0.755
4 d     0.628
5 e     0.615
6 f     0.651


Notice all the means are fairly close together. Given the total number of observations, and the number of observations in each group, and most critically, the number of groups (only 6) it's hard to say variation exists. I recommend looking into rstarnarm or brms to fit the analogous model from a fully Bayesian point of view. (However, this is my opinion, and there is by no means consensus in the statistical community on how to deal with this problem, see ?isSingular)

library(rstanarm)
fit2 <- stan_glmer(y ~ 1 + (1 | f), family = binomial, data = dat)
summary(fit2)

> fit2
stan_glmer
family:       binomial [logit]
formula:      y ~ 1 + (1 | f)
observations: 364
------
(Intercept) 0.6    0.1

Error terms:
Groups Name        Std.Dev.
f      (Intercept) 0.24
Num. levels: f 6

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg


• Thanks for the links. Still there are several things I do not understand : 1) Ok there is no significant differences between the groups. Still, given that we do observe differences, would it be not more likely that the sigma differs from 0 ? 2) I don't know exactly what are the calculations involved in the full Bayesian models, but if we have flat priors on the parameters, wouldn't the outputs of that kind of model be exactly the same as the model that I used, considering that the latter is supposed to find the values of the parameters by maximizing the likelihood ? Jan 13 '20 at 15:07
• 1. Yes, it's more reasonable that $\sigma_f$ is non-zero. Unfortunately, in your glmm, the estimate $\hat\sigma_f$ that jointy maximizes the likelihood is equal to zero according to glmer. 2. To your second question, not exactly. This depends on the parameterization of $\sigma_f$ on which the prior is flat (could be $\sigma_f$, $\sigma_f^2$, $\log(\sigma_f)$, etc..). In one case, only the MAPs will match the output of glmer, but I wouldn't take the MAPs as parameter estimates. Rather, consider the means or medians of your posterior distributions.
– JTH
Jan 14 '20 at 1:20
• Lastly, not that I used the default priors for stan_glmer above. You will want to consider these defaults before adopting them.
– JTH
Jan 14 '20 at 1:25

As noted in the comments of the other responders you have a quite small dataset, which makes fitting a mixed model tricky.

In general, you could give a try to different optimization algorithms, and altering the defaults. For example, the simple random intercepts model seems to converge with GLMMadaptive when you increase the number of EM iterations, i.e.,

library("GLMMadaptive")

dat <- data.frame(f = rep(letters[1:6], c(26 + 40, 24 + 36, 12 + 37,
35 + 59, 20 + 32, 15 + 28)),
y = rep(c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1),
c(26, 40, 24, 36, 12, 37, 35, 59, 20, 32, 15, 28)))

gm <- mixed_model(y ~ 1, random = ~ 1 | f, family = binomial(),
data = dat, iter_EM = 100)

summary(gm)