Statistics 244 – Poisson log-linear models, and related issues

Count models

• Log-linear models for Poisson counts

• Models with offsets

• Zero-inflated and hurdle Poisson models

Poisson regression

For observation i, let xi1 , . . . , xip be p predictor variables, and assume a response variable Yi where

Yi ∼ Po(µi ).

The probability mass function (pmf) is given by

e−µi µyi i
pi (yi ) =
yi !

for yi = 0, 1, 2, . . ..

Note that µi = E(Yi ). Also, Var(Yi ) = µi .

Example: Species richness of ants

Study was conducted (Gottelli and Ellison, 2002) that counted the number of ant species in 44
locations in CT, MA and VT. Sampling occurred in 64 square meter grids.

In addition to the number of species per location (the response), the following predictor informa-
tion was recorded:

Habitat: Whether the location was a forest or a bog

Latitude: The latitude (in decimal degrees) for the site
Elevation: Elevation of the site in meters above sea level

Want to model the mean number of species per location as a function of these three predictors.

Numerical summaries of ant data

Srich Habitat Latitude Elevation

Min. : 2.000 Bog :22 Min. :41.97 Min. : 1.0
1st Qu.: 4.000 Forest:22 1st Qu.:42.17 1st Qu.: 95.0

1
Median : 6.000 Median :42.56 Median :223.0
Mean : 7.023 Mean :43.02 Mean :232.7
3rd Qu.: 8.250 3rd Qu.:44.29 3rd Qu.:353.0
Max. :18.000 Max. :44.95 Max. :543.0

2
3
4
5
As we have seen previously, a Poisson pmf is a member of the EDF, where
 
yi θi − b(θi )
pi (yi ) = exp + c(yi , ϕ)
ϕ
and where θi = log µi , and

ϕ = 1
b(θi ) = exp(θi )
c(yi , ϕ) = − log y!

Can show V (µ) = µ.

Because µi > 0, ideally want to choose a link function that maps (0, ∞) to (−∞, ∞).

Conventional choice of link function: log link

Choosing η = g(µ) = log µ is the standard choice for Poisson models because

• The log link is the canonical link function.

• The log link maps (0, ∞) to (−∞, ∞).
• The log function is monotone increasing and differentiable, and has a simple inverse func-
tion (µ = exp(η)).
• Difficult to come up with a simple alternative.

Using a log link is so common, Poisson GLMs are often referred to as “log-linear models.”

Poisson log-linear model Full log-linear model specification:

Yi ∼ Po(µi )

with
ηi = log µi = xi β.

Statistical inference

Most of the same issues carry over from binary response models to Poisson models.

Maximum likelihood estimation

Use Fisher scoring to estimate β.

Estimated covariance matrix for β̂:

var(β̂) = (X T W X)−1

where W is the diagonal matrix with entries

1
wii = = µ̂i = exp(xi β̂).
(∂ηi /∂µi )2µi =µ̂i ϕV (µ̂i )

6
The standard error of β̂j is given by
q
s.e.(β̂j ) = var(β̂)jj .

Wald tests and confidence intervals

The test statistic as usual is

β̂j
z= .
s.e.(β̂j )
Under Ho , z ∼ N(0, 1) approximately.

Similarly, a 100(1 − α) confidence interval for βj is given by

∗
β̂j ± z1−α/2 s.e.(β̂j )
∗
where z1−α/2 is the (1 − α/2) critical value on N(0, 1).

Better to construct profile likelihood-based confidence intervals.

Log-linear model for ants data

> [Link] = glm(Srich ˜ Habitat + Latitude + Elevation,

family=poisson, data=ants)
> summary([Link])

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 11.9368121 2.6214970 4.553 5.28e-06 ***
HabitatForest 0.6354389 0.1195664 5.315 1.07e-07 ***
Latitude -0.2357930 0.0616638 -3.824 0.000131 ***
Elevation -0.0011411 0.0003749 -3.044 0.002337 **
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 102.763 on 43 degrees of freedom

Residual deviance: 40.690 on 40 degrees of freedom
AIC: 209.04

Coefficient interpretation

Value of βj (or its estimate) is the change in log µ to an increase of 1.0 in xj , holding the other
variables constant.

