GLMM. Why are my confidence intervals so wide when my p-values are so small?

Question

I am trying to find out whether there is a significant effect of Treatment (factor with 3 levels: GR, BC and WF) on the number of moths captured per night in a field experiment. I have 15 Blocks that each contain all three Treatment levels. I sampled four Blocks per night over 65 nights. I have run a generalised linear mixed effects model with a negative binomial error structure and a partially-crossed random effects structure using lme4 in R. I'm pretty sure I have the random effects structure right, I asked a question about it here. I did all my residual diagnostic checks and the model fit seems fine. A likelihood ratio test comparing the full model with the null model (no Treatment effect) suggested that Treatment was highly significant. This was corroborated by a parametric bootstrap test as well as a comparison of the AICcs. I also looked at the confidence intervals of the parameters using confint(Model) and saw that the 95% CIs for the Treatment "WF" did not overlap zero in relation to the baseline Treatment “GR”. However, when I plot my model predictions with 95% confidence intervals, the CIs are strongly overlapping.

I know that overlapping CIs do not always imply statistical non-significance, but the size of the discrepancy concerns me. I.e., the p-values suggest a highly significant effect (p = 1.987e-10) but the CIs of the predicted values overlap the predicted means of each other. I would not be comfortable showing this plot in a paper and claiming the significance of Treatment effect was p < 0.0001, as it just seems implausible from the plot! To get the CIs, I used Ben Bolker’s function I found here and also tried using bootpredictlme4::bootMer with almost identical results.

I far as I can tell, there are four possibilities. (1) the model is specified wrong and the p-values are unreliable, (2) I have used an inappropriate way to determine the significance of the Treatment effect, (3) the model predictions that I've plotted are wrong, or (4) I don't understand the connection between statistical significance and confidence intervals in model predictions. Or possibly all four.

The R code I used is here...

# Starting with data frame "Moths"

# Relevel Treatment so GR comes first
Moths <- within(Moths, Treatment <- relevel(Treatment, ref = "GR"))
levels(Moths$Treatment) # "GR" "BC" "WF"

# Look at data
str(Moths)
#'data.frame':  711 obs. of  4 variables:
#  $ Abundance: int  16 4 29 47 16 0 12 2 3 0 ...
#$ Treatment: Factor w/ 3 levels "GR","BC","WF": 2 2 2 2 2 2 2 2 2 2 ...
#$ Night_re : Factor w/ 65 levels "Y2018_N06","Y2018_N07",..: 3 7 11 15 19 23 27 31 35 39 ...
#$ Block_re : Factor w/ 15 levels "M_1","M_10","M_11",..: 2 2 2 2 2 2 2 2 2 2 ...

head(Moths, 10)
# Abundance   Treatment  Night_re Block_re
#1         16        BC Y2018_N08     M_10
#2          4        BC Y2018_N12     M_10
#3         29        BC Y2018_N16     M_10
#4         47        BC Y2018_N20     M_10
#5         16        BC Y2018_N24     M_10
#6          0        BC Y2018_N28     M_10
#7         12        BC Y2018_N32     M_10
#8          2        BC Y2018_N36     M_10
#9          3        BC Y2019_N03     M_10
#10         0        BC Y2019_N07     M_10

# Model Abundance as function of Treatment with a random intercept for Night 
# and a random intercept for Block
Mod_1 <- glmer.nb(Abundance ~ Treatment + (1|Night_re) + (1|Block_re), data = Moths)

# Have a look at model
summary(Mod_1)
#Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#Family: Negative Binomial(5.6851)  ( log )
#Formula: Abundance ~ Treatment + (1 | Night_re) + (1 | Block_re)
#Data: Moths
#
#AIC      BIC   logLik deviance df.resid 
#4803.6   4831.0  -2395.8   4791.6      705 #
#
#Scaled residuals: 
#  Min      1Q  Median      3Q     Max 
#-1.9733 -0.6972 -0.1803  0.4955  4.4589 #

#Random effects:
#  Groups   Name        Variance Std.Dev.
#Night_re (Intercept) 1.33753  1.1565  
#Block_re (Intercept) 0.04765  0.2183  
#Number of obs: 711, groups:  Night_re, 65; Block_re, 15

