# How to perform model comparison based on multinom( ) function of nnet package in R?

My independent variables are gender and sequence, and the dependent variable is intervention (including 3 intervention methods). I established a multinomial logistic regression model to examine the impact of gender and sequence on intervention. The main steps are:

I first use the multinom() function to establish a null model, then establish a full model, and then use the anova() function to compare the models.

The code is as follows:

model_null<-  multinom(intervention~ 1,data = test_data)
model_full<- multinom(intervention~ 1+ gender+sequence,data = test_data)
anova(model_full,model_null,test = "Chisq")


The result output is as follows:

I have two questions:

1. I am not sure if this analysis method is correct, because this is my first time using multinomial logistic regression. Since linear regression generally uses this method for model comparison, I adopted this method;
2. If it is correct, I don’t know how to get the chi-square value for model comparison. The current output result does not have a chi-square value, only a P value.
• The LR stat. is your statistic of interest, you can confirm this via pchisq(5.364442, df=4, lower.tail=FALSE). Commented Jul 10 at 15:57
• I disagree with the closure and have voted to reopen in case someone has a better answer than mine. The answer to this question is statistical, which is why I posted my (currently accepted) answer about the underlying statistics.
– Dave
Commented Jul 11 at 2:42

This makes sense, as it is just a likelihood ratio "chunk" test of your gender and sequence variables.

The anova function does not make this explicit, but the LR stat is the $$\chi^2$$ statistic. You can see this in the lmtest package by explicitly calling a likelihood ratio test and getting the same value for the Chisq as you get for LR stat.

library(nnet)
library(lmtest)
set.seed(2024)

N <- 1000
n <- 5
p <- runif(n, 0, 1); p <- p/sum(p) # Some probabilities that sum to 1
x <- runif(N)
y <- t(rmultinom(N, 1, p))

L0 <- nnet::multinom(y ~ 1)
L1 <- nnet::multinom(y ~ x)

A <- anova(L1, L0, test = "Chisq")
L <- lmtest::lrtest(L0, L1)
A
L


The math is is that a function of the multinomial likelihoods of the two models is (asymptotically) $$\chi^2$$-distributed, discussed on pages 131-134 of Agresti (2015). You can work through the outputs of anova(L1, L0, test = "Chisq") and lmtest::lrtest(L0, L1) to show them to have equivalent representations.

L$LogLik[1] * -2 == A$Resid. Dev[1] # TRUE
L$LogLik[2] * -2 == A$Resid. Dev[2] # TRUE


You should get similar results if you work through this with your data and your models. You do not have to use the lmtest package to do this testing, though I do because I remember the lmtest::lrtest better than I remember anova(..., test = "Chisq"). Some skeptics may be more convinced if you use a package that explicitly calls the statistic the $$\chi^2$$ statistic instead of having to go through the mathematics of the deviance values and log-likelihoods.

REFERENCES

Agresti, Alan. Foundations of Linear and Generalized Linear Models. John Wiley & Sons, 2015.

• Wow! I did consider whether LR stat is the chi-square value, but GPT gave me a negative answer. . . . . . Thank you very much for your timely help and solving my urgent problem! Have a nice day! Commented Jul 10 at 13:46
• ChatGPT is not known for being a great source for statistics, particularly earlier (free) versions.
– Dave
Commented Jul 10 at 13:47