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?