4
$\begingroup$

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)

enter image description here

enter image description here

> 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)

enter image description here

enter image description here

> 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

  1. Which method should be preferred?
  2. 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))
$\endgroup$
3
  • $\begingroup$ It looks like glmmTMB might not be handling weights 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$
    – amoeba
    Commented Apr 30, 2019 at 11:18
  • $\begingroup$ Thanks. I'm not very familiar with nonlinear regression - any suggestions for packages/resources to get started on this? $\endgroup$
    – Joris
    Commented Apr 30, 2019 at 11:33
  • $\begingroup$ It's still linear. You can start with glm and include sin(month/12 * pi) and cos(month/12 * pi) as fixed effects. More fancily one would use cyclic splines and gam. But I am not an expert in any of that. $\endgroup$
    – amoeba
    Commented Apr 30, 2019 at 13:49

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.