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