Variable selection with bayesian linear mixed models (the brms package)

Question

I am fitting a bayesian linear mixed model in R with 6 variables and 2 random effects. Inclusion of all 6 variables is motivated by a well-founded hypothesis. Does it make sense to do variable selection or can I report whether or not each fixed effect has an effect on the dependent variable (CI does not overlap 0)? And then, perhaps, compare the final model to an intercept only model using LOOIC to evaluate model fit? Or run the full model, select only the significant variables, and compare?

If I should do stepwise variable selection is there a way to automate it as with dredge() in MuMIn?

I am trying to figure out the most correct way to test if each covariate has an effect on y without running ~100 models (which I would have to do to examine model fit for all combinations of covariates).

eg,

library(brms)
library(loo)
model.full <- brm(data = data,
             family = gaussian,
             formula = y1 ~ x1 + x2 + x3 + x4 + x5 + x6 +
               (1|r1) + (1|r2),
             prior = c(prior(normal(0, 10), class = Intercept),
                       prior(cauchy(0, 1), class = sd),
                       prior(cauchy(0, 1), class = sigma)),
             iter = 2000, warmup = 1000, chains = 4, cores = 2,
             control = list(adapt_delta = .975, max_treedepth = 20),
             seed = 12261996)

model.intercept <- brm(data = data,
             family = gaussian,
             formula = y1 ~ 1 +
               (1|r1) + (1|r2),
             prior = c(prior(normal(0, 10), class = Intercept),
                       prior(cauchy(0, 1), class = sd),
                       prior(cauchy(0, 1), class = sigma)),
             iter = 2000, warmup = 1000, chains = 4, cores = 2,
             control = list(adapt_delta = .975, max_treedepth = 20),
             seed = 12261996)

loo(model.full, model.intercept)

Answer 1

Given you have well-founded hypotheses about which variables should be included as predictors of your outcome, there is little reason to worry about variable selection.

As you point out, the model comparison you are running with loo is evaluating the fit of a model with only intercepts and no predictors to a model with all the predictors plus the intercepts. Unless most of your predictors are not at all helping to explain the outcome (or predict to new observations), I would be shocked if the results indicated the random intercept only model was the superior model. Thus I am not sure it makes sense to run the model comparison using these two models.

Is it the case that one or more of the predictors are less well-founded in the literature? If so, you might instead remove those two predictors, run that model, and compare it to the model with all predictors included. Or, if there is some reason to believe that the effect of one or more of the predictors on the outcome vary across groups (i.e., a random slope), you might estimate such a model and compare it to the model without random slopes.

Given Bayesian estimation and loo takes some time to run, depending on sample size and model complexity, I would be more judicious about model comparison. Or, I would consider doing some of these comparisons using maximum likelihood in lmer or perhaps rstanarm because the models are pre-compiled. Perhaps also consider using variational inference with either brms or rstanarm with the (algorithm="meanfield") option just to speed things up for model testing.