如何使用 R 中的 rms 包进行负二项式回归？答案

【问题标题】：How to do negative binomial regression with the rms package in R?如何使用 R 中的 rms 包进行负二项式回归？
【发布时间】：2021-10-11 14:37:08
【问题描述】：

如何使用 R 中的 rms 包来执行负二项式回归？（I originally posted this question on Statistics SE，但它显然已关闭，因为它更适合这里。）

使用MASS 包，我使用glm.nb 函数，但我试图切换到rms 包，因为在使用glm.nb 和其他一些函数引导时有时会出现奇怪的错误。但我无法弄清楚如何使用 rms 包进行负二项式回归。

这是我想做的示例代码（复制自rms::Glm 函数文档）：

library(rms)
## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
f <- Glm(counts ~ outcome + treatment, family=poisson())

f
anova(f)
summary(f, outcome=c('1','2','3'), treatment=c('1','2','3'))

所以，我不想使用family=poisson()，而是使用family=negative.binomial() 之类的东西，但我不知道该怎么做。

在family {stats} 的文档中，我在“另请参阅”部分找到了这条注释：

对于二项式系数，选择；二项式和负二项式分布、Binomial 和 NegBinomial。

但即使点击了?NegBinomial 的链接，我也无法理解这一点。

对于如何使用 R 中的 rms 包执行负二项式回归，我将不胜感激。

【问题讨论】：

查看@Tripartio下面的答案

标签： r regression non-linear-regression

【解决方案1】：

基于this，以下似乎可行：

library(rms)
library(MASS)

counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)

Glm(counts ~ outcome + treatment, family = negative.binomial(theta = 1))
General Linear Model
 
 rms::Glm(formula = counts ~ outcome + treatment, family = negative.binomial(theta = 1))
 
                    Model Likelihood    
                          Ratio Test    
    Obs       9    LR chi2      0.31    
 Residual d.f.4    d.f.            4    
    g 0.2383063    Pr(> chi2) 0.9892    
 
             Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept    3.0756 0.2121 14.50  <0.0001 
 outcome=2   -0.4598 0.2333 -1.97  0.0487  
 outcome=3   -0.2962 0.2327 -1.27  0.2030  
 treatment=2 -0.0347 0.2333 -0.15  0.8819  
 treatment=3 -0.0503 0.2333 -0.22  0.8293

【讨论】：

抱歉延迟回复；我正在度假。我很确定 theta == 1 不是我想要的。正如您链接的文章所解释的那样，这是几何模型的简单假设。 MASS::glm.nb 能够自动估计 theta。我的问题本质上是如何在事先不知道 theta 的情况下使用 rms 包运行负二项式（因为事先几乎不知道它）。

【解决方案2】：

预先提出意见您最好发布（作为一个单独的问题）一个可重现的示例，说明您的引导尝试中的“奇怪错误”，看看人们是否有解决这些问题的想法。当数据分散均匀或分散不足时，NB 拟合程序抛出警告或错误是相当常见的，因为在这种情况下分散参数的估计值变得无限......

@coffeinjunky 是正确的，使用family = negative.binomial(theta=VALUE) 将起作用（其中VALUE 是一个数字常量，例如theta=1 用于几何分布[NB 的一种特殊情况]）。但是：您将无法（无需大量工作）拟合一般的 NB 模型，即在拟合过程中估计分散参数 (theta) 的模型。这就是 MASS::glm.nb 所做的，而 AFAICS 在 rms 包中没有类似物。

除了 MASS::glm.nb 之外，还有一些其他的包/函数适合负二项式模型，包括（至少）bbmle 和 glmmTMB — 可能还有其他的，例如 gamlss。

## Dobson (1990) Page 93: Randomized Controlled Trial :
dd < data.frame(
   counts = c(18,17,15,20,10,20,25,13,12)
   outcome = gl(3,1,9),
   treatment = gl(3,3))

MASS::glm.nb

library(MASS)
m1 <- glm.nb(counts ~ outcome + treatment, data = dd)
## "iteration limit reached" warning

glmmTMB

library(glmmTMB)
m2 <- glmmTMB(counts ~ outcome + treatment, family = nbinom2, data = dd)
## "false convergence" warning

bbmle

library(bbmle)
m3 <- mle2(counts ~ dnbinom(mu = exp(logmu), size = exp(logtheta)),
     parameters = list(logmu ~outcome + treatment),
     data = dd,
     start = list(logmu = 0, logtheta = 0)
)

signif(cbind(MASS=coef(m1), glmmTMB=fixef(m2)$cond, bbmle=coef(m3)[1:5]), 5)
                   MASS     glmmTMB       bbmle
(Intercept)  3.0445e+00  3.04540000  3.0445e+00
outcome2    -4.5426e-01 -0.45397000 -4.5417e-01
outcome3    -2.9299e-01 -0.29253000 -2.9293e-01
treatment2  -1.1114e-06  0.00032174  8.1631e-06
treatment3  -1.9209e-06  0.00032823  6.5817e-06

这些都非常一致（至少对于拦截/结果参数）。这个例子对于 NB 模型来说是相当困难的（5 个参数 + 9 个观测值的离散度，数据是泊松而不是 NB）。

【讨论】：