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
