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

Question

如何使用 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} 的文档中，我在“另请参阅”部分找到了此注释:

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

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

如果您能提供有关如何使用 R 中的 rms 包执行负二项式回归的帮助，我将不胜感激。

Answer 1

预先提出意见您可能最好发布(作为一个单独的问题)您的引导尝试中的“奇怪错误”的可重现示例，并看看人们是否有解决这些问题的想法。当数据分散或分散不足时，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))

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