How to interpret the summary of a multinomial model?

Question

I have a model with 5 multinomial levels coded 0 to 4:

library(mgcv)

gam = gam(
  list(
    level ~ generation * gender + s(long, lat) + s(id, bs = "re"),
    ~ generation * gender + s(long, lat) + s(id, bs = "re"),
    ~ generation * gender + s(long, lat) + s(id, bs = "re"),
    ~ generation * gender + s(long, lat) + s(id, bs = "re")
  ),
  data = data,
  family = multinom(K = 4)
)

summary(gam)

How can interpret the summary?:

Family: multinom 
Link function: 

Formula:
level ~ generation * gender + s(long, lat) + s(id, bs = "re")
~generation * gender + s(long, lat) + s(id, bs = "re")
~generation * gender + s(long, lat) + s(id, bs = "re")
~generation * gender + s(long, lat) + s(id, bs = "re")

Parametric coefficients:
                                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)                       -2.2018     0.2435  -9.042  < 2e-16 ***
generationYounger                 -0.2989     0.2890  -1.034  0.30096    
genderFemale                      -0.0975     0.2870  -0.340  0.73403    
generationYounger:genderFemale     0.3714     0.4097   0.907  0.36460    
(Intercept).1                     -3.5076     0.3136 -11.185  < 2e-16 ***
generationYounger.1               -0.2449     0.4032  -0.607  0.54353    
genderFemale.1                     0.3558     0.3952   0.900  0.36793    
generationYounger:genderFemale.1  -1.0224     0.6592  -1.551  0.12088    
[...]
generationYounger:genderFemale.3  -0.7765     1.1310  -0.687   0.4924    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

[...]

A multinomial model compares every level from 1 to x to the reference level 0 and has therefore x formulae. Do the coefficients with no extention in the summary correspond to level 1 compared to level 0 and .1 to level 2 compared to level 0 etc.? So is following relation correct?:

level number	coding for model	extention in summary
1	0	[not in summary]
2	1	[none]
3	2	.1
4	3	.2
5	4	.3

And does generationYounger:genderFemale.1 correspond to the effect of Young and Female to manifest as level 2 compared to Older, Male, and level 0 (globally; to the intercept of the reference level) or compared to Older, Male, and level 2 (its own intercept)?

Answer 1

With a large number of coefficients like this, it's safest to rely on predict() functions to get estimates and standard errors for particular scenarios of interest. Those calculations need to use information about the covariances of coefficient estimates that aren't part of the standard summary but are kept within the model. If mgcv allows for calculations in the probability ("response") scale from a multinomial model, that would provide a more useful summary of your results for your audience.

That said, the "parametric coefficients" you show do have straightforward interpretations, similar to coefficients in logistic regression, in the context of the smooths performed in the generalized additive model. The difference from logistic regression is that the coefficients in these log-odds scales are for each group relative to the single reference group. That is, if $y$ represents the (numbered but not ordered) category and $y=0$ is the reference, the coefficients are for $\log(\pi_{y=i}/\pi_{y=0})$, where $\pi_{y=i}$ is the probability of being in outcome class $i$. That's more complicated than the binary outcome situation, in which you can simply substitute $\pi_{y=0}= 1-\pi_{y=1}$.

The numbering of the multiple models is a bit awkward, but you seem to have them identified properly. Thus the model with extension .1 has coefficients related to $\log(\pi_{y=2}/\pi_{y=0})$, in the coding of levels provided to the model.

I'm not completely clear about how you think the coefficients for each model should be interpreted, however. With the default predictor coding in R, each is a difference (in log-odds scale here) from a lower-level situation. For the model with extension .1, that means:

The (Intercept).1 of -3.5076 is the log-odds, $\log(\pi_{y=2}/\pi_{y=0})$, when both predictors are at reference levels (presumably Older,Male).

The generationYounger.1 coefficient of -0.2449 is the difference in log odds from the Intercept for generationYounger and Male. The log-odds, $\log(\pi_{y=2}/\pi_{y=0})$, for that scenario is -3.5076-0.2449.

The genderFemale.1 coefficient of 0.3558 is the difference in log odds from the Intercept when generationOlder holds but gender is Female. The log-odds, $\log(\pi_{y=2}/\pi_{y=0})$, for that scenario is -3.5076+0.3558.

The generationYounger:genderFemale.1 interaction coefficient of -1.0224 is the extra difference from the above individual coefficients when you have both generationYounger and genderFemale. The log-odds, $\log(\pi_{y=2}/\pi_{y=0})$, for that scenario is -3.5076-0.2449+0.3558-1.0224.

To get standard errors in the log-odds scale you then have to use the covariances of the coefficient estimates, and to translate to probabilities you have to take all the coefficients together with the constraint imposed by the model that the probabilities across classes for any set of predictor values adds up to 1. Hence the suggestion to rely on predict() functions rather than on model coefficient summaries.