Example The coefficient estimate −0.2357 of Latitude means that an increase in 1 degree in
latitude corresponds to a drop of 0.2357 of the mean species richness of ants on the log scale.

7
But there is a more natural interpretation on the scale of the mean.

Multiplicative interpretation

The log link implies

µ = exp(β1 x1 + . . . + βp xp )
so that increasing xj by 1.0 corresponds the a multiplicative increase in µ by exp(βj ).

Example The coefficient estimate −0.2357 of Latitude means that an increase in 1 degree in
latitude corresponds to multiplying µ by exp(−0.2357) = 0.79. In other words, a 1 degree latitude
increase corresponds to a 1 − 0.79 = 21% drop in the mean number of ant species.

Multiplicative interpretation for factors

For treatment contrast coded factors, it is straightforward to interpret the coefficient by exponen-
tiating.

Example The coefficient estimate 0.6354 for HabitatForest implies that the mean species rich-
ness of ants in forest locations is exp(0.6354) = 1.8877 times as large as the mean species richness
of ants in bog locations, holding the other variables constant.

Profile likelihood confidence intervals

> confint([Link])

2.5 % 97.5 %
(Intercept) 6.907572250 17.1971550450
HabitatForest 0.403636184 0.8728861025
Latitude -0.359658887 -0.1176287695
Elevation -0.001885055 -0.0004146158

Estimating µ

For log-linear models, once β̂ is obtained, then because

log µ = xβ,

we can make inferences about µ from

µ̂ = exp(xβ̂).
using either the Delta method for standard errors, or obtaining confidence intervals for η = xβ
and then transforming the interval endpoints.

Example Ant species richness, again

New data:

> ants1new = [Link](Habitat=c("Bog","Forest"),

8
 Latitude=c(43,43),Elevation=c(232,232))
> ants1new

Habitat Latitude Elevation

1 Bog 43 232
2 Forest 43 232

Predictions from these two observations, with standard errors:

> predict([Link], newdata=ants1new,

type="response",[Link]=T)
$fit
1 2
4.631969 8.744464

$[Link]
1 2
0.4556655 0.6355438

Can see that the expected count for “Forest” is greater than for “Bog” with the given latitude and
elevation.

Moderately large standard errors reflect small sample size.

Deviance function, and LRT

For a Poisson GLM, the deviance function is

n    
X yi
D(µ̂ | y) = 2 yi log − (yi − µ̂i ) .
µ̂i
i=1

Because ϕ = 1 in Poisson models, the scaled deviance is equal to the ordinary deviance. That is,
D∗ (µ̂ |y) = D(µ̂ |y).

Likelihood ratio test

Want to test model M0 with p0 coefficients nested within model M1 with p1 coefficients.

Letting D∗ (µ̂0 | y) be the deviance for model M0 and letting D∗ (µ̂1 | y) be the deviance for the
model M1 , then compute
D∗ (µ̂0 | y) − D∗ (µ̂1 | y).

Determine p-value based on a χ2p1 −p0 distribution.

As with binary response models, it is usually best to use the LRT to compare models rather than to
evaluate single models (relative to the saturated model).

Ants data, continued

9
Fit models that interacts habitat (bog/forest) with latitude, elevation, and both.

To keep the nested structure of the models, first ignore the interaction between habitat and eleva-
tion:

Analysis of Deviance Table

Model 1: Srich ˜ Habitat + Latitude + Elevation

Model 2: Srich ˜ Habitat * Latitude + Elevation
Model 3: Srich ˜ Habitat * (Latitude + Elevation)
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 40 40.69
2 39 40.65 1 0.04019 0.8411
3 38 40.27 1 0.38049 0.5373

Now ignore the interaction between habitat and latitude:

Analysis of Deviance Table

Model 1: Srich ˜ Habitat + Latitude + Elevation

Model 2: Srich ˜ Habitat * Elevation + Latitude
Model 3: Srich ˜ Habitat * (Latitude + Elevation)
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 40 40.690
2 39 40.273 1 0.41771 0.5181
3 38 40.270 1 0.00297 0.9566

Either way, the data does not provide strong evidence of interactions with habitat.

Diagnostics: Nothing new here

Can compute (as with binary response models)

• Fitted values

• Deviance residuals

