在负二项式回归中，色散参数的 MLE 与回归系数的 MLE 渐近不相关（http://pointer.esalq.usp.br/departamentos/lce/arquivos/aulas/2011/LCE5868/OverdispersionBook.pdf） .
R 中 MASS 包中的 glm.nb 函数拟合负二项式回归，并给出了离散参数的标准误差和回归系数的方差协方差矩阵，但没有给出这些之间协方差的估计，大概是因为前面提到的渐近零相关。
负二项式回归的其他实现，例如 SAS 的 PROC GENMOD 或 Stata 的 nbreg，确实报告了离散度和回归系数估计之间的协方差。
此外，glm.nb 采用的方法的结果是我认为回归系数的标准误差与 SAS 或 Stata 的标准误差不匹配，因为前者是在假设回归系数和分散参数之间独立的情况下计算的。
问题：有谁知道 R 中的另一个负二项式回归实现确实允许这种相关性，因此给出的标准误差与 SAS 或 Stata 给出的标准误差相匹配？
我还没有找到另一个可以执行此操作的 R 包，但我编写了代码，该代码基于装有 glm.nb 的模型的最大似然估计，使用观察到的信息矩阵计算完整的方差协方差矩阵。
与 SAS 的值相比，这似乎匹配，但如果有人发现错误或发现它与 SAS 或 Stata 的方差协方差矩阵不匹配，请在此答案中添加评论。
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)
}