How to construct group-level and family summaries in brms?

Interfaces brms
1

I’m trying to understand how are group-level effects constructed from brms object in summary.brmsfit. For example, with the following output from summary:

summary(fm)
    Formula: Y ~ 1 + (1 | subject)
       Data: dat (Number of observations: 2604) 
    Samples: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
             total post-warmup samples = 2000

Group-Level Effects:
~subject (Number of levels: 124) 
                Estimate  Est.Error   l-95% CI   u-95% CI Eff.Sample       Rhat
sd(Intercept) 0.07724968 0.00629353 0.06620370 0.09049358        540 1.00791419

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample    Rhat
sigma  0.15308   0.00212  0.14910  0.15729       3493 0.99895

I couldn’t figure out how to obtain those numbers above for ‘subject’ such as Estimate, Est.Error, l-95% CI and u-95% CI for the group-level effects of ‘subject’.

With

bb ← ranef(fm, summary = FALSE)
cc ← bb[[‘subject’]][,‘Intercept’]
mean(apply(cc, 1, sd))

the last line gave me 0.07627706, which is slightly different from the summary results above. I don’t know how to get Est.Error and the quantile intervals either. For example,

sd(apply(cc, 1, sd))
[1] 0.003229776
quantile(apply(cc, 1, sd), c(0.025, 0.975))
      2.5%      97.5% 
0.06999685 0.08267143

Similarly, how to obtain the estimate, error and quantile interval for the residuals (sigma)? Again, the following

quantile(apply(residuals(fm, summary=FALSE), 1, sd), c(0.5, 0.025, 0.975))
      50%      2.5%     97.5% 
0.1530026 0.1520525 0.1541562

gives me different results from summary(fm).

2

The values may be slightly different since the value in the summary refers to the actual SD parameter whereas the SD based on the ranefs is kind of the empirical SD of the coefficients. In any case, they will usually be very similar as long as you have enough group levels.

To extract the values from the summary, you simply use posterior_summary. You can check the exact names of the parameters via parnames.

3

Thanks, Paul!

You are right, the “Estimate” from summary() is very close to the empirical results from ranef(). However, I’m still puzzled by the dramatically different results of the other 3 numbers from summary(), “Est.Error”, “l-95% CI”, and “u-95% CI”, which are are very different from what I derived from ranef(). Could you shed some light on these large discrepancies?

4

The same applies. You are computing a similar concept using two different ways. One is directly via a parameter and the other is empirically. Different values may occur here and given the estimation uncertainty are probably not too different at all.

5

Thanks again, Paul!

Do you mind showing how those variances and their uncertainties are derived through their corresponding parameters? If I want to obtain the posterior distribution of each variance, is the empirical data through ranef() the only way?

6

Please reread my Initial comment. There is an SD parameter for this purpose whose posterior samples you can transform and summarise to your needs. See eg posterior_samples, posterior_summary and parnames.

1 Like