• Pearson residuals, and their Jackknifed/studentized variation

• Cooks distances

Examine plots for deviation from model assumptions. No need to bin observations this time.

10
11
12
13
Poisson models with different exposures

• Poisson counts are usually measured over a specified time period, or a specified region, or
some other exposure measure.

• Different observations may be associated with different exposures.

• In such situations, we should model the mean per unit exposure (e.g., per unit time).

Basic idea

Suppose Yi is the number of mosquitoes zapped in a bug zapper on the i-th observation. Assume

Y1 = Number of mosquitoes zapped in a 1-minute period

Y2 = Number of mosquitoes zapped in a 2-minute period

Reasonable to assume Y1 ∼ Po(µ1 ) and Y2 ∼ Po(µ2 ).

Because the exposure for Y1 is half that of Y2 , should model 2µ1 = µ2 .

Let λ be the mean number of mosquitoes per minute that are zapped. Then rewrite model as

Y1 ∼ Po(λ)
Y2 ∼ Po(2λ).

In general, suppose for observation i,

Ti = Known length of time for unit i over which counts are measured
Yi = Number of events observed for unit i in Ti units of time

and let λi be the mean number of counts per unit time for observation i. Then model µi = Ti λi so
that
Yi ∼ Po(Ti λi ).

Example: Ship damage data

Example from McCullagh and Nelder, section 6.3.2.

Number of damage incidents for 34 ships. In addition to damage counts per ship, the following
was recorded:

• Ship type (A, B, C, D and E)

• Year of construction (1960-64, 1964-69, 1970-74, 1975-79)

• Period of operation (1960-74, 1975-79)

• Total months of service

14
Want to examine the effects of the first three factors on ship damage counts per month.

Offsets

With a Poisson GLM, we really want the predictors (which are assumed not to depend on expo-
sure) to model the mean counts per unit time, not the overall mean counts.

• Let λi be the mean rate of counts per unit time.

• Set log λi = xi β to link λi with the predictors for unit i.

• The observations are Poisson with mean µi = Ti λi over Ti units of time.

The linear predictor can be expressed as

log λi = xi β.

Adding log Ti to both sides results in

log Ti + log λi = log Ti λi = log Ti + xi β

which then implies

log µi = log Ti + xi β
Notice that log Ti is not multiplied by an unknown coefficient. This quantity (log Ti ) is called an
“offset” in the log-linear model.

Straightforward to adjust Fisher scoring algorithm to include offsets. (homework)

Ship damage modeling: Additive model

> [Link] = glm(Y ˜ offset(log(Months)) + Type + Constr

+ OpPer, family=poisson, data=ships)
> summary([Link])

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.40590 0.21744 -29.460 < 2e-16 ***
TypeB -0.54334 0.17759 -3.060 0.00222 **
TypeC -0.68740 0.32904 -2.089 0.03670 *
TypeD -0.07596 0.29058 -0.261 0.79377
TypeE 0.32558 0.23588 1.380 0.16750
Constr65 0.69714 0.14964 4.659 3.18e-06 ***
Constr70 0.81843 0.16977 4.821 1.43e-06 ***
Constr75 0.45343 0.23317 1.945 0.05182 .
OpPer75 0.38447 0.11827 3.251 0.00115 **
---
Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

15
Ship damage modeling Interaction between year of construction and period of operation:

Analysis of Deviance Table

Model 1: Y ˜ offset(log(Months)) + Type + Constr + OpPer

Model 2: Y ˜ offset(log(Months)) + Type + Constr * OpPer
Resid. Df Resid. Dev Df Deviance P(>|Chi|)
1 25 38.695
2 23 36.908 2 1.7875 0.4091

Zero-inflated count models

In many random phenomenon when counting occurrences over a fixed exposure, one might ex-
pect an excess of zero counts (compared to Poisson variation):

• Randomly sample college students in the Boston area during flu season. Count the number
of days a student suffered with the flu. (some fraction may have never have contracted the
flu)

• Randomly sample a couple that has been married at least 10 years. Count the number of
children born to the couple. (some fraction may not be able to have children even if they
wanted them, other “structural” reasons)

• Random sample college students. Count the number of days in the last year that the student
has engaged in a physical recreational activity. (some fraction of people never engage in
physical activity)

