Wald Test for negative binomial regression

[StarryCatte]

I'm trying to run a negative binomial regression using fenegbin() function. Next, the function linearHypothesis in car package is used to test whether there is significant difference between two coefficients. I recieved the error information: "Error in L %*% V : non-conformable arguments". However, such error didn't occur when using fepois(). How should I do to deal with it?

library(fixest)
library(car)

data(mtcars)

## it works well with fepois
mod_pois <- fepois(hp ~ drat+wt, data=mtcars)
linearHypothesis(mod_pois, 'drat=wt')

# Linear hypothesis test

# Hypothesis:
# drat - wt = 0

# Model 1: restricted model
# Model 2: hp ~ drat + wt

#   Df  Chisq Pr(>Chisq)    
# 1                         
# 2  1 96.621  < 2.2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## using fenegbin and an error ocurred
mod_negbin <- fenegbin(hp ~ drat+wt, data=mtcars)
linearHypothesis(mod_negbin, 'drat=wt')

# Error in L %*% V : non-conformable arguments

[kylebutts]

Hi @StarryCatte,

library(fixest)
library(car)
#> Loading required package: carData
library(marginaleffects)

mod_pois <- fepois(hp ~ drat + wt, data = mtcars)
mod_negbin <- fenegbin(hp ~ drat+wt, data=mtcars)
vcov_negbin <- vcov(mod_negbin)
vcov_negbin
#>              (Intercept)         drat            wt        .theta
#> (Intercept)  0.686233453 -0.132369379 -0.0641033566 -0.0070288964
#> drat        -0.132369379  0.027630576  0.0102607140  0.0028403353
#> wt          -0.064103357  0.010260714  0.0084307030 -0.0008987525
#> .theta      -0.007028896  0.002840335 -0.0008987525  6.7745504686

The error comes from the vcov(mod_negbin) including .theta. You can fix that by dropping that row/column:

linearHypothesis(mod_negbin, 'drat=wt', vcov. =vcov_negbin[1:3, 1:3])
#> Linear hypothesis test
#> 
#> Hypothesis:
#> drat - wt = 0
#> 
#> Model 1: restricted model
#> Model 2: hp ~ drat + wt
#> 
#> Note: Coefficient covariance matrix supplied.
#> 
#>   Df  Chisq Pr(>Chisq)   
#> 1                        
#> 2  1 8.8224   0.002975 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Additionally, the function hypotheses from the marginaleffects package works without needing to do that:

hypotheses(mod_pois, "drat=wt")
#> 
#>     Term Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
#>  drat=wt   -0.306     0.0311 -9.83   <0.001 73.3 -0.366 -0.245
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

hypotheses(mod_negbin, 'drat=wt')
#> Warning in names(model[["coefficients"]]) == names(coefs): longer object length
#> is not a multiple of shorter object length

#> Warning in names(model[["coefficients"]]) == names(coefs): longer object length
#> is not a multiple of shorter object length

#> Warning in names(model[["coefficients"]]) == names(coefs): longer object length
#> is not a multiple of shorter object length

#> Warning in names(model[["coefficients"]]) == names(coefs): longer object length
#> is not a multiple of shorter object length
#> 
#>     Term Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
#>  drat=wt    -0.37      0.125 -2.97  0.00298 8.4 -0.615 -0.126
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high

^{Created on 2023-12-11 with reprex v2.0.2}

[StarryCatte]

Yes, i knew hypothesesfunction could do that, and just want to know how to use linearHypothesis function to get a chisq statistics. It works for me, thanks!

[grantmcdermott]

@StarryCatte Please close this issue if you think your problem has been addressed.