'Interpreting DHARMa residuals for a glmer.nb regression using count data
1
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))
Solution 1:[1]
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.
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 | kjetil b halvorsen |