0

I asked this question before but never got an answer, so I am trying again and providing a sample data set so someone can tell me why I'm getting the errors I'm getting when I try implementing the Wald Test from aod and lmtest packages.

Sample data:

marital <- sample(1:5, 64614, replace = T)
race <- sample(1:3, 64614, replace = T)
educ <- sample(1:20, 64614, replace = T)
test <- data.frame(educ, marital, race)

test$marital <- as.factor(test$marital)
test$race <- as.factor(test$race)

test$marital <- relevel(test$marital, ref = "3")

require(nnet)
require(aod)
require(lmtest)

testmod <- multinom(marital ~ race*educ, data = test)
testnull <- multinom(marital ~ 1, data = test) #null model for the global test

waldtest(testnull, testmod)
wald.test(b = coef(testmod), Sigma = vcov(testmod), Terms = 1:24) #testing all terms for the global test

As you can see, when I use the waldtest function from lmtest package I get the following error: 

Error in solve.default(vc[ovar, ovar]) : 'a' is 0-diml

When I use the wald.test function from aod, I get the following error:

Error in L %*% b : non-conformable arguments

I assume these are related errors as they both seem to do with the variance matrix. I'm not sure why I'm getting these errors though as the data set has no missing values.

Michael
  • 111
  • 9
  • Those functions (`waldtest` from lmtest and `wald.test` from aod) appear to only accept linear or general linear models. The `multinom` function from the nnet package performs multinomial regression models via neural networks. You may have to perform the wald test(s) manually. – Edward Feb 22 '20 at 02:18
  • Oh... That makes infinitely more sense. Thank you. – Michael Feb 22 '20 at 04:48
  • And this should give you the p-values from Wald's test for each term: `z <- summary(testmod)$coefficients/summary(testmod)$standard.errors; (1 - pnorm(abs(z), 0, 1)) * 2` – Edward Feb 22 '20 at 04:54

1 Answers1

0

Just as heads up when using nnet package with multinom: You can also use the package broom to tidy things a bit by doing this:

tidy(multinom_model, conf.int= True, conf.level = 0.95, exponentiate = T)

This will return a tibble with the coefficients exponentiated, confidence intervals (similiar to confint used in lm), as well as the Z-scores, Standard errors and the respective p-value for the Wald Z test (essentially doing z = summary(multinom_model)$coefficients/summary(multinom_model)$standard.errors and the round((1 - pnorm(abs(z), 0, 1)) * 2,digits=5) already