Example: Boating trips in Texas

Data were obtained on the number of recreational boating trips to Lake Somerville, Texas, in 1980,
based on 659 respondents to a survey of registered leisure boat owners in 23 counties in eastern
Texas.

• trips: Number of recreational boating trips

• quality: Facility’s quality ranking (1 to 5)

• ski: Did the respondent water ski at the lake (no/yes)?

• income: Household income in $1000s.

• userfee: Did the individual pay an annual user fee at the lake (no/yes)?

• costC, costS, costH: Expenses (in $) when visiting Lake Conroe, Somerville, and Hous-
ton, respectively

Data summary

16
> summary(RecreationDemand)

trips quality ski income userfee

Min. : 0.000 Min. :0.000 no :417 Min. :1.000 no :646
1st Qu.: 0.000 1st Qu.:0.000 yes:242 1st Qu.:3.000 yes: 13
Median : 0.000 Median :0.000 Median :3.000
Mean : 2.244 Mean :1.419 Mean :3.853
3rd Qu.: 2.000 3rd Qu.:3.000 3rd Qu.:5.000
Max. :88.000 Max. :5.000 Max. :9.000

costC costS costH

Min. : 4.34 Min. : 4.767 Min. : 5.70
1st Qu.: 28.24 1st Qu.: 33.312 1st Qu.: 28.96
Median : 41.19 Median : 47.000 Median : 42.38
Mean : 55.42 Mean : 59.928 Mean : 55.99
3rd Qu.: 69.67 3rd Qu.: 72.573 3rd Qu.: 68.56
Max. :493.77 Max. :491.547 Max. :491.05

Log-linear model fit

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.264993 0.093722 2.827 0.00469 **

17
quality 0.471726 0.017091 27.602 < 2e-16 ***
skiyes 0.418214 0.057190 7.313 2.62e-13 ***
income -0.111323 0.019588 -5.683 1.32e-08 ***
userfeeyes 0.898165 0.078985 11.371 < 2e-16 ***
costC -0.003430 0.003118 -1.100 0.27131
costS -0.042536 0.001670 -25.467 < 2e-16 ***
costH 0.036134 0.002710 13.335 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 4849.7 on 658 degrees of freedom

Residual deviance: 2305.8 on 651 degrees of freedom

Compare observed and expected counts for the model

To estimate the expected counts:

• We have n = 659 observations.

• An estimated probability of Pr(Y = y) for y = 0, 1, 2, . . . is

n
1 X e−µ̂i µ̂yi
p̂y =
n y!
i=1

where µ̂i = exp(xi β̂).

The expected counts are then np̂y .

In R:

rbind(obs = table(RecreationDemand$trips)[1:10],
exp = round(sapply(0:9, function(x) sum(dpois(x, fitted([Link])))))))

0 1 2 3 4 5 6 7 8 9
obs 417 68 38 34 17 13 11 2 8 1
exp 277 146 68 41 30 23 17 13 10 7

The observed number of zero counts is quite a bit higher than explained by the log-linear model.

If some of the respondents were of the type who would never really go on boat trips, then we
might expect an excess zero count.

Zero-inflated Poisson [ZIP] model (Lambert, 1992)

This model for excess zero counts assumes that the count data is a mixture of two types:

18
 • Observations coming from a process that always produces zero counts
• Observations coming from a process that produces Poisson counts

The mixture probability is assumed unknown, but can be inferred from the model fitting.

More specifically, assume

