6.11 The Zero inflated model
The purpose of this model (Mullahy 1986) is the same as the previous model in handling large number of zero’s in the data. The difference is that this model combines a mixtures of two types of zero’s, the ones that belong to the hurdle zero model, and those that belong to the standard count model. This model also can be used with Poisson or negative binomial model. It is called inflated because, first we compute the probability of obtaining zero count using the binomial distribution \(\pi=f_{zero}(y=0|Z,\gamma)\), and then add it to the model as follows:
\[\begin{equation} f_{ZI}(y|X,Z,\beta, \gamma)=\pi I(y=0)+(1-\pi) f_{count}(y|X, \beta) \tag{6.81} \end{equation}\]
For instance, if we use the Poisson model, the regression equation of \(\lambda_i\) will be:
\[\begin{equation} \lambda_i=(1-\pi_i)exp(x_i\beta) \tag{6.82} \end{equation}\]
Example 6.5 The data that we use to illustrate the count data models is DoctorVisits data from the R package AER
In R:
suppressPackageStartupMessages(library(AER))
data("DoctorVisits")
counts_R <- DoctorVisits
glimpse(counts_R)[out] Rows: 5,190
[out] Columns: 12
[out] $ visits <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1,…
[out] $ gender <fct> female, female, male, male, male, female, female, female, female, male, male, ma…
[out] $ age <dbl> 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.…
[out] $ income <dbl> 0.55, 0.45, 0.90, 0.15, 0.45, 0.35, 0.55, 0.15, 0.65, 0.15, 0.45, 0.25, 0.55, 0.…
[out] $ illness <dbl> 1, 1, 3, 1, 2, 5, 4, 3, 2, 1, 1, 2, 3, 4, 3, 2, 1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 2,…
[out] $ reduced <dbl> 4, 2, 0, 0, 5, 1, 0, 0, 0, 0, 0, 0, 13, 7, 1, 0, 0, 1, 0, 0, 3, 0, 2, 14, 14, 2,…
[out] $ health <dbl> 1, 1, 0, 0, 1, 9, 2, 6, 5, 0, 0, 2, 1, 6, 0, 7, 5, 0, 0, 0, 0, 1, 0, 2, 11, 6, 6…
[out] $ private <fct> yes, yes, no, no, no, no, no, no, yes, yes, no, no, no, no, yes, no, yes, no, ye…
[out] $ freepoor <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, …
[out] $ freerepat <fct> no, no, no, no, no, no, no, no, no, no, no, yes, no, no, no, no, no, no, no, no,…
[out] $ nchronic <fct> no, no, no, no, yes, yes, no, no, no, no, no, no, yes, yes, yes, no, no, yes, no…
[out] $ lchronic <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, …This data come from Australian health survey where visits is the number of doctor visits in past two weeks with \(11\) features listed above. For simplicity, we restrict our analysis to only three regressors: gender, age, health.
Let us first get the summary of all the data.
[out] visits gender age health
[out] Min. :0.0000 male :2488 Min. :0.1900 Min. : 0.000
[out] 1st Qu.:0.0000 female:2702 1st Qu.:0.2200 1st Qu.: 0.000
[out] Median :0.0000 Median :0.3200 Median : 0.000
[out] Mean :0.3017 Mean :0.4064 Mean : 1.218
[out] 3rd Qu.:0.0000 3rd Qu.:0.6200 3rd Qu.: 2.000
[out] Max. :9.0000 Max. :0.7200 Max. :12.000As we see we do not have missing values and the visits values ranges from 0 to 9 but it should be of integer type rather than double.
[out]
[out] 0 1 2 3 4 5 6 7 8 9
[out] 4141 782 174 30 24 9 12 12 5 1We intentionally used this data that has excess of zero to justify the use of the hurdle and the zero-inflated models discussed earlier.
The Poisson model:
This model belongs to the generalized linear model families, so in the function glm we set the argument family to poisson. In practice this model is sufficient with a wide range of count data.
[out]
[out] Call:
[out] glm(formula = visits ~ ., family = "poisson", data = counts_R)
[out]
[out] Coefficients:
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -2.182264 0.066961 -32.590 < 2e-16 ***
[out] genderfemale 0.208666 0.054611 3.821 0.000133 ***
[out] age 1.377451 0.126268 10.909 < 2e-16 ***
[out] health 0.152990 0.008051 19.002 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] (Dispersion parameter for poisson family taken to be 1)
[out]
[out] Null deviance: 5634.8 on 5189 degrees of freedom
[out] Residual deviance: 5160.9 on 5186 degrees of freedom
[out] AIC: 7500.4
[out]
[out] Number of Fisher Scoring iterations: 6For the interpretation of the coefficient estimates, we should exponentiate these values to get the marginal effect since the poisson model uses the log link function to preclude negative values. For continuous predictor, say age, if this predictor increases by one year, ceteris-paribus, we expect the doctor visits will be $exp(1.377)=$3.9629948 times larger. whereas, for categorical predictor, say gender, the female has $exp(0.208)=$1.2312132 larger doctor visits than male.
To check the overdispersion, we use the function dispersiontest from the AER package. The argument trafo defines the dispersion relationships such that:
- \(trafo=0\): \(Var(y)=E(y)\) The default.
- \(trafo=1\): \(Var(y)=E(y)+\alpha E(y)\).
- \(trafo=2\): \(Var(y)=E(y)+\alpha E(y)^2\).
[out]
[out] Overdispersion test
[out]
[out] data: poisson_R
[out] z = 7.0729, p-value = 7.585e-13
[out] alternative hypothesis: true dispersion is greater than 1
[out] sample estimates:
[out] dispersion
[out] 1.917024[out]
[out] Overdispersion test
[out]
[out] data: poisson_R
[out] z = 7.0729, p-value = 7.585e-13
[out] alternative hypothesis: true alpha is greater than 0
[out] sample estimates:
[out] alpha
[out] 0.9170243As we see, the test detects overdispersion problem so that we should specify another model.
The quasi Poisson model:
We have to change only the argument family as follows:
[out]
[out] Call:
[out] glm(formula = visits ~ ., family = "quasipoisson", data = counts_R)
[out]
[out] Coefficients:
[out] Estimate Std. Error t value Pr(>|t|)
[out] (Intercept) -2.18226 0.09253 -23.584 < 2e-16 ***
[out] genderfemale 0.20867 0.07547 2.765 0.00571 **
[out] age 1.37745 0.17449 7.894 3.54e-15 ***
[out] health 0.15299 0.01113 13.751 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] (Dispersion parameter for quasipoisson family taken to be 1.909559)
[out]
[out] Null deviance: 5634.8 on 5189 degrees of freedom
[out] Residual deviance: 5160.9 on 5186 degrees of freedom
[out] AIC: NA
[out]
[out] Number of Fisher Scoring iterations: 6This model uses the quasi maximum likelihood that gives identical coefficient estimates but with different (corrected) standard errors. Since all the variables are similarly significant , we see that the models are the same except the correction of the standard errors that are now more larger. In other words, the poisson distribution under overdispersion underestimates the standard errors, and hence the t test would be biased towards the rejection of the null hypothesis. To better understand what is going on with quasi-poisson model, we put the estimates and the standard errors of both models into one table, then we add a column that resulted from dividing the second standard errors vector by the first one.
table1 <- tidy(poisson_R)
colnames(table1) <- NULL
table2 <- tidy(quasi_poisson_R)
colnames(table2) <- NULL
tibble(term=table1[[1]], estimate1=table1[[2]], std1=table1[[3]],
estimate2=table2[[2]], std2=table2[[3]], dispersion= std2/std1)[out] # A tibble: 4 × 6
[out] term estimate1 std1 estimate2 std2 dispersion
[out] <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
[out] 1 (Intercept) -2.18 0.0670 -2.18 0.0925 1.38
[out] 2 genderfemale 0.209 0.0546 0.209 0.0755 1.38
[out] 3 age 1.38 0.126 1.38 0.174 1.38
[out] 4 health 0.153 0.00805 0.153 0.0111 1.38The first two columns contain the coefficients of both models. Not surprisingly that the result of the last column is constant since this is exactly what the quasi maximum likelihood does, it computes the corrected standard errors from the original ones by \(std2=dispersion\times std1\) , with the dispersion value being estimated as \(1.381868\). if you want to know where this value came from, the answer is simple. this model computes the sigma of the standardized residuals resulted from the original model. we can thus get this value by specifying the argument type to pear then computing sigma by hand as follows:
[out] [1] 1.381867The negative binomial model:
[out]
[out] Call:
[out] glm.nb(formula = visits ~ ., data = counts_R, init.theta = 0.5146961771,
[out] link = log)
[out]
[out] Coefficients:
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -2.20101 0.08327 -26.433 < 2e-16 ***
[out] genderfemale 0.26650 0.06987 3.814 0.000137 ***
[out] age 1.30315 0.16367 7.962 1.69e-15 ***
[out] health 0.16368 0.01259 13.001 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] (Dispersion parameter for Negative Binomial(0.5147) family taken to be 1)
[out]
[out] Null deviance: 3313.4 on 5189 degrees of freedom
[out] Residual deviance: 3030.7 on 5186 degrees of freedom
[out] AIC: 6917.2
[out]
[out] Number of Fisher Scoring iterations: 1
[out]
[out]
[out] Theta: 0.5147
[out] Std. Err.: 0.0401
[out]
[out] 2 x log-likelihood: -6907.1890The hurdle model with poisson distribution:
To perform this model we use the function hurdle from the package pscl.
suppressPackageStartupMessages(library(pscl))
hurdle_poisson_R <- hurdle(visits~., data=counts_R, dist = "poisson")
summary(hurdle_poisson_R)[out]
[out] Call:
[out] hurdle(formula = visits ~ ., data = counts_R, dist = "poisson")
[out]
[out] Pearson residuals:
[out] Min 1Q Median 3Q Max
[out] -1.1176 -0.4569 -0.3591 -0.3040 15.3653
[out]
[out] Count model coefficients (truncated poisson with log link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -0.43707 0.11542 -3.787 0.000153 ***
[out] genderfemale -0.07514 0.08684 -0.865 0.386852
[out] age 0.26540 0.20287 1.308 0.190788
[out] health 0.08512 0.01182 7.203 5.9e-13 ***
[out] Zero hurdle model coefficients (binomial with logit link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -2.52720 0.09191 -27.498 < 2e-16 ***
[out] genderfemale 0.33242 0.07563 4.396 1.1e-05 ***
[out] age 1.65830 0.17705 9.366 < 2e-16 ***
[out] health 0.17730 0.01453 12.201 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Number of iterations in BFGS optimization: 17
[out] Log-likelihood: -3544 on 8 DfAs we see, this output has two tables. The above is for the Poisson model performed only on the truncated positive values. The other one below is the result of the logistic regression with only two classes (zero or positive value). For instance, the variable age is not significant in the first model, so it is not a good predictor for the count variable. But it could be a good one to predict whether the predicted value is zero or positive since it is highly significant (from the second model). That is why we said, in our earlier discussion about this model, that we could use different regressors for each part. So in our case, we can remove the variable age from the Poisson model and leave it in the zero hurdle model.
hurdle model with negative binomial distribution:
[out]
[out] Call:
[out] hurdle(formula = visits ~ ., data = counts_R, dist = "negbin")
[out]
[out] Pearson residuals:
[out] Min 1Q Median 3Q Max
[out] -0.7795 -0.4317 -0.3424 -0.2894 14.6020
[out]
[out] Count model coefficients (truncated negbin with log link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -11.69311 52.65821 -0.222 0.824
[out] genderfemale -0.06821 0.14388 -0.474 0.635
[out] age 0.14608 0.33676 0.434 0.664
[out] health 0.12518 0.02553 4.904 9.38e-07 ***
[out] Log(theta) -11.48561 52.65855 -0.218 0.827
[out] Zero hurdle model coefficients (binomial with logit link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -2.52720 0.09191 -27.498 < 2e-16 ***
[out] genderfemale 0.33242 0.07563 4.396 1.1e-05 ***
[out] age 1.65830 0.17705 9.366 < 2e-16 ***
[out] health 0.17730 0.01453 12.201 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Theta: count = 0
[out] Number of iterations in BFGS optimization: 85
[out] Log-likelihood: -3410 on 9 DfUsing the negative binomial model, all the regressors are not significant except for the variable health in the Poisson model. But all of which are significant in the zero hurdle model.
The zero inflated model with poisson distribution:
We use the function zeroinfl from the same package pscl
[out]
[out] Call:
[out] zeroinfl(formula = visits ~ ., data = counts_R, dist = "poisson")
[out]
[out] Pearson residuals:
[out] Min 1Q Median 3Q Max
[out] -1.0886 -0.4650 -0.3626 -0.2975 15.3397
[out]
[out] Count model coefficients (poisson with log link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -0.42148 0.11415 -3.692 0.000222 ***
[out] genderfemale -0.10284 0.08529 -1.206 0.227884
[out] age 0.26123 0.19684 1.327 0.184481
[out] health 0.08535 0.01202 7.102 1.23e-12 ***
[out]
[out] Zero-inflation model coefficients (binomial with logit link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) 1.8107 0.1593 11.365 < 2e-16 ***
[out] genderfemale -0.5566 0.1306 -4.260 2.04e-05 ***
[out] age -1.8813 0.3173 -5.929 3.05e-09 ***
[out] health -0.1411 0.0250 -5.645 1.65e-08 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Number of iterations in BFGS optimization: 17
[out] Log-likelihood: -3543 on 8 DfThe result here is similar to the hurdle poisson model.
The zero inflated model with negative binomial distribution:
[out]
[out] Call:
[out] zeroinfl(formula = visits ~ ., data = counts_R, dist = "negbin")
[out]
[out] Pearson residuals:
[out] Min 1Q Median 3Q Max
[out] -0.6729 -0.4547 -0.3885 -0.2690 14.5419
[out]
[out] Count model coefficients (negbin with log link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) -1.53588 0.18618 -8.249 < 2e-16 ***
[out] genderfemale -0.10357 0.10811 -0.958 0.338052
[out] age 0.88147 0.23152 3.807 0.000141 ***
[out] health 0.13085 0.01657 7.897 2.86e-15 ***
[out] Log(theta) -0.44462 0.11817 -3.763 0.000168 ***
[out]
[out] Zero-inflation model coefficients (binomial with logit link):
[out] Estimate Std. Error z value Pr(>|z|)
[out] (Intercept) 1.0308 0.3565 2.891 0.00384 **
[out] genderfemale -2.6596 1.6623 -1.600 0.10961
[out] age -3.4206 0.9967 -3.432 0.00060 ***
[out] health -0.3606 0.2230 -1.617 0.10585
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Theta = 0.6411
[out] Number of iterations in BFGS optimization: 31
[out] Log-likelihood: -3436 on 9 DfCompare all the models using RMSE:
we will compute the root mean squared errors for all the above models using the ModelMetrics package, then choose the best one based on this metric.
suppressPackageStartupMessages(library(ModelMetrics))
# poisson model
pred_poisson_R<- predict(poisson_R,counts_R, type = "response")
rmse_poisson_R<-rmse(counts_R$visits,round(pred_poisson_R))
#Quasi poisson model
pred_quasi_poisson_R<- predict(quasi_poisson_R,counts_R, type = "response")
rmse_quasi_poisson_R<-rmse(counts_R$visits,round(pred_quasi_poisson_R))
# negative binomial model
pred_nbinom_R<- predict(nbinom_R,counts_R, type = "response")
rmse_nbinom_R<-rmse(counts_R$visits,round(pred_nbinom_R))
# hurdle poisson model
pred_hurdle_poisson_R<- predict(hurdle_poisson_R,counts_R, type = "response")
rmse_hurdle_poisson_R<-rmse(counts_R$visits,round(pred_hurdle_poisson_R))
# hurdle negative binomial model
pred_hurdle_nbinom_R<- predict(hurdle_nbinom_R,counts_R, type = "response")
rmse_hurdle_nbinom_R<-rmse(counts_R$visits,round(pred_hurdle_nbinom_R))
# zero poisson model
pred_zero_poisson_R<- predict(zero_poisson_R, counts_R, type = "response")
rmse_zero_poisson_R<-rmse(counts_R$visits,round(pred_zero_poisson_R))
# zero negative binomial model
pred_zero_nbinom_R<- predict(zero_nbinom_R, counts_R, type = "response")
rmse_zero_nbinom_R<-rmse(counts_R$visits,round(pred_zero_nbinom_R))
# put everything in a data frame
rmses <- c(rmse_poisson_R, rmse_quasi_poisson_R, rmse_nbinom_R,
rmse_hurdle_poisson_R, rmse_hurdle_nbinom_R, rmse_zero_poisson_R,
rmse_zero_nbinom_R)
models <- c("poisson_R", "quasi_poisson_R", "nbinom_R",
"hurdle_poisson_R", "hurdle_nbinom_R", "zero_poisson_R",
"zero_nbinom_R")
tibble(models,rmses)%>%
arrange(rmses)[out] # A tibble: 7 × 2
[out] models rmses
[out] <chr> <dbl>
[out] 1 zero_nbinom_R 0.817
[out] 2 zero_poisson_R 0.820
[out] 3 hurdle_poisson_R 0.822
[out] 4 hurdle_nbinom_R 0.823
[out] 5 poisson_R 0.825
[out] 6 quasi_poisson_R 0.825
[out] 7 nbinom_R 0.827As we see, the zero inflated model is the best one. However, we should be careful when using the whole data to evaluate models. in machine learning methodology we keep some data aside as a testing data to get a reliable evaluation.
We first move the data to python
In Python we will use the statsmodels package.
The Poisson model:
import statsmodels.formula.api as smf
poisson_py = smf.poisson(f"visits~{'+'.join(counts_py.columns[1:])}", counts_py).fit()[out] Optimization terminated successfully.
[out] Current function value: 0.721815
[out] Iterations 6[out] <class 'statsmodels.iolib.summary.Summary'>
[out] """
[out] Poisson Regression Results
[out] ==============================================================================
[out] Dep. Variable: visits No. Observations: 5190
[out] Model: Poisson Df Residuals: 5186
[out] Method: MLE Df Model: 3
[out] Date: Thu, 17 Sept 2026 Pseudo R-squ.: 0.05949
[out] Time: 11:14:28 Log-Likelihood: -3746.2
[out] converged: True LL-Null: -3983.2
[out] Covariance Type: nonrobust LLR p-value: 2.103e-102
[out] ====================================================================================
[out] coef std err z P>|z| [0.025 0.975]
[out] ------------------------------------------------------------------------------------
[out] Intercept -2.1823 0.067 -32.590 0.000 -2.314 -2.051
[out] gender[T.female] 0.2087 0.055 3.821 0.000 0.102 0.316
[out] age 1.3775 0.126 10.909 0.000 1.130 1.625
[out] health 0.1530 0.008 19.002 0.000 0.137 0.169
[out] ====================================================================================
[out] """As we see, The results are the same.
The negative binomial model:
import statsmodels.formula.api as smf
nbinom_py = smf.negativebinomial(f"visits~{'+'.join(counts_py.columns[1:])}", counts_py).fit()[out] Optimization terminated successfully.
[out] Current function value: 0.665432
[out] Iterations: 23
[out] Function evaluations: 25
[out] Gradient evaluations: 25[out] <class 'statsmodels.iolib.summary.Summary'>
[out] """
[out] NegativeBinomial Regression Results
[out] ==============================================================================
[out] Dep. Variable: visits No. Observations: 5190
[out] Model: NegativeBinomial Df Residuals: 5186
[out] Method: MLE Df Model: 3
[out] Date: Thu, 17 Sept 2026 Pseudo R-squ.: 0.03692
[out] Time: 11:14:29 Log-Likelihood: -3453.6
[out] converged: True LL-Null: -3586.0
[out] Covariance Type: nonrobust LLR p-value: 4.128e-57
[out] ====================================================================================
[out] coef std err z P>|z| [0.025 0.975]
[out] ------------------------------------------------------------------------------------
[out] Intercept -2.2009 0.084 -26.282 0.000 -2.365 -2.037
[out] gender[T.female] 0.2665 0.069 3.866 0.000 0.131 0.402
[out] age 1.3028 0.162 8.049 0.000 0.986 1.620
[out] health 0.1637 0.013 12.346 0.000 0.138 0.190
[out] alpha 1.9429 0.152 12.814 0.000 1.646 2.240
[out] ====================================================================================
[out] """The zero inflated model with poisson distribution:
The statsmodels.formula.api does not have an attribute for the zero inflated model, so we should call from the statsmodels.api:
import statsmodels.formula.api as smf
import statsmodels.api as sm
X = counts_py.iloc[:,1:]
X = sm.add_constant(X)
# convert categorical variables to dammies
# pandas >= 2.0 returns boolean dummies by default; statsmodels needs numeric.
X = pd.get_dummies(X, drop_first=True, dtype=float)
Y = counts_py.iloc[:, 0]
zero_poisson_py = sm.ZeroInflatedPoisson(Y, X).fit()[out] Optimization terminated successfully.
[out] Current function value: 0.692347
[out] Iterations: 26
[out] Function evaluations: 28
[out] Gradient evaluations: 28[out] <class 'statsmodels.iolib.summary.Summary'>
[out] """
[out] ZeroInflatedPoisson Regression Results
[out] ===============================================================================
[out] Dep. Variable: visits No. Observations: 5190
[out] Model: ZeroInflatedPoisson Df Residuals: 5186
[out] Method: MLE Df Model: 3
[out] Date: Thu, 17 Sept 2026 Pseudo R-squ.: 0.03256
[out] Time: 11:14:29 Log-Likelihood: -3593.3
[out] converged: True LL-Null: -3714.2
[out] Covariance Type: nonrobust LLR p-value: 3.686e-52
[out] =================================================================================
[out] coef std err z P>|z| [0.025 0.975]
[out] ---------------------------------------------------------------------------------
[out] inflate_const 0.1744 0.050 3.482 0.000 0.076 0.273
[out] const -1.2580 0.074 -16.928 0.000 -1.404 -1.112
[out] age 1.1328 0.141 8.030 0.000 0.856 1.409
[out] health 0.1253 0.009 14.176 0.000 0.108 0.143
[out] gender_female 0.2066 0.061 3.392 0.001 0.087 0.326
[out] =================================================================================
[out] """