Factor included based on AIC from anova, yet no significant comparisons using PostHoc

Question

Using step-wise model reduction I have reduced my full model down (Based on AIC and model comparisons using the anova() function in R. This resulted in the following model:

m1_DD4 <- glm(GotoPB ~ Pb_type + ZF_Pb + ZF_NotPb, data = DF_DD_5, family = binomial(link = "logit"))

Pb_type is a factor with 4 levels. ZF_Pb and ZF_NotPb are numerical. The summary of this model looks as follows:

Call:
glm(formula = GotoPB ~ Pb_type + ZF_Pb + ZF_NotPb, family = binomial(link = "logit"), 
    data = DF_DD_5)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.15642  -1.07365  -0.00006   1.02744   1.56573  

Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.013410   0.136495   0.098  0.92174    
Pb_typeNG        0.152678   0.163990   0.931  0.35184    
Pb_typeSilence -20.152316 452.945271  -0.044  0.96451    
Pb_typeSong      0.221061   0.156019   1.417  0.15652    
ZF_Pb            0.079849   0.013634   5.857 4.72e-09 ***
ZF_NotPb        -0.027021   0.007097  -3.807  0.00014 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2047.0  on 1559  degrees of freedom
Residual deviance: 1385.9  on 1554  degrees of freedom
AIC: 1397.9

Number of Fisher Scoring iterations: 18

From the model reduction/AIC I conclude that the model fit improves when Pb_type is included, and it should therefore be included. Yet here all are shown as non-significant. From what I understand/read, this is because the p-values reported here are from Wald tests, which tests if each coefficient is different from 0. Whereas the anova() function uses likelihood ratio tests, which indicates which model has a statistically better fit. Please correct me if I am wrong here.

I have also been told to not worry about this too much, and to rely on the outcome of anova(). But I would like to undestand why.

Under the assumption that Pb_type indeed matters, I now wish to determine how the different levels of the Pb_type differ when scored against each other using a PostHoc pair-wise comparison. I do this using the emmeans package, via the following formula:

emmeans::emmeans(m1_DD4, pairwise ~ Pb_type, adjust = "Tukey")

Which results in the following output:

$emmeans
 Pb_type  emmean      SE  df asymp.LCL asymp.UCL
 DC        0.128   0.120 Inf -1.72e-01     0.427
 NG        0.280   0.115 Inf -5.30e-03     0.566
 Silence -20.025 452.945 Inf -1.15e+03  1108.223
 Song      0.349   0.104 Inf  8.96e-02     0.608

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 

$contrasts
 contrast       estimate      SE  df z.ratio p.value
 DC - NG         -0.1527   0.164 Inf -0.931  0.7882 
 DC - Silence    20.1523 452.945 Inf  0.044  1.0000 
 DC - Song       -0.2211   0.156 Inf -1.417  0.4887 
 NG - Silence    20.3050 452.945 Inf  0.045  1.0000 
 NG - Song       -0.0684   0.152 Inf -0.449  0.9698 
 Silence - Song -20.3734 452.945 Inf -0.045  1.0000 

Results are given on the log odds ratio (not the response) scale. 
P value adjustment: tukey method for comparing a family of 4 estimates 

Here a few things stand out:

1. My df = INF. From here, I undestand that this presents no issue and that this is how emmeans labels asymptotic results.

2. In both my original model summary and here, it can be seen that for Silence there is a huge Standard error.

3. None of the values are significant.

This leaves me with 2 questions:

1. Why is my standard error so high, and does this present a problem? To answer this it might be necessary to know more about my data. So in short: All Pb_type treatments are different playbacks, and the silence served as a control treatment. The response variable GotoPB is either Y/N and represents whether individuals approached the playback (Y) or not (N). In the case of silence, as there never truly was playback, all values for GotoPB are N.
2. How can none of the comparisons of Pb_type be significant, yet including it in the model results in a significantly better fit (with a much lower AIC). Is something wrong here? And if not, how should I interpret this?

Kind regards.

Answer 1

First, the standard error for the Silence condition is huge because you are trying to estimate $-\infty$; i.e., the logit of zero. As you explained, the response is defined to always be N under silence. Therefore, Silence is not an experimental condition, at least for this response, because you don't observe a response under that condition.

You should omit or subset-out all the Silence results from the dataset and re-fit the model without Silence as a level of Pb_type. Your AIC results, etc. may well be quite different, because all those zeros have destabilized everything.

Second, just because a term is selected in a model does not guarantee that any comparison of its levels is statistically significant. Keep in mind that the primary goal of fitting a statistical model is to understand the patterns. I suggest plotting the estimates, e.g. using emmip(), and generally being more descriptive about your findings. Consider including interactions among these factors. Put less emphasis on how many asterisks you get, and more on the patterns.