Why don't lines from my conditional effects plots pass through the data points?

General
, ,
1

I’m fitting a model of population counts over time in three different locations with some random effects structure. The model fits fine using brms but when I plot the conditional effects the lines of the resultant graph don’t pass through any of the data points. The posterior predictive plot looks OK so I’m wondering if this is a consequence of the choice of priors?

conditional HV

Here’s the model:


# specify the priors
newprior_hv <-
  c(
    set_prior("normal(0, .2)", class = "b", coef = "Intercept"),
    set_prior("normal(0, .2)", class = "b", coef = "NPRuaha"),
    set_prior("normal(0, .2)", class = "b", coef = "NPSelous")
  )


#' fit the model 
nbglm_norm_prior_tight_hv <-
  brm(
    HV ~  0 + Intercept + # this allows control of prior on the intercept
      ndate * NP + Season + CarcassPres +
      (1 | NP / StandardTransect) +
      offset(log(Tlength)),
    family = negbinomial,
    data = mydata,
    chains = 4,
    control = list(adapt_delta = 0.999, max_treedepth = 15),
    prior = newprior_hv,
    save_pars = save_pars(all = TRUE)
  )


#' check the models
plot(nbglm_norm_prior_tight_hv)
plot(conditional_effects(nbglm_norm_prior_tight_hv, prob = 0), ask = FALSE, points = T)

Here’s the data:

mydata <- structure(list(HV = c(0, 1, 21, 2, 1, 5, 8, 3, 1, 7, 2, 0, 1, 
0, 2, 0, 0, 0, 0, 12, 2, 0, 4, 0, 3, 5, 2, 0, 0, 3, 7, 0, 6, 
8, 8, 0, 0, 2, 2, 3, 2, 3, 3, 0, 0, 0, 0, 4, 1, 1, 0, 8, 1, 0, 
3, 11, 0, 2, 1, 0, 5, 0, 4, 0, 2, 0, 5, 4, 2, 9, 5, 0, 0, 16, 
7, 2, 1, 0, 0, 2, 0, 1, 2, 5, 1, 7, 0, 1, 0, 2, 0, 6, 16, 0, 
8, 4, 1, 0, 7, 6, 6, 0, 1, 1, 3, 2, 4, 2), 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")
2

I think this is because of the offset in your model. The documentation on this is a few layers down because conditional_effects is passing arguments in ... down to posterior_epred, which in turn passes them down to the prepare_predictions helper function. Here you’ll see that it defaults to offset = TRUE. Toggle that off in your last line of code, and you should be all set:

plot(conditional_effects(nbglm_norm_prior_tight_hv, prob = 0, offset = FALSE), ask = FALSE, points = T)

Rplot

2 Likes
3

Thanks for this. Strangely for another data set on a separate species the plot fits lines to the points when the offset is true and is completely off when it’s false i.e. the opposite to the initial one.

Here’s the case for true
AWBV TRUE.jpeg

and here’s the case for false
AWBV FALSE.jpeg

And for completeness here’s the code and data:

newprior_awb <-
  c(
    set_prior("normal(0,.7)", class = "b", coef = "Intercept"),
    set_prior("normal(0,.7)", class = "b", coef = "NPRuaha"),
    set_prior("normal(0,.7)", class = "b", coef = "NPSelous")
  )

#' fit the model with tighter normal priors
nbglm_norm_prior_tight <-
  brm(
    AWB ~  0 + Intercept + # this allows control of prior on the intercept
      ndate * NP + Season + CarcassPres +
      (1 | NP / StandardTransect) +
      offset(log(Tlength)),
    family = negbinomial,
    data = mydata,
    chains = 4,
    #control = list(adapt_delta = 0.999, max_treedepth = 15),
    prior = newprior_awb,
    save_pars = save_pars(all = TRUE)
  )

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")