# Building hierarchical multinomial models brms

**URL:** <https://discourse.mc-stan.org/t/building-hierarchical-multinomial-models-brms/34764>\
**Category:** Modeling\
**Tags:** brms, specification\
**Created:** [April 11, 2024, 2:36pm UTC](https://discourse.mc-stan.org/t/building-hierarchical-multinomial-models-brms/34764 "2024-04-11T14:36:56Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![E\_michels](https://avatars.discourse-cdn.com/v4/letter/e/b487fb/32.png) [@E\_michels](https://discourse.mc-stan.org/u/E_michels)\
**Post date:** [April 11, 2024, 2:36pm UTC](https://discourse.mc-stan.org/t/building-hierarchical-multinomial-models-brms/34764/1 "2024-04-11T14:36:56Z")

</div>

Hi Folks, apologies if there’s already a post covering this topic.

I’m following the example presented in Solomon Kurz’s: [11 God Spiked the Integers | Statistical rethinking with brms, ggplot2, and the tidyverse: Second edition](https://bookdown.org/content/4857/god-spiked-the-integers.html) CH11.3.2, which I should first say, is great. But I’m not sure how one would go about adding random effects to such a model.

# adapted data generation.

n ← 1000

#adapted f\_income to better fit the data I’m actually working on

f\_income ← rep(-0.5:0.5, time = 500) # the treatment which I’m interested in  
culture\_random ← rep(runif(5, -0.1, 0.1), time = 200) # random effect.

#unique co-efficient for each event type

b ← c(1, 0.7, 1.3) # for fixed effect  
r ← c(1, 0.8, 0.7) # for random effect

#simulate career choices  
career ← rep(NA, n)

for( i in 1:n) {  
score ← 0.5 \* (1:3) + b \* f\_income[i] + r \* culture\_random[i]  
p ← softmax(score[1], score[2], score[3])  
career[i] ← sample( 1:3, size = 1, prob = p)  
}

d ← tibble(  
career = career,  
f\_income = f\_income,  
cr = rep(1:5, time = 200)  
)

# models

multi1 ← # this model runs fine, as expected, its almost the same as the example model  
brm(data = d,  
family = categorical(link = logit, refcat = 3),  
bf(career ~ 1,  
nlf(mu1 ~ a1 + b1 \* f\_income),  
nlf(mu2 ~ a2 + b2 \* f\_income),  
a1 + a2 + b1 + b2 ~ 1),  
prior = c(prior(normal(0, 1.5), class = b, nlpar = a1),  
prior(normal(0, 1.5), class = b, nlpar = a2),  
prior(normal(0, 1), class = b, nlpar = b1),  
prior(normal(0, 1), class = b, nlpar = b2)),  
iter = 2000, warmup = 1000, cores = 4, chains = 4)

multi2 ← # how I would imagine adding random effects doesn’t work  
brm(data = d,  
family = categorical(link = logit, refcat = 3),  
bf(career ~ 1,  
nlf(mu1 ~ a1 + b1 \* f\_income + (a1|cr)),  
nlf(mu2 ~ a2 + b2 \* f\_income + (a2|cr)),  
a1 + a2 + b1 + b2 ~ 1),  
prior = c(prior(normal(0, 1.5), class = b, nlpar = a1),  
prior(normal(0, 1.5), class = b, nlpar = a2),  
prior(normal(0, 1), class = b, nlpar = b1),  
prior(normal(0, 1), class = b, nlpar = b2)),  
iter = 2000, warmup = 1000, cores = 4, chains = 4)

Any help, advice or signposting would be greatly appreciated,  
Many thanks

---

<div class="post-metadata">

**Author:** ![wpetry](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wpetry/32/4852_2.png) [@wpetry](https://discourse.mc-stan.org/u/wpetry)\
**Post date:** [April 11, 2024, 3:55pm UTC](https://discourse.mc-stan.org/t/building-hierarchical-multinomial-models-brms/34764/2 "2024-04-11T15:55:48Z")

</div>

I don’t think `nlf()` allows the typical `(gterms | group)` syntax. It might be easier to re-specify your model in the “terse” or “verbose” syntax (following Solomon’s nomenclature). Alternatively, you could estimate the group-level effects on a separate line and simply add them into the non-linear formula. I haven’t thought carefully about the right way to do this, but it would look something like:

```no-highlight
multi2 <- brm(data = d,
              family = categorical(link = logit, refcat = 3),
              bf(career ~ 1,
                 c1 ~ 0 + (1 | cr),
                 c2 ~ 0 + (1 | cr),
                 nlf(mu1 ~ a1 + b1 * f_income + c1),
                 nlf(mu2 ~ a2 + b2 * f_income + c2),
                 a1 + a2 + b1 + b2 ~ 1),
              prior = c(prior(normal(0, 1.5), class = b, nlpar = a1),
                        prior(normal(0, 1.5), class = b, nlpar = a2),
                        prior(normal(0, 1), class = b, nlpar = b1),
                        prior(normal(0, 1), class = b, nlpar = b2)),
              iter = 2000, warmup = 1000, cores = 4, chains = 4)

```

---

<div class="post-metadata">

**Author:** ![E\_michels](https://avatars.discourse-cdn.com/v4/letter/e/b487fb/32.png) [@E\_michels](https://discourse.mc-stan.org/u/E_michels)\
**Post date:** [April 11, 2024, 4:32pm UTC](https://discourse.mc-stan.org/t/building-hierarchical-multinomial-models-brms/34764/3 "2024-04-11T16:32:27Z")

</div>

Thankyou, your solution did indeed work.