Summary output of binomial GLMM shows significant effects, but graph shows overlapping CI error bars?

Question

I have run this binomial GLMM:

m17 <- glmer(Attendance ~ Day + Sex + (1|AgeClass) + (1|Year) + (1|Plastic), data = test2, family = binomial(link = "logit"))

This is the summary output of the model:

Generalized linear mixed model fit by maximum likelihood (Laplace  Approximation)
 [glmerMod]
 Family: binomial  ( logit )
Formula: Attendance ~ Day + Sex + (1 | AgeClass) + (1 | Year) + (1 | Plastic)
   Data: test2

     AIC      BIC   logLik deviance df.resid 
 16138.4  16182.7  -8063.2  16126.4    11868 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.5394 -0.9688 -0.5740  0.9656  1.8762 

Random effects:
 Groups   Name        Variance  Std.Dev.
 Plastic  (Intercept) 0.1374817 0.3708  
 Year     (Intercept) 0.0531973 0.2306  
 AgeClass (Intercept) 0.0009734 0.0312  
Number of obs: 11874, groups:  Plastic, 237; Year, 4; AgeClass, 3

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.135928   0.137669   0.987  0.32347    
Day         -0.005972   0.002061  -2.897  0.00377 ** 
SexFemales  -0.248756   0.063937  -3.891    1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
           (Intr) Day   
Day        -0.378       
SexFemales -0.240  0.004

As evidenced by the p-values, both Day and Sex have a significant effect on attendance. However, when I use plot_model to graph it, I get this result:

As you can see, the error bars are nearly overlapping, or overlapping, the predicted probability values for each sex. I haven't had this issue before with other graphs I've made with plot_model, so I'm confused as to why I'm getting these contradictory results. Does it have to do with the conversion of log-odds beta coefficients to predicted probability values? Or, is there something wrong with my code below?

plot_model(m17, type = "pred" ,
                    terms = c("Day [all]", "Sex"), 
                    colors = c("blue", "red"), 
                    ci.lvl = .95,
                    legend.title = NA,
                    line.size = 1,
                    show.legend = FALSE) +
  labs(x = "Days Relative to Clutch Initiation", y = "P(Attend at night)", title = NULL) +
  scale_x_continuous(expand = c(0, 0), breaks = breaks, labels = labels) + 
  scale_y_continuous(limit = c(0, 1), expand = c(0, 0)) + 
  theme_classic() + 
  guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
  theme(axis.text=element_text(size=22, color = "black", family = "serif"), 
        axis.title=element_text(size=22, color = "black", family = "serif"), 
        axis.title.x = element_text(vjust=-0.5),
        axis.title.y = element_text(vjust=2.5),
        axis.ticks.length=unit(0.1,"inch"), 
        legend.key.size = unit(1, 'cm'), 
        legend.position = c(0.15, 0.92),
        legend.text = element_text(colour="black", size=20, family = "serif"),
        plot.margin = margin(1,.5,1,1, "cm"))

All help is appreciated, and please let me know if there is another place I should post this question!

Answer 1

Interpreting this type of plot takes some care, for a few reasons.

First, the simple relationship between 95% confidence intervals (CI) and $p < 0.05$ only holds for comparisons of the CI against fixed values. Non-overlap of 95% CI for two estimates from a model is much more stringent than $p < 0.05$ as a significance criterion for the difference between the estimates. This answer shows that such non-overlap is approximately equivalent to $p < 0.005$ in the context of a t-test. It's similarly too stringent in any comparison between two statistical estimates.

Second, the 95% CI shown on your plot include the joint error estimates for all of the coefficients involved, not just the error in the SexFemales coefficient. The calculations involve the variance-covariance matrix of all the coefficients involved in the modeled estimates, together with the formula for the variance of a weighted sum of variables.

To get a sense of what's going on, try constructing not only the point estimates but also the 95% CI at Day 11 using only the standard error of the SexFemales coefficient. Work first in the (log-odds) coefficient scale where that coefficient estimate has an asymptotic normal distribution, then transform the CI back to the probability scale. You will probably see that the CI you get are narrower than what's displayed on the graph. The plot involves uncertainties beyond that in the SexFemales coefficient.

Third, in some implementations, this type of plot includes a correction for multiple comparisons based on the number of values of Day that were used for the estimates. That correction will widen the 95% CI. A quick glance at the documentation for plot_model() (from the sjPlot package) suggests that might not be the case here, but you need to check the documentation carefully to know how the software you use handles multiple comparisons.