How do I add temporal autocorrelation structure to a brms mixed effects model?

Interfaces brms
, ,
1

I have a data set of counts of various bird species in three national parks. My aim is to see if they’re declining in any of those areas. The data is at a monthly resolution. I’m using brms to fit a hierarchical model (I give some more detail of the model below*). I fit this model already but never included an autoregressive term which I’d like to include to see how it changes the estimates. I’m just not sure how to write it for my regression. Hope you can help.

brm(
    AWB ~  
      ndate * NP + Season + CarcassPres +
      (1 | NP / StandardTransect) +
      offset(log(Tlength)),
    family = negbinomial,
    data = mydata,
    chains = 4)

Here’s my data:


mydata <- structure(list(AWB = c(7, 66, 15, 44, 22, 60, 45, 32, 30, 33, 
14, 0, 45, 39, 39, 24, 37, 66, 37, 60, 18, 3, 25, 13, 34, 38, 
58, 0, 12, 6, 33, 2, 34, 18, 75, 20, 4, 9, 15, 4, 0, 21, 50, 
24, 21, 9, 5, 87, 13, 43, 1, 19, 13, 1, 28, 56, 18, 42, 13, 2, 
53, 16, 37, 51, 79, 5, 49, 11, 34, 91, 30, 2, 0, 15, 3, 57, 5, 
18, 31, 14, 56, 72, 35, 94, 10, 45, 8, 29, 33, 34, 8, 53, 54, 
24, 5, 21, 11, 27, 83, 23, 24, 4, 10, 13, 17, 11, 51, 6), ndate = c(0, 
0, 0, 0, 0, 14, 14, 14, 14, 14, 14, 19, 19, 19, 19, 20, 20, 25, 
25, 25, 25, 25, 25, 32, 32, 33, 33, 33, 33, 33, 38, 38, 38, 38, 
38, 38, 39, 39, 39, 39, 39, 41, 41, 42, 42, 42, 42, 43, 43, 43, 
44, 46, 46, 46, 46, 51, 51, 51, 51, 51, 51, 54, 54, 54, 55, 56, 
56, 56, 58, 58, 61, 61, 61, 61, 61, 63, 63, 66, 66, 67, 67, 67, 
68, 68, 68, 68, 68, 71, 72, 73, 78, 78, 79, 79, 89, 89, 89, 90, 
90, 91, 91, 91, 96, 96, 96, 96, 97, 97), nyear = c(2013, 2013, 
2013, 2013, 2013, 2014, 2014, 2014, 2014, 2014, 2014, 2015, 2015, 
2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2015, 2016, 
2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 
2016, 2016, 2016, 2016, 2016, 2016, 2017, 2017, 2017, 2017, 2017, 
2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 
2017, 2017, 2017, 2017, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 
2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2019, 2019, 
2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 
2020, 2020, 2020, 2020, 2021, 2021, 2021, 2021, 2021, 2021, 2021, 
2021, 2021, 2021, 2021, 2021, 2021, 2021), NP = structure(c(1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 3L, 3L, 
2L, 2L, 2L, 3L, 3L, 2L, 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, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 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", "Dry", "Dry", "Dry", "Wet", "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", "Wet", "Wet", 
"Wet", "Dry", "Dry", "Dry", "Dry", "Dry", "Dry"), CarcassPres = structure(c(1L, 
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 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, 
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, 2L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L), .Label = c("0", 
"1"), class = "factor"), StandardTransect = structure(c(4L, 3L, 
1L, 5L, 7L, 5L, 1L, 8L, 6L, 4L, 3L, 8L, 5L, 6L, 1L, 4L, 3L, 4L, 
3L, 5L, 1L, 8L, 6L, 8L, 6L, 5L, 1L, 8L, 6L, 1L, 5L, 5L, 1L, 6L, 
4L, 3L, 8L, 6L, 5L, 1L, 8L, 4L, 3L, 5L, 1L, 8L, 6L, 5L, 1L, 6L, 
8L, 5L, 1L, 8L, 6L, 4L, 3L, 5L, 1L, 8L, 6L, 3L, 4L, 10L, 2L, 
1L, 5L, 6L, 2L, 10L, 5L, 8L, 8L, 1L, 6L, 3L, 4L, 4L, 3L, 9L, 
2L, 10L, 5L, 1L, 7L, 5L, 7L, 2L, 10L, 9L, 4L, 3L, 10L, 2L, 5L, 
1L, 7L, 4L, 3L, 2L, 10L, 9L, 5L, 7L, 1L, 2L, 9L, 4L), .Label = c("Jongomero", 
"Kidai", "LakeChada", "LakeKatavi", "Lunda", "Magangwe", "Mbagi-Mdonya", 
"Mpululu", "Msolwa", "Mtemere"), class = "factor"), Tlength = c(35.2, 
86.7, 93, 75, 27.2, 74.4, 93, 10.3, 45.8, 35.2, 78.2, 10.3, 71, 
45.8, 93, 35.2, 63.9, 35.2, 77.9, 86.6, 93, 10.3, 45.8, 10.3, 
45.8, 68.9, 93, 10.3, 45.8, 93, 86.7, 90.5, 93, 45.8, 35.2, 81.6, 
10.3, 45.8, 88.2, 93, 10.3, 35.2, 64.6, 82.3, 93, 10.3, 45.8, 
77.9, 93, 45.8, 10.3, 90.3, 93, 10.3, 45.8, 35.2, 77.4, 87.5, 
93, 10.3, 45.8, 66, 35.2, 71.2, 85.7, 93, 87.5, 45.8, 85.5, 69.6, 
97.8, 10.3, 10.3, 93, 45.8, 86.6, 35.2, 35.2, 71.9, 77.9, 88.5, 
80, 85.2, 93, 56.1, 85.5, 56.1, 97.6, 81.8, 79.7, 35.2, 71.1, 
81.9, 53.8, 86.5, 68.7, 60.7, 78.7, 56.6, 66.9, 71.8, 79.2, 82.6, 
71.8, 92.4, 85.6, 78.5, 77.3)), row.names = c(NA, -108L), class = "data.frame")

*I modelled the counts as a function of time (in months) interacting with park, the presence of carrion, the season (wet or dry), and a random effect of the transect nested within the park. We also included the log length of the transect as an offset to account for unequal sampling effort.

2

Your time variable is month, ndate, and the data is grouped by transect. However, the observations are not regularly spaced and there are a few instances of two observations for the same (transect, month) pair. Since you don’t have exactly one observation for each time (month) and group (transect), an autoregressive term is not really applicable. Why do you want to add one?

Here is a plot to illustrate the irregular sampling in time: x-axis = time (in months), y-axis = year, panel = transect, shape indicates the season (triangle = wet, circle = dry).

mc-stan-33362
mc-stan-333624500×1500 144 KB

There seem to be some inconsistencies in the season labels. For example, here is Mpululu in 2018:

  AWB ndate nyear    NP Season CarcassPres StandardTransect Tlength
    2    61  2018 Ruaha    Wet           0          Mpululu    10.3
    0    61  2018 Ruaha    Dry           0          Mpululu    10.3

How does the same ndate fall in bot the wet and dry season?

1 Like
3

Thanks for spotting that mistake! I’ll deal with that first. Nice figure to help see that.