Strange output for pairwise comparisons on glm with quasi-binomial distribution

Question

I'm new to CrossValidated - I've read up on how to ask questions properly but sorry if I do anything slightly wrong.

My data is showing whether microplastics were present or absent in the gut of fish larvae.

> dput(data) structure(list(age = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("5", "20"), class = "factor"), concentration = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L), .Label = c("20", "200", "2000", "20000"), class = "factor"), replicate = c(1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L), present = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 0L, 2L, 2L, 7L, 7L, 6L, 7L, 0L, 1L, 0L, 0L, 2L, 0L, 0L, 2L, 7L, 3L, 8L, 0L, 11L, 16L, 17L, 19L), absent = c(20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 16L, 20L, 18L, 18L, 13L, 13L, 14L, 13L, 20L, 19L, 20L, 20L, 18L, 20L, 20L, 18L, 13L, 17L, 12L, 20L, 9L, 4L, 3L, 1L)), row.names = c(NA, -32L), class = "data.frame")

There are two main effects: concentration (the concentration of plastic the fish were exposed to) and age (of fish)

I am running a glm with quasibinomial distribution

model <- glm(cbind(present, absent) ~ age + concentration + age:concentration,
          family = quasibinomial(link = logit), data = data)

Firstly, I don't understand why the coefficients are listed like that in those strange combinations? I.e. missing some of the factor levels and combining some

glm(formula = cbind(present, absent) ~ age + concentration + 
age:concentration, family = quasibinomial(link = logit), 
data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.1931  -0.7093  -0.0001   0.2812   2.0579  

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)
(Intercept)              -2.175e+01  4.269e+03  -0.005    0.996
age20                     1.738e+01  4.269e+03   0.004    0.997
concentration200         -1.746e-10  6.037e+03   0.000    1.000
concentration2000         1.955e+01  4.269e+03   0.005    0.996
concentration20000        2.108e+01  4.269e+03   0.005    0.996
age20:concentration200    1.425e+00  6.037e+03   0.000    1.000
age20:concentration2000  -1.642e+01  4.269e+03  -0.004    0.997
age20:concentration20000 -1.540e+01  4.269e+03  -0.004    0.997

(Dispersion parameter for quasibinomial family taken to be 1.416656)

Null deviance: 296.406  on 31  degrees of freedom
Residual deviance:  40.676  on 24  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 18

I then run an ANOVA using the following code

anova(model, test = "F")

Which makes sense

Analysis of Deviance Table

Model: quasibinomial, link: logit

Response: cbind(present, absent)

Terms added sequentially (first to last)


                  Df Deviance Resid. Df Resid. Dev       F    Pr(>F)    
NULL                                 31    296.406                      
age                1   27.202        30    269.204 19.2016    0.0002 ***
concentration      3  223.765        27     45.439 52.6510 1.047e-10 ***
age:concentration  3    4.763        24     40.676  1.1207    0.3603    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The main problem I am having with this data is when I go to run post hoc pairwise comparisons.

summary(glht(model, mcp(concentration="Tukey")))
summary(glht(model, mcp(age="Tukey")))

These give me values which make no sense as there are very clear differences between the groups when the data is plotted. I have tried with emmeans and get the same strange results.

Does anyone know if I am going wrong somewhere? Thanks in advance! And please let me know if I haven't asked this question correctly.

Edit: These are the results of glht

For concentration:

summary(glht(model, mcp(concentration="Tukey")))


Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = cbind(present, absent) ~ age + concentration + 
    age:concentration, family = quasibinomial(link = logit), 
    data = data)

Linear Hypotheses:
                    Estimate Std. Error z value Pr(>|z|)  
200 - 20 == 0     -1.746e-10  6.037e+03   0.000   1.0000  
2000 - 20 == 0     1.955e+01  4.269e+03   0.005   1.0000  
20000 - 20 == 0    2.108e+01  4.269e+03   0.005   1.0000  
2000 - 200 == 0    1.955e+01  4.269e+03   0.005   1.0000  
20000 - 200 == 0   2.108e+01  4.269e+03   0.005   1.0000  
20000 - 2000 == 0  1.523e+00  5.253e-01   2.899   0.0135 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

And for age:

summary(glht(model, mcp(age="Tukey")))

Multiple Comparisons of Means: Tukey Contrasts


Fit: glm(formula = cbind(present, absent) ~ age + concentration + 
    age:concentration, family = quasibinomial(link = logit), 
    data = data)

Linear Hypotheses:
            Estimate Std. Error z value Pr(>|z|)
20 - 5 == 0    17.38    4268.50   0.004    0.997
(Adjusted p values reported -- single-step method)

The reason I find these results "strange" is that when plotted there are obvious significant differences (very small standard error) yet these do not show up using glht. When I see it plotted I just can't get my head around how they could NOT be significant?

Answer 1

@RussLenth and @GordonSmyth is spot on in pointing out you should not do a posthoc on the main effects when you fitted a model with the interaction. If you look at the first model with the interaction terms, all of them have very small effects and huge errors. Most likely you don't need them, and you are overfitting the model. We can look at the data:

library(ggplot2)
data$ratio = data$present/(data$present+data$absent)
ggplot(data,aes(x=age,y=ratio))+
geom_point(position=position_jitter(width=0.1,height=0)) + 
facet_wrap(~concentration,scale="free_y")

Take note the y-axis are on different scales, so you can see that the overall presence probability in 20 and 200 are quite low to start with (0-0.1) compared to 2000,20000. Indicating there is a strong concentration effect on the presence/absence. However it's quite hard to judge whether the effect of age is different in these compared to the other two concentrations.

So for your research question the following model should work:

fit = glm(cbind(present, absent) ~ age + concentration, family = quasibinomial, 
data = data)

anova(fit,test="F")
              Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
NULL                             31    296.406                     
age            1   27.202        30    269.204 18.380 0.0002063 ***
concentration  3  223.765        27     45.439 50.397  3.41e-11 ***

And if we look at the coefficients, I think it makes sense given what we see, in this model, concentration20 and age 5 are taken as a reference levels and the coefficients gives the log-odds ratio compared to these references. For example, concentration20000 gives a coefficient of 5.6 which indicates the log-odds of presence is 5.6882 more in this group, compared to concentration 20.. which is roughly what we saw in the plot (you can convert this to probability).

summary(model)

                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -6.2375     1.2564  -4.965 3.35e-05 ***
age20                1.6964     0.3454   4.911 3.87e-05 ***
concentration200     1.4145     1.3703   1.032   0.3111    
concentration2000    3.5130     1.2521   2.806   0.0092 ** 
concentration20000   5.6882     1.2480   4.558 1.00e-04 ***

You can now do the posthoc using glht to make the pairwise comparisons

library(multcomp)
summary(glht(model, mcp(concentration="Tukey")))
summary(glht(age, mcp(age="Tukey")))