Building hierarchical multinomial models brms

Modeling
,
1

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 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

2

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:

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)
3

Thankyou, your solution did indeed work.