`zeroinfl` из пакетов `pscl` и `countreg` дают очень разные результаты. Почему?
Я экспериментирую с mdvis набором данных из COUNT пакета R в учебных целях. Я подогнал отрицательно-биномиальную модель с нулевым раздуванием, используя функцию zeroinfl из пакетов pscl и countreg. Однако результаты zeroinfl из пакета pscl и из пакета countreg сильно отличаются.
Модели и результаты представлены ниже.
zeroinfl из pscl:
mdvisit_zeroinf<-pscl::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F)
summary(mdvisit_zeroinf)
Call:
pscl::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc |
reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin")
Pearson residuals:
Min 1Q Median 3Q Max
-0.9667 -0.8037 -0.3597 0.3154 9.4878
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.30316 0.56565 -0.536 0.5920
reform -0.09916 0.05543 -1.789 0.0736 .
badh 1.11555 0.07517 14.841 <2e-16 ***
agegrp2 0.11376 0.06990 1.628 0.1036
agegrp3 0.21126 0.08620 2.451 0.0143 *
educ2 0.07605 0.07209 1.055 0.2914
educ3 -0.03979 0.07295 -0.545 0.5854
loginc 0.13209 0.07499 1.761 0.0782 .
Log(theta) 0.04747 0.05815 0.816 0.4144
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -30.7110 772.2930 -0.040 0.968
reform 14.3274 763.7283 0.019 0.985
badh -1.9503 3.5392 -0.551 0.582
agegrp2 10.9236 113.6292 0.096 0.923
agegrp3 9.7723 113.6558 0.086 0.931
educ2 -0.6632 1.9878 -0.334 0.739
educ3 -0.8743 1.3501 -0.648 0.517
loginc 0.5437 1.9718 0.276 0.783
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 1.0486
Number of iterations in BFGS optimization: 106
Log-likelihood: -4557 on 17 Df`
zeroinfl из countreg:
mdvisit_zeroinf2<-countreg::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F))
summary(mdvisit_zeroinf2)
Call:
countreg::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc |
reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin",
control = zeroinfl.control(method = "BFGS", EM = F))
Pearson residuals:
Min 1Q Median 3Q Max
-0.9523 -0.8057 -0.3615 0.3106 9.6300
Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.44473 0.53752 -0.827 0.40803
reform -0.12962 0.05106 -2.538 0.01113 *
badh 1.13558 0.07413 15.319 < 2e-16 ***
agegrp2 0.05851 0.06359 0.920 0.35756
agegrp3 0.18678 0.07176 2.603 0.00925 **
educ2 0.06398 0.06621 0.966 0.33385
educ3 -0.04825 0.07087 -0.681 0.49599
loginc 0.15338 0.07056 2.174 0.02973 *
Log(theta) 0.01232 0.04778 0.258 0.79652
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -119.137 95.259 -1.251 0.2111
reform 4.821 2.646 1.822 0.0684 .
badh -21.461 4614.989 -0.005 0.9963
agegrp2 9.261 28.341 0.327 0.7438
agegrp3 3.502 27.984 0.125 0.9004
educ2 -19.790 203.802 -0.097 0.9226
educ3 -7.894 4.758 -1.659 0.0971 .
loginc 13.010 10.463 1.243 0.2137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Theta = 1.0124
Number of iterations in BFGS optimization: 197
Log-likelihood: -4554 on 17 Df
В чем может быть причина такой разницы? Результаты Stata очень похожи на результаты zeroinfl из пакета pscl (особенно результаты для счетной части модели).
Я подогнал модели точно так же, указав одинаковые параметры с помощью zeroinfl.control(). Я также попытался найти в Интернете, чтобы узнать, сообщали ли другие о подобных проблемах ранее, и если ответы уже доступны. Я искал ответы также на stackoverflow и CV. Но я не мог получить ни одного.
@BenBolker Я изучил данные и подозреваю, что может быть полное разделение из-за переменной loginc, которая в основном содержит уникальные значения для каждого случая. Кросс-таблица дает для каждого уникального значения loginc количество либо zero, либо non-zero. Но мне интересно, как одна и та же функция только из двух пакетов приводит к разным результатам как для счетной, так и для бинарной частей модели.
Я пишу разработчикам по электронной почте, когда возникают подобные проблемы. Я всегда получал ответы от разработчиков, с которыми я связывался, которые отвечали на мой вопрос. Разработчики, как правило, весьма полезны на этом фронте.
@LeroyTyrone Хорошо. Попробую и это.
Спасибо, что указали мне на это, я добавил ответ сейчас. Дополнительное замечание: функции в обоих пакетах изначально были написаны мной, а countreg поддерживается мной, но содержит мою первоначальную реализацию. Несколько лет назад сопровождающий pscl интегрировал патч в zeroinfl(), который иногда может привести к немного лучшей конвергенции, а иногда нет (как здесь). Поэтому результаты могут немного отличаться. Но я не видел случаев, когда различия действительно имеют практическое значение (включая ваш пример). Здесь они просто выглядят иначе, но на самом деле их не так много.
Ответы 1
Источник проблемы:
Проблема в том, что нет нулевой инфляции, но на самом деле меньше нулей, чем ожидалось от отрицательной биномиальной модели, особенно в некоторых подгруппах по отношению к reform и agegrp.
Симптомы:
Из-за отсутствия нулевой инфляции бинарная часть с нулевой инфляцией изо всех сил пытается получить предсказанные вероятности нулевой инфляции, которые очень близки к нулю для определенных подгрупп. Обратите особое внимание на перехват в той части модели, которая чрезвычайно мала с большой стандартной ошибкой.
В целом это похоже на квазиполное разделение в моделях бинарной регрессии. Вероятность становится очень плоской и зависит от настроек оптимизатора, где именно она останавливается (и эти настройки различаются между pscl и countreg).
Поскольку два компонента модели нулевой инфляции (количество и двоичная часть) не могут быть оценены по отдельности, проблемы с оценкой одного компонента могут привести к проблемам/различиям в оценке другого компонента.
Альтернатива:
Модель барьеров имеет здесь несколько преимуществ: (1) Она может иметь дело не только с нулевой инфляцией, но и с отсутствием нулей. (2) Два компонента модели могут быть оценены по отдельности, и, следовательно, проблемы не могут распространяться друг на друга.
Иллюстрация:
Давайте сравним базовую отрицательную биномиальную модель с ее аналогами с нулевой инфляцией и барьером:
data("mdvis", package = "COUNT")
mdvis <- transform(mdvis,
reform = factor(reform),
badh = factor(badh),
agegrp = factor(agegrp),
educ = factor(educ)
)
library("countreg")
f <- numvisit ~ reform + badh + agegrp + educ + loginc
m <- glm.nb(f, data = mdvis)
m2 <- zeroinfl(f, dist = "negbin", data = mdvis)
m3 <- hurdle(f, dist = "negbin", data = mdvis)
С точки зрения логарифмической вероятности модель с нулевой инфляцией лишь незначительно улучшает базовую модель, несмотря на то, что ей требуется почти в два раза больше параметров. Модель барьеров улучшается немного больше, но тоже ненамного:
c(logLik(m), logLik(m2), logLik(m3))
## [1] -4560.910 -4554.146 -4550.520
С точки зрения как AIC, так и BIC модель с нулевой инфляцией является наихудшей из этих трех моделей. AIC немного предпочитает барьерную модель, в то время как BIC несколько более явно предпочитает базовую модель.
AIC(m, m2, m3)
## df AIC
## m 9 9139.819
## m2 17 9142.293
## m3 17 9135.040
BIC(m, m2, m3)
## df BIC
## m 9 9191.195
## m2 17 9239.336
## m3 17 9232.083
Глядя на кореньограмму как на диагностический дисплей для базовой модели, мы видим, что в целом наблюдаемых нулей меньше, чем ожидалось от отрицательной биномиальной модели. Но в целом базовая модель уже подходит достаточно хорошо.
rootogram(m)
(Примечание. Версия рутограммы с доверительным интервалом на самом деле создается другим разрабатываемым пакетом. Версия в countreg не имеет таких доверительных интервалов.)
Большое спасибо. Это превосходно отвечает на мой вопрос, развеивает мои сомнения и дает мне новое понимание данных и методов.
Единственное, что бросается в глаза, это явное полное разделение компонента нулевой инфляции (любой коэффициент с абсолютным значением больше 10 является здесь красным флажком). Не знаю, должно ли/почему это испортить оценку компонента счетчика модели, но кажется разумным, что это может быть.