Missing standard errors and confidence for mixed model and ggeffects

Question

I'm trying to use ggeffects::ggpredict to make some effects plots for my model. I find that the standard errors and confidence limits are missing for many of the results. I can reproduce the problem with some simulated data. It seems specifically for observations where the standard error puts the predicted probability close to 0 or 1.

I tried to get predictions on the link scale to diagnose if it's a problem with the translation from link to response, but I don't believe this is supported by the package.

Any ideas how to address this? Many thanks.

library(tidyverse)
library(lme4)
library(ggeffects)

# number of simulated observations
n <- 1000

# simulated data with a numerical predictor x, factor predictor f, response y
# the simulated effects of x and f are somewhat weak compared to the noise, so expect high standard errors
df <- tibble(
  x = seq(-0.1, 0.1, length.out = n),
  g = floor(runif(n) * 3),
  f = letters[1 + g] %>% as.factor(),
  y = pracma::sigmoid(x + (runif(n) - 0.5) + 0.1 * (g - mean(g))),
  z = if_else(y > 0.5, "high", "low") %>% as.factor()
)

# glmer model
model <- glmer(z ~ x + (1 | f), data = df, family = binomial)
print(summary(model))
#> Generalized linear mixed model fit by maximum likelihood (Laplace
#>   Approximation) [glmerMod]
#>  Family: binomial  ( logit )
#> Formula: z ~ x + (1 | f)
#>    Data: df
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>   1373.0   1387.8   -683.5   1367.0      997 
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -1.3858 -0.9928  0.7317  0.9534  1.3600 
#> 
#> Random effects:
#>  Groups Name        Variance Std.Dev.
#>  f      (Intercept) 0.0337   0.1836  
#> Number of obs: 1000, groups:  f, 3
#> 
#> Fixed effects:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.02737    0.12380   0.221    0.825    
#> x           -4.48012    1.12066  -3.998 6.39e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Correlation of Fixed Effects:
#>   (Intr)
#> x -0.001

# missing standard errors
ggpredict(model, c("x", "f")) %>% print()
#> Data were 'prettified'. Consider using `terms="x [all]"` to get smooth plots.
#> # Predicted probabilities of z
#> 
#> # f = a
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.54, 0.69]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |             
#> 
#> # f = b
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.56, 0.67]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |             
#> 
#> # f = c
#> 
#>     x | Predicted |       95% CI
#> --------------------------------
#> -0.10 |      0.62 | [0.54, 0.69]
#>  0.00 |      0.51 |             
#>  0.10 |      0.40 |
ggpredict(model, c("x", "f")) %>% as_tibble() %>% print(n = 20)
#> Data were 'prettified'. Consider using `terms="x [all]"` to get smooth plots.
#> # A tibble: 9 x 6
#>       x predicted std.error conf.low conf.high group
#>   <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>
#> 1  -0.1     0.617     0.167    0.537     0.691 a    
#> 2  -0.1     0.617     0.124    0.558     0.672 b    
#> 3  -0.1     0.617     0.167    0.537     0.691 c    
#> 4   0       0.507    NA       NA        NA     a    
#> 5   0       0.507    NA       NA        NA     b    
#> 6   0       0.507    NA       NA        NA     c    
#> 7   0.1     0.396    NA       NA        NA     a    
#> 8   0.1     0.396    NA       NA        NA     b    
#> 9   0.1     0.396    NA       NA        NA     c

^{Created on 2022-04-12 by the reprex package (v2.0.1)}

Answer 1

I think this may be due to the singular model fit.

I dug down into the guts of the code as far as here, where there appears to be a mismatch between the dimensions of the covariance matrix of the predictions (3x3) and the number of predicted values (15).

I further suspect that the problem may happen here:

rows_to_keep <- as.numeric(rownames(unique(model_matrix_data[
             intersect(colnames(model_matrix_data), terms)])))

Perhaps the function is getting confused because the conditional modes/BLUPs for every group are the same (which will only be true, generically, when the random effects variance is zero) ... ?

This seems worth opening an issue on the ggeffects issues list ?