(
0 with probability 1 − πi
Yi ∼
Po(µi ) with probability πi .

From the law of total probability, this implies


 (1 − πi ) + πi exp(−µi ) for y = 0
Pr(Yi = y) =
exp(−µi )µyi
 πi y! for y = 1, 2, . . ..

Mean and variance of Yi for ZIP model

Trick for latent class models: Assume Zi ∼ Bern(πi ) is an indicator of observation i being in the
Poisson group.

Then
E(Yi ) = E(E(Yi |Zi )) = E(Zi µi ) = µi E(Zi ) = µi πi
Var(Yi ) = Var(E(Yi |Zi )) + E(Var(Yi |Zi )) = Var(Zi µi ) + E(Zi µi )
= µ2i Var(Zi ) + µi πi = µ2i πi (1 − πi ) + µi πi
= µi πi (µi (1 − πi ) + 1)
Notice that the variance is larger than the mean. This means that the ZIP model accounts for
“over-dispersion” relative to the Poisson model.

We can now create two link functions; one to link πi to xi , and a second to link µi to xi . A common
choice (others are possible)
g1 (πi ) = logit πi = xi β 1
g2 (µi ) = log µi = xi β 2
where some of the elements of β 1 and β 2 may be set to 0 by design.

The log-likelihood can be written (up to a constant additive term) as

ℓ(β 1 , β 2 |y) =
X n
X
log[1 + exp(xi β 1 ) exp(− exp(xi β 2 ))] − log(1 + exp(xi β 1 ))
i:yi =0 i=1
X
+ [xi β 1 + yi xi β 2 − exp(xi β 2 )]
i:yi ≥1

Example: Boating trips

19
> library("pscl")
> [Link] = zeroinfl(trips ˜quality+ski+income+
userfee+costC+costS+costH |
quality+income, dist="pois", data=RecreationDemand)
> print(summary([Link]))
Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.099163 0.111397 18.844 < 2e-16 ***
quality 0.033833 0.023914 1.415 0.157
skiyes 0.471691 0.058187 8.106 5.21e-16 ***
income -0.099780 0.020779 -4.802 1.57e-06 ***
userfeeyes 0.610488 0.079435 7.685 1.53e-14 ***
costC 0.002369 0.003818 0.620 0.535
costS -0.037600 0.002038 -18.454 < 2e-16 ***
costH 0.025234 0.003355 7.522 5.40e-14 ***

Boating trips, cont’d

Zero-inflation model coefficients (binomial with logit link):

Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.29191 0.51608 6.379 1.79e-10 ***
quality -1.91407 0.20619 -9.283 < 2e-16 ***
income -0.04502 0.10797 -0.417 0.677
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

Number of iterations in BFGS optimization: 23

Log-likelihood: -1181 on 11 Df

Important For the binomial part of the ZIP model, the pscl library models the probability of the
zero component, not the probability of the Poisson component.

Improvement over log-linear model?

> print(lrtest([Link],[Link]))
Likelihood ratio test

Model 1: trips ˜ quality + ski + income + userfee +

costC + costS + costH
Model 2: trips ˜ quality + ski + income + userfee +
costC + costS + costH | quality + income

#Df LogLik Df Chisq Pr(>Chisq)

1 8 -1529.4
2 11 -1180.8 3 697.27 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

20
Yes-sir!

Actually, need to be careful about this interpretation because the null hypothesis is on the bound-
ary of the parameter space (πi = 1).

Comparison of observed and expected ZIP counts

> rbind(obs = table(RecreationDemand$trips)[1:10],

exp = round(colSums(
predict([Link], type=’prob’)[, 1:10]))))

0 1 2 3 4 5 6 7 8 9
obs 417 68 38 34 17 13 11 2 8 1
exp 414 25 31 32 30 26 22 18 14 11

Much better consistency with observed counts compared to log-linear model.

Comments about ZIP model

• Can use Newton-Raphson algorithm for model fitting. Other popular alternatives are the
EM algorithm (a natural choice when working with latent class models), or the BFGS algo-
rithm (Broyden–Fletcher–Goldfarb–Shanno) which is a quasi-Newton method.

• Sometimes difficult to disentangle β 1 and β 2 when both link functions have overlapping
covariates. That is because E(yi ) = µi πi , each factor involving β 1 and β 2 separately.

• For small to moderate sized data sets, need to be aware that the number of parameters is
doubled (if all covariates are used) so less precision may result.

Poisson hurdle models (Cragg, 1971)

Another possible way to address an excess of zero counts is to assume a different type of two-part
process:

• An observation either has a zero count or a positive count.

• If the observation has a positive count (“crossed the hurdle”), then the count follows a dis-
tribution over the positive integers.

Unlike the ZIP model, Pr(Y = 0) is not a mixture of two processes.

Poisson hurdle model If Yi follows a Poisson hurdle with parameters µi and πi , then

 (1 − πi ) for y = 0
Pr(Yi = y) = y
 π exp(−µi )µi /y! for y = 1, 2, . . ..
