Negative binomial regression in R allowing for correlation between dispersion & regression coefficients

In negative binomial regression, the MLE of the dispersion parameter is asymptotically uncorrelated with the MLEs of the regression coefficients (http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf).

The glm.nb function in the MASS package in R fits negative binomial regression, and gives standard errors for the dispersion parameter and a variance covariance matrix for the regression coefficients, but does not give an estimate of the covariance between these, presumably because of the aforementioned asymptotic zero correlation.

Other implementations of negative binomial regression, e.g. SAS's PROC GENMOD or Stata's nbreg, do report the covariance between the dispersion and regression coefficient estimates.

Moreover, a consequence of the approach taken by glm.nb is I think that the standard errors for the regression coefficients do not match those from SAS or Stata, because the former are calculated assuming independence between the regression coefficients and the dispersion parameter.

QUESTION: does anyone know of another negative binomial regression implementation in R that does allow for this correlation, and therefore gives standard errors that match those given by SAS or Stata?

• (1) I suspect PROC GENMOD & nbreg may calculate standard errors from the expected Fisher information rather than the observed. (2) See ?glm.nb- the algorithm alternates between optimizing the estimate of the dispersion keeping coefficients fixed & optimizing the estimates of coefficients keeping dispersion fixed - there's no implicit assumption that they're uncorrelated involved. – Scortchi - Reinstate Monica Jul 1 '16 at 10:35
• I don't think they do - SAS and Stata give you a covariance between the dispersion and regression coefficient parameters. This covariance has (as per my post) expectation zero, and so if they were giving the expected Fisher information these terms would I think be zero? – Jonathan Bartlett Jul 1 '16 at 10:47
• I think you're right. Some potentially useful refs for answerers: Lawless (1987), "Negative binomial & mixed Poisson regression", Canadian J. Statistics, 15, 3, gives the observed & expected information matrices, & the parametrization used in glm.nb is explained here. – Scortchi - Reinstate Monica Jul 1 '16 at 11:24

I haven't found another R package which does this, but I have written code which, based on the maximum likelihood estimates of a model fitted with glm.nb, calculates the full variance covariance matrix using the observed information matrix.

Comparing to values from SAS this appears to match, but if anyone spots an error or finds that it does not match the variance covariance matrix from SAS or Stata, please add a comment to this answer.

glm.nb.cov <- function(mod) {
#given a model fitted by glm.nb in MASS, this function returns a variance covariance matrix for the
#regression coefficients and dispersion parameter, without assuming independence between these
#note that the model must have been fitted with x=TRUE argument so that design matrix is available

#formulae based on p23-p24 of http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf
#and http://www.math.mcgill.ca/~dstephens/523/Papers/Lawless-1987-CJS.pdf

k <- mod$theta #p is number of regression coefficients p <- dim(vcov(mod))[1] #construct observed information matrix obsInfo <- array(0, dim=c(p+1, p+1)) #first calculate top left part for regression coefficients for (i in 1:p) { for (j in 1:p) { obsInfo[i,j] <- sum( (1+mod$y/mod$theta)*mod$fitted.values*mod$x[,i]*mod$x[,j] / (1+mod$fitted.values/mod$theta)^2  )
}
}

#information for dispersion parameter
obsInfo[(p+1),(p+1)] <- -sum(trigamma(mod$theta+mod$y) - trigamma(mod$theta) - 1/(mod$fitted.values+mod$theta) + (mod$theta+mod$y)/(mod$theta+mod$fitted.values)^2 - 1/(mod$fitted.values+mod$theta) + 1/mod$theta)

#covariance between regression coefficients and dispersion
for (i in 1:p) {
obsInfo[(p+1),i] <- -sum(((mod$y-mod$fitted.values) * mod$fitted.values / ( (mod$theta+mod$fitted.values)^2 )) * mod$x[,i] )
obsInfo[i,(p+1)] <- obsInfo[(p+1),i]
}

#return variance covariance matrix
solve(obsInfo)
}

• (+1) Looks like I had it the wrong way round, & that glm.nb uses expected info. for std errors of the coefficients (for a dispersion parameter fixed at the MLE) & observed info. just for the std error of the dispersion parameter, whereas proc genmod uses observed info. throughout. You could check your code's results against a numerical evaluation of the Hessian (numDeriv::hessian). The negative binomial isn't the only model commonly used with a non-canonical link so this issue is of wider interest. – Scortchi - Reinstate Monica Jul 4 '16 at 13:23
• Great and very useful. For anyone using it, remember to specify x=TRUE in the glm.nb call. Also, you can use the VCOV variable in the EMMEANS call to get LS-Means type estimates (just remember to omit the part of the covariance matrix that refers to $\hat{\theta}$ - i.e. the last row and column, because EMMEANS only wants to the variance-covariance matrix for the fixed effects of the model). – Björn Jan 23 '18 at 14:21