我有一些实际数据,恐怕有些讨厌。
它本质上是一个正负二项分布(没有任何零计数(。 但是,有一些异常值似乎会导致一些错误的计算发生(可能是下溢或 NaN? 前 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
但是,如果我重新运行这个切断异常值,我会得到一个 posnegbinomial 的解决方案
> 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 包。 这向我表明,也许包编写者需要更好地预测计算可能会失败。
我认为您的近端问题是您的极值的负二项式的预测均值非常接近于零,以至于它们以包作者没有预料到/保护的方式溢出到零。(关于非线性优化/拟合,要意识到的一件事是,总是可以通过提供极端数据来破坏拟合方法......
我无法让它在VGAM
中工作,但我会提供一些其他建议。
plot(log(counts)~t,data=dd)
并目测数据以获得参数值的初始估计值(至少对于平均模型(:
m0 <- lm(log(counts)~t,data=subset(dd,t<10))
我以为我可以通过设置起始值来让vglm()
工作,但这实际上并没有成功,即使我从其他平台获得了相当好的值(见下文(。
格拉姆ADMB
glmmADMB
包可以通过family="truncnbinom"
处理正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
,简单的截断lm
适合...
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()
相同的问题......