Regresión binomial negativa en R que permite la correlación entre la dispersión y los coeficientes de regresión

Question

En la regresión binomial negativa, la MLE del parámetro de dispersión no está correlacionada asintóticamente con las MLE de los coeficientes de regresión ( http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf ).

La función glm.nb en el paquete MASS en R ajusta la regresión binomial negativa, y da errores estándar para el parámetro de dispersión y una matriz de covarianza de varianza para los coeficientes de regresión, pero no da una estimación de la covarianza entre estos, presumiblemente debido a la correlación cero asintótica antes mencionada.

Otras implementaciones de la regresión binomial negativa, por ejemplo PROC GENMOD de SAS o nbreg de Stata, sí informan de la covarianza entre la dispersión y las estimaciones del coeficiente de regresión.

Además, una consecuencia del enfoque adoptado por glm.nb es que creo que los errores estándar de los coeficientes de regresión no coinciden con los de SAS o Stata, porque los primeros se calculan asumiendo la independencia entre los coeficientes de regresión y el parámetro de dispersión.

PREGUNTA: ¿conoce alguien otra implementación de regresión binomial negativa en R que tenga en cuenta esta correlación y, por lo tanto, proporcione errores estándar que coincidan con los de SAS o Stata?

Answer 1

No he encontrado otro paquete de R que haga esto, pero he escrito un código que, basándose en las estimaciones de máxima verosimilitud de un modelo ajustado con glm.nb, calcula la matriz de covarianzas de varianza completa utilizando la matriz de información observada.

Comparando con los valores de SAS parece coincidir, pero si alguien detecta un error o encuentra que no coincide con la matriz de covarianza de varianza de SAS o Stata, por favor, añada un comentario a esta respuesta.

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