Model selection using ELPD LOOIC

Hi!

I want to assess the effects of first- and second-order temporal trends (t) and survey timing (doy) on my binomial response variable (occ = territory occupancy). The goal is to select the “best” baseline model to use in subsequent analyses. I am using ELPD from LOOIC to compare nested models and would like to verify that I am not selecting the “wrong” model.

The most complex model is:

to_base8 <- brm(occ ~ z_t + I(z_t^2) + z_doy + I(z_doy^2) + (1|year) + (1|terrID), 
                data = to_df,
                family = bernoulli(link="logit"),
                iter = 6000) 

The rest of the models are simplified versions of the most complex model, with the simplest being:

to_base0 <- brm(occ ~ 1 + (1|year) + (1|terrID), 
                data = to_df,
                family = bernoulli(link="logit"),
                iter = 6000) 

All models include year and terrID as random effects.

Below is my output using loo_compare()

         elpd_diff se_diff elpd_loo se_elpd_loo p_loo   se_p_loo looic   se_looic
to_base6     0.0       0.0 -1504.5     29.8       105.9     3.3   3009.0    59.6 
to_base8    -1.2       1.0 -1505.7     29.9       106.9     3.3   3011.4    59.8 
to_base3    -2.7       3.1 -1507.2     29.9       113.2     3.4   3014.4    59.8 
to_base7    -4.2       3.2 -1508.7     29.9       113.7     3.5   3017.4    59.8 
to_base2    -4.4       3.8 -1508.9     29.7       117.2     3.5   3017.8    59.4 
to_base5    -5.7       3.9 -1510.2     29.7       118.3     3.5   3020.3    59.5 
to_base4   -30.8       7.7 -1535.3     29.8        97.7     3.1   3070.5    59.6 
to_base1   -40.4       8.3 -1544.9     29.8       105.7     3.3   3089.8    59.7 
to_base0   -46.1       8.0 -1550.6     29.6       116.8     3.5   3101.2    59.2 

My understanding is that there is a negligible difference in the predictive performance between base6 and base8/base3 because Δelpd < 4. Similarily, base7, base2, and base5 are not significantly different from base6 because their Δelpd < SE*2.

In these scenarios, Liu et al., 2025 suggests selecting the simplest model among those not significantly worse than the top model to avoid selection bias.

In my case, this would be base2 which is

to_base2 <- brm(occ ~ z_doy + (1|year) + (1|terrID),
data = to_df,
family = bernoulli(link="logit"),
iter = 6000)

compared to the top model which is

to_base6 <- brm(occ ~ z_t + I(z_t^2) + z_doy + (1|year) + (1|terrID), 
                    data = to_df,
                    family = bernoulli(link="logit"),
                    iter = 6000) 

Is this correct?

Thank you in advance for any insight!

Katie

It has been common advice to choose the simplest model among the best with similar performance, and your interpretation of the values is correct. However, in [2606.22850] To select or not to select: predictively consistent priors instead of model selection, we argue it is better to choose the biggest model unless it has clearly worse predictive performance.

Thank you so much for providing a link to the preprint! To confirm, I would select the biggest (i.e., most complex) model among those with Δelpd ≤ 4 (i.e., base6, base8, base3), or would this also include those with Δelpd < SE*2 (i.e., base6, base8, base3, base7, base2, base5)? All my models use flat priors. Either way, base8 (the second “best” model) is the biggest in this case.

Thanks again!

As the number of models is small, it’s enough to compare the biggest model to the best model. I recommend using R2D2 and R2D2M2 priors (supported by brms), but if you have big data, then there is not much difference to using flat priors on coefficients

My datasets span 33 years with 1173-3118 observations so I imagine the priors you recommended won’t change the results much. Based on that, I’ll go with the biggest model rather than the best model. Thank you so much!!