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.