2
$\begingroup$

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 ?

$\endgroup$
  • $\begingroup$ Please can you post the output from glm() (ie without random effects) $\endgroup$ – Robert Long Jan 12 at 19:56
  • $\begingroup$ please see below the comment of JTH $\endgroup$ – M. Vande Velde Jan 13 at 14:58
2
$\begingroup$

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
------
            Median MAD_SD
(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

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ 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 ? $\endgroup$ – M. Vande Velde Jan 13 at 15:07
  • 1
    $\begingroup$ 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. $\endgroup$ – JTH Jan 14 at 1:20
  • $\begingroup$ Lastly, not that I used the default priors for stan_glmer above. You will want to consider these defaults before adopting them. $\endgroup$ – JTH Jan 14 at 1:25
1
$\begingroup$

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)
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.