- Email: [email protected]
Zero-inflated Poisson modeling to evaluate occupational safety interventions
Safety Science 41 (2003) 53–63 www.elsevier.com/locate/ssci
Zero-inflated Poisson modeling to evaluate occupational safety interventions Philip J.W. Carrivicka,*, Andy H. Leeb, Kelvin K.W. Yauc a
Department of Occupational Health, Sir Charles Gairdner Hospital, Hospital Avenue, Nedlands, WA 6009, Australia b School of Public Health, Curtin University of Technology, Perth, Australia c Department of Management Sciences, City University of Hong Kong, Hong Kong
Abstract Much of the data collected on occupational safety involves accident or injury counts. Traditional approaches to analyzing this type of population-based count data, including the Poisson model and the Negative Binomial model, do not take into consideration the relatively few events and hence, many observed zeros. In this paper, we propose a Zero-Inflated Poisson (ZIP) model to evaluate the effectiveness of a consultative manual handling workplace risk assessment team (WRATS) in reducing the risk of occupational injury among cleaners within a 600-bed hospital. We modify the ZIP model (which adjusts for the excess zeros) to take account of the extent of exposure in terms of hours worked during the study period. The findings highlight that the modified ZIP approach provided a satisfactory fit to the data, and that the manual handling WRATS intervention was associated with a reduction in the proportion of cleaners injured at work. In the presence of extra zeros, the ZIP model can be an alternative to the Poisson and Negative Binomial models, when appropriate according to statistical criteria. # 2002 Elsevier Science Ltd. All rights reserved. Keywords: Excess zeros; Exposure; Lost-time injury; Participatory ergonomics teams; Poisson distribution; Population-based count data; Safety intervention
1. Introduction Traditionally the Poisson distribution is used to model accident counts in occupational health (Altree-Williams, 1990; Bailer, 1997; Glazner et al., 1999). However, analysis of such discrete events may be hindered by their infrequent occurrence over extended time-periods (Quintana and Pawlowitz, 1999). This phenomenon is even * Corresponding author. Tel.: +61-8-9346-3414; fax: +61-8-9346-3100. E-mail address: [email protected] (P.J.W. Carrivick). 0925-7535/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0925-7535(01)00057-1
54
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
more pronounced when analyzing population-based data where there are few incidents and hence, many observed zeros. If the number of observed zeros far exceeds the expected number of zeros under the Poisson assumption, the implicit mean– variance relationship of the Poisson structure will be violated, and over-dispersion occurs (variance is greater than the mean). The general problem of over-dispersion has been considered in the statistical literature. An excellent review of approaches to analyze over-dispersed count data can be found in McLachlan (1997). In the current context, the extra-Poisson variation induced by the excess zeros can be accommodated through a compound probability model for the event, namely the Zero-Inflated Poisson (ZIP) distribution. This model has received much attention in the literature recently. For example, van den Broek (1995) demonstrated its usefulness in modeling the frequency of urinary tract infection of HIV-infected men, where the majority of patients have no episodes of infection. Similarly, Bo¨hning et al. (1999) applied the ZIP model to study caries prevention in dental epidemiology. In many applications, the accident or injury count data are typically observed in conjunction with the extent of exposure such as the number of hours worked over a specified period. In this paper the ZIP model is further modified to accommodate individual exposure information. Maximum likelihood estimation and asymptotic inference on the parameters are provided. Based on evidence drawn from an empirical study, detailed in Section 2, we recommend adjusting the Poisson model for both the extra zeros and the extent of exposure.
2. Empirical study This paper utilizes data from a quasi-experimental 7-year retrospective study into the effectiveness of a participatory ergonomics team. The study group was a population of cleaners employed by a 600-bed adult teaching hospital in Western Australia that had been experiencing historically high rates of injury relative to most other workers. The cleaners were exposed to risk control strategies recommended by an eight-member group of peers who had, with assistance from an ergonomist and support by management, formed a workplace risk assessment team (WRATS). WRATS developed and utilized a risk assessment tool, based on a manual handling code of practice, and met on 32 occasions over a 3-year period. Implemented hazard controls included changes to cleaners’ equipment, duties and working environment. Since it was impractical to randomly allocate cleaners to WRATS or not, the population of orderlies from the same hospital, with similar size and rates of injury to cleaners but not exposed to WRATS intervention, formed the comparison group. The orderlies were selected to help control for intra-hospital injury risk factors. Apart from their immediate supervisors, both groups shared the same line management, were employed under the same industrial award, and are categorized as Labourer under the Australian Bureau of Statistics occupational classification system. Orderly duties, like the cleaners, were predominantly manual handling, but the
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
55
nature of work was different. Although orderlies were not exposed to WRATS, they did continue with their usual induction training and had ongoing access to safety and health representatives, and to resolution of safety issues processes. WRATS commenced activity on 1 November 1992 and the groups were followed through until 31 October 1995. Pre-intervention data were collected for the 52 months from 1 July 1988 to 31 October 1992. All cleaners and all orderlies ever employed during the study period were identified. The average number of cleaners (approximately 145) and orderlies (140) working during the study period remained relatively constant; although the cleaners experienced a higher turnover (507 people ever employed) compared with the orderlies (279). The outcome of interest was any post-intervention change in the likelihood of an individual sustaining a lost-time injury (a new workers’ compensation injury resulting in the loss of a shift or more off work). Lost-time injuries (LTI) were categorized by date of occurrence and person. Total hours actually worked (i.e. not including leave) were determined pre- and post-intervention for each subject. This was to account for the difference in pre–post period duration, the variation in hours worked by each subject and the fact that not all subjects worked in both the pre- and post-intervention periods.
3. Methods 3.1. Negative binomial distribution A number of studies suggested the Negative Binomial (NB) as an alternative to the Poisson distribution when the count data are over- or under-dispersed in relation to the mean (Miaou, 1994; Poch and Mannering, 1996). The NB distribution allows for different levels of risk for different individuals. It explicitly accounts for the heterogeneity by modeling the Poisson mean as a Gamma random variable and introducing an extra dispersion parameter . Details of its derivation can be found in Poch and Mannering (1996). For a discrete random variable Y with a NB distribution PrðY ¼ yÞ ¼
yþk1 y
k þk
k
y ; þk
y ¼ 0; 1 . . .
where k>0, and =E(Y ). It can be shown that VarðY Þ ¼ þ 2 ; where =1/k is referred to as the dispersion parameter. Therefore, unlike the Poisson distribution which constrained the variance to be equal to the mean, Var(Y ) is allowed to exceed . The Poisson model can be regarded as a special case of the NB by letting tends to zero. The inappropriateness of the Poisson relative to the NB model is thus determined by the statistical significance of the parameter.
56
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
3.2. ZIP distribution ZIP is a model for count data with excess zeros. Consider a discrete random variable Y with a ZIP distribution (Johnson et al., 1992, p. 314): PrðY ¼ 0Þ ¼ þ ð1 Þe PrðY ¼ yÞ ¼ ð1 Þ
e ðÞy y ¼ 1; 2; . . . y!
where 0 < < 1 so that it incorporates extra zeros than those allowed by the Poisson (=0). The ZIP distribution may be regarded as a mixture of Poisson () and a degenerate component whereby, all of its mass is at zero. A plausible interpretation in the context of the cleaners population is its (unobserved) two-point heterogeneity: a sub-population of inherently ‘safe’ cleaners (contributing to the excess zeros), and another sub-population whose members are susceptible and may incur injury several times (due to their job nature and possibly their personal characteristics). It can be shown that EðY Þ ¼ ð1 Þ
and
VarðY Þ ¼ EðY Þ þ EðY Þ ð EðY ÞÞ
so that the ZIP model incorporates extra variation than the Poisson. Further properties and inference, including maximum likelihood estimation of and , can be found in Gupta et al. (1996). 3.3. Tests for over-dispersion and zero-inflation If the Poisson assumption is valid, the mean and variance should coincide. Bo¨hning et al. (1997, Eq. (2)) proposed an over-dispersion test that compares the sample mean and the sample variance. Their statistic has an approximate normal distribution. Provided that the over-dispersion test is significant, it is then important to determine whether the apparent over-dispersion is induced by the extra number of zeros. Discrimination between the Poisson and ZIP distributions amounts to testing the hypothesis H0:=0 against the alternative H1:>0. A score test for zeroinflation is proposed by van den Broek (1995). The score test has the advantage that one does not need to fit the ZIP but just the null Poisson distribution. The score statistic, as given by van den Broek (1995, Eq. (5)), has an approximate chi-squared distribution with one degree of freedom under H0, for large sample size n. Alternatively, the Wald statistic (^ over its standard error) and likelihood ratio test may also be used. The adequacy of the postulated model can then be assessed through the Pearson’s 2 goodness-of-fit statistic. 3.4. Conditional ZIP model incorporating exposure information To facilitate the interpretation of the ZIP model, let p=(1) (1e). The ZIP distribution can be equivalently expressed as:
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
57
PrðY ¼ 0Þ ¼ 1 p PrðY ¼ yÞ ¼ p
e y =ð1 e Þ y!
y ¼ 1; 2; . . .
In the above equation p is the probability of observing at least one injury and, conditional on at least one injury, the accident count is described by the truncated Poisson distribution. This conditional model has the advantage of an orthogonal parameterization, which is simple to interpret and readily fitted than its original formulation. In the context of our empirical study, individual exposure information are available. It is then necessary to adjust the ZIP model accounting for such information so that the mean parameter in the Poisson part can be interpreted as the mean injury rate. The ZIP model is now generalized to take into account of individual workplace exposure t: PrðY ¼ 0Þ ¼ 1 p PrðY ¼ yÞ ¼ p
et ðtÞy =ð1 et Þ y ¼ 1; 2; . . . y!
Suppose a sample of n observations consists of distinct injury counts y ð04y4kÞ with associated exposure tyi ; i ¼ 1; :::; ny ; where ny is the observed frequency of y. Maximum likelihood estimation and inference for parameters p and are provided in the Appendix. To evaluate whether the WRATS intervention is effective in reducing the injury incidence, one can test the hypothesis H0: p(pre-WRATS)=p(post-WRATS). Furthermore, for those that are injured, comparison of the pre- and post-mean injury rates may be based on the truncated Poisson mean ! = /(1e- ) of the conditional ZIP distribution, i.e. testing H0: !(pre-WRATS)=!(post-WRATS). It can be shown, using the delta method, that 2 1 e ð1 þ Þ Varð!^ Þ ffi Varð^ Þ ð1 e Þ2 and assuming the asymptotic normality of the maximum likelihood estimator ^ ; inference on ! can be drawn. A non-significant difference between !(pre-WRATS) and !(post-WRATS) will confirm the effectiveness of the intervention in the sense that injuries are not merely ‘shifted’ to those who are prone to accidents.
4. Results The ZIP models are fitted to data using macros written in S-Plus (MathSoft, 1999). The program codes are available from the authors. The NB model is also
58
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
fitted using the STATA statistical software (Stata Corporation, 1999). The empirical distributions of LTI counts before and after the WRATS intervention are given in the first two columns of Tables 1–4. It is evident that a peak of extra zeros (representing injury-free individuals) exists in each population. Such a pattern is typical of population-based occupational health data. Table 5 provides the summary statistics. Firstly, the over-dispersion test (Bo¨hning et al., 1997) delivers highly significant values, especially for the Cleaner population pre-WRATS, indicating strong overdispersion. The statistical significance of the NB dispersion parameter , as reflected by the likelihood ratio test, implies that the Poisson distribution is inappropriate for these data. Moreover, the score, Wald, and likelihood ratio statistics for testing the hypothesis H0:=0 are all significant (for brevity, only the score test results are reported). The tests indicated that the observed frequency of zeros far exceeds that expected under the Poisson assumption. In other words, the spurious overdispersion is merely an artifact of the zero-inflation problem. It is clear from Tables 1–4 that, according to the Pearson goodness-of-fit criterion, the homogeneous Poisson distribution fits badly to these data whereas the standard ZIP distribution does not indicate any lack of fit. To take account of the length of pre- and post-intervention periods and the extent of individual’s exposure to the Table 1 Observed and fitted frequency distributions for Cleaner pre-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5 6 7 8
267 52 16 4 2 0 0 0 1
246.49 80.72 13.22 1.44 0.12 0.01 0 0 0
267 47.21 20.39 5.87 1.27 0.22 0.03 0 0
267 49.38 17.30 5.98 1.77 0.45 0.10 0.02 0
270.71 44.56 14.98 6.19 2.81 1.35 0.67 0.34 0.18
^ =0.327 (0.031)
^ =0.864 (0.134) p^= 0.219 (0.022) ^ =0.621 (0.054)
^ =0.602 (0.093) p^= 0.219 (0.022)
^ =0.421 (0.051) ^=1.317 (0.431)
7.39 D=1.45 P-value=0.228
8.32 D=0.45 P-value=0.504
11.75 D=3.28 P-value=0.070
2+
23
3+
7
a b c
14.79 D=16.48 P-value<0.001
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
59
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63 Table 2 Observed and fitted frequency distributions for Cleaner post-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2+
194 17 4
191.40 22.26 1.34 ^ =0.116 (0.023)
194 17.45 3.55 ^=0.360 (0.175) p^=0.098 (0.020) ^ =0.677 (0.147)
194 17.6 3.4 ^ =0.373 (0.180) p^= 0.098 (0.020)
194.79 13.86 6.35 ^ =0.163 (0.040) ^=4.572 (2.831)
D=6.56 P-value=0.011
D=0.07
D=0.13
D=1.59
a b c
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
Table 3 Observed and fitted frequency distributions for Orderly pre-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5
138 40 16 4 3 2
120.43 62.88 16.42 2.86 0.37 0.04 ^ = 0.522 (0.051)
138 36.24 19.45 6.96 1.87 0.40 ^ =1.073 (0.157) p^=0.320 (0.033) ^ =0.513 (0.064)
138 37.29 18.00 6.96 2.10 0.52 ^=0.530 (0.077) p^=0.320 (0.033)
139.39 37.75 14.60 6.22 2.75 1.24 ^ =0.319 (0.041) ^=1.184 (0.406)
3+
9
3.27 D=20.94 P-value<0.001
9.31 D=1.01 P-value=0.314
8.71 D=0.43 P-value=0.513
11.27 D=0.74 P-value=0.390
a b c
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
workplace, the modified ZIP model of Section 3.4 is next fitted. We found that the modified ZIP model utilizing exposure information tyi provides a marginal improvement in overall fit, as evident from the smaller Pearson’s 2 goodness-of-fit statistic D values, when compared with the ZIP and NB models. To evaluate the effectiveness of the WRATS intervention, the hypotheses H0: p(pre-WRATS)=p(post-WRATS) and H0:!(pre-WRATS)=!(post-WRATS) are tested separately for the Cleaner and Orderly populations, in view of their different work duties. The results are graphically presented in Figs. 1 and 2. The intervention is
60
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
Table 4 Observed and fitted frequency distributions for Orderly post-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5
120 41 12 6 0 1
110.40 53.97 13.19 2.15 0.26 0.03
120 38.68 15.90 4.36 0.90 0.15
120 38.93 15.53 4.39 0.96 0.17
119.82 40.88 13.39 4.17 1.24 0.36
^ =0.489 (0.052)
^=0.822 (0.147) p^=0.333 (0.035) ^ =0.405 (0.093)
^ =0.547 (0.098) p^= 0.333 (0.035)
^ =0.413 (0.049) ^=0.447 (0.281)
5.42 D=1.56 P-value=0.212
5.54 D=1.297 P-value=0.255
5.91 D=0.35 P-value=0.554
2+
19
3+
7
a b c
15.63 D=4.68 P-value=0.031
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
Table 5 Summary statistics, over-dispersion test, score test for zero-inflation, and likelihood ratio test of =0 Population
n
LTIa
Total hours worked
Over-dispersion test
Score test for zero-inflation
Likelihood ratio test of =0
Cleaner pre-WRATSb
342
112
839,946
11.652 P-value<0.001
11.967 P-value<0.001
26.12 P-value<0.001
Cleaner post-WRATS
215
25
559,563
2.165 P-value=0.015
5.649 P-value=0.017
6.61 P-value=0.01
Orderly pre-WRATS
203
106
1,006,500
7.353 P-value<0.001
26.438 P-value<0.001
22.07 P-value<0.001
Orderly post-WRATS
180
88
621,356
4.050 P-value<0.001
9.621 P-value=0.002
4.16 P-value=0.041
a b
LTI, lost-time injuries. WRATS, workplace risk assessment team.
found to be significant in terms of reducing the injury incidence of the Cleaner population but not the Orderly population. The p^ estimates show that approximately 22% of the Cleaner population has been involved in at least one accident during the pre-intervention period, as contrasted to the 10% post-intervention. No significant difference between pre- and post-! is detected. This confirms that the
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
61
Fig. 1. Comparison of pre- and post-workplace risk assessment team (WRATS) injury incidence and mean injury rate for cleaners, with standard errors of estimates within parentheses.
Fig. 2. Comparison of pre- and post-workplace risk assessment team (WRATS) injury incidence and mean injury rate for orderlies, with standard errors of estimates within parentheses.
decrease in proportion of injured cleaners has not been achieved by penalising those already at higher risk of injury. However, it is of interest to remark that, if the standard ZIP model without adjusting for exposure were adopted, estimates of the Cleaners pre- and post-! are 1.493 (S.E.=0.086) and 1.191 (S.E.=0.098), respectively. The test statistic then becomes 2.324 which has a P-value of 0.01, implying a spurious significance of the intervention in reducing the Cleaner mean injury rates for those that are injured. For the improperly specified Poisson model, testing the mean injury pre- and post- WRATS outputs a test statistic of 5.466 with P-value < 0.01 leading again to the erroneous conclusion. This highlights a major finding, namely, the importance of adjusting the count model for both the extra zeros and the extent of exposure, in order to avoid spurious outcomes.
5. Discussion In evaluating intervention outcomes it is crucial that appropriate statistical methods have been utilized in determining results (Shannon et al., 1999; Rushton, 2000). In this regard, a potential for drawing misleading inferences exists if an inappropriate distributional assumption is made, particularly when excess zeros are observed due to the infrequent event. Ironically this is more likely where a successful intervention reduces injury count. To avoid this, when analyzing population-based accident counts, it is essential to assess the empirical distribution of the data carefully prior to adopting the homogeneous Poisson assumption. In our study, both the NB
62
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
and ZIP models provided adequate fits to the data, despite the different motivation behind the models, the former simply assumes the accident risk for individuals to be heterogeneous. In the presence of concomitant information, a Poisson distribution with log-linear modeling to include covariates is the standard method of analysis. However, a ZIP regression model (Lambert, 1992) is recommended as an alternative to Poisson regression once zero-inflation is confirmed. In our context, a modified ZIP regression accounting for the extent of individual workplace exposure appears worthwhile to further analyze factors affecting injury risk. Preferably, after allowing for exposure, covariates (such as age, gender and work experience) would be utilized for ZIP regression. In this empirical study, these data were not available.
Appendix The log-likelihood function for the conditional ZIP model incorporating individual exposure tyi is given by ny ny k X k k X X X X tyi þ ðlnÞ yny lnð1 etyi Þ ‘ ¼ n0 lnð1 pÞ þ ðn n0 Þlnp y¼1 i¼1
y¼1
y¼1 i¼1
The score equations, upon simplifying, are then @‘ n0 n n0 ¼ þ ¼0 @p 1p p ny k k X X tyi @‘ 1 X ¼ yny ¼0 @ y¼1 1 etyi y¼1 i¼1 The first equation yields p^ ¼ ðn n0 Þ=n whereas the second equation yields ^ ¼
ny k k X X X yny = y¼1
y¼1 i¼1
tyi 1 e^tyi
which can be solved iteratively for the maximum likelihood estimator ^ : The second derivatives of the log-likelihood are @2 ‘ n0 n n0 ¼ 2 2 @p p2 ð1 pÞ ny k k X X t2yi etyi @2 ‘ 1X ¼ yn þ y tyi Þ2 @2 2 y¼1 y¼1 i¼1 ð1 e
@2 ‘ ¼0 @p@
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
63
Based on the observed information, 2 1 @ ‘ p^ð1 p^Þ ¼ Varðp^Þ ¼ 2 @p p^ n !1 2 1 ny k k X X X t2yi e^tyi @ ‘ 1 Varð^ Þ ¼ 2 yny ¼ ^tyi Þ2 @ ^ ^ 2 y¼1 y¼1 i¼1 ð1 e Covðp^; ^ Þ ¼ 0 which can be used to construct asymptotic confidence intervals and make inference for the parameters p and . Furthermore, estimates for the standard ZIP model are found by setting tyi =1. References Altree-Williams, S., 1990. The appropriate measure for work injury rate. Journal of Occupational Health and Safety Australia and New Zealand 6, 199–204. Bailer, A.J., 1997. Modeling fatal injury rates using poisson regression: a case study of workers in agriculture, forestry and fishing. Journal of Safety Research 28, 177–186. Bo¨hning, D., Dietz, E., Schlattmann, P., 1997. Zero-inflated count models and their applications in public health and social science. In: Rost, J., Langeheine, R. (Eds.), Applications of Latent Trait and Latent Class Models in the Social Sciences. Waxmann, Mu¨nster, pp. 333–344. Bo¨hning, D., Dietz, E., Schlattmann, P., Mendonc¸a, L., Kirchner, U., 1999. The zero-inflated Poisson model and the decayed, missing and filled teeth index in dental epidemiology. Journal of Royal Statistical Society A 162, 195–209. Glazner, J.E., Borgerding, J., Bondy, J., et al., 1999. Contractor safety practices and injury rates in construction of the Denver international airport. American Journal of Industrial Medicine 35, 175–185. Gupta, P.L., Gupta, R.C., Tripathi, R.C., 1996. Analysis of zero-adjusted count data. Computational Statistics & Data Analysis 23, 207–218. Johnson, N.L., Kotz, S., Kemp, A.W., 1992. Univariate Discrete Distributions, 2nd Edition. Wiley, New York. Lambert, D., 1992. Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34, 1–14. MathSoft, 1999. S-PLUS 2000: Modern Statistics and Advanced Graphics. MathSoft Inc, Seattle. McLachlan, G.J., 1997. On the EM algorithm for overdispersed count data. Statistical Methods in Medical Research 6, 76–98. Miaou, S.P., 1994. The relationship between truck accidents and geometric design of road sections: poisson versus negative binomial regressions. Accident Analysis and Prevention 26, 471–482. Poch, M., Mannering, F., 1996. Negative binomial analysis of intersection-accident frequencies. Journal of Transportation Engineering 122, 105–113. Quintana, R., Pawlowitz, I., 1999. A poisson model for work-related musculoskeletal disorder cost estimation. Safety Science 32, 19–31. Rushton, L., 2000. Reporting of occupational and environmental research: use and misuse of statistical and epidemiological methods. Occupational and Environmental medicine 57, 1–9. Shannon, H.S., Robson, L.S., Guastello, S.J., 1999. Methodological criteria for evaluating occupational safety intervention research. Safety Science 31, 161–179. Stata Corporation, 1999. STATA Release 6 User’s Guide. Stata Press, College Station, Texas. Van den Broek, J., 1995. A score test for zero inflation in a poisson distribution. Biometrics 51, 738–743.
Zero-inflated Poisson modeling to evaluate occupational safety interventions Philip J.W. Carrivicka,*, Andy H. Leeb, Kelvin K.W. Yauc a
Department of Occupational Health, Sir Charles Gairdner Hospital, Hospital Avenue, Nedlands, WA 6009, Australia b School of Public Health, Curtin University of Technology, Perth, Australia c Department of Management Sciences, City University of Hong Kong, Hong Kong
Abstract Much of the data collected on occupational safety involves accident or injury counts. Traditional approaches to analyzing this type of population-based count data, including the Poisson model and the Negative Binomial model, do not take into consideration the relatively few events and hence, many observed zeros. In this paper, we propose a Zero-Inflated Poisson (ZIP) model to evaluate the effectiveness of a consultative manual handling workplace risk assessment team (WRATS) in reducing the risk of occupational injury among cleaners within a 600-bed hospital. We modify the ZIP model (which adjusts for the excess zeros) to take account of the extent of exposure in terms of hours worked during the study period. The findings highlight that the modified ZIP approach provided a satisfactory fit to the data, and that the manual handling WRATS intervention was associated with a reduction in the proportion of cleaners injured at work. In the presence of extra zeros, the ZIP model can be an alternative to the Poisson and Negative Binomial models, when appropriate according to statistical criteria. # 2002 Elsevier Science Ltd. All rights reserved. Keywords: Excess zeros; Exposure; Lost-time injury; Participatory ergonomics teams; Poisson distribution; Population-based count data; Safety intervention
1. Introduction Traditionally the Poisson distribution is used to model accident counts in occupational health (Altree-Williams, 1990; Bailer, 1997; Glazner et al., 1999). However, analysis of such discrete events may be hindered by their infrequent occurrence over extended time-periods (Quintana and Pawlowitz, 1999). This phenomenon is even * Corresponding author. Tel.: +61-8-9346-3414; fax: +61-8-9346-3100. E-mail address: [email protected] (P.J.W. Carrivick). 0925-7535/02/$ - see front matter # 2002 Elsevier Science Ltd. All rights reserved. PII: S0925-7535(01)00057-1
54
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
more pronounced when analyzing population-based data where there are few incidents and hence, many observed zeros. If the number of observed zeros far exceeds the expected number of zeros under the Poisson assumption, the implicit mean– variance relationship of the Poisson structure will be violated, and over-dispersion occurs (variance is greater than the mean). The general problem of over-dispersion has been considered in the statistical literature. An excellent review of approaches to analyze over-dispersed count data can be found in McLachlan (1997). In the current context, the extra-Poisson variation induced by the excess zeros can be accommodated through a compound probability model for the event, namely the Zero-Inflated Poisson (ZIP) distribution. This model has received much attention in the literature recently. For example, van den Broek (1995) demonstrated its usefulness in modeling the frequency of urinary tract infection of HIV-infected men, where the majority of patients have no episodes of infection. Similarly, Bo¨hning et al. (1999) applied the ZIP model to study caries prevention in dental epidemiology. In many applications, the accident or injury count data are typically observed in conjunction with the extent of exposure such as the number of hours worked over a specified period. In this paper the ZIP model is further modified to accommodate individual exposure information. Maximum likelihood estimation and asymptotic inference on the parameters are provided. Based on evidence drawn from an empirical study, detailed in Section 2, we recommend adjusting the Poisson model for both the extra zeros and the extent of exposure.
2. Empirical study This paper utilizes data from a quasi-experimental 7-year retrospective study into the effectiveness of a participatory ergonomics team. The study group was a population of cleaners employed by a 600-bed adult teaching hospital in Western Australia that had been experiencing historically high rates of injury relative to most other workers. The cleaners were exposed to risk control strategies recommended by an eight-member group of peers who had, with assistance from an ergonomist and support by management, formed a workplace risk assessment team (WRATS). WRATS developed and utilized a risk assessment tool, based on a manual handling code of practice, and met on 32 occasions over a 3-year period. Implemented hazard controls included changes to cleaners’ equipment, duties and working environment. Since it was impractical to randomly allocate cleaners to WRATS or not, the population of orderlies from the same hospital, with similar size and rates of injury to cleaners but not exposed to WRATS intervention, formed the comparison group. The orderlies were selected to help control for intra-hospital injury risk factors. Apart from their immediate supervisors, both groups shared the same line management, were employed under the same industrial award, and are categorized as Labourer under the Australian Bureau of Statistics occupational classification system. Orderly duties, like the cleaners, were predominantly manual handling, but the
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
55
nature of work was different. Although orderlies were not exposed to WRATS, they did continue with their usual induction training and had ongoing access to safety and health representatives, and to resolution of safety issues processes. WRATS commenced activity on 1 November 1992 and the groups were followed through until 31 October 1995. Pre-intervention data were collected for the 52 months from 1 July 1988 to 31 October 1992. All cleaners and all orderlies ever employed during the study period were identified. The average number of cleaners (approximately 145) and orderlies (140) working during the study period remained relatively constant; although the cleaners experienced a higher turnover (507 people ever employed) compared with the orderlies (279). The outcome of interest was any post-intervention change in the likelihood of an individual sustaining a lost-time injury (a new workers’ compensation injury resulting in the loss of a shift or more off work). Lost-time injuries (LTI) were categorized by date of occurrence and person. Total hours actually worked (i.e. not including leave) were determined pre- and post-intervention for each subject. This was to account for the difference in pre–post period duration, the variation in hours worked by each subject and the fact that not all subjects worked in both the pre- and post-intervention periods.
3. Methods 3.1. Negative binomial distribution A number of studies suggested the Negative Binomial (NB) as an alternative to the Poisson distribution when the count data are over- or under-dispersed in relation to the mean (Miaou, 1994; Poch and Mannering, 1996). The NB distribution allows for different levels of risk for different individuals. It explicitly accounts for the heterogeneity by modeling the Poisson mean as a Gamma random variable and introducing an extra dispersion parameter . Details of its derivation can be found in Poch and Mannering (1996). For a discrete random variable Y with a NB distribution PrðY ¼ yÞ ¼
yþk1 y
k þk
k
y ; þk
y ¼ 0; 1 . . .
where k>0, and =E(Y ). It can be shown that VarðY Þ ¼ þ 2 ; where =1/k is referred to as the dispersion parameter. Therefore, unlike the Poisson distribution which constrained the variance to be equal to the mean, Var(Y ) is allowed to exceed . The Poisson model can be regarded as a special case of the NB by letting tends to zero. The inappropriateness of the Poisson relative to the NB model is thus determined by the statistical significance of the parameter.
56
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
3.2. ZIP distribution ZIP is a model for count data with excess zeros. Consider a discrete random variable Y with a ZIP distribution (Johnson et al., 1992, p. 314): PrðY ¼ 0Þ ¼ þ ð1 Þe PrðY ¼ yÞ ¼ ð1 Þ
e ðÞy y ¼ 1; 2; . . . y!
where 0 < < 1 so that it incorporates extra zeros than those allowed by the Poisson (=0). The ZIP distribution may be regarded as a mixture of Poisson () and a degenerate component whereby, all of its mass is at zero. A plausible interpretation in the context of the cleaners population is its (unobserved) two-point heterogeneity: a sub-population of inherently ‘safe’ cleaners (contributing to the excess zeros), and another sub-population whose members are susceptible and may incur injury several times (due to their job nature and possibly their personal characteristics). It can be shown that EðY Þ ¼ ð1 Þ
and
VarðY Þ ¼ EðY Þ þ EðY Þ ð EðY ÞÞ
so that the ZIP model incorporates extra variation than the Poisson. Further properties and inference, including maximum likelihood estimation of and , can be found in Gupta et al. (1996). 3.3. Tests for over-dispersion and zero-inflation If the Poisson assumption is valid, the mean and variance should coincide. Bo¨hning et al. (1997, Eq. (2)) proposed an over-dispersion test that compares the sample mean and the sample variance. Their statistic has an approximate normal distribution. Provided that the over-dispersion test is significant, it is then important to determine whether the apparent over-dispersion is induced by the extra number of zeros. Discrimination between the Poisson and ZIP distributions amounts to testing the hypothesis H0:=0 against the alternative H1:>0. A score test for zeroinflation is proposed by van den Broek (1995). The score test has the advantage that one does not need to fit the ZIP but just the null Poisson distribution. The score statistic, as given by van den Broek (1995, Eq. (5)), has an approximate chi-squared distribution with one degree of freedom under H0, for large sample size n. Alternatively, the Wald statistic (^ over its standard error) and likelihood ratio test may also be used. The adequacy of the postulated model can then be assessed through the Pearson’s 2 goodness-of-fit statistic. 3.4. Conditional ZIP model incorporating exposure information To facilitate the interpretation of the ZIP model, let p=(1) (1e). The ZIP distribution can be equivalently expressed as:
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
57
PrðY ¼ 0Þ ¼ 1 p PrðY ¼ yÞ ¼ p
e y =ð1 e Þ y!
y ¼ 1; 2; . . .
In the above equation p is the probability of observing at least one injury and, conditional on at least one injury, the accident count is described by the truncated Poisson distribution. This conditional model has the advantage of an orthogonal parameterization, which is simple to interpret and readily fitted than its original formulation. In the context of our empirical study, individual exposure information are available. It is then necessary to adjust the ZIP model accounting for such information so that the mean parameter in the Poisson part can be interpreted as the mean injury rate. The ZIP model is now generalized to take into account of individual workplace exposure t: PrðY ¼ 0Þ ¼ 1 p PrðY ¼ yÞ ¼ p
et ðtÞy =ð1 et Þ y ¼ 1; 2; . . . y!
Suppose a sample of n observations consists of distinct injury counts y ð04y4kÞ with associated exposure tyi ; i ¼ 1; :::; ny ; where ny is the observed frequency of y. Maximum likelihood estimation and inference for parameters p and are provided in the Appendix. To evaluate whether the WRATS intervention is effective in reducing the injury incidence, one can test the hypothesis H0: p(pre-WRATS)=p(post-WRATS). Furthermore, for those that are injured, comparison of the pre- and post-mean injury rates may be based on the truncated Poisson mean ! = /(1e- ) of the conditional ZIP distribution, i.e. testing H0: !(pre-WRATS)=!(post-WRATS). It can be shown, using the delta method, that 2 1 e ð1 þ Þ Varð!^ Þ ffi Varð^ Þ ð1 e Þ2 and assuming the asymptotic normality of the maximum likelihood estimator ^ ; inference on ! can be drawn. A non-significant difference between !(pre-WRATS) and !(post-WRATS) will confirm the effectiveness of the intervention in the sense that injuries are not merely ‘shifted’ to those who are prone to accidents.
4. Results The ZIP models are fitted to data using macros written in S-Plus (MathSoft, 1999). The program codes are available from the authors. The NB model is also
58
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
fitted using the STATA statistical software (Stata Corporation, 1999). The empirical distributions of LTI counts before and after the WRATS intervention are given in the first two columns of Tables 1–4. It is evident that a peak of extra zeros (representing injury-free individuals) exists in each population. Such a pattern is typical of population-based occupational health data. Table 5 provides the summary statistics. Firstly, the over-dispersion test (Bo¨hning et al., 1997) delivers highly significant values, especially for the Cleaner population pre-WRATS, indicating strong overdispersion. The statistical significance of the NB dispersion parameter , as reflected by the likelihood ratio test, implies that the Poisson distribution is inappropriate for these data. Moreover, the score, Wald, and likelihood ratio statistics for testing the hypothesis H0:=0 are all significant (for brevity, only the score test results are reported). The tests indicated that the observed frequency of zeros far exceeds that expected under the Poisson assumption. In other words, the spurious overdispersion is merely an artifact of the zero-inflation problem. It is clear from Tables 1–4 that, according to the Pearson goodness-of-fit criterion, the homogeneous Poisson distribution fits badly to these data whereas the standard ZIP distribution does not indicate any lack of fit. To take account of the length of pre- and post-intervention periods and the extent of individual’s exposure to the Table 1 Observed and fitted frequency distributions for Cleaner pre-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5 6 7 8
267 52 16 4 2 0 0 0 1
246.49 80.72 13.22 1.44 0.12 0.01 0 0 0
267 47.21 20.39 5.87 1.27 0.22 0.03 0 0
267 49.38 17.30 5.98 1.77 0.45 0.10 0.02 0
270.71 44.56 14.98 6.19 2.81 1.35 0.67 0.34 0.18
^ =0.327 (0.031)
^ =0.864 (0.134) p^= 0.219 (0.022) ^ =0.621 (0.054)
^ =0.602 (0.093) p^= 0.219 (0.022)
^ =0.421 (0.051) ^=1.317 (0.431)
7.39 D=1.45 P-value=0.228
8.32 D=0.45 P-value=0.504
11.75 D=3.28 P-value=0.070
2+
23
3+
7
a b c
14.79 D=16.48 P-value<0.001
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
59
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63 Table 2 Observed and fitted frequency distributions for Cleaner post-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2+
194 17 4
191.40 22.26 1.34 ^ =0.116 (0.023)
194 17.45 3.55 ^=0.360 (0.175) p^=0.098 (0.020) ^ =0.677 (0.147)
194 17.6 3.4 ^ =0.373 (0.180) p^= 0.098 (0.020)
194.79 13.86 6.35 ^ =0.163 (0.040) ^=4.572 (2.831)
D=6.56 P-value=0.011
D=0.07
D=0.13
D=1.59
a b c
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
Table 3 Observed and fitted frequency distributions for Orderly pre-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5
138 40 16 4 3 2
120.43 62.88 16.42 2.86 0.37 0.04 ^ = 0.522 (0.051)
138 36.24 19.45 6.96 1.87 0.40 ^ =1.073 (0.157) p^=0.320 (0.033) ^ =0.513 (0.064)
138 37.29 18.00 6.96 2.10 0.52 ^=0.530 (0.077) p^=0.320 (0.033)
139.39 37.75 14.60 6.22 2.75 1.24 ^ =0.319 (0.041) ^=1.184 (0.406)
3+
9
3.27 D=20.94 P-value<0.001
9.31 D=1.01 P-value=0.314
8.71 D=0.43 P-value=0.513
11.27 D=0.74 P-value=0.390
a b c
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
workplace, the modified ZIP model of Section 3.4 is next fitted. We found that the modified ZIP model utilizing exposure information tyi provides a marginal improvement in overall fit, as evident from the smaller Pearson’s 2 goodness-of-fit statistic D values, when compared with the ZIP and NB models. To evaluate the effectiveness of the WRATS intervention, the hypotheses H0: p(pre-WRATS)=p(post-WRATS) and H0:!(pre-WRATS)=!(post-WRATS) are tested separately for the Cleaner and Orderly populations, in view of their different work duties. The results are graphically presented in Figs. 1 and 2. The intervention is
60
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
Table 4 Observed and fitted frequency distributions for Orderly post-WRATSa Injury count
Observed frequency
Poisson
ZIPb
ZIP (with t)
NBc (with t)
0 1 2 3 4 5
120 41 12 6 0 1
110.40 53.97 13.19 2.15 0.26 0.03
120 38.68 15.90 4.36 0.90 0.15
120 38.93 15.53 4.39 0.96 0.17
119.82 40.88 13.39 4.17 1.24 0.36
^ =0.489 (0.052)
^=0.822 (0.147) p^=0.333 (0.035) ^ =0.405 (0.093)
^ =0.547 (0.098) p^= 0.333 (0.035)
^ =0.413 (0.049) ^=0.447 (0.281)
5.42 D=1.56 P-value=0.212
5.54 D=1.297 P-value=0.255
5.91 D=0.35 P-value=0.554
2+
19
3+
7
a b c
15.63 D=4.68 P-value=0.031
WRATS, workplace risk assessment team. ZIP, zero-inflated poisson. NB, negavitive binominal.
Table 5 Summary statistics, over-dispersion test, score test for zero-inflation, and likelihood ratio test of =0 Population
n
LTIa
Total hours worked
Over-dispersion test
Score test for zero-inflation
Likelihood ratio test of =0
Cleaner pre-WRATSb
342
112
839,946
11.652 P-value<0.001
11.967 P-value<0.001
26.12 P-value<0.001
Cleaner post-WRATS
215
25
559,563
2.165 P-value=0.015
5.649 P-value=0.017
6.61 P-value=0.01
Orderly pre-WRATS
203
106
1,006,500
7.353 P-value<0.001
26.438 P-value<0.001
22.07 P-value<0.001
Orderly post-WRATS
180
88
621,356
4.050 P-value<0.001
9.621 P-value=0.002
4.16 P-value=0.041
a b
LTI, lost-time injuries. WRATS, workplace risk assessment team.
found to be significant in terms of reducing the injury incidence of the Cleaner population but not the Orderly population. The p^ estimates show that approximately 22% of the Cleaner population has been involved in at least one accident during the pre-intervention period, as contrasted to the 10% post-intervention. No significant difference between pre- and post-! is detected. This confirms that the
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
61
Fig. 1. Comparison of pre- and post-workplace risk assessment team (WRATS) injury incidence and mean injury rate for cleaners, with standard errors of estimates within parentheses.
Fig. 2. Comparison of pre- and post-workplace risk assessment team (WRATS) injury incidence and mean injury rate for orderlies, with standard errors of estimates within parentheses.
decrease in proportion of injured cleaners has not been achieved by penalising those already at higher risk of injury. However, it is of interest to remark that, if the standard ZIP model without adjusting for exposure were adopted, estimates of the Cleaners pre- and post-! are 1.493 (S.E.=0.086) and 1.191 (S.E.=0.098), respectively. The test statistic then becomes 2.324 which has a P-value of 0.01, implying a spurious significance of the intervention in reducing the Cleaner mean injury rates for those that are injured. For the improperly specified Poisson model, testing the mean injury pre- and post- WRATS outputs a test statistic of 5.466 with P-value < 0.01 leading again to the erroneous conclusion. This highlights a major finding, namely, the importance of adjusting the count model for both the extra zeros and the extent of exposure, in order to avoid spurious outcomes.
5. Discussion In evaluating intervention outcomes it is crucial that appropriate statistical methods have been utilized in determining results (Shannon et al., 1999; Rushton, 2000). In this regard, a potential for drawing misleading inferences exists if an inappropriate distributional assumption is made, particularly when excess zeros are observed due to the infrequent event. Ironically this is more likely where a successful intervention reduces injury count. To avoid this, when analyzing population-based accident counts, it is essential to assess the empirical distribution of the data carefully prior to adopting the homogeneous Poisson assumption. In our study, both the NB
62
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
and ZIP models provided adequate fits to the data, despite the different motivation behind the models, the former simply assumes the accident risk for individuals to be heterogeneous. In the presence of concomitant information, a Poisson distribution with log-linear modeling to include covariates is the standard method of analysis. However, a ZIP regression model (Lambert, 1992) is recommended as an alternative to Poisson regression once zero-inflation is confirmed. In our context, a modified ZIP regression accounting for the extent of individual workplace exposure appears worthwhile to further analyze factors affecting injury risk. Preferably, after allowing for exposure, covariates (such as age, gender and work experience) would be utilized for ZIP regression. In this empirical study, these data were not available.
Appendix The log-likelihood function for the conditional ZIP model incorporating individual exposure tyi is given by ny ny k X k k X X X X tyi þ ðlnÞ yny lnð1 etyi Þ ‘ ¼ n0 lnð1 pÞ þ ðn n0 Þlnp y¼1 i¼1
y¼1
y¼1 i¼1
The score equations, upon simplifying, are then @‘ n0 n n0 ¼ þ ¼0 @p 1p p ny k k X X tyi @‘ 1 X ¼ yny ¼0 @ y¼1 1 etyi y¼1 i¼1 The first equation yields p^ ¼ ðn n0 Þ=n whereas the second equation yields ^ ¼
ny k k X X X yny = y¼1
y¼1 i¼1
tyi 1 e^tyi
which can be solved iteratively for the maximum likelihood estimator ^ : The second derivatives of the log-likelihood are @2 ‘ n0 n n0 ¼ 2 2 @p p2 ð1 pÞ ny k k X X t2yi etyi @2 ‘ 1X ¼ yn þ y tyi Þ2 @2 2 y¼1 y¼1 i¼1 ð1 e
@2 ‘ ¼0 @p@
P.J.W. Carrivick et al. / Safety Science 41 (2003) 53–63
63
Based on the observed information, 2 1 @ ‘ p^ð1 p^Þ ¼ Varðp^Þ ¼ 2 @p p^ n !1 2 1 ny k k X X X t2yi e^tyi @ ‘ 1 Varð^ Þ ¼ 2 yny ¼ ^tyi Þ2 @ ^ ^ 2 y¼1 y¼1 i¼1 ð1 e Covðp^; ^ Þ ¼ 0 which can be used to construct asymptotic confidence intervals and make inference for the parameters p and . Furthermore, estimates for the standard ZIP model are found by setting tyi =1. References Altree-Williams, S., 1990. The appropriate measure for work injury rate. Journal of Occupational Health and Safety Australia and New Zealand 6, 199–204. Bailer, A.J., 1997. Modeling fatal injury rates using poisson regression: a case study of workers in agriculture, forestry and fishing. Journal of Safety Research 28, 177–186. Bo¨hning, D., Dietz, E., Schlattmann, P., 1997. Zero-inflated count models and their applications in public health and social science. In: Rost, J., Langeheine, R. (Eds.), Applications of Latent Trait and Latent Class Models in the Social Sciences. Waxmann, Mu¨nster, pp. 333–344. Bo¨hning, D., Dietz, E., Schlattmann, P., Mendonc¸a, L., Kirchner, U., 1999. The zero-inflated Poisson model and the decayed, missing and filled teeth index in dental epidemiology. Journal of Royal Statistical Society A 162, 195–209. Glazner, J.E., Borgerding, J., Bondy, J., et al., 1999. Contractor safety practices and injury rates in construction of the Denver international airport. American Journal of Industrial Medicine 35, 175–185. Gupta, P.L., Gupta, R.C., Tripathi, R.C., 1996. Analysis of zero-adjusted count data. Computational Statistics & Data Analysis 23, 207–218. Johnson, N.L., Kotz, S., Kemp, A.W., 1992. Univariate Discrete Distributions, 2nd Edition. Wiley, New York. Lambert, D., 1992. Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34, 1–14. MathSoft, 1999. S-PLUS 2000: Modern Statistics and Advanced Graphics. MathSoft Inc, Seattle. McLachlan, G.J., 1997. On the EM algorithm for overdispersed count data. Statistical Methods in Medical Research 6, 76–98. Miaou, S.P., 1994. The relationship between truck accidents and geometric design of road sections: poisson versus negative binomial regressions. Accident Analysis and Prevention 26, 471–482. Poch, M., Mannering, F., 1996. Negative binomial analysis of intersection-accident frequencies. Journal of Transportation Engineering 122, 105–113. Quintana, R., Pawlowitz, I., 1999. A poisson model for work-related musculoskeletal disorder cost estimation. Safety Science 32, 19–31. Rushton, L., 2000. Reporting of occupational and environmental research: use and misuse of statistical and epidemiological methods. Occupational and Environmental medicine 57, 1–9. Shannon, H.S., Robson, L.S., Guastello, S.J., 1999. Methodological criteria for evaluating occupational safety intervention research. Safety Science 31, 161–179. Stata Corporation, 1999. STATA Release 6 User’s Guide. Stata Press, College Station, Texas. Van den Broek, J., 1995. A score test for zero inflation in a poisson distribution. Biometrics 51, 738–743.