I want to build a model to assess whether a species is declining in three different national parks.
My dependent variable is count data of the species and I have date, park, season and food as predictors. I also have the transect site nested within park as a random effect. Finally, I include the length of the transect as an offset.
This is all done in a Bayesian framework using brms. My problem is that there seems to be some temporal auto-correlation present despite including the date which you can see in the plot below (each row represents a chain in the sampler).
This response suggests that I'm capturing trend with my date variable but not the non-trend autocorrelation.
How would I capture this remaining correlation?
Full data and model here for completeness. Note that I converted my dates to numeric values and then z-scored them.
Data:
df <- structure(list(AWB = c(15, 66, 7, 44, 22, 45, 14, 33, 60, 30,
32, 39, 45, 39, 0, 37, 24, 18, 37, 66, 60, 25, 3, 34, 13, 58,
6, 38, 33, 12, 0, 34, 20, 75, 2, 18, 4, 15, 9, 4, 0, 50, 21,
21, 24, 5, 9, 13, 87, 43, 1, 13, 19, 28, 1, 13, 18, 56, 42, 53,
2, 16, 37, 51, 79, 5, 49, 11, 2, 34, 91, 15, 30, 3, 0, 57, 5,
31, 18, 56, 14, 72, 94, 35, 45, 10, 8, 29, 33, 34, 53, 8, 24,
54, 21, 5, 11, 83, 27), ndate = c(-2.03388968672155, -2.04126857410741,
-2.04274435158459, -2.03241390924438, -2.0309381317672, -1.42144203369491,
-1.40668425892318, -1.40963581387753, -1.42291781117208, -1.41849047874056,
-1.41996625621774, -1.18384185987012, -1.18679341482446, -1.18531763734729,
-1.18826919230163, -1.15580208780384, -1.15727786528101, -0.916726136501869,
-0.927056578842078, -0.92853235631925, -0.918201913979042, -0.904919916684488,
-0.906395694161661, -0.583200426660855, -0.584676204138028, -0.58024887170651,
-0.553684877117402, -0.581724649183682, -0.347076030313234, -0.556636432071748,
-0.55811220954892, -0.345600252836062, -0.317560480769782, -0.319036258246954,
-0.347076030313234, -0.344124475358889, -0.295423818612192, -0.296899596089365,
-0.298375373566537, -0.29985115104371, -0.29394804113502, -0.177361620438382,
-0.178837397915554, -0.141942960986239, -0.143418738463411, -0.138991406031893,
-0.140467183509066, -0.118330521351477, -0.119806298828649, -0.116854743874304,
-0.0592994222645715, 0.0233441164570958, 0.0218683389799232,
0.0262956714114411, 0.0248198939342684, 0.257992735327544, 0.241759183078645,
0.240283405601472, 0.256516957850371, 0.260944290281889, 0.259468512804717,
0.396715818181771, 0.398191595658944, 0.404094705567634, 0.40704626052198,
0.476407801949093, 0.477883579426266, 0.479359356903438, 0.712532198296714,
0.573809115442487, 0.575284892919659, 0.714007975773887, 0.711056420819541,
0.715483753251059, 0.712532198296714, 0.770087519906446, 0.771563297383619,
0.942753484735644, 0.941277707258471, 0.961938591938888, 0.958987036984543,
0.963414369416061, 0.992929918959514, 0.991454141482341, 1.00768769373124,
0.994405696436686, 1.01063924868559, 1.16707166126588, 1.17887788108326,
1.22610276035279, 1.47698493147214, 1.47255759904062, 1.50207314858407,
1.49764581615255, 1.96841883137062, 1.96694305389345, 1.96989460884779,
1.98907971605104, 1.98760393857387), NP = structure(c(2L, 1L,
1L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 3L, 3L, 2L,
2L, 2L, 2L, 3L, 3L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 3L, 3L, 3L,
2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 1L,
1L), .Label = c("Katavi", "Ruaha", "Selous"), class = "factor"),
Season = c("Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry",
"Dry", "Dry", "Dry", "Dry", "Wet", "Wet", "Wet", "Wet", "Wet",
"Wet", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Wet", "Wet",
"Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Dry", "Dry", "Dry",
"Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Wet", "Wet",
"Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Dry",
"Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry",
"Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Dry",
"Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry", "Wet", "Wet",
"Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Dry",
"Dry", "Dry", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet", "Wet",
"Wet", "Wet"), CarcassPres = structure(c(2L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L), .Label = c("0", "1"), class = "factor"),
StandardTransect = structure(c(1L, 3L, 4L, 5L, 7L, 1L, 3L,
4L, 5L, 6L, 8L, 1L, 5L, 6L, 8L, 3L, 4L, 1L, 3L, 4L, 5L, 6L,
8L, 6L, 8L, 1L, 1L, 5L, 5L, 6L, 8L, 1L, 3L, 4L, 5L, 6L, 1L,
5L, 6L, 8L, 8L, 3L, 4L, 1L, 5L, 6L, 8L, 1L, 5L, 6L, 8L, 1L,
5L, 6L, 8L, 1L, 3L, 4L, 5L, 6L, 8L, 3L, 4L, 10L, 2L, 1L,
5L, 6L, 8L, 2L, 10L, 1L, 5L, 6L, 8L, 3L, 4L, 3L, 4L, 2L,
9L, 10L, 1L, 5L, 5L, 7L, 7L, 2L, 10L, 9L, 3L, 4L, 2L, 10L,
1L, 5L, 7L, 3L, 4L), .Label = c("Jongomero", "Kidai", "LakeChada",
"LakeKatavi", "Lunda", "Magangwe", "Mbagi-Mdonya", "Mpululu",
"Msolwa", "Mtemere"), class = "factor"), Tlength = c(93,
86.7, 35.2, 75, 27.2, 93, 78.2, 35.2, 74.4, 45.8, 10.3, 93,
71, 45.8, 10.3, 63.9, 35.2, 93, 77.9, 35.2, 86.6, 45.8, 10.3,
45.8, 10.3, 93, 93, 68.9, 86.7, 45.8, 10.3, 93, 81.6, 35.2,
90.5, 45.8, 93, 88.2, 45.8, 10.3, 10.3, 64.6, 35.2, 93, 82.3,
45.8, 10.3, 93, 77.9, 45.8, 10.3, 93, 90.3, 45.8, 10.3, 93,
77.4, 35.2, 87.5, 45.8, 10.3, 66, 35.2, 71.2, 85.7, 93, 87.5,
45.8, 10.3, 85.5, 69.6, 93, 97.8, 45.8, 10.3, 86.6, 35.2,
71.9, 35.2, 88.5, 77.9, 80, 93, 85.2, 85.5, 56.1, 56.1, 97.6,
81.8, 79.7, 71.1, 35.2, 53.8, 81.9, 68.7, 86.5, 60.7, 56.6,
78.7)), row.names = c(NA, 99L), class = "data.frame")
Model:
library(tidyverse)
nbglm <-
brm(
AWB ~ ndate * NP + Season + CarcassPres
+ (1 | NP/StandardTransect) +
offset(log(Tlength)),
family = negbinomial,
data = df,
chains = 4,
control = list(adapt_delta = 0.99, max_treedepth = 15) # to correct for divergent transitions
)
Make the plot:
mcmc_plot(nbglm,
type = "acf")
