Zero-inflated counts — R cross-check (pscl::zeroinfl)¶

Companion to zip_python.ipynb / zip_pymc.ipynb. The pscl package fits zero-inflated models by maximum likelihood — the standard frequentist benchmark for the fish data (this is the UCLA worked example). zeroinfl(count ~ count_covariates | zero_covariates, dist=...): the part before | is the count model (log link), the part after is the structural-zero model (logit link), matching the from-scratch $\beta$ and $\gamma$.

In [1]:
ok <- suppressWarnings(suppressMessages(require(pscl)))
if (!ok) { install.packages('pscl', repos='https://cloud.r-project.org'); library(pscl) }
d <- read.csv('fish.csv')
zip  <- zeroinfl(count ~ child + camper | persons, data = d, dist = 'poisson')
zinb <- zeroinfl(count ~ child + camper | persons, data = d, dist = 'negbin')
cat('================ ZIP ================\n');  print(summary(zip)$coefficients)
cat('\n================ ZINB ===============\n'); print(summary(zinb)$coefficients)
cat(sprintf('\nZINB theta (= r, dispersion) = %.3f\n', zinb$theta))

Installing package into 'C:/Users/user/R/win-library/4.6'
(as 'lib' is unspecified)


package 'pscl' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\user\AppData\Local\Temp\RtmpEBnoIj\downloaded_packages

Classes and Methods for R originally developed in the
Political Science Computational Laboratory
Department of Political Science
Stanford University (2002-2015),
by and under the direction of Simon Jackman.
hurdle and zeroinfl functions by Achim Zeileis.


================ ZIP ================

$count
              Estimate Std. Error    z value     Pr(>|z|)
(Intercept)  1.5978883 0.08553816  18.680413 7.146347e-78
child       -1.0428391 0.09998825 -10.429616 1.816209e-25
camper       0.8340236 0.09362679   8.907959 5.198064e-19

$zero
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  1.2974403  0.3738523  3.470463 0.0005195626
persons     -0.5643474  0.1629638 -3.463023 0.0005341420


================ ZINB ===============

$count
              Estimate Std. Error   z value     Pr(>|z|)
(Intercept)  1.3710482  0.2561131  5.353291 8.636874e-08
child       -1.5152547  0.1955913 -7.747047 9.405440e-15
camper       0.8790513  0.2692732  3.264534 1.096444e-03
Log(theta)  -0.9853536  0.1759500 -5.600191 2.141163e-08

$zero
             Estimate Std. Error   z value   Pr(>|z|)
(Intercept)  1.603106  0.8365073  1.916428 0.05531064
persons     -1.666566  0.6792845 -2.453413 0.01415076


ZINB theta (= r, dispersion) = 0.373
In [2]:
# compact comparison table (count = before |, zero = after |)
cz <- coef(zip); cn <- coef(zinb)
cat(sprintf('%-26s %10s %10s\n','parameter','ZIP','ZINB'))
cat(sprintf('%-26s %10.3f %10.3f\n','count: intercept', cz['count_(Intercept)'], cn['count_(Intercept)']))
cat(sprintf('%-26s %10.3f %10.3f\n','count: child',     cz['count_child'],       cn['count_child']))
cat(sprintf('%-26s %10.3f %10.3f\n','count: camper',    cz['count_camper'],      cn['count_camper']))
cat(sprintf('%-26s %10.3f %10.3f\n','zero:  intercept', cz['zero_(Intercept)'],  cn['zero_(Intercept)']))
cat(sprintf('%-26s %10.3f %10.3f\n','zero:  persons',   cz['zero_persons'],      cn['zero_persons']))
cat(sprintf('%-26s %10s %10.3f\n','dispersion theta(r)','-', zinb$theta))

# vuong-style fit: predicted vs observed P(0)
cat(sprintf('\nObserved P(0) = %.3f   ZIP P(0) = %.3f   ZINB P(0) = %.3f\n',
    mean(d$count==0), mean(predict(zip,type='prob')[,1]), mean(predict(zinb,type='prob')[,1])))

parameter                         ZIP       ZINB

count: intercept                1.598      1.371

count: child                   -1.043     -1.515

count: camper                   0.834      0.879

zero:  intercept                1.297      1.603

zero:  persons                 -0.564     -1.667

dispersion theta(r)                 -      0.373

Observed P(0) = 0.568   ZIP P(0) = 0.539   ZINB P(0) = 0.587

Results¶

  • pscl::zeroinfl reproduces the from-scratch and PyMC fits: ZIP child ≈ −1.04, camper ≈ +0.83, zero-persons ≈ −0.56; ZINB child ≈ −1.5, camper ≈ +0.9, θ = r ≈ 0.37. (These are the canonical UCLA numbers.)
  • The ZINB zero-model coefficients (intercept, persons) are the most uncertain part — once the NB dispersion explains much of the excess zeros, the structural-zero population is only weakly identified, so the three engines agree most tightly on the count model and r, and more loosely on the zero model. The Bayesian notebooks show this directly as wide posteriors there.
  • Three engines, one mixture: from-scratch augmentation Gibbs ≈ PyMC (marginalised NUTS) ≈ R pscl (ML). Both Bayesian routes additionally quantify the zero-vs-overdispersion trade-off through the joint posterior.