Zero-augmented models

\[ \newcommand{\E}{\mathsf{E}} \newcommand{\VAR}{\mathsf{VAR}} \newcommand{\COV}{\mathsf{COV}} \newcommand{\Prob}{\mathsf{P}} \]

1 Hurdle models

In addition to over-dispersion, many empirical count data sets exhibit more zero observations than would be allowed for by the Poisson model. One model class capable of capturing both properties is the hurdle model, originally proposed by Mullahy (1986) in the econometrics literature (see Cameron and Trivedi 1998; Cameron and Trivedi 2005 for an overview). They are two-component models: A truncated count component, such as Poisson, geometric or negative binomial, is employed for positive counts, and a hurdle component models zero vs. larger counts. For the latter, either a binomial model or a censored count distribution can be employed.

More formally, the hurdle model combines a count data model \(f_\mathrm{count}(y; x, \beta)\) (that is left-truncated at \(y = 1\)) and a zero hurdle model \(f_\mathrm{zero}(y; z, \gamma)\) (right-censored at \(y = 1\)): \[ f_\mathrm{hurdle}(y; x, z, \beta, \gamma) = \left\{ \begin{array}{ll} f_\mathrm{zero}(0; z, \gamma) & \mbox{if } y = 0, \\ (1 - f_\mathrm{zero}(0; z, \gamma)) \cdot f_\mathrm{count}(y; x, \beta)/(1 - f_\mathrm{count}(0; x, \beta)) & \mbox{if } y > 0 \end{array} \right. \tag{1}\] The model parameters \(\beta\), \(\gamma\), and potentially one or two additional dispersion parameters \(\theta\) (if \(f_\mathrm{count}\) or \(f_\mathrm{zero}\) or both are negative binomial densities) are estimated by ML, where the specification of the likelihood has the advantage that the count and the hurdle component can be maximized separately. The corresponding mean regression relationship is given by \[ \log(\mu_i) \quad = \quad x_i^\top \beta + \log(1 - f_\mathrm{zero}(0; z_i, \gamma)) - \log(1 - f_\mathrm{count}(0; x_i, \beta)), \tag{2}\] again using the canonical log link. For interpreting the zero model as a hurdle, a binomial GLM is probably the most intuitive specification¹. Another useful interpretation arises if the same regressors \(x_i = z_i\) are used in the same count model in both components \(f_\mathrm{count} = f_\mathrm{zero}\): A test of the hypothesis \(\beta = \gamma\) then tests whether the hurdle is needed or not.

In R, hurdle count data models can be fitted with the hurdle() function from the countreg package (Zeileis and Kleiber 2013). Both its fitting function and the returned model objects of class hurdle are modelled after the corresponding GLM functionality in R. The arguments of hurdle() are given by

hurdle(formula, data, subset, na.action, weights, offset,
  dist = "poisson", zero.dist = "binomial", link = "logit",
  control = hurdle.control(...),
  model = TRUE, y = TRUE, x = FALSE, ...)

where the first line contains the standard model-frame specifications, the second and third lines have the arguments specific to hurdle models and the arguments in the last line control some components of the return value.

If a formula of type y ~ x1 + x2 is supplied, it not only describes the count regression relationship of \(y_i\) and \(x_i\) but also implies that the same set of regressors is used for the zero hurdle component \(z_i = x_i\). This is could be made more explicit by equivalently writing the formula as y ~ x1 + x2 | x1 + x2. Of course, a different set of regressors could be specified for the zero hurdle component, e.g., y ~ x1 + x2 | z1 + z2 + z3, giving the count data model y ~ x1 + x2 conditional on (|) the zero hurdle model y ~ z1 + z2 + z3.

The model likelihood can be specified by the dist, zero.dist and link arguments. The count data distribution dist is "poisson" by default (it can also be set to "negbin" or "geometric"), for which the canonical log link is always used. The distribution for the zero hurdle model can be specified via zero.dist. The default is a binomial model with link (defaulting to "logit", but all link functions of the binomial() family are also supported), alternatively a right-censored count distribution (Poisson, negative binomial or geometric, all with log link) could be specified.

ML estimation of all parameters employing analytical gradients is carried out using R’s optim() with control options set in hurdle.control(). Starting values can be user-supplied, otherwise they are estimated by glm.fit() (the default). The covariance matrix estimate is derived numerically using the Hessian matrix returned by optim(). See Section 4 for further technical details.

The returned fitted-model object of class hurdle is a list similar to glm objects. Some of its elements—such as coefficients or terms—are lists with a zero and count component, respectively. For details see Section 4.

A set of standard extractor functions for fitted model objects is available for objects of class hurdle, including the usual summary() method that provides partial Wald tests for all coefficients. No anova() method is provided, but the general coeftest(), waldtest() from lmtest, and linearHypothesis() from car can be used for Wald tests and lrtest() from lmtest for LR tests of nested models. The function hurdletest() is a convenience interface to linearHypothesis() for testing for the presence of a hurdle (which is only applicable if the same regressors and the same count distribution are used in both components).

2 Zero-inflated models

Zero-inflated models (Mullahy 1986; Lambert 1992) are another model class capable of dealing with excess zero counts (see Cameron and Trivedi 1998; 2005 for an overview). They are two-component mixture models combining a point mass at zero with a count distribution such as Poisson, geometric or negative binomial. Thus, there are two sources of zeros: zeros may come from both the point mass and from the count component. For modeling the unobserved state (zero vs. count), a binary model is used: in the simplest case only with an intercept but potentially containing regressors.

Formally, the zero-inflated density is a mixture of a point mass at zero \(I_{\{0\}}(y)\) and a count distribution \(f_\mathrm{count}(y; x, \beta)\). The probability of observing a zero count is inflated with probability \(\pi = f_\mathrm{zero}(0; z, \gamma)\): \[ f_\mathrm{zeroinfl}(y; x, z, \beta, \gamma) \quad = \quad f_\mathrm{zero}(0; z, \gamma) \cdot I_{\{0\}}(y) \; + \; (1 - f_\mathrm{zero}(0; z, \gamma)) \cdot f_\mathrm{count}(y; x, \beta), \tag{3}\] where \(I(\cdot)\) is the indicator function and the unobserved probability \(\pi\) of belonging to the point mass component is modelled by a binomial GLM \(\pi = g^{-1}(z^\top \gamma)\). The corresponding regression equation for the mean is \[ \mu_i \quad = \quad \pi_i \cdot 0 \; + \; (1 - \pi_i) \cdot \exp(x_i^\top \beta), \tag{4}\] using the canonical log link. The vector of regressors in the zero-inflation model \(z_i\) and the regressors in the count component \(x_i\) need not to be distinct—in the simplest case, \(z_i = 1\) is just an intercept. The default link function \(g(\pi)\) in binomial GLMs is the logit link, but other links such as the probit are also available. The full set of parameters of \(\beta\), \(\gamma\), and potentially the dispersion parameter \(\theta\) (if a negative binomial count model is used) can be estimated by ML. Inference is typically performed for \(\beta\) and \(\gamma\), while \(\theta\) is treated as a nuisance parameter even if a negative binomial model is used.

In R, zero-inflated count data models can be fitted with the zeroinfl() function from the countreg package. Both the fitting function interface and the returned model objects of class zeroinfl are almost identical to the corresponding hurdle() functionality and again modelled after the corresponding GLM functionality in R. The arguments of zeroinfl() are given by

zeroinfl(formula, data, subset, na.action, weights, offset,
  dist = "poisson", link = "logit", control = zeroinfl.control(...),
  model = TRUE, y = TRUE, x = FALSE, ...)

where all arguments have almost the same meaning as for hurdle(). The main difference is that there is no zero.dist argument: a binomial model is always used for distribution in the zero-inflation component.

Again, ML estimates of all parameters are obtained from optim(), with control options set in zeroinfl.control() and employing analytical gradients. Starting values can be user-supplied, estimated by the expectation maximization (EM) algorithm, or by glm.fit() (the default). The covariance matrix estimate is derived numerically using the Hessian matrix returned by optim(). Using EM estimation for deriving starting values is typically slower but can be numerically more stable. It already maximizes the likelihood, but a single optim() iteration is used for determining the covariance matrix estimate. See Section 5 for further technical details.

The returned fitted model object is of class zeroinfl whose structure is virtually identical to that of hurdle models. As above, a set of standard extractor functions for fitted model objects is available for objects of class zeroinfl, including the usual summary() method that provides partial Wald tests for all coefficients. Again, no anova() method is provided, but the general functions coeftest() and waldtest() from lmtest, as well as linearHypothesis() from car can be used for Wald tests, and lrtest() from lmtest for LR tests of nested models.

3 Illustrations

In the following, we illustrate hurdle and zero-inflated models by applying them to the cross-sectional data set used in Deb and Trivedi (1997). It is based on the US National Medical Expenditure Survey (NMES) for 1987/88 and is available from the data archive of the Journal of Applied Econometrics at https://journaldata.zbw.eu/dataset/demand-for-medical-care-by-the-elderly-a-finite-mixture-approach. It was prepared for the R package AER accompanying Kleiber and Zeileis (2008) and is also available as DebTrivedi.rda in the Journal of Statistical Software together with Zeileis (2006). The same data set is used to illustrate some basic count data models in the article “Basics”, which contains a more detailed description and a basic exploratory analysis of the data.

The objective is to model the number of physician office visits visits using the health status variables health (self-perceived health status), chronic (number of chronic conditions), as well as the socio-economic variables gender, school (number of years of education), and insurance (private insurance indicator) as regressors. For convenience, we select the variables used from the full data set:

data("NMES1988", package = "AER")
dt <- NMES1988[, c(1, 7, 8, 13, 15, 18)]

At the end, we will provide a brief comparison between the zero-adjusted models and the basic count data models illustrated in “Basics”.

3.1 Hurdle regression

The exploratory analysis in “Basics” conveyed the impression that there might be more zero observations than explained by the basic count data distributions, hence a negative binomial hurdle model is fitted via

fm_hurdle0 <- hurdle(visits ~ ., data = dt, dist = "negbin")

This uses the same type of count data model as in the preceeding section but it is now truncated for visits < 1 and has an additional hurdle component modeling zero vs. count observations. By default, the hurdle component is a binomial GLM with logit link which contains all regressors used in the count model. The associated coefficient estimates and partial Wald tests for both model components are displayed via

summary(fm_hurdle0)

Call:
hurdle(formula = visits ~ ., data = dt, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.1362 -0.7069 -0.2790  0.3066 17.1420 

Count model coefficients (truncated negbin with log link):
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      1.194286   0.060182  19.844  < 2e-16 ***
healthpoor       0.382645   0.048855   7.832 4.79e-15 ***
healthexcellent -0.368340   0.067493  -5.457 4.83e-08 ***
chronic          0.147537   0.012654  11.659  < 2e-16 ***
gendermale      -0.056464   0.033144  -1.704  0.08846 .  
school           0.021187   0.004632   4.574 4.79e-06 ***
insuranceyes     0.130454   0.043511   2.998  0.00272 ** 
Log(theta)       0.269981   0.043070   6.268 3.65e-10 ***
Zero hurdle model coefficients (binomial with logit link):
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.07459    0.13919   0.536   0.5920    
healthpoor       0.04995    0.15969   0.313   0.7544    
healthexcellent -0.30883    0.14256  -2.166   0.0303 *  
chronic          0.55842    0.04505  12.395  < 2e-16 ***
gendermale      -0.40145    0.08737  -4.595 4.33e-06 ***
school           0.05821    0.01196   4.867 1.13e-06 ***
insuranceyes     0.74670    0.10059   7.424 1.14e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.3099
Number of iterations in BFGS optimization: 14 
Log-likelihood: -1.215e+04 on 15 Df

The coefficients in the count component resemble those from the previous models, but the increase in the log-likelihood (see also Table 1) conveys that the model has improved by including the hurdle component. However, it might be possible to omit the health variable from the hurdle model. To test this hypothesis, the reduced model is fitted via

fm_hurdle <- hurdle(visits ~ . | chronic + insurance + school + gender,
  data = dt, dist = "negbin")

and can then be compared to the full model in a Wald test

waldtest(fm_hurdle0, fm_hurdle)

Wald test

Model 1: visits ~ .
Model 2: visits ~ . | chronic + insurance + school + gender
  Res.Df Df  Chisq Pr(>Chisq)  
1   4391                       
2   4393 -2 4.8785    0.08723 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

or an LR test

lrtest(fm_hurdle0, fm_hurdle)

which leads to virtually identical results.

3.2 Zero-inflated regression

A different way of augmenting the negative binomial count model fm_nbin with additional probability weight for zero counts is a zero-inflated negative binomial (ZINB) regression. The default model is fitted via

fm_zinb0 <- zeroinfl(visits ~ ., data = dt, dist = "negbin")

As for the hurdle model above, all regressors from the count model are also used in the zero-inflation model. Again, we can modify the regressors in the zero-inflation part, e.g., by fitting a second model

fm_zinb <- zeroinfl(visits ~ . | chronic + insurance + school + gender,
  data = dt, dist = "negbin")

that has the same variables in the zero-inflation part as the hurdle component in fm_hurdle. By omitting the health variable, the fit does not change significantly which can again be brought out by a Wald test

waldtest(fm_zinb0, fm_zinb)

Wald test

Model 1: visits ~ .
Model 2: visits ~ . | chronic + insurance + school + gender
  Res.Df Df  Chisq Pr(>Chisq)
1   4391                     
2   4393 -2 0.0937     0.9542

or an LR test lrtest(fm_zinb0, fm_zinb) that produces virtually identical results. The chosen fitted model can again be inspected via

summary(fm_zinb)

See Table 1 for a more concise summary.

Table 1: Summary of fitted count regression models for NMES data (including the models from “Basics”): coefficient estimates from count model, zero-inflation model (both with standard errors in parantheses) , number of estimated parameters, maximized log-likelihood, AIC, BIC and expected number of zeros (sum of fitted densities evaluated at zero). The observed number of zeros is 683 in 4406 observations.

	ML-Pois	Adj-Pois	Quasi-Pois	NB	NB-Hurdle	ZINB
(Intercept)	1.035	1.035	1.035	0.940	1.194	1.198
	(0.024)	(0.065)	(0.063)	(0.055)	(0.060)	(0.057)
healthpoor	0.318	0.318	0.318	0.368	0.383	0.349
	(0.017)	(0.056)	(0.046)	(0.049)	(0.049)	(0.046)
healthexcellent	-0.379	-0.379	-0.379	-0.374	-0.368	-0.354
	(0.030)	(0.078)	(0.080)	(0.062)	(0.067)	(0.061)
chronic	0.169	0.169	0.169	0.196	0.148	0.150
	(0.004)	(0.012)	(0.012)	(0.012)	(0.013)	(0.012)
gendermale	-0.108	-0.108	-0.108	-0.115	-0.056	-0.072
	(0.013)	(0.036)	(0.034)	(0.032)	(0.033)	(0.032)
school	0.026	0.026	0.026	0.027	0.021	0.022
	(0.002)	(0.005)	(0.005)	(0.004)	(0.005)	(0.004)
insuranceyes	0.216	0.216	0.216	0.250	0.130	0.151
	(0.017)	(0.043)	(0.045)	(0.040)	(0.044)	(0.042)
(Intercept)					0.054	-0.098
					(0.137)	(0.274)
chronic					0.577	-1.302
					(0.043)	(0.191)
insuranceyes					0.742	-1.193
					(0.100)	(0.228)
school					0.056	-0.087
					(0.012)	(0.027)
gendermale					-0.406	0.583
					(0.087)	(0.204)
AIC	36597.0	36597.0		24469.9	24331.1	24336.9
BIC	36641.7	36641.7		24521.0	24414.2	24419.9
Log.Lik.	-18291.494	-18291.494		-12226.953	-12152.560	-12155.431
no. of parameters	7	7	8	8	13	13
$\sum_{i} \hat{f}_{i}(0)$	46			618	683	712

3.3 Comparison

Having fitted hurdle and zero-inflated models to the demand for medical care in the NMES data, it is, of course, of interest to understand what these models have in common and what their differences are, especially in comparison to the basic count data models described in “Basics”. In this section, we show how to compute the components of Table 1 and provide some further comments and interpretations.

Before we start the analysis, we refit the models from “Basics”:

fm_pois <- glm(visits ~ ., data = dt, family = poisson)
fm_qpois <- glm(visits ~ ., data = dt, family = quasipoisson)
fm_nbin <- glm.nb(visits ~ ., data = dt)

As a first comparison, it is of natural interest to inspect the estimated regression coefficients in the count data model

fm <- list("ML-Pois" = fm_pois, "Quasi-Pois" = fm_qpois, "NB" = fm_nbin,
  "Hurdle-NB" = fm_hurdle, "ZINB" = fm_zinb)
sapply(fm, function(x) coef(x)[1:7])

The result (see Table 1) shows that there are some small differences, especially between the GLMs and the zero-augmented models. However, the zero-augmented models have to be interpreted slightly differently: While the GLMs all have the same mean function (\(g(\mu_i) = x_i^\top \beta\), see “Basics”), the zero-augmentation also enters the mean function, see (4) and (2). Nevertheless, the overall impression is that the estimated mean functions are rather similar. Moreover, the associated estimated standard errors are very similar as well (see Table 1):

cbind("ML-Pois" = sqrt(diag(vcov(fm_pois))),
  "Adj-Pois" = sqrt(diag(sandwich(fm_pois))),
  sapply(fm[-1], function(x) sqrt(diag(vcov(x)))[1:7]))

The only exception are the model-based standard errors for the Poisson model, when treated as a fully specified model, which is obviously not appropriate for this data set.

In summary, the models are not too different with respect to their fitted mean functions. The differences become obvious if not only the mean but the full likelihood is considered:

rbind(logLik = sapply(fm, function(x) round(logLik(x), digits = 0)),
  Df = sapply(fm, function(x) attr(logLik(x), "df")))

       ML-Pois Quasi-Pois     NB Hurdle-NB   ZINB
logLik  -18291         NA -12227    -12153 -12155
Df           7          8      8        13     13

The ML Poisson model is clearly inferior to all other fits. The quasi-Poisson model and the sandwich-adjusted Poisson model are not associated with a fitted likelihood. The negative binomial already improves the fit dramatically but can in turn be improved by the hurdle and zero-inflated models which give almost identical fits. This also reflects that the over-dispersion in the data is captured better by the negative-binomial-based models than the plain Poisson model. Additionally, it is of interest how the zero counts are captured by the various models. Therefore, the observed zero counts are compared to the expected number of zero counts for the likelihood-based models:

round(c("Obs" = sum(dt$visits < 1),
  "ML-Pois" = sum(dpois(0, fitted(fm_pois))),
  "NB" = sum(dnbinom(0, mu = fitted(fm_nbin), size = fm_nbin$theta)),
  "NB-Hurdle" = sum(predict(fm_hurdle, type = "density", at = 0)),
  "ZINB" = sum(predict(fm_zinb, type = "density", at = 0))))

      Obs   ML-Pois        NB NB-Hurdle      ZINB 
      683        46       618       683       712 

Thus, the ML Poisson model is again not appropriate whereas the negative-binomial-based models are much better in modeling the zero counts. By construction, the expected number of zero counts in the hurdle model matches the observed number.

In summary, the hurdle and zero-inflation models lead to the best results (in terms of likelihood) on this data set. Above, their mean function for the count component was already shown to be very similar, below we take a look at the fitted zero components:

t(sapply(fm[4:5], function(x) round(x$coefficients$zero, digits = 3)))

          (Intercept) chronic insuranceyes school gendermale
Hurdle-NB       0.054   0.577        0.742  0.056     -0.406
ZINB           -0.098  -1.302       -1.193 -0.087      0.583

This shows that the absolute values are rather different—which is not surprising as they pertain to slightly different ways of modeling zero counts—but the signs of the coefficients match, i.e., are just inversed. For the hurdle model, the zero hurdle component describes the probability of observing a positive count whereas, for the ZINB model, the zero-inflation component predicts the probability of observing a zero count from the point mass component. Overall, both models lead to the same qualitative results and very similar model fits. Perhaps the hurdle model is slightly preferable because it has the nicer interpretation: there is one process that controls whether a patient sees a physician or not, and a second process that determines how many office visits are made.

4 Technical details for hurdle models

The fitting of hurdle models via ML in hurdle() is controlled by the arguments in the hurdle.control() wrapper function:

hurdle.control(method = "BFGS", maxit = 10000, trace = FALSE,
  separate = TRUE, start = NULL, ...)

This modifies some default arguments passed on to the optimizer optim(), such as method, maxit and trace. The latter is also used within hurdle() and can be set to produce more verbose output concerning the fitting process. The argument separate controls whether the two components of the model are optimized separately (the default) or not. This is possible because there are no mixed sources for the zeros in the data (unlike in zero-inflation models). The argument start controls the choice of starting values for calling optim(), all remaining arguments passed through ... are directly passed on to optim().

By default, starting values are estimated by calling glm.fit() for both components of the model separately, once for the counts and once for zero vs. non-zero counts. If starting values are supplied, start needs to be set to a named list with the parameters for the $count and $zero part of the model (and potentially a $theta dispersion parameter if a negative binomial distribution is used).

The fitted model object of class hurdle is similar to glm objects and contains sufficient information on all aspects of the fitting process. In particular, the estimated parameters and associated covariances are included as well as the result from the optim() call. Furthermore, the call, formula, terms structure etc. is contained, potentially also the model frame, dependent variable and regressor matrices.

Following glm.nb(), the \(\theta\) parameter of the negative binomial distribution is treated as a nuisance parameter. Thus, the $coefficients component of the fitted model object just contains estimates of \(\beta\) and \(\gamma\) while the estimate of \(\theta\) and its standard deviation (on a log scale) are kept in extra list elements $theta and $SE.logtheta.

5 Technical details for zero-inflated models

Both the interface of the zeroinfl() function as well as its fitted model objects are virtually identical to the corresponding hurdle functionality. Hence, we only provide some additional information for those aspects that differ from those discussed above. The details of the ML optimization are again provided by a zeroinfl.control() wrapper:

zeroinfl.control(method = "BFGS", maxit = 10000, trace = FALSE,
  EM = FALSE, start = NULL, ...)

The only new argument here is the argument EM which allows for EM estimation of the starting values. Instead of calling glm.fit() only once for both components of the model, this process can be iterated until convergence of the parameters to the ML estimates. The optimizer is still called subsequently (for a single iteration) to obtain the Hessian matrix from which the estimated covariance matrix can be computed.

6 Methods for fitted zero-inflated and hurdle models

Users typically should not need to compute on the internal structure of hurdle or zeroinfl objects because a set of standard extractor functions is provided, an overview is given in Table 2. This includes methods to the generic functions print() and summary() which print the estimated coefficients along with further information. The summary() in particular supplies partial Wald tests based on the coefficients and the covariance matrix. As usual, the summary() method returns an object of class summary.hurdle or summary.zeroinfl, respectively, containing the relevant summary statistics which can subsequently be printed using the associated print() method.

The methods for coef() and vcov() by default return a single vector of coefficients and their associated covariance matrix, respectively, i.e., all coefficients are concatenated. By setting their model argument, the estimates for a single component can be extracted. Concatenating the parameters by default and providing a matching covariance matrix estimate (that does not contain the covariances of further nuisance parameters) facilitates the application of generic inference functions such as coeftest(), waldtest(), and linearHypothesis(). All of these compute Wald tests for which coefficient estimates and associated covariances is essentially all information required and can therefore be queried in an object-oriented way with the coef() and vcov() methods.

Similarly, the terms() and model.matrix() extractors can be used to extract the relevant information for either component of the model. A logLik() method is provided, hence AIC() can be called to compute information criteria and lrtest() for conducting LR tests of nested models.

The predict() method computes predicted means (default) or probabilities (i.e., likelihood contributions) for observed or new data. Additionally, the means from the count and zero component, respectively, can be predicted. For the count component, this is the predicted count mean (without hurdle/inflation): \(\exp(x_i^\top \beta)\). For the zero component, this is the the ratio of probabilities \((1 - f_\mathrm{zero}(0; z_i, \gamma))/(1 - f_\mathrm{count}(0; x_i, \beta))\) of observing non-zero counts in hurdle models. In zero-inflation models, it is the probability \(f_\mathrm{zero}(0; z_i, \gamma)\) of observing a zero from the point mass component in zero-inflated models

Predicted means for the observed data can also be obtained by the fitted() method. Deviations between observed counts \(y_i\) and predicted means \(\hat \mu_i\) can be obtained by the residuals() method returning either raw residuals \(y_i - \hat \mu_i\) or the Pearson residuals (raw residuals standardized by square root of the variance function) with the latter being the default.

Table 2: Functions and methods for objects of class zeroinfl and hurdle. The first ten rows refer to methods, the remaining rows contain generic functions whose default methods work because of the information supplied by the methods above.

Function	Description
print()	simple printed display with coefficient estimates
summary()	standard regression output (coefficient estimates, standard errors, partial Wald tests); returns an object of class “summary.class” containing the relevant summary statistics (which has a print() method)
coef()	extract coefficients of model (full or components), a single vector of all coefficients by default
vcov()	associated covariance matrix (with matching names)
predict()	predictions (means or probabilities) for new data
fitted()	fitted means for observed data
residuals()	extract residuals (response or Pearson)
terms()	extract terms of model components
model.matrix()	extract model matrix of model components
logLik()	extract fitted log-likelihood
coeftest()	partial Wald tests of coefficients
waldtest()	Wald tests of nested models
linearHypothesis()	Wald tests of linear hypotheses
lrtest()	likelihood ratio tests of nested models
AIC()	compute information criteria (AIC, BIC, …)

Replication of textbook results

Cameron and Trivedi (1998, 204) use a somewhat extended version of the model employed above. Because not all variables in that extended model are significant, a reduced set of variables was used throughout the main text. Here, however, we use the full model to show that the tools in countreg reproduce the results of Cameron and Trivedi (1998).

After omitting the responses other than visits and setting "other" as the reference category for region using

dt2 <- NMES1988[, -(2:6)]
dt2$region <- relevel(dt2$region, "other")

we fit a model that contains all explanatory variables, both in the count model and the zero hurdle model:

fm_hurdle2 <- hurdle(visits ~ ., data = dt2, dist = "negbin")

The resulting coefficient estimates are virtually identical to those published in Cameron and Trivedi (1998, 204). The associated Wald statistics are also very similar provided that sandwich standard errors are used (which is not stated explicitely in Cameron and Trivedi 1998).

cfz <- coef(fm_hurdle2, model = "zero")
cfc <- coef(fm_hurdle2, model = "count")
se <- sqrt(diag(sandwich(fm_hurdle2)))
round(cbind(zero = cfz, zero_t = cfz/se[-seq(along = cfc)], 
  count = cfc, count_t = cfc/se[seq(along = cfc)]),
  digits = 3)[c(3, 2, 4, 5, 7, 6, 8, 9:17, 1),]

                  zero zero_t  count count_t
healthexcellent -0.329 -2.310 -0.378  -4.312
healthpoor       0.071  0.420  0.333   5.863
chronic          0.557 10.547  0.143  10.520
adllimited      -0.188 -1.448  0.129   2.504
regionmidwest    0.101  0.880 -0.016  -0.344
regionnortheast  0.129  1.033  0.104   1.974
regionwest       0.202  1.509  0.123   2.444
age              0.190  2.348 -0.075  -2.339
afamyes         -0.327 -2.450  0.002   0.023
gendermale      -0.464 -4.715  0.004   0.098
marriedyes       0.247  2.379 -0.092  -2.110
school           0.054  4.109  0.022   3.824
income           0.007  0.365 -0.002  -0.380
employedyes     -0.012 -0.085  0.030   0.401
insuranceyes     0.762  6.501  0.227   4.007
medicaidyes      0.554  3.055  0.185   2.777
(Intercept)     -1.475 -2.283  1.631   6.017

logLik(fm_hurdle2)

'log Lik.' -12110.49 (df=35)

1/fm_hurdle2$theta

    count 
0.7437966 

There are some small and very few larger deviations in the Wald statistics which are probably explicable by different approximations to the gradient of \(\theta\) (or \(1/\theta\) or \(\log(\theta)\)) and the usage of different non-linear optimizers (and at least ten years of software development).

More replication exercises are performed in the example sections of AER (Zeileis and Kleiber 2008), the software package accompanying Kleiber and Zeileis (2008).

References

Cameron, A. Colin, and Pravin K. Trivedi. 1998. Regression Analysis of Count Data. Cambridge: Cambridge University Press.

———. 2005. Microeconometrics: Methods and Applications. Cambridge: Cambridge University Press.

Deb, Partha, and Pravin K. Trivedi. 1997. “Demand for Medical Care by the Elderly: A Finite Mixture Approach.” Journal of Applied Econometrics 12: 313–36.

Kleiber, Christian, and Achim Zeileis. 2008. Applied Econometrics with R. New York: Springer-Verlag.

Lambert, Diane. 1992. “Zero-Inflated Poisson Regression, with an Application to Defects in Manufacturing.” Technometrics 34: 1–14.

Mullahy, J. 1986. “Specification and Testing of Some Modified Count Data Models.” Journal of Econometrics 33: 341–65.

Zeileis, Achim. 2006. “Object-Oriented Computation of Sandwich Estimators.” Journal of Statistical Software 16 (9): 1–16. http://www.jstatsoft.org/v16/i09/.

Zeileis, Achim, and Christian Kleiber. 2008. AER: Applied Econometrics with R. http://CRAN.R-project.org/package=AER.

———. 2013. countreg: Count Data Regression. http://R-Forge.R-project.org/projects/countreg/.

Footnotes

1. Note that binomial logit and censored geometric models as the hurdle part both lead to the same likelihood function and thus to the same coefficient estimates (Mullahy 1986).