paul-buerkner / brms

brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan
https://paul-buerkner.github.io/brms/
GNU General Public License v2.0
1.28k stars 184 forks source link

Allow group-level estimates in hypothesis() #327

Closed mvuorre closed 6 years ago

mvuorre commented 6 years ago

Among other things, hypothesis() can be used to easily calculate posteriors of parameter transformations. For example, if I'm interested in the posterior distribution of b0/b1, I can easily get a summary of it with

hypothesis(fit, "b0/b1 = 0")

Background

This is a really powerful function for what is sometimes known as "general linear hypothesis testing". However, sometimes (apart from hypothesis testing) we would like to also obtain the transformed parameters for varying effects in multilevel models. So, for example

fit <- brm(y ~ 1 + x + (1 + x | id), data = d)

If we use

hypothesis(fit, "Intercept / x = 0")

The hypothesis is by default evaluated for the population-level Intercept and x (b_Intercept, b_x). We can evaluate it for the standard deviations with:

hypothesis(fit, "Intercept / x = 0", class = "sd", group = "id")

Question

Would it be possible to allow calculating this quantity for the group-level parameters as well?

I'm imagining a syntax something like this:

hypothesis(fit, "Intercept / x = 0", class = "r", group = "id")

I can already calculate this for one group at a time by using the default method (for example for id = 1):

hypothesis(as.data.frame(fit), "r_id.1.Intercept. / r_id.1.x. = 0")

I could also create strings representing that hypothesis for all ids, and enter them to hypothesis(), but I think the proposed syntax would make the process a lot neater.

This would be really useful, for example, when one wants to visualize multiple transformed estimates at the group- and population-levels simultaneously. Thoughts?

Thanks!

paul-buerkner commented 6 years ago

This sounds like a good idea. Let's see if I can get this implemented. As a side note, you can also use the brmsfit method of hypothesis for group-level effects. For instance, the following code works:

fit1 <- brm(count ~ log_Age_c + log_Base4_c * Trt + (1|patient),
            data = epilepsy, family = poisson())
hypothesis(fit1, "r_patient[1,Intercept] = 0", class = "")
paul-buerkner commented 6 years ago

It should be working now. Try out argument scope in the new github version.

mvuorre commented 6 years ago

Wonderful, this works very well. Thanks!