我正在运行一个模型,它的响应变量是计数数据和非正常(非常右尾)有很多零。
我正在使用 glmmTMB 包并运行带有泊松分布的初始 glmm,其中 id(代码)和月份(montyear)作为随机因素。
m0 <- glmmTMB(deg~ SE_score + species + season + TL +
(1|code) + (1|monthyear), family=poisson(), data=node_dat)
我使用了check_overdispersion()
并且毫不奇怪地发现我在模型中存在过度分散。我使用负二项分布 (nbinom2) 运行了第二个模型,并使用 DHARMa 包绘制了 QQ 图和残差,但我仍然有很多分散。
m1 <- glmmTMB(deg~ SE_score + species + season + TL +
(1|code) + (1|monthyear), family=nbinom2(), data=node_dat)
res1 <- simulateResiduals(m1)
plot(res1)
我检查了零通货膨胀并且没有任何通货膨胀,并尝试了 nbinom1 系列,但在所有测试(KS、分散和异常值)中给出了更高的 AIC 和显着偏差。
我不确定下一个处理这个问题的方法是什么,也不确定为什么残差看起来很奇怪。
我的数据的一个子集如下
structure(list(station = structure(c(25L, 22L, 53L, 45L, 24L,
65L, 12L, 64L, 27L, 28L, 2L, 15L, 46L, 53L, 10L, 16L, 11L, 59L,
19L, 19L, 19L, 8L, 25L, 46L, 45L, 36L, 39L, 67L, 66L, 52L, 10L,
10L, 3L, 67L, 41L, 12L), .Label = c("AR01", "BE01", "BE02", "BEUWM01",
"BL01", "BL02", "GCB01", "GCB02", "GCB03", "GCB04", "NI01", "NI01b",
"NI02", "NI03", "PB01", "PB02", "PB03", "PB04", "PB05", "PB06",
"PB07", "PB08", "PB09", "PB10", "PB11", "PB12", "PB13", "PB14",
"PB15", "PB16", "PB17", "PB18", "PB19", "PB20", "PB21", "PB22",
"PB23", "PB24", "PB25", "PB26", "PB27", "PB28", "PB29", "PB30",
"PB4G01", "PB4G02", "PBUWM01", "PBUWM02", "SA01", "SA02", "SA02b",
"SA03", "SA04", "SA05", "SA06", "SA07", "SA11", "SA4G01", "SAUWM01",
"SB01", "SB02/AR02", "SB03/AR05", "SB04/AR06", "VB01", "VB02",
"VB03", "VB04"), class = "factor"), monthyear = structure(c(56L,
27L, 52L, 57L, 40L, 57L, 7L, 45L, 51L, 46L, 3L, 46L, 4L, 11L,
43L, 34L, 57L, 37L, 33L, 42L, 47L, 21L, 56L, 11L, 51L, 30L, 54L,
18L, 57L, 55L, 12L, 26L, 27L, 7L, 2L, 25L), .Label = c("2014/01",
"2014/02", "2014/03", "2014/04", "2014/05", "2014/06", "2014/07",
"2014/08", "2014/09", "2014/10", "2014/11", "2014/12", "2015/01",
"2015/02", "2015/03", "2015/04", "2015/05", "2015/06", "2015/07",
"2015/08", "2015/09", "2015/10", "2015/11", "2015/12", "2016/01",
"2016/02", "2016/03", "2016/04", "2016/05", "2016/06", "2016/07",
"2016/08", "2016/09", "2016/10", "2016/11", "2016/12", "2017/01",
"2017/02", "2017/03", "2017/04", "2017/05", "2017/06", "2017/07",
"2017/08", "2017/09", "2017/10", "2017/11", "2017/12", "2018/01",
"2018/02", "2018/03", "2018/04", "2018/05", "2018/06", "2018/07",
"2018/08", "2018/09", "2018/10", "2018/11", "2018/12"), class = "factor"),
deg = c(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),
gs = c(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),
btw = c(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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
code = structure(c(169L, 120L, 87L, 197L, 121L, 65L, 173L,
183L, 11L, 191L, 22L, 191L, 141L, 5L, 85L, 130L, 90L, 188L,
197L, 82L, 99L, 202L, 193L, 19L, 69L, 142L, 195L, 117L, 118L,
89L, 5L, 132L, 121L, 4L, 172L, 15L), .Label = c("2390", "13573",
"13574", "13575", "13576", "19318", "19319", "19321", "19322",
"19506", "19514", "19519", "19520", "19524", "25537", "25540",
"25541", "25543", "25546", "25549", "25552", "25553", "27583",
"27585", "27586", "27591", "27592", "27593", "27594", "27595",
"27596", "27597", "27600", "27601", "27605", "27607", "27608",
"27613", "27614", "27617", "27619", "27620", "27621", "27626",
"27627", "27629", "27630", "27631", "27632", "28608", "28611",
"28612", "28618", "28625", "28628", "28629", "28631", "28632",
"28633", "28638", "28641", "28644", "28662", "28672", "28674",
"52978", "54815", "54846", "54852", "54860", "54863", "54865",
"54866", "54868", "54877", "54882", "54883", "54884", "54886",
"54890", "54892", "54895", "54896", "54901", "54904", "54914",
"54919", "54920", "54922", "54925", "54931", "54932", "54938",
"54952", "54954", "54955", "54958", "54959", "54962", "59950",
"59953", "59954", "59955", "59957", "59958", "59959", "59960",
"59961", "59962", "59964", "59966", "59969", "59970", "59971",
"59972", "59973", "59975", "59976", "59979", "59981", "59988",
"2388", "12950", "12952", "12956", "12958", "12960", "12962",
"12964", "12966", "12968", "13577", "14203", "19320", "19523",
"25534", "25535", "25536", "25539", "25542", "25544", "25545",
"25547", "25548", "25550", "27584", "27588", "27589", "27590",
"27598", "27599", "27602", "27603", "27604", "27606", "27609",
"27610", "27611", "27615", "27616", "27618", "27622", "27624",
"27625", "28624", "28627", "28637", "28639", "28642", "28660",
"28670", "34176", "34177", "34178", "34179", "52975", "52977",
"54817", "54821", "54822", "54825", "54845", "54849", "54880",
"54887", "54889", "54893", "54898", "54899", "54905", "54911",
"54912", "54915", "54933", "54947", "54957", "54961", "59951",
"59963", "59968", "59978", "59991", "59992", "59993", "59994",
"59995"), class = "factor"), species = structure(c(2L, 1L,
1L, 2L, 1L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 1L, 2L, 1L,
2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 2L, 1L), .Label = c("Grey Reef Shark", "Silvertip Shark"
), class = "factor"), SE_score = c(0.13, 0.29, 0.08, 0.1,
0.14, 0.07, 0.18, 0.01, 0.12, 0.09, 0.09, 0.1, 0.18, 0.06,
0.08, 0.21, 0.09, 0.03, 0.16, 0.11, 0.15, 0.15, 0.13, 0.14,
0.22, 0.21, 0.08, 0.14, 0, 0.03, 0.1, 0.09, 0.09, 0, 0.31,
0.31), region = structure(c(6L, 6L, 7L, 6L, 6L, 9L, 5L, 9L,
6L, 6L, 2L, 6L, 6L, 7L, 4L, 6L, 5L, 7L, 6L, 6L, 6L, 4L, 6L,
6L, 6L, 6L, 6L, 9L, 9L, 7L, 4L, 4L, 2L, 9L, 6L, 5L), .Label = c("Acoustic Release",
"Benares", "Blenheim", "Grand Chagos Bank", "Nelsons Island",
"Peros Banhos", "Saloman", "Speakers Bank", "Victory Bank"
), class = "factor"), date = structure(c(1533078000, 1456790400,
1522537200, 1535756400, 1491001200, 1535756400, 1404169200,
1504220400, 1519862400, 1506812400, 1393632000, 1506812400,
1396306800, 1414800000, 1498863600, 1475276400, 1535756400,
1483228800, 1472684400, 1496271600, 1509494400, 1441062000,
1533078000, 1414800000, 1519862400, 1464735600, 1527807600,
1433113200, 1535756400, 1530399600, 1417392000, 1454284800,
1456790400, 1404169200, 1391212800, 1451606400), tzone = "", class = c("POSIXct",
"POSIXt")), month = c(8, 3, 4, 9, 4, 9, 7, 9, 3, 10, 3, 10,
4, 11, 7, 10, 9, 1, 9, 6, 11, 9, 8, 11, 3, 6, 6, 6, 9, 7,
12, 2, 3, 7, 2, 1), season = structure(c(1L, 2L, 1L, 1L,
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L,
2L, 2L), .Label = c("dry.season", "wet.season"), class = "factor"),
sex = structure(c(3L, 3L, 2L, 3L, 2L, 2L, 3L, 3L, 3L, 2L,
3L, 2L, 3L, 2L, 2L, 3L, 2L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L,
2L, 2L, 2L, 3L, 2L, 2L, 3L, 2L, 2L, 2L, 2L), .Label = c("",
"F", "M", "U"), class = "factor"), TL = c(103, 131, 118,
86, 112, 123, 115, 137, 70, 104, 131, 104, 134, 100, 122,
117, 135, 151, 86, 106, 117, 132, 139, 125, 114, 108, 144,
114, 106, 132, 100, 119, 112, 94, 140, 88)), row.names = c(NA,
-36L), na.action = structure(c(`107971` = 107971L), class = "omit"), class = "data.frame")