16

glm.nb在某些输入上引发异常错误。虽然有多种值会导致此错误，但即使是非常轻微地更改输入也可以防止错误发生。

一个可重现的例子：

set.seed(11)
pop <- rnbinom(n=1000,size=1,mu=0.05)
glm.nb(pop~1,maxit=1000)

运行此代码会引发错误：

Error in while ((it <- it + 1) < limit && abs(del) > eps) { : 
  missing value where TRUE/FALSE needed

起初我认为这与算法不收敛有关。然而，我惊讶地发现即使是非常轻微地改变输入也可以防止错误。例如：

pop[1000] <- pop[1000] + 1
glm.nb(pop~1,maxit=1000)

我发现它会在 1 到 500 之间的 19.4% 的种子上引发此错误：

fit.with.seed = function(s) {
    set.seed(s)
    pop <- rnbinom(n=1000, size=1, mu=0.05)
    m = glm.nb(pop~1, maxit=1000)
}

errors = sapply(1:500, function(s) {
    is.null(tryCatch(fit.with.seed(s), error=function(e) NULL))
})

mean(errors)

我发现在任何地方都只提到了这个错误，在一个没有响应的线程上。

什么可能导致此错误，以及如何修复它（除了每次glm.nb引发错误时随机排列输入？）

ETA：设置control=glm.control(maxit=200,trace = 3)发现theta.ml算法通过变得非常大而中断，然后变为-Inf，然后变为NaN：

theta.ml: iter67 theta =5.77203e+15
theta.ml: iter68 theta =5.28327e+15
theta.ml: iter69 theta =1.41103e+16
theta.ml: iter70 theta =-Inf
theta.ml: iter71 theta =NaN
4

1 回答 1

8

这有点粗糙，但在过去，我已经能够glm.nb通过采用直接最大似然估计来解决问题（即没有像 中使用的聪明的迭代估计算法glm.nb）

一些探索/分析表明 theta 参数的 MLE 实际上是无限的。我决定把它放在反比例上，这样我就可以把一个边界放在 0 处（一个更高级的版本会设置一个对数似然函数，它会在 theta=0 处恢复为泊松，但这会取消尝试想出一个快速的罐头解决方案）。

对于上面给出的两个不好的例子，这工作得相当好，尽管它确实警告参数拟合在边界上......

library(bbmle)
m1 <- mle2(Y~dnbinom(mu=exp(logmu),size=1/invk),
           data=d1,
           parameters=list(logmu~X1+X2+offset(X3)),
           start=list(logmu=0,invk=1),
           method="L-BFGS-B",
           lower=c(rep(-Inf,12),1e-8))

第二个例子实际上更有趣，因为它在数值上证明了 theta 的 MLE 本质上是无限的，即使我们有一个完全由负二项式偏差生成的大数据集（或者我对某些事情感到困惑......）

set.seed(11);pop <- rnbinom(n=1000,size=1,mu=0.05);glm.nb(pop~1,maxit=1000)
m2 <- mle2(pop~dnbinom(mu=exp(logmu),size=1/invk),
           data=data.frame(pop),
           start=list(logmu=0,invk=1),
           method="L-BFGS-B",
           lower=c(-Inf,1e-8))
于 2012-08-01T02:20:48.127 回答