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)

応答のないスレッドで、このエラーに関する言及が 1 つだけ見つかりました。

glm.nb何がこのエラーを引き起こしている可能性があり、どのように修正できますか (エラーをスローするたびに入力をランダムに並べ替える以外に?)

ETA: 設定は、アルゴリズムが非常に大きくなり、次に になり、次に になることによって壊れてcontrol=glm.control(maxit=200,trace = 3)いることを発見しました:theta.ml-InfNaN

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）

いくつかの突っ込み/プロファイリングは、シータパラメータのMLEが事実上無限であることを示しています。境界を0に設定できるように、逆スケールで近似することにしました（より洗練されたバージョンでは、theta = zeroでポアソンに戻る対数尤度関数が設定されますが、これにより、迅速な缶詰の解決策を考え出す）。

上記の2つの悪い例では、これはかなりうまく機能しますが、パラメーターの適合が境界上にあることを警告しています...

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))

2番目の例は、負の二項偏差から正確に生成された適切なサイズのデータ​​セットがあるにもかかわらず、シータの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 に答える