#Fixed effects:
#  Estimate Std. Error z value Pr(>|z|)    
#(Intercept)  2.29096    0.15856  14.449  < 2e-16 ***
#TreatmentBC  0.05655    0.04897   1.155    0.248    
#TreatmentWF  0.30265    0.04858   6.229 4.68e-10 ***

# Check model assumptions
# First, with DHARMa
library(DHARMa)
sim <- simulateResiduals(fittedModel = Mod_1, n = 1000)
plot(sim) # Looks fine

# Check residuals against fitted values and covariates
modfort <- fortify(Mod_1)

# Residual vs fitted
ggplot(data = modfort, aes(x = .fitted, y= .scresid)) +
  geom_point() # Good

# Residuals vs covariate
ggplot(data = modfort, aes(x = Treatment, y= .scresid)) +
  geom_boxplot() +
  geom_point() # Good

# Check for overdispersion (using Ben Bolker's function)
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun(Mod_1) # = 0.93 Fine

# Testing significance of parameters

# Test whether Treatment is significant with a likelihood ratio test
drop1(Mod_1, test = "Chisq")
# Model:
# Abundance ~ Treatment + (1 | Night_re) + (1 | Block_re)
# Df    AIC    LRT   Pr(Chi)    
# <none>       4803.6                     
# Treatment  2 4844.3 44.679 1.987e-10 ***

# Test with a parametric bootstrap

# Make a null model
Null_mod_1 <- update(Mod_1,~.-Treatment)

# Compare with parametric bootstrapping
pbkrtest::PBmodcomp.merMod(Mod_1, Null_mod_1, nsim = 1000)
#Parametric bootstrap test; time: 1240.38 sec; samples: 1000 extremes: 0;
#large : Abundance ~ Treatment + (1 | Night_re) + (1 | Block_re)
#small : Abundance ~ (1 | Night_re) + (1 | Block_re)
#stat df   p.value    
#LRT    44.679  2 1.987e-10 ***
#PBtest 44.679     0.000999 ***

# Could be smaller than 0.000999 if nsim was larger

# Compare with AICc
AICc(Mod_1, Null_mod_1)
#           df   AICc
#Mod_1       6 4803.768
#Null_mod_1  4 4844.384

# What are the confidence intervals of the parameters?
confint(Mod_1, oldNames = FALSE)
#Computing profile confidence intervals ...
#                            2.5 %    97.5 %
#sd_(Intercept)|Night_re  0.97646643 1.3952734
#sd_(Intercept)|Block_re  0.13976365 0.3636530
#(Intercept)              1.97519988 2.6046183
#TreatmentBC             -0.03940728 0.1525319
#TreatmentWF              0.20746858 0.3979320

# Everything suggests treatment effect is highly significant

# Plot model predictions

# Using Ben Bolker's function for calculating fixed effect CIs for mixed models

easyPredCI <- function(model,newdata=NULL,alpha=0.05) {
  ## baseline prediction, on the linear predictor (logit) scale:
  pred0 <- predict(model,re.form=NA,newdata=newdata)
  ## fixed-effects model matrix for new data
  X <- model.matrix(formula(model,fixed.only=TRUE)[-2],newdata)
  beta <- fixef(model) ## fixed-effects coefficients
  V <- vcov(model)     ## variance-covariance matrix of beta
  pred.se <- sqrt(diag(X %*% V %*% t(X))) ## std errors of predictions
  ## inverse-link function
  linkinv <- family(model)$linkinv
  ## construct 95% Normal CIs on the link scale and
  ##  transform back to the response (probability) scale:
  crit <- -qnorm(alpha/2)
  linkinv(cbind(conf.low=pred0-crit*pred.se,
                conf.high=pred0+crit*pred.se))
}

# First, make a new data frame with the three Treatment levels
New_X <- expand.grid(Treatment = c("GR", "BC", "WF"))

# Get the predicted responses.
New_Y <- predict(Mod_1, newdata = New_X, type = "response", re.form = NA)

# Add together into a single data frame. 
New_df <- data.frame(New_X, New_Y)

# Get the 95% confidence intervals
cpred1.CI <- easyPredCI(model = Mod_1, newdata = New_X)

# Add it to the data frame:
New_df <- data.frame(New_df, cpred1.CI)

# Plot it
ggplot(data = New_df, aes(x = Treatment, y = New_Y)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2)

# Not very convincing!