Thanks to the helpful folks in this thread (especially Paul), I’ve now got a multinomial GLMM to analyze, with varying intercepts for one group-level factor. There are 4 response categories c_1...c_4, with c_1 the arbitrarily chosen baseline category. The random effects are ignored at this stage of the analysis, which focuses on quantifying the population-level effects.

As most readers will know, multinomial models are notoriously hard to interpret, because unlike in binomial or linear regression, interpretation of a given coefficient for any c is non-monotonous and also depends on that same predictor’s coefficients for all the other c's.

Let x_1 be a binary predictor of interest in a multinomial model containing many predictors. Fox&Hong and Long, who are the only authors I’ve found giving detailed advice on the interpretation of multinomial models, suggest analyzing x_1's effect by fixing all the other x's at their means (if quantitative) or proportions (if factor levels), then calculating the difference between the fitted probabilities of each c when x_1 = 1 versus when x_1 = 0. The result is called the marginal effect at the mean (MEM).

The average marginal effect (AME), on the other hand, is supposedly (Long 2012: 25) computed by going through the i_1...i_n data points, once for each c, and

- calculating the fitted probability of c for subject i when x_1 = x_{1i}
- calculating the hypothetical fitted probability of c when x_{1i} is changed one unit (from 0 to 1 or vice versa)
- dividing the difference by the sign of the change in x_{1i}

The AME of x_1 on a given response category is the average of this difference over the entire dataset.

Does it make sense (or is it feasible) to calculate a confidence interval for the MEM or the AME? That would look real good in an effect display plot. But I haven’t seen such things anywhere. What comes closest is the `effects`

package in R, which computes confidence intervals for the estimated category probabilities at x_1 = 1 and x_1 = 0, holding the other x's at their means. But that’s not a confidence interval for the MEM itself, is it? It’s just a pair of two confidence intervals from whose possible overlap we may infer whether the MEM differs “significantly” from zero. In particular, using these two confidence intervals, rather than a single one for the MEM itself, forces us to plot two points (namely \hat{\mu}(y = c|x_1 = 1) and \hat{\mu}(y = c|x_1 = 0)) per category per predictor, which takes a lot of plot space when there are very many predictors and categories.

I’ve written an R function that calculates a variant of AME – one that bases inference only on those pairs of predictor combinations that A) actually occur in the data and B) contain information about x_1, so nothing is extrapolated to counterfactual data. But I’m not sure how (or even **if**) the confidence intervals of these AMEs should be computed. ATM, it is done by setting x_1's coefficient to the upper(lower) boundary of its CI, then going through every informative pair of predictor combinations in the data and calculating each \hat{\mu}(y = c|x_1 = 1) using that upper(lower) boundary value for x_1's coefficient. There are no explicit coefficients for effects on c_1, so the confidence interval for x_1's effect on c_1 is calculated by setting the coefficients of c_2, c_3 and c_4 simultaneously to the upper limit of their 70.7% CI’s (1-\sqrt[3]{\frac{.05}{2}} = .707). But this is likely wrong, because the resulting intervals look funky:

The baseline category always has the smallest confidence interval, which hardly makes sense unless there’s an error.

**Does anyone have an idea of how the CIs could be calculated correctly?**

Since my own approach doesn’t seem to work so far, I also wonder if the `effects()`

function could be made to work or replicated for `brms`

models. The function uses something called the *delta method* to calculate SEs for fitted multinomial probabilities computed from a multinomial logit model fit by `nnet`

's `multinom()`

function. How could we calculate them for a `brms`

multinomial GLMM? Page 22 in the aforementioned article provides the details, but as a non-STEM student I don’t understand all the mathematical notation.

As expected, the `effects()`

function currently fails when applied to a `brms`

model, complaining:

`Error in terms.default(model) : no terms component nor attribute`

From what I do understand of the cited article, the function requires access to at least the estimated beta coefficients and their covariance matrix in order to work, so if I could either:

A) figure out how to extract the necessary components from a `brms`

multinomial model to give `effects()`

what it needs to work, or

B) figure out how exactly the SEs of fitted probabilities from a multinomial model are calculated,

then I could probably write a function to do all that automatically. That would be an alternate solution.

**References:**

Fox, J. and J. Hong (2009). *Effect displays in R for multinomial and proportional-odds logit models: Extensions to the effects package*. Journal of Statistical Software 32:1, 1–24, http://www.jstatsoft.org/v32/i01/.

Long, S. (2012) *Regression Models for Nominal and Ordinal Outcomes*. Forthcoming in Best and Wolf (editors), Regression Models , Sage Publications