Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Plot method for default priors #205

Draft
wants to merge 13 commits into
base: develop
Choose a base branch
from

Conversation

venpopov
Copy link
Owner

@venpopov venpopov commented Apr 6, 2024

Here's my current prototype. It depends on ggdist and distributional. ggdist has a nice parse_prior function to work with brms priors, but I noticed it gave errors for a bunch of cases, so I wrote a function prep_brmsprior which does some preprocessing of a brmsprior object to ensure no errors happen and to add the family/model links to the prior object.

I also added a change to the bmm default_prior method to add a column to the prior with the link form.

Then we have a few functions for plotting distributions, but the main one is plot.brmsprior, which is a method for plot. So now you can do this:

library(bmm)
#> Loading required package: brms
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.21.1). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
# get the default prior for a bmm model given a formula
priors <- default_prior(bmf(c ~ 0 + set_size, 
                            kappa ~ set_size + (set_size | ID),
                            mu ~ 1),
                        oberauer_lin_2017,
                        sdm('dev_rad'))

# see the prior object
priors
#>                     prior     class      coef group resp  dpar nlpar   lb   ub
#>                    lkj(1)       cor                                           
#>                    lkj(1)       cor              ID                           
#>     student_t(5, 2, 0.75)         b set_size1                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size2                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size3                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size4                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size5                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size6                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size7                c       <NA> <NA>
#>     student_t(5, 2, 0.75)         b set_size8                c       <NA> <NA>
#>              normal(0, 1)         b set_size2            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size3            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size4            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size5            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size6            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size7            kappa       <NA> <NA>
#>              normal(0, 1)         b set_size8            kappa       <NA> <NA>
#>      student_t(3, 0, 2.5)        sd                      kappa          0     
#>      student_t(3, 0, 2.5)        sd              ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd Intercept    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size2    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size3    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size4    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size5    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size6    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size7    ID      kappa          0     
#>      student_t(3, 0, 2.5)        sd set_size8    ID      kappa          0     
#>        student_t(1, 0, 1) Intercept                         mu       <NA> <NA>
#>     student_t(5, 2, 0.75)         b                          c       <NA> <NA>
#>              normal(0, 1)         b                      kappa       <NA> <NA>
#>  student_t(5, 1.75, 0.75) Intercept                      kappa       <NA> <NA>
#>        source     link
#>       default identity
#>  (vectorized) identity
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>       default      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>  (vectorized)      log
#>          user tan_half
#>          user      log
#>          user      log
#>          user      log

# plot the priors on the sampling scale
plot(priors)

# plot the priors on the native scale by applying the inverse link
plot(priors, transform = TRUE)

# same, but manual control of max limits
plot(priors, transform = TRUE, stat_slab_control = list(limits = c(-4, 30)))

Created on 2024-04-06 with reprex v2.1.0

A few problems:

  • for intercepts or ~0+factor coding the transformation is meaningful. But when we have an intercept + effects, the prior for the effects coefficients could be confusing on the native scale, because this is now (in the log-link case), a prior on the multiplicative coefficient. So I see two options - not to display priors for coefficients on the transformed scale, or to give a warning that these should be interpreted carefully.
  • a similar problem is for the sd of the random effects. For those I'm not even sure how to interpret them on the native scale. So maybe skip those completely if plotted on the transformed scale.
  • as you can see above, sometimes the scale of the plots can be quite wide, so it is sometimes necessary to manualyl control the limits
  • for now I've opted for printing the label of the facet to be a combination of the parameter class and parameter name, followed by the prior distribution. For transformed distributions, the default from distributional is to wrap the distribution in t(), but maybe it would be better to put the inverse link function there. But I'm also not sure how to communicate the parameter name that it is clear this is for the transformed parameter (e.g. currently for the transformed plot it says intercept_kappa, but this is exp(intercept_kappa)

Any thoughts or feedback, on the above things or anything else, are welcomed

@venpopov venpopov linked an issue Apr 6, 2024 that may be closed by this pull request
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add plot method for default_priors
1 participant