1
$\begingroup$

enter image description here1

I am modeling overdispersed count data (detection of species) in a GLMM to account for changes in the number of detections of the individual (response variable) to covid period, area (rural vs urban) and year. I modeled animal id a random effect and after running a few different models this regressions seemed like the best fit:

mod_6 <- glmer.nb(num_detections ~ year + 
    covid_period*area + (1|animal_id),   
    data=nursesharks_3)

when I evaluated the model using DHARMa I got the following graphs --- the QQ plot looks fine but the residuals look weird to me. I am wondering how to interpret this? I have seen that poisson/interger GLMMs often do not have good residual plots but I thought DHARMa was supposed to account for that.

Here is my data set:

nursesharks <- structure(list(animal_id = c("A69-1602-28094", 
    "A69-1602-28094", "A69-1602-28094", "A69-1602-28094", 
    "A69-1602-28095", "A69-1602-28095", "A69-1602-28095", 
    "A69-1602-28095", "A69-1602-28096", "A69-1602-28096", 
    "A69-1602-28096", "A69-1602-28096", "A69-1602-28096",  
    "A69-1602-28096", "A69-1602-28097""A69-1602-28097","A69-1602-28097", 
    "A69-1602-28098", "A69-1602-28098", "A69-1602-28098", 
    "A69-1602-28098", "A69-1602-28098", "A69-1602-28099", 
    "A69-1602-28099", "A69-1602-28099", "A69-1602-28099", 
    "A69-1602-28100", "A69-1602-28100", "A69-1602-28101", 
    "A69-1602-28101", "A69-1602-28101", "A69-1602-28101",  
    "A69-1602-28101", "A69-1602-28102", "A69-1602-28102", 
    "A69-1602-28102", "A69-1602-28102", "A69-1602-28102", 
    "A69-1602-28102", "A69-1602-28102", "A69-1602-28103", 
    "A69-1602-28103", "A69-1602-28103", "A69-1602-28103", 
    "A69-1602-28103", "A69-1602-28103", "A69-1602-28104", 
    "A69-1602-28104", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16327", 
    "A69-9001-16327", "A69-9001-16327", "A69-9001-16327", 
    "A69-9001-16327", "A69-9001-16327", "A69-9001-16327", 
    "A69-9001-16330", "A69-9001-16330", "A69-9001-16330", 
    "A69-9001-16330", "A69-9001-16330", "A69-9001-16331", 
    "A69-9001-18403", "A69-9001-18405", "A69-9001-18405", 
    "A69-9001-18405", "A69-9001-18405", "A69-9001-18405", 
    "A69-9001-18405", "A69-9001-18405", "A69-9001-18405", 
    "A69-9001-18405", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18416", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18416", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18416", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18416", "A69-9001-18416", "A69-9001-18420", 
    "A69-9001-18420", "A69-9001-18420", "A69-9001-18420", 
    "A69-9001-18420", "A69-9001-18420", "A69-9001-18420", 
    "A69-9001-18420", "A69-9001-18420", "A69-9001-18420", 
    "A69-9001-18422", "A69-9001-18422", "A69-9001-18422",  
    "A69-9001-18422", "A69-9001-18422", "A69-9001-18422", 
    "A69-9001-18422", "A69-9001-18422", "A69-9001-18422", 
    "A69-9001-18422", "A69-9001-18423", "A69-9001-18423", 
    "A69-9001-18423", "A69-9001-18423", "A69-9001-18423", 
    "A69-9001-18423", "A69-9001-18425", "A69-9001-18425", 
    "A69-9001-18425", "A69-9001-18425", "A69-9001-18425",  
    "A69-9001-18425", "A69-9001-18425", "A69-9001-18425", 
    "A69-9001-18425", "A69-9001-18425", "A69-9001-18425", 
    "A69-9001-20772", "A69-9001-20772"), 
    covid_period = structure(c(3L, 3L, 1L, 1L, 
     3L, 3L, 1L, 1L, 3L, 3L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 3L, 3L, 1L, 
     1L, 2L, 3L, 3L, 3L, 1L, 3L, 2L, 3L, 3L, 1L, 1L, 2L, 3L, 3L, 3L, 
     3L, 1L, 1L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 
     3L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 1L, 1L, 1L, 
     2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 3L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 
     1L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
     3L, 3L, 3L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 
     1L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 1L, 
     1L, 1L, 2L, 2L, 2L, 3L, 2L), 
     .Label = c("before", "during", "after"), 
     class = "factor"), area = c("rural", "rural", "rural", "rural", 
     "rural", "rural", "rural", "rural", "rural", "rural", "rural",    
     "rural", "rural", "rural", "rural", "rural", "rural", "rural", 
     "rural", "rural", "rural", "rural", "rural", "rural", "urban", 
     "rural", "urban", "urban", "rural", "rural", "rural", "rural", 
     "rural", "rural", "rural", "urban", "urban", "rural", "urban", 
     "urban", "rural", "rural", "rural", "rural", "rural", "rural", 
     "rural", "rural", "rural", "rural", "urban", "urban", "urban", 
     "rural", "rural", "urban", "urban", "urban", "urban", "rural", 
     "urban", "urban", "urban", "urban", "urban", "rural", "urban", 
     "urban", "rural", "urban", "urban", "rural", "urban", "urban", 
     "urban", "urban", "urban", "urban", "rural", "rural", "urban", 
     "urban", "urban", "urban", "urban", "urban", "urban", "rural", 
     "urban", "urban", "urban", "rural", "urban", "urban", "urban", 
     "rural", "rural", "urban", "urban", "urban", "urban", "urban", 
     "urban", "rural", "urban", "urban", "urban", "urban", "urban", 
     "urban", "rural", "urban", "urban", "rural", "urban", "urban", 
     "urban", "rural", "urban", "urban", "rural", "urban", "urban", 
     "urban", "urban", "urban", "rural", "urban", "urban", "urban", 
     "urban", "urban", "urban", "urban", "urban", "urban", "urban", 
     "rural", "urban"), 
     year = c(2, 3, 3, 4, 2, 3, 3, 4, 2, 3, 4, 
     2, 3, 4, 3, 2, 3, 2, 3, 3, 4, 3, 2, 3, 3, 4, 2, 3, 2, 3, 3, 4, 
     3, 2, 3, 2, 3, 4, 3, 3, 1, 2, 3, 2, 3, 4, 2, 3, 2, 3, 1, 2, 3, 
     2, 4, 1, 2, 3, 4, 2, 1, 2, 3, 4, 1, 1, 1, 2, 2, 1, 2, 2, 2, 3, 
     1, 3, 1, 1, 1, 2, 1, 2, 1, 2, 3, 2, 3, 2, 1, 2, 3, 4, 1, 2, 4, 
     2, 4, 1, 2, 3, 1, 2, 3, 3, 1, 2, 3, 1, 2, 3, 3, 1, 3, 3, 2, 3, 
     4, 3, 1, 3, 2, 1, 2, 3, 1, 4, 2, 1, 2, 3, 1, 2, 3, 4, 1, 2, 3, 
     1, 1), num_detections = c(1L, 4L, 1L, 8L, 5L, 5L, 1L, 1L, 29L, 
     22L, 2L, 2L, 16L, 5L, 1L, 2L, 9L, 64L, 19L, 6L, 7L, 17L, 3L, 
     25L, 1L, 3L, 2L, 3L, 19L, 39L, 8L, 13L, 6L, 4L, 22L, 20L, 6L, 
     3L, 14L, 2L, 4L, 9L, 13L, 2L, 6L, 4L, 1L, 1L, 12L, 6L, 214L, 
     165L, 251L, 5L, 2L, 16L, 116L, 115L, 12L, 48L, 55L, 36L, 44L, 
     2L, 94L, 1L, 7L, 7L, 3L, 8L, 2L, 2L, 2L, 2L, 1L, 1L, 56L, 1L, 
     1L, 3L, 57L, 63L, 10L, 13L, 3L, 6L, 2L, 9L, 55L, 41L, 1L, 2L, 
     5L, 37L, 10L, 5L, 1L, 12L, 18L, 11L, 17L, 22L, 3L, 1L, 7L, 12L, 
     11L, 1L, 15L, 22L, 3L, 26L, 14L, 2L, 8L, 3L, 4L, 2L, 1L, 3L, 
     1L, 1L, 3L, 1L, 1L, 1L, 3L, 54L, 50L, 31L, 11L, 10L, 22L, 9L, 
     6L, 15L, 28L, 1L, 3L)), 
     class = c("grouped_df", "tbl_df", "tbl", "data.frame"), 
     row.names = c(NA, -139L), groups = structure(list(
    animal_id = c("A69-1602-28094", "A69-1602-28094", "A69-1602-28095", 
    "A69-1602-28095", "A69-1602-28096", "A69-1602-28096",  
    "A69-1602-28096", 
    "A69-1602-28097", "A69-1602-28097", "A69-1602-28098",  
    "A69-1602-28098", 
    "A69-1602-28098", "A69-1602-28099", "A69-1602-28099", 
    "A69-1602-28099", 
    "A69-1602-28100", "A69-1602-28100", "A69-1602-28101", 
    "A69-1602-28101", 
    "A69-1602-28101", "A69-1602-28102", "A69-1602-28102", 
    "A69-1602-28102", 
    "A69-1602-28102", "A69-1602-28102", "A69-1602-28103", 
    "A69-1602-28103", 
    "A69-1602-28104", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16326", 
    "A69-9001-16326", "A69-9001-16326", "A69-9001-16326", 
    "A69-9001-16327", 
    "A69-9001-16327", "A69-9001-16327", "A69-9001-16327", 
    "A69-9001-16327", 
    "A69-9001-16330", "A69-9001-16330", "A69-9001-16330", 
    "A69-9001-16331", 
    "A69-9001-18403", "A69-9001-18405", "A69-9001-18405", 
    "A69-9001-18405", 
    "A69-9001-18405", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18416", 
    "A69-9001-18416", "A69-9001-18416", "A69-9001-18416", 
    "A69-9001-18420", 
    "A69-9001-18420", "A69-9001-18420", "A69-9001-18420", 
    "A69-9001-18422", 
    "A69-9001-18422", "A69-9001-18422", "A69-9001-18422", 
    "A69-9001-18422", 
    "A69-9001-18422", "A69-9001-18423", "A69-9001-18423", 
    "A69-9001-18423", 
    "A69-9001-18425", "A69-9001-18425", "A69-9001-18425", 
    "A69-9001-18425", 
    "A69-9001-20772", "A69-9001-20772"), covid_period = structure(c(1L, 
    3L, 1L, 3L, 1L, 2L, 3L, 1L, 2L, 1L, 2L, 3L, 1L, 3L, 3L, 2L, 
    3L, 1L, 2L, 3L, 1L, 1L, 2L, 3L, 3L, 1L, 3L, 3L, 1L, 1L, 2L, 
    2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 1L, 3L, 3L, 3L, 1L, 1L, 2L, 
    3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 3L, 1L, 1L, 2L, 
    2L, 3L, 3L, 1L, 3L, 3L, 1L, 2L, 3L, 3L, 2L, 3L), 
    .Label = c("before", 
    "during", "after"), class = "factor"), area = c("rural", 
    "rural", "rural", "rural", "rural", "rural", "rural", "rural", 
    "rural", "rural", "rural", "rural", "rural", "rural", "urban", 
    "urban", "urban", "rural", "rural", "rural", "rural", "urban", 
    "urban", "rural", "urban", "rural", "rural", "rural", "rural", 
    "urban", "rural", "urban", "rural", "urban", "rural", "urban", 
    "rural", "urban", "urban", "urban", "rural", "urban", "urban", 
    "urban", "urban", "urban", "rural", "urban", "rural", "urban", 
    "rural", "urban", "rural", "urban", "rural", "urban", "urban", 
    "urban", "rural", "urban", "rural", "urban", "rural", "urban", 
    "urban", "rural", "urban", "urban", "urban", "rural", "urban", 
    "urban", "rural"), .rows = structure(list(3:4, 1:2, 7:8, 
        5:6, 11L, 12:14, 9:10, 15L, 16:17, 20:21, 22L, 18:19, 
        26L, 23:24, 25L, 28L, 27L, 31:32, 33L, 29:30, 38L, 39L, 
        40L, 34:35, 36:37, 44:46, 41:43, 47:48, 54:55, 56:59, 
        60L, 61:64, 49:50, 51:53, 66L, 67:68, 69L, 70:71, 65L, 
        75:76, 72L, 73:74, 77L, 78L, 83:85, 86:87, 79:80, 81:82, 
        92L, 93:95, 96:97, 98:100, 88L, 89:91, 104L, 105:107, 
        108:110, 101:103, 114L, 115:117, 118L, 119:120, 111L, 
        112:113, 125:126, 121L, 122:124, 131:134, 135:137, 127L, 
        128:130, 139L, 138L), ptype = integer(0), 
        class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"), 
     row.names = c(NA, -73L), .drop = TRUE))
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

The right plot shows heteroskedasticity.

You can account for that by specifying the same NB model with glmmTMB, but adding a ~ disp term to model the dispersion of the NB as a function of the predictors (see help of glmmTMB). Creating DHARMa plots with residuals against predictors will help deciding on the appropriate form the the ~ disp model.

Improve this answer
$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.