Multilevel regression for 3-way categorical outcome

Interfaces brms
1

This is probably an easy-to-answer question. I am analyzing a 3-way categorical outcome (baseline: “Absolute” vs. “Relative” and “Other”):

brm(Response3way ~ 1 + Context + (1 | ObjectType), family = categorical(link=logit), ...)

Yielding the following output:

Group-Level Effects: 
...
Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
muRelative_Intercept      -0.41      1.68    -4.10     3.52       5826 1.00
muOther_Intercept         -0.87      1.99    -5.14     3.71       6647 1.00
muRelative_ContextWith     0.75      0.23     0.31     1.20      14885 1.00
muOther_ContextWith        1.27      0.30     0.73     1.88      15075 1.00

I understand that the estimates for the two Context parameters describe the effect of Context (sum-coded) on the respective relative odds of a) a “Relative” response compared to an “Absolute” response and b) an “Other” response compared to an “Absolute” response. Is there any way, e.g., through a hypothesis() query to ask whether Context has an effect on, for instance, the log-odds of a “Relative” responses, compared to both “Other” and “Absolute” responses? (I know I could run a separate logistic regression to address this question, but I’m wondering whether one can also obtain the answer to such questions from the model described above.)

2

You could use the refcat argument of the categorical family function to estimate parameters for all categories at once and then transform things according to your needs. Be aware though that without proper priors, this model will not be idenfitified and in any case will take much more estimation time.

3

Thank you. I looked into refcat argument but am still a bit confused, as to how I would test from such a model whether there are effects on the absolute odds of a particular level. Here’s an example with a three-way outcome prog:

library("brms")
library("foreign")
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")

summary(ml)
contrasts(ml$female) = cbind("yes" = c(-1, 1))

my.priors = c(
  prior(normal(0, 5), class = b), 
  prior(normal(0, 5), class = Intercept)
)

m = brm(
  prog ~ 1 + read * female, 
  family = categorical(link=logit, refcat = NA), 
  prior = my.priors,
  data = ml,
)

This yields a model with estimates for all three outcomes:

Family: categorical 
  Links: mugeneral = logit; muacademic = logit; muvocation = logit 
Formula: prog ~ 1 + read * female 
   Data: ml (Number of observations: 200) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
mugeneral_Intercept          27.21    169.02  -291.18   389.49         36 1.06
muacademic_Intercept         24.03    168.99  -294.35   384.63         36 1.06
muvocation_Intercept         29.61    169.00  -289.11   390.29         36 1.06
mugeneral_read               -0.54      3.23    -7.46     5.55         37 1.06
mugeneral_femaleyes           1.61      2.78    -4.03     6.65         72 1.03
mugeneral_read:femaleyes      0.12      2.86    -5.98     5.37         51 1.09
muacademic_read              -0.46      3.23    -7.35     5.61         37 1.06
muacademic_femaleyes         -0.38      2.78    -6.10     4.66         72 1.03
muacademic_read:femaleyes     0.16      2.86    -5.93     5.43         51 1.09
muvocation_read              -0.59      3.23    -7.48     5.52         37 1.06
muvocation_femaleyes         -0.39      2.80    -6.17     4.65         73 1.03
muvocation_read:femaleyes     0.16      2.86    -5.93     5.42         51 1.09

How do I interpret this output? are the effects given in each row the effects on the absolute log-odds of that outcome? I.e., mugeneral_read describes the effect of read on the log-odds of prog == “general”? I didn’t find any further documentation on this.

1 Like
4

The usual multivariate logit link looks as follows for category i

exp(eta_i) / (exp(eta_1) + exp(eta_2) + ...)

For identification we usually set one of the etas to zero, but in your case above you literally
estimate all of the etas at once.