VGAM有Bug吗？vglm family=posnegbinomial =>“if (take.half.step) {：需要TRUE/FALSE时缺少值”

Question

我有一些实际的数据，我担心这些数据有点令人讨厌。

它本质上是一个正负二项分布(没有任何零计数)。然而，有一些异常值似乎会导致一些错误的计算发生(可能是下溢或NaNs?)前8个条目是合理的，但我猜最后几个条目造成了一些拟合问题。

以下是数据：

> df
   counts  t
1    1968  1
2     217  2
3      55  3
4      26  4
5      11  5
6       5  6
7       8  7
8       3  8
9       1 10
10      1 11
11      1 12
12      1 13
13      1 15
14      1 18
15      1 26
16      1 59

此命令运行一段时间，然后显示错误消息

> vglm(counts ~ t, data=df, family = posnegbinomial)
Error in if (take.half.step) { : missing value where TRUE/FALSE needed

但是，如果我重新运行这个截断离群值的方法，我就得到了一个关于solution二项式的解决方案

> vglm(counts ~ t, data=df[1:9,], family = posnegbinomial)
Call:
vglm(formula = counts ~ t, family = posnegbinomial, data = df[1:9,])

Coefficients:
(Intercept):1 (Intercept):2             t 
    7.7487404     0.7983811    -0.9427189 

Degrees of Freedom: 18 Total; 15 Residual
Log-likelihood: -36.21064 

如果我尝试家族pospoisson (正泊松:没有零值)，我得到一个类似的错误“参数不能解释为逻辑”。

我确实注意到在Stackoverflow中有许多类似的问题，关于需要TRUE/FALSE的缺失值，但对于其他R包。这向我表明，也许包编写者需要更好地预测计算可能会失败。

Answer 1

我认为你最近的问题是，你的极值的负二项式的预测均值非常接近于零，以至于它们下溢到零，以一种包作者没有预料到/防止的方式。(关于非线性优化/拟合，要认识到的一件事是，总是有可能通过提供极端数据来打破拟合方法……)

我不能让这个在VGAM中工作，但我会提供一些其他的建议。

plot(log(counts)~t,data=dd)

并观察数据以获得参数值的初始估计(至少对于均值模型)：

m0 <- lm(log(counts)~t,data=subset(dd,t<10))

我认为我可以通过设置初始值来让vglm()工作，但这实际上并不起作用，即使我从其他平台获得了相当好的值(见下文)。

glmmADMB

通过family="truncnbinom"，glmmADMB包可以处理正NB

library(glmmADMB)
m1 <- glmmadmb(counts~t, data=dd, family="truncnbinom")

(有一些警告消息...)

bbmle::mle2()

这需要更多的工作:它在标准模型中失败了，但如果我在预测平均值上设置下限，则可以工作……

library(VGAM)  ## for dposnegbin
library(bbmle)
m2 <- mle2(counts~dposnegbin(size=exp(logk),
                         munb=pmax(exp(logeta),1e-7)),
           parameters=list(logeta~t),
           data=dd,
           start=list(logk=0,logeta=0))

再次警告消息。

比较glmmADMB，mle2，simple truncated lm fit ...

cc <- cbind(coef(m2),
  c(log(m1$alpha),coef(m1)),
  c(NA,coef(m0)))
dimnames(cc) <- list(c("log_k","log_int","slope"),

                 c("mle2","glmmADMB","lm"))

##               mle2   glmmADMB         lm
## log_k    0.8094678  0.8094625         NA
## log_int  7.7670604  7.7670637  7.1747551
## slope   -0.9491796 -0.9491778 -0.8328487

原则上，使用glmmTMB也可以做到这一点，但它遇到了与vglm()相同的问题……