Interpretability of Population-Level Effects with smooth terms in brms

Modeling
,
1

Hello all,

I am working with seabird count data, producing models in brms using the mgcv smooths available:

  m1 <- brms::brm(CountIndivs ~ s(Date, k = 6) +
                    s(DayofYear, bs = 'cc', k = 6) +
                    s(daily_wind_speed, k = 4) +
                    s(wind_30dav, k = 4) +
                    s(upwelling_30dav, k = 4) +
                    s(SOI, k = 4) +
                    s(SAM, k = 4) +
                    (1|Port),
                  data = speciesdf,
                  family = zero_inflated_negbinomial(link = "log"),
                  cores = 4,
                  seed = 17,
                  iter = 4000,
                  warmup = 1000,
                  thin = 10,
                  refresh = 1000,
                  control = list(adapt_delta = 0.99))

I am running into issues when trying to interpret the output of the summary function:

summary(m1)

 Family: zero_inflated_negbinomial 
  Links: mu = log; shape = identity; zi = identity 
   Data: speciesdat (Number of observations: 397) 
Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 10;
         total post-warmup samples = 1200

Smooth Terms: 
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(sDate_1)                 1.83      1.78     0.06     6.60 1.00     1164     1171
sds(sDayofYear_1)            2.37      1.19     0.84     5.14 1.00     1192     1215
sds(sdaily_wind_speed_1)     2.55      2.43     0.09     8.96 1.00     1089     1162
sds(swind_30dav_1)           2.84      2.35     0.11     8.74 1.00     1239     1135
sds(supwelling_30dav_1)      1.87      1.69     0.06     6.56 1.00     1206     1202
sds(sSOI_1)                  4.27      4.55     0.13    16.13 1.00     1098      884
sds(sSAM_1)                  2.81      2.24     0.15     8.41 1.00     1185     1101

Group-Level Effects: 
~Port (Number of levels: 4) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.66      0.65     0.02     2.33 1.00     1233     1210

Population-Level Effects: 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              -2.40      0.71    -3.91    -1.12 1.00      832     1163
sDate_1                -4.02      2.98   -10.51    -1.61 1.00     1195     1178
sdaily_wind_speed_1    -0.52      5.93   -14.26     9.32 1.00      978     1030
swind_30dav_1          -0.97      4.79    -8.98    10.03 1.00     1211     1160
supwelling_30dav_1     -1.74      4.84   -11.48     7.29 1.00     1308     1067
sSOI_1                  2.39      8.98   -11.28    27.01 1.00     1141     1156
sSAM_1                  4.54      4.01    -3.85    12.20 1.00     1116     1174

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.36      0.20     0.15     0.87 1.00     1166     1032
zi        0.29      0.17     0.02     0.58 1.00     1060     1238

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

As I understand it, the ‘Smooth Terms’ (e.g. sds(sDate_1) represent the wiggliness of each term, and where the CI cross zero, there is little usefulness in including a smooth term on the given variable.

My question concerns how I interpret the ‘Population-Level Effects’ section of the summary output. Would terms with Population-Level Effects CI’s crossing zero be of little significance in explaining variability in the data used (i.e. the term is not a good predictor of the response variable)? In this example, ‘sDate_1’ is the only term where the CI’s do not cross zero.

I ask as I am running models for each species in my dataset (67 species) with each one containing the same predictors. This is producing a significant number of results and smooth plots, and I am trying to determine which terms should be included for each species in the publication, and which are of little significance and resigned to supplementary.

I realise this is quite a basic question but I have yet to be able to find a definitive answer.
Thank you very much.

2 Likes
2

@nick1 Maybe doing a quick summary plot of all the effects might help to better understand the results and help clarify any relevant interpretations:

require(patchwork)

p1 = plot(conditional_effects(m1), plot=FALSE)

wrap_plots(p1) + 
  plot_annotation(title="seabird species: albatross", 
                  theme=theme(plot.title=element_text(hjust=0.5)))
1 Like
3

No, I don’t think this is correct. See here and the link to Tristan Mahr’s explanation in my first response to my own question. In Tristan’s article, there is some explanation of what the ‘population-level effects’ in brms are as relates to these smooth terms, but I think they are the coefficient for one of the ‘natural’ parameterized basis functions. Thus, they really have no direct interpretation.

Are you trying to determine which terms should be included at all? or which terms don’t really need a smooth? You could use LOO to compare the models. I think if you are trying to decide which out of a bunch of terms to include in the models based on predictive performance, then maybe projection predictive variable selection might be the way to go, but I am not really familiar with this.

As @MilaniC said, making plots from predictions is the way to go to understand the results from models with smooth terms (at least that is what I do). If you want to understand the slopes of the smooths, you can take the first derivative and plot this with uncertainty intervals. See here for the concept, and for brms I think you would need to use posterior_smooths() for reasons described here.

4

Thank you for the suggestions. All of the models I am running contain the same suite of variables, however the variable show differing predictive performance between the models. I am not currently looking to perform variable selection across each model as I do not have access to the computing power required.

My biggest issue is currently in determining whether a predictor is having an appreciable effect on my response. I can eyeball this by examining my pots, however I am looking for a more concrete way of determining importance.

As an example the below plot shows very little meaningful information:

Rplot
Rplot581×595 6.62 KB

However a much more meaningful trend can be observed in this plot:

Rplot01
Rplot01581×595 8.15 KB

Is there any way of drawing a line between these two in brms? If I was using mgcv I would rely on the p-values reported for each term, however I am unsure of an equivalent way in brms.

Thanks again.

5

@nick1 Ok not really helpful to give plots from mgcv rather than your brms model results.

Anyway, the following might be helpful to you to better understand where you might be going here:

6

Thanks MilaniC, those plots are from brms models, the y-axis titles are slightly misleading as the brms terms have k values set using the edf from similar mgcv GAMs, but those plots are definitely brms.

Thanks for linking the Gavin Simpson article, I have read through it a few times previously, but I cannot find a solution there.

Cheers

7

I’ve reached out to Gavin myself when having questions. Don’t be shy to contact him and ask him. He’s a friendly person :)

2 Likes
8

I suspect the sds for the smooth for upwelling in the first plot is very close to zero. The sds can tell you about the ‘wiggliness’ and whether a smooth might be needed. You can remove the smooth for that variable and compare the model via LOO. You can also look at p_loo and the se, which should be similar to the EDF in gam. Doing these things will help you determine if the smooth is needed at all.

In terms of effect on the response, estimating the first derivative and plotting so that you can see the value of the slope over a range of the predictor values can give information on the relationship between predictor and response.

Let’s take your example plots above and say that you implemented the suggestions above. I suspect you would find this: For the first plot you would likely see the sds close to zero and LOO would likely show better performance for the model without the smooth. When plotting the estimate of the first derivative with 95% UI’s, you would likely see that the UI’s broadly cover both sides of zero across the entire range of plausible predictor values, suggesting that the slope of the smooth at any given value of the predictor is inconclusive. For the predictor in the second plot, you would probably see that the sds is well away from zero and LOO would likely show better performance for the model with the smooth. Plotting the estimated first derivative would likely show slopes with UI’s well away from zero on a range of days of year, and those slopes that were around zero would have much narrower UI’s than the predictor shown in the first plot.

I don’t think there is a way to simply look at the ‘Population-level effects’ like it seems like you want. But if you reach out to Gavin Simpson and find out otherwise, I would like to hear about it:-)

1 Like