i 1−exp(−µi )

21
The term in the second line
exp(−µi )µyi /y!
1 − exp(−µi )
is the Poisson probability for Y = y, but normalized so that the sum over y ≥ 1 is 1.

That is, this is Pr(Y = y|Y ≥ 1) if Y ∼ Po(µ).

Actually, without any further discussion, the hurdle model is identical to the ZIP model.

Let πi = πi∗ (1 − exp(−µi )), Then the hurdle model becomes

(1 − [πi∗ (1 − exp(−µi ))) for y = 0



Pr(Yi = y) = y
 π ∗ (1 − exp(−µ )) exp(−µi )µi /y! for y = 1, 2, . . .
i i 1−exp(−µi )

or after simplification

(1 − πi∗ ) + πi∗ exp(−µi ) for y = 0

(
Pr(Yi = y) =
πi∗ exp(−µi )µyi /y! for y = 1, 2, . . ..

which is the ZIP model with parameters µi and πi∗ .

The way in which the models are distinguished is the application of link functions.

As with the ZIP model, let

g1 (πi ) = logit πi = xi β 1
g2 (µi ) = log µi = xi β 2

where some of the elements of β 1 and β 2 may be set to 0 by design.

Now the ZIP and hurdle Poisson models are distinct.

The log-likelihood conveniently separates into separate additive terms that can be maximized
separately:

ℓ(β 1 , β 2 |y) = ℓ1 (β 1 |y) + ℓ2 (β 2 |y) =

 
X Xn 
xi β 1 − log(1 + exp(xi β 1 ))
 
i:yi >0 i=1
  
X 
+ yi xi β 2 − exp(xi β 2 ) − log[1 − exp(− exp(xi β 2 ))]
 
i:yi >0

Note that the first bracketed term above is equivalent to the log-likelihood of a logistic regression
with response variable 1y>0 (binary indicator of y > 0).

The second bracketed term has no special convenient form.

Example: Boating trips again

22
> [Link] = hurdle(trips ˜quality+ski+income+userfee+
costC+costS+costH | quality+income,
dist="pois", data = RecreationDemand)
> summary([Link])

Count model coefficients (truncated poisson with log link):

Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.150401 0.111784 19.237 < 2e-16 ***
quality 0.044261 0.023851 1.856 0.0635 .
skiyes 0.467022 0.058789 7.944 1.96e-15 ***
income -0.097695 0.020575 -4.748 2.05e-06 ***
userfeeyes 0.600686 0.079524 7.554 4.24e-14 ***
costC 0.001409 0.003959 0.356 0.7219
costS -0.036606 0.002043 -17.922 < 2e-16 ***
costH 0.023876 0.003466 6.889 5.62e-12 ***

Zero hurdle model coefficients (binomial with logit link):

Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.76631 0.36230 -7.635 2.25e-14 ***
quality 1.50291 0.10032 14.981 < 2e-16 ***
income -0.04467 0.07853 -0.569 0.569
---
Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

Number of iterations in BFGS optimization: 16

Log-likelihood: -1188 on 11 Df

Comparison to log-linear model

• Cannot perform a likelihood ratio test like for zero-inflated Poisson models.

• The problem is that an ordinary log-linear model is not nested within the Poisson hurdle
model.

• Instead one can only perform likelihood ratio tests of nested hurdle models.

Observed and expected counts

> rbind(obs = table(RecreationDemand$trips)[1:10],

exp = round(colSums(
predict([Link], type=’prob’)[, 1:10]))))

0 1 2 3 4 5 6 7 8 9
obs 417 68 38 34 17 13 11 2 8 1
exp 417 28 31 31 29 25 21 17 13 10

23
Reasonably good (but not great) consistency, and certainly good at y = 0.

Comments

• Generally it is more difficult to imagine a sampling process that corresponds to a hurdle

models than a zero-inflated model, the latter being a mixture model.

• Given that Poisson hurdle and ZIP models are nearly identical except for the way in which
covariates are linked, it is not surprising that fitted values will be similar.

• While both “excess zero” models can address a particular deficiency with the ordinary Pois-
son model, they should not be considered the go-to model when a lack of fit exists in log-
linear models.

• We will see other count model alternatives to Poisson shortly.

24