Multinomial logit regression in R: how to estimate confidence intervals for a NULL model and from model averaging? The main question relates to the nnet package and how to assess the predicted uncertainties of the category-specific probabilities obtained with the multinom() function, including those from a NULL model or from a model averaging approach.
But first, the context:
A length-based fishing regulation ("harvest slot") was implemented in 2011 for Walleye (Sander vitreus) in the St. Lawrence River, Québec, Canada, and relying on annual surveys conducted between 2002 and 2021, we are interested to model the probability to sample a small (SM: Total Length < 381 mm), medium (MD: 381 to 545 mm) or large (LG: > 545 mm) Walleye over this period at different sampling sites. The harvest slot allows to keep only MD walleye, whereas SM and LG ones are catch and release only since 2011.
We are especially interested to know whether the probability to sample a MD or LG walleye has increased over time during this period, or if we can see a change in these probabilities in recent years that can cause the long-term trend to move upward.
The multinom() function of the nnet package was used to perform a multinomial logit regression, analyzing all 3 classes in a same model and using SM walleye as the reference category. I'm providing an example below for the Lac Saint-Pierre (LSPI) from the NORTH shore of the St. Lawrence River were used.
The variables used are:
SIZE: response categorical variable with 3 categories (SM, MD, LG) from individual walleye data
YEAR0: independent continuous variable "centered" on 2000, such that YEAR = 2021 is YEAR0 = 21
The R script used is:
library(nnet)
summary(m.LSPI_NORTH_SIZE.multinom<-multinom(SIZE~YEAR0,data=SPECIMENS_SAVI_LSPI_NORTH,Hess=TRUE))
The predictions of this first model including YEAR0 as the sole predictor (AIC = 1959.745) are shown in the figure below.
We are also considering the NULL version of this multinomial model, i.e. testing whether an apparent stability is more likely instead by using SIZE ~ 1 (not shown, but this model has AIC = 1960.432).
Thus, the delta AIC between these two models is relatively low at 0.687, suggesting that both alternatives have some statistical support. Based on second-order Akaike weight, SIZE ~ YEAR0 gets 58.2% of support, whereas SIZE ~ 1 obtains 41.8%.
Here are the observed annual proprotions for each SIZE class (filled circles) together with the size-specific (matching colors) predicted probabilites from the SIZE ~ YEAR (solid line) and SIZE ~ 1 (dashed line) multinomial logit models:

Model averaging would clearly provide (from my point of view) a more likely statistical description of the biological process being examined.
If I use the MNLpred package with its mnl_pred_ova() function, I can obtain through a simulation process the lower and upper 95%CI of the model including YEAR0, but not for the NULL model nor a compromise obtained with the model.sel() and model.avg() functions of the MuMIn package.
If the NULL model would have obtained 100% of the second-order Akaike weight, I would not be able to plot the central tendency together with its uncertainty.
As the NULL model gets nearly equivalent statistical support in this case, the compromise I get with MuMIn allows me to get the coefficients to then draw the central tendency for each SIZE class (not shown), but again without its associated uncertainty.
So, my main question in the end is:
Is there a way to obtain the predicted uncertainty (95%CI) from a multinomial model for each considered category when the NULL version is retained or when we opt instead for model averaging (either with MuMIn or other package) due to support found for both alternatives as it's the case here?
As a second question:
Would the mgcv package with its gam() function and family=multinomial be an interesting alternative for such case and if so, is there an available example other than that from the mgcv package document that could be accessible to illustrate how to obtain category-specific probabilities and their respective uncertainties?
And as a last one:
What is the best available goodness-of-fit test that we can apply to a multinomial logit model? To be more specific: to make sure that the predicted values reflect the observed ones in a satisfactory manner (i.e., model adequacy).
Thanks for your help in advance.
Julien
 A: It looks like the best option may be to "decompose" the multinomial logit regression into separate logistic regressions (see Hosmer and Lemeshow 2000) and this link:

*

*probability of SM|(MD+LG)

*probability of MD|(SM+LG)

*probability of LG|(SM+MD)

which will produce nearly (if not) identical predicted probabilities as those obtained from the multinomial logit model above.
In a second step, one may use the predict() function with the arguments type="link" and se.fit=TRUE to then calculate the 95% CI on the link scale, simply by multiplying the SE by + or - 1.96 and then back-transform the values obtained on the response scale with plogis().
The three 95% CI's generated this way would however be somewhat independent, which is not ideal, but with this approach, one could get an approximate idea of the associated uncertainty for each size category through time - which may be better than nothing.
This approach would thus also allow to estimate the 95% CI of an alternative NULL model as well as any model averaging performed in the MuMIn package with model.avg(). In the last case, using full=TRUE with the predict function() would be required.
I'm pretty sure that using the same approach with the mgcv package, the gam() function and the argument family=binomial(link=logit) without a smoothing function applied to YEAR0 and using method="ML" instead of "REML" would give the same results and also allow to compute an 95% CI for each size category, again with the predict() function as above.
Regarding model adequacy, I would look with LogisticDx package and its gof() function provided each of the 3 logistic regressions obtained with the simpler glm() function offers a sufficiently adequate fit.
To judge if the fit is adequate, and given that the independent variable is continuous, use the standardized Pearson Chi-squared test of Osius and Rojek (1992). Furthermore, make sure that the logit link is correctly specified with the Stukel (1988) score test, which is also simultaneously performed by the gof() function (with other tests too, being a "wrapper" of adequacy tests).
In both cases a $p$-value < 0.05 is indicative of a lack-of-fit or misspecified link. If all three separate logistic regressions are adequate, the multinomial logit regression is likely to be adequate too. If not, one may want to be more careful with the interpretation of the (overall) multinomial model.
In the example I've originally provided, using s() for YEAR0 with a number of knots varying from a minimum of $k = 3$ to a maximum of $k = 8$ (with 8 sampling years) with the gam() function of mgcv may offer a better fit than the original multinomial model obtained with the nnet package, but I'm really not sure of how the adequacy could be checked (i.e., predicted vs. observed values) for each binomial regressions that are allowed to "wiggle".
Any suggestion or comment to help better quantify the uncertainty of different multinomial logit models is welcome!
Cheers,
Julien
