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

Question

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?

Answer 1

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)
}