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

Question

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

对于 MASS 包，我使用 glm.nb 函数，但我正在尝试切换到 rms 包，因为在使用 [=13 引导时有时会出现奇怪的错误=] 和一些其他功能。但我无法弄清楚如何使用 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} 的文档中，我在“另请参阅”部分找到了这条注释：

For binomial coefficients, choose; the binomial and negative binomial distributions, Binomial, and NegBinomial.

但即使在点击 ?NegBinomial 的 link 之后，我也无法理解这一点。

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

Answer 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

Answer 2

前面的意见 你最好张贴（作为一个单独的问题）你的 bootstrap 尝试中的“奇怪错误”的可重现示例，看看是否人们有解决这些问题的想法。当数据分散或分散不足时，NB 拟合程序会发出警告或错误是相当常见的，因为在这种情况下分散参数的估计值变为无限...

@coffeinjunky 是正确的，使用 family = negative.binomial(theta=VALUE) 会起作用（其中 VALUE 是一个数字常数，例如 theta=1 用于几何分布 [NB 的特例]）。但是：您将无法（没有更多的工作）能够拟合一般的 NB 模型，即色散参数 (theta) 被估计为拟合的一部分的模型程序。这就是 MASS::glm.nb 所做的，AFAICS rms 包中没有类似物。

除了MASS::glm.nb之外，还有其他一些packages/functions符合负二项式模型，包括（至少）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))

质量::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

这些都非常吻合（至少对于 intercept/outcome 参数而言）。这个例子对于 NB 模型来说相当困难（5 个参数 + 9 个观测值的分散，数据是泊松而不是 NB）。

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

How to do negative binomial regression with the rms package in R?

r

regression

non-linear-regression

质量::glm.nb

glmmTMB

bbmle