Same model with different default priors in rstanarm and brms?

Interfaces rstanarm
1

In a simple glm with one predictor it seems as if rstanarm is using different default priors compared to brms.

Now for the rstanarm package the default priors (and adustment) are clear to me, but does anyone know what brms does in order to set priors scale and location?

The model code for rstanarm:

p_1_ec_g_brm_rstan <- stan_glm(ec ~ g, data = p_1_36ec, chains = 4, iter = 10000, warmup = 1000)

The model code for the brms:

p_1_ec_g_brm_n <- brm(ec ~ g, data = p_1_36ec, chains = 4, iter = 10000, warmup = 1000)

After modelling the following priors are used:

Priors for model 'p_1_ec_g_brm_rstan' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 53, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 53, scale = 16)

Coefficients
  Specified prior:
    ~ normal(location = 0, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 0, scale = 40)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.16)
------

While the brms shows the following in the summary:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: ec ~ g 
   Data: p_1_36ec (Number of observations: 1057) 
Samples: 4 chains, each with iter = 10000; warmup = 1000; thin = 1;
         total post-warmup samples = 36000

Priors: 
Intercept ~ student_t(3, 55, 7.4)
sigma ~ student_t(3, 0, 7.4)

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    51.10      0.43    50.25    51.96 1.00    39083    27776
gv            2.72      0.48     1.77     3.67 1.00    39362    27428

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     6.18      0.14     5.92     6.46 1.00    40910    27733
  • Operating System: mac os 11.3
  • rstanarm Version: 2.21.2
  • brms Version 2.15.0
1 Like
2

I know you can see the brms priors with get_prior get_prior: Overview on Priors for 'brms' Models in brms: Bayesian Regression Models using 'Stan' . As for how it picks those, Iā€™d need to dig around in the manual and get back to you.

1 Like
3

Thanks for the tips. The code with the priors for the brms was obtained with summary(ā€¦ priors = TRUE) which the gave, obviously, the same values as with the get_prior function; the priors for the rstanarm-model was obtained with the prior_summary function. The problem is, I donā€™t know how the brms package ā€˜calculatesā€™ the priorā€™s distribution, scale and location.

I couldnā€™t find anything in the manual or the website info. But possibly I overlooked somethingā€¦

1 Like
4

I donā€™t think brms's approach is completely documented anywhere beyond the code of the package and has in fact evolved over time. Some relevant discussion and pointers at: How does brms use the given data to adjust the default priors?

I think the philosohpy of brms is that the moment you start to be curious about priors, you should be defining your own anyway. The defaults are not some kind of ā€œbest practceā€, they are just meant to avoid the model breaking badly for most use cases - nothing more, nothing less.

Best of luck with your model!

2 Likes
5

For the intercept, the location and scale is given by (as also mentioned in the post Martin linked):

location = round(median(data$outcome), 1)
scale = max(round(mad(data$outcome), 1), 2.5)

The prior for sigma is calculated using the same approach for the scale, but the location is defaulted to 0.

You can can verify this using some of the example data:

> get_prior(count ~ zAge, data = epilepsy, family = gaussian())
                prior     class coef group resp dpar nlpar bound       source
               (flat)         b                                       default
               (flat)         b zAge                             (vectorized)
 student_t(3, 4, 4.4) Intercept                                       default
 student_t(3, 0, 4.4)     sigma                                       default
> round(median(epilepsy$count), 1)
[1] 4
> max(round(mad(epilepsy$count), 1), 2.5)
[1] 4.4
3 Likes