Thanks for your reply, @jsocolar ! The explanation and worked example were extremely helpful. One thing I neglected to mention in my original post was that I had already successfully fit the model I was considering (with some strengthening of the default priors, see below). However, the results from sampling the posterior (spoiler alert: all my coefficients are belong to zero [1]), together with my lack of experience, led me to doubt whether I was specifying the model appropriately.
In case anyone is interested, I’ve reproduced my workflow and thought process for fitting the model below. I’d be more than happy to hear any thoughts or critiques you may have regarding the choices I’ve made.
Data
birds.csv (67.6 KB)
I have a decent sample—905 individuals of 30 species. For a few species, sample size is >100; many are between 30 and 100. Several species were poorly sampled, with <10 samples/species.
library(brms)
library(bayesplot)
br <- read.csv("tmp/birds.csv") # import data and set factor references
br$age_sex <- factor(br$age_sex, levels=c("HY","AHY-F","AHY-M"))
br$stratum <- factor(br$stratum, levels=c("G","B","LC","HC"))
hist(table(br$species))
In a maximum likelihood framework I’d be inclined to drop these samples, considering it’s difficult to estimate species-level prevalence with so few trials. However, my understanding of the Bayesian approach is that it shouldn’t be detrimental to keep those species; they’ll help to shape the posterior for population effects, and I simply won’t try to make any strong inferences about prevalence in those particular species. Ultimately, due to missing predictor data, I end up losing a few of the poorly-sampled species anyway, dropping to 894 samples from 25 species.
Priors
# define model structure and priors
model_form <- bf(y ~ age_sex + stratum + group_size + nest + centrality + longevity + migratory + (1|species))
model_priors <- c(set_prior("normal(0,1.5)","b"), set_prior("normal(0,1.5)","Intercept"))
I defined the model structure as described in my earlier post. My first attempt to fit this with the default (flat) priors that brms places on the \beta coefficients led to divergent transitions and warnings about effective sample sizes. Following a suggestion made in another thread, I tried placing somewhat tighter priors (\text{Normal}(0,1.5)) on both the \beta coefficients and the intercept. This seemed to look better than the defaults when folded back into probability space, whereas the default prior on the intercept seems to place a lot of density near 0 and 1 (though perhaps I’m thinking about this incorrectly?).
hist(plogis(rnorm(1e3,0,1.5))) # quick prior check
hist(plogis(rstudent_t(1e3,3,0,2.5))) # default prior for the intercept and scale parameters
Sampling from the prior still showed a lot of density around 0 and 1, but there was at least some coverage of all reasonably-expected values, and I proceeded with sampling from the posterior.
# draw samples from prior
prior_samples <- brm(model_form, data=br, family=bernoulli(), chains=1, cores=1, prior=model_priors, seed=1138, sample_prior="only")
br_y <- get_y(prior_samples) # prior predictive checks
br_yrep_naive <- posterior_predict(prior_samples)
br_grps <- prior_samples$data$species
prevalence <- function(y) mean(y==1) # prevalence computation for ppc_stat()
pp_check(prior_samples)
ppc_stat(br_y, br_yrep_naive, stat="prevalence")
ppc_stat_grouped(br_y, br_yrep_naive, br_grps, stat="prevalence")
Posterior
Using the more informed priors allowed the chains to converge without issue. There were no warnings about divergences or effective sample sizes as when using the naïve priors. Additionally, it seems that the likelihood was able to dominate the posterior, which leads me to believe that the priors were not inappropriate (or if they were, then at least they were not too strong).
# sample posterior
post_samples <- update(prior_samples, chains=4, cores=4, sample_prior="no")
br_yrep <- posterior_predict(post_samples) # posterior predictive checks
pp_check(post_samples)
ppc_stat(br_y, br_yrep, stat="prevalence")
ppc_stat_grouped(br_y, br_yrep, br_grps, stat="prevalence")
Conclusion
summary(post_samples)
All of this was simply a prelude to finding that there is a lot of posterior density around 0 for all of the group-level effects. This wasn’t totally unexpected, particularly considering that prevalence was extremely low in the majority of species sampled (including several of the best-sampled species). I’m happy to accept that result and conclude a lack of evidence for any strong effects from the predictors—provided that I am confident that I’ve specified the model appropriately to address that question. But, I didn’t want to draw that conclusion prematurely if it turned that the reason for it was a model mis-specification. I think, now, that it wasn’t, and that the results can be accepted and interpreted accordingly.
Thanks again for your help, Jacob!