Calculating individual specific intercepts in brms


#1

I have a question about the calculation of the individual specific intercepts in a multilevel model using brms. I make a multilevel model using a panel dataset using the following command:
result_brms = brm(Activations ~ avg_CATALOGUE_PRICE + (1 | ID), data = aggregated_data, chains = 4, iter = 10000, warmup = 5000, save_model = "test_code_stan.txt", cores = getOption("mc.cores", 4L)),

so in my model I only have a individual-specific intercept. To get these intercepts I use the command

ind_specific_coefficient = coef(result_brms) and
ind_specific_coefficient = as.data.frame(ind_specific_coefficient$ID).

I also tried do construct these individual-specific-intercepts, \alpha_i, by myself using the stan file in the stan function generated by the brm function. To get the same estimates as with brm I just take the mean of the individual specific draws r_1_1[i] and add the mean of the b_intercept. However, if the individual intercepts are determined by:

\alpha = \mu + u, where u \sim N(0,\sigma)

I would expect that I can construct the \alpha_i's by drawing from a normal distribution with mean equal to mean(temp_intercept) (which equals the \mu I think) and standard deviation sd(r_1_1[i]) (which equals \sigma). However, in this case I get different \alpha_i's than from the brm function, but I do not get why this is wrong. Could someone explain this to me? Also, when I make a histogram of the intercepts from the brm function they look heavily censored while this is not the case in my approach. Any help is greatly appreciated.


#2

Your intuition about sd(r_1_1[i]) is incorrect. The SD you are computing this way does not equal sigma. Instead, you are estimating the SD of the posterior distribution of each of the r_ coefficients, which tells you something about the precision with with each of these coefficients was estimated, not the varation between coefficients, which is what \sigma represents.


#3

So the mean of the r_1_1[i] draws are a measure for \sigma? Or is it just coincidence that in this case the mean of the r_1_1[i] draws added on top of the b_intercept give approximately the same individual specific intercepts as using coef(result_brms)?


#4

You have to compute stuff on the levels of posterior samples. Compute

ran <- ranef(<model>, summary = FALSE)[[<group name>]]
str(ran)
sigma <- apply(ran, 1, sd)
str(sigma)

to get the posterior samples of sigma.

Similarily, if you add ran to the vector of posterior samples of b_intercept, you get the same posterior samples as returnd by coef(<model>, summary = FALSE).


#5

Ah thank you so much for your explanation!