Hi,
are you interested in the fixed effects model or the multilevel model?
If it’s the multilevel model, the trick is to include group means of x
into your model. The tidyverse way, that would be:
library(dplyr)
d <- d %>%
mutate(x_mean = mean(x),
.by = firm_id)
And then fit the model:
library(brms)
m1 <- brm(y ~ x + x_mean + (1|firm_id),
data = d,
backend = "cmdstanr",
cores = 4,
seed = 123)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#Intercept 0.06 0.19 -0.32 0.44 1.00 4058 3137
#x 1.03 0.04 0.94 1.11 1.00 3346 2973
#x_mean 0.65 0.05 0.56 0.75 1.00 3096 2926
x
corresponds to the within-group-variation (a1
parameter in the blog post) and x_mean
is the between-group-variation (b
in the blogpost).
In econometrics, this is called Mundlak model. Snijders and Bosker explain it in detail in chapter 4 of their book Multilevel analysis: An introduction to Basic and Advanced Multilevel Modeling.
If you are asking about the fixed effects model, then for a single fixed effect, getting the point estimate is easy. You just need to mean center x
.
d$x_centered <- d$x - d$x_mean
m2 <- brm(y ~ x_centered,
data = d,
backend = "cmdstanr",
cores = 4,
seed = 123)
# Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#Intercept 5.58 0.41 4.78 6.39 1.00 4003 2919
#x_centered 1.03 0.11 0.82 1.24 1.00 3608 2950
Unfortunately, this model makes errors nonindependent (because observation from the same group are correlated with each other). With frequentist/maximum likelihood models, this issue can be solved using clustered standard errors. But I don’t think clustered/robust standard error exist in bayesian context? So there is no easy way to replicate the results of fixest package. Also, with more than one fixed effect, the centering becomes much more complex (using so-called alternating projections) and I don’t think there is a way to do it in brms. So the multilevel approach is your best best.