(brms version is 2.19.0, cmdstan version 2.32.0, cmdstanr version 0.5.3.)
For a logit coefficient of interest, whose posterior has been sampled by brms, I need to test the null hypothesis that the parameter equals 0. I’ve been told that this can be accomplished by “evaluating the posterior probability of the region of values having smaller density than at \theta_j = 0.” (Agresti 2015: 339). To this end I’ve written a convenience function that does the following (key lines only):
post <- as.data.frame(as_draws_df(brmsmodel))
...
densities <- sapply(post[,CoefName], function(x) density_ratio(post[,CoefName], point = x))
cutoff <- density_ratio(post[,CoefName], point = 0)
p.val <- mean(densities < cutoff)
I presume that this is better than simply dividing the coefficient by its SE and calculating 2 \times the p-value of values as or more extreme as the posterior mean, assuming a normal distribution with \mu = 0 and \sigma = \text{SE}, particularly because not all the relevant parameter posteriors look symmetrically bell-shaped.
But here’s the question: I read in the helpfile of density_ratio() that “you may need an effective sample size of 10,000 or more to reliably estimate the densities.” How do I find this effective sample size? The brms summary output only lists bulk_ess and tail_ess, and I don’t know how to infer the crucial quantity from those numbers.
Reference
Agresti, A., 2015. Foundations of linear and generalized linear models, Wiley series in probability and statistics. John Wiley & Sons Inc, Hoboken, New Jersey.
UPDATE: Since I’m interested in posterior tail probabilities here, could I use use ess_tail instead of the unreported “effective sample size”? I.e. could I just increase sampling iterations until summary(brmsmodel) reports ess_tail at \geq10000, and trust that all manner of things shall be well?
Second question: should I be concerned that posterior::ess_bulk() and posterior::ess_tail() return values which Differ Slightly from those reported by summary(brmsmodel) ?
hypothesis() is indeed a useful function. It doesn’t report two-tailed p-values for null point hypotheses though — only whether they’re rejected. That’s dichotomous and thus not quite as useful as a p-value (Andrew Gelman for one has advocated that p-values be interpreted continuously) — but still useful.
Anyway, I’m very curious about something. What exactly is the Post.Prob that hypothesis() produces for point hypotheses of the form Parameter = 0? It’s not a density ratio (that would be the evidence ratio), and it’s not a probability density (those can exceed one). What is it? I suppose it must be the area under the curve of some interval around 0, but how wide an interval exactly, and what determines the width? Is there some conventional value for it? Inquiring minds want to learn more statistics.
For an two-sided (point) hypothesis, the evidence ratio is a Bayes factor between the hypothesis and its alternative computed via the Savage-Dickey density ratio method. That is the posterior density at the point of interest divided by the prior density at that point. Values greater than one indicate that evidence in favor of the point hypothesis has increased after seeing the data. In order to calculate this Bayes factor, all parameters related to the hypothesis must have proper priors and argument sample_prior of function brm must be set to "yes". Otherwise Evid.Ratio (and Post.Prob) will be NA. Please note that, for technical reasons, we cannot sample from priors of certain parameters classes. Most notably, these include overall intercept parameters (prior class "Intercept") as well as group-level coefficients. When interpreting Bayes factors, make sure that your priors are reasonable and carefully chosen, as the result will depend heavily on the priors. In particular, avoid using default priors.