I am fitting a GLMM to test if parasite prevalence in snails (positive snails divided by total snails) differs between different sites (site_type
). Sites were sampled repeatedly over several months, so month
was included as a random effect.
I have noticed that using either lme4
or glmmTMB
provides different model fits (visualized and tested with the DHARMa
package).
glmmTMB
glmmTMB_model <- glmmTMB::glmmTMB(BT_pos_tot/BT_tot ~ (1|month) + site_type,
weights = BT_tot,
data= df,
family= binomial)
sim_residuals_glmmTMB <- DHARMa::simulateResiduals(glmmTMB_model, 1000)
plot(sim_residuals_glmmTMB)
DHARMa::testDispersion(sim_residuals_glmmTMB)
> DHARMa::testDispersion(sim_residuals_glmmTMB)
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
ratioObsSim = 0.042, p-value <0.0000000000000002
alternative hypothesis: two.sided
lme4
lme4_model <- lme4::glmer(BT_pos_tot/BT_tot ~ (1|month) + site_type,
weights = BT_tot,
data= df,
family= binomial)
sim_residuals_lme4 <- DHARMa::simulateResiduals(lme4_model, 1000)
plot(sim_residuals_lme4)
DHARMa::testDispersion(sim_residuals_lme4)
> DHARMa::testDispersion(sim_residuals_lme4)
DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
data: simulationOutput
ratioObsSim = 1.9, p-value = 0.05
alternative hypothesis: two.sided
Questions
- Which method should be preferred?
- Why is the model fit so different?
The data
> dput(df)
structure(list(BT_pos_tot = c(0, 0, 2, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 2, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0, 0, 5, 0, 2, 0, 0, 0, 0,
4, 9, 0, 0, 0, 0, 4, 0, 0, 0, 0, 1, 0, 0, 0, 5, 11, 0, 0, 1,
0, 0, 2, 0, 0, 0, 0, 4, 0, 0, 0, 25, 0, 1, 0, 2, 0, 0, 0, 4,
0, 1, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 21,
34, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 8, 0, 0, 0, 19, 0, 0, 0, 0,
0, 0, 0, 0, 0, 2, 0, 1, 2, 5, 9, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 5, 2, 0, 2, 2, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 17, 9, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 7, 4, 0,
0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 4, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,
1, 1, 7, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 6, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
3, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0), BT_tot = c(9, 129, 367,
0, 0, 20, 61, 21, 21, 0, 0, 0, 16, 17, 250, 27, 33, 0, 0, 0,
0, 0, 101, 93, 0, 213, 1123, 0, 1, 48, 0, 80, 0, 1, 837, 13,
5, 0, 13, 0, 36, 8, 0, 0, 0, 0, 105, 95, 157, 11, 0, 6, 13, 37,
154, 0, 0, 0, 23, 116, 89, 65, 2, 223, 62, 63, 0, 0, 297, 120,
0, 566, 909, 12, 0, 170, 68, 12, 25, 3, 1260, 22, 29, 2, 85,
76, 80, 45, 3, 0, 0, 0, 71, 481, 696, 0, 12, 0, 11, 76, 1, 0,
0, 30, 2, 230, 140, 0, 592, 95, 0, 0, 0, 1456, 351, 62, 401,
1699, 0, 0, 297, 49, 47, 0, 6, 1141, 67, 59, 2, 97, 47, 6, 0,
8, 8, 4, 1, 6, 166, 103, 0, 86, 64, 108, 116, 103, 847, 31, 95,
28, 0, 16, 64, 1, 0, 0, 103, 167, 0, 0, 0, 0, 3, 0, 1, 0, 0,
0, 0, 203, 221, 0, 127, 261, 0, 127, 0, 759, 44, 0, 109, 1262,
0, 0, 1, 0, 0, 0, 13, 818, 86, 34, 0, 32, 60, 0, 0, 38, 11, 1,
0, 0, 462, 635, 1, 4, 0, 0, 3, 19, 1, 0, 0, 0, 0, 150, 4, 0,
229, 252, 0, 0, 1, 703, 119, 0, 10, 886, 0, 0, 0, 0, 3, 0, 11,
948, 103, 17, 1, 1, 0, 0, 0, 21, 0, 0, 0, 3, 240, 600, 1, 0,
0, 0, 3, 14, 0, 0, 0, 21, 1, 37, 0, 21, 263, 78, 0, 0, 0, 799,
82, 42, 4, 53, 0, 0, 110, 0, 0, 0, 0, 679, 136, 6, 0, 0, 0, 0,
0, 2, 0, 0, 0, 85, 137, 34, 0, 9, 0, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0, 16, 13, 0, 0, 0, 550, 37, 0, 82, 33, 0, 1, 26, 0, 14,
0, 0, 412, 48, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 28, 196, 493, 0,
0, 0, 0, 0, 3, 2, 0, 0, 18, 2, 51, 7, 1, 6, 77, 0, 0, 0, 85,
4, 0, 480, 250, 0, 0, 2, 1, 6, 0, 0, 219, 38, 0, 0, 0, 0, 0,
7, 34, 0, 0, 0, 64, 202, 223, 7, 7, 0, 3, 14, 18, 0, 0, 0, 4,
64, 30, 2, 0, 20, 113, 5, 0, 0, 375, 42, 11, 160, 487, 0, 0,
25, 1, 7, 0, 0, 325, 44, 31, 0, 0, 10, 9, 5, 0, 0, 0, 0, 8, 78,
219, 0, 0, 0, 17, 1, 18, 0, 0, 0, 2, 21, 20, 9, 7, 59, 2, 0,
0, 0, 208, 7, 0, 187, 747, 0, 0, 3, 0, 3, 1, 0, 853, 19, 256,
0, 1, 13, 0, 1, 0, 0, 0, 0, 40, 463, 149, 2, 0, 7, 21, 11, 58,
0, 0, 0, 15, 36, 173, 0, 11, 46, 57, 0, 0, 0, 114, 56, 1, 79,
363, 0, 0, 9, 0, 1, 0, 0, 736, 30, 369, 29, 3, 0, 25, 2, 0, 0,
0, 0), month = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 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, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L,
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L,
10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L,
11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L), .Label = c("1", "2", "3", "4", "5", "6", "7", "8",
"9", "10", "11", "12"), class = "factor"), site_type = structure(c(6L,
1L, 2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L,
1L, 2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L,
7L, 1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L,
1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L,
3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L,
2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L, 4L,
5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L,
1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L, 1L, 2L, 1L, 2L,
1L, 2L, 7L, 1L, 1L, 2L, 3L, 4L, 2L, 4L, 6L, 1L, 2L, 1L, 2L, 4L,
5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L,
1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L,
3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L, 2L, 4L, 5L, 6L,
3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L, 1L, 2L,
4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L, 3L, 4L,
5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L,
2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L,
2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L,
1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L,
5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L,
2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L,
4L, 7L, 6L, 1L, 2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L,
1L, 2L, 3L, 1L, 2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L,
5L, 1L, 5L, 7L, 1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L,
6L, 1L, 2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L,
3L, 1L, 2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L,
5L, 7L, 1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L,
2L, 1L, 2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L,
2L, 3L, 4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L,
1L, 2L, 1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L, 6L, 1L, 2L, 1L,
2L, 4L, 5L, 6L, 3L, 1L, 2L, 4L, 5L, 5L, 1L, 2L, 3L, 1L, 2L, 3L,
4L, 7L, 1L, 2L, 4L, 1L, 2L, 1L, 2L, 3L, 5L, 1L, 5L, 7L, 1L, 2L,
1L, 2L, 3L, 4L, 5L, 6L, 1L, 2L, 4L, 7L), .Label = c("can.2",
"canal.3", "pond", "rice.p", "river", "rivulet", "spillway",
"stream"), class = "factor")), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -523L))
glmmTMB
might not be handlingweights
correctly. However, I'm not sure that modeling month as a random effect is very meaningful. I think it should rather be fixed, ideally treated as a circular variable using polynomial terms or splines. $\endgroup$glm
and includesin(month/12 * pi)
andcos(month/12 * pi)
as fixed effects. More fancily one would use cyclic splines andgam
. But I am not an expert in any of that. $\endgroup$