Ошибка с VGAM? 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

НО, если я перезапущу эту отсечку, я получу решение для поснегиномиальной

> 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 (Positive Poisson: нет нулевых значений), я получу похожую ошибку "аргумент не интерпретируется как логический".

Я заметил, что в Stackru есть ряд похожих вопросов о пропущенных значениях, где требуется TRUE/FALSE, но с другими пакетами R. Это указывает на то, что, возможно, разработчикам пакетов лучше предвидеть, что вычисления могут потерпеть неудачу.

Answer 1

Я думаю, что ваша проксимальная проблема заключается в том, что предсказанные средние значения отрицательного бинома для ваших экстремальных значений настолько близки к нулю, что они снижаются до нуля, что не было предусмотрено / защищено авторами пакетов. (Одна вещь, которую нужно понять о нелинейной оптимизации / подгонке, это то, что всегда можно сломать метод подгонки, дав ему экстремальные данные...)

Я не мог заставить это работать в VGAM, но я предложу пару других предложений.

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

И просматривая данные, чтобы получить начальную оценку значений параметров (по крайней мере, для средней модели):

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

Я думал, что смогу получить vglm() работать, устанавливая начальные значения, но это на самом деле не удавалось, даже когда у меня есть довольно хорошие значения от других платформ (см. ниже).

glmmADMB

glmmADMB Пакет может обрабатывать положительные NB, через family="truncnbinom":

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