A comparison between Poisson and zero-inflated Poisson regression models with an application to number of black spots in Corriedale sheep

Dark spots in the fleece area are often associated with dark fibres in wool, which limits its competitiveness with other textile fibres. Field data from a sheep experiment in Uruguay revealed an excess number of zeros for dark spots. We compared the performance of four Poisson and zero-inflated Poisson (ZIP) models under four simulation scenarios. All models performed reasonably well under the same scenario for which the data were simulated. The deviance information criterion favoured a Poisson model with residual, while the ZIP model with a residual gave estimates closer to their true values under all simulation scenarios. Both Poisson and ZIP models with an error term at the regression level performed better than their counterparts without such an error. Field data from Corriedale sheep were analysed with Poisson and ZIP models with residuals. Parameter estimates were similar for both models. Although the posterior distribution of the sire variance was skewed due to a small number of rams in the dataset, the median of this variance suggested a scope for genetic selection. The main environmental factor was the age of the sheep at shearing. In summary, age related processes seem to drive the number of dark spots in this breed of sheep.


INTRODUCTION
The presence of black-brown fibres in wool from Corriedale sheep is recognised as a fault [13,20]. This issue limits the competitiveness of wool with other textile fibres and reduces its value by 15-18% when the number exceeds 300 fibresAEkg À1 top (Frank Racket, 1997, personal communication). In Uruguayan wool, this value can be as large as 5000 fibresAEkg À1 top, with most of the dark fibres having an environmental origin, e.g. faeces and urine dyeing [3,18]. With appropriate clip preparation, values ranging from 800 to 1000 fibres have been found, and these probably have a genetic background. Skin spots with blackbrown fibres and isolated pigmented fibres are the probable origin of these fibres [2,9,12,20].
With the aim of investigating factors involved in the development of pigmented fibres, an experiment was carried out in which fleeces of animals from two experimental flocks were sampled yearly at shearing for laboratory analysis. Each animal was inspected, and the number of black spots, their diameter and the estimated percentage of dark fibres in each spot were recorded. While genetic selection should focus on reducing the number of dark fibres, it is expensive and cumbersome to record such a value for each animal on a routine basis. Laboratory techniques are labour intensive and slow.
In this context, the number of dark spots in the fleece area of animals may be a useful indicator trait, for several reasons. First, our empirical observations suggest that dark fibres are associated with dark spots, hinting a positive correlation between the two variables. Second, spots can be assessed easily and quickly, and scoring is less subjective than for other candidate measures such as the percentage of spot area with dark fibres [1,10,11]. Third, we have observed that in spots without or with dark fibres in young animals, the presence of dark fibres increases with age. Hence, the presence of spots indicates dark fibres in adult animals. If laboratory analyses confirm that black spots are positively correlated with the number of dark fibres, recording on a nation-wide basis would be straightforward.
Previous studies in Romney sheep [6] have addressed the occurrence of black wool spots at weaning (BWS w ) and at yearling age (BWS y ). Enns and Nicoll [6] used a threshold model for a binary response variable (the presence or absence of pigmented spots), and their largest heritability estimates were 0.070 (0.018) and 0.072 (0.014) for BWS w and BWS y , respectively. In contrast, in our research, focus has been on modelling the number of dark spots in each animal, irrespective of the presence of dark fibres. As a count variable, the number of spots could plausibly follow a Poisson distribution. However, as shown in Figure 1, there is an excess of zeros in the empirical distribution for field records, 380 relative to their expected value under Poisson sampling with a homogeneous parameter. If Y follows a Poisson distribution, then E(Y) = Var(Y), where E(AE) and Var(AE) represent the mean and variance, respectively. In a Poisson distribution, the variance-to-mean-ratio (VTMR) is 1. In the observed data in Figure 1, VTMR was 6.8. A zero-inflated Poisson model (ZIP) [17], may provide a better description of the data. This model assumes that observations come from one of two different components, a ''perfect'' state which produces only zeros with probability h, and an ''imperfect'' one that follows a Poisson distribution, with probability (1 À h) and Poisson parameter k. It can be shown that the mean and variance of a ZIP variate are and respectively, which accounts for VTMR > 1 provided that overdispersion arises from an excess of zeros. Zero-inflated models for count data in animal breeding have been discussed by Gianola [15] and used by Rodrigues-Motta and collaborators [27] in an analysis of the number of mastitis cases in dairy cattle. Figure 1. Distribution of the number of black spots in field data (n = 497). The solid line represents the best fit of a Poisson distribution to the observed data, fitted with package ''gnlm'' (http://popgen.unimaas.nl/~jlindsey/rcode.html) of R [26].

Poisson and ZIP models for spots in sheep
From previous exploratory analysis [16,24], the age of animals appears to be a main source of variability of the number of spots, with flock and year having marginal effects. Modelling can proceed along the lines of generalised linear models [21] or generalised linear mixed models [25] provided that the link function used is appropriate. However when mixture distributions are assumed, as in the ZIP model, estimation is more involved, since indicator variables (e.g., from which of the two states a zero originates) are not observed. However, imputation of non-observed parameters given the data fits naturally in the Bayesian framework [19]. In recent years, Bayesian Markov chain Monte Carlo (MCMC) methods have become widely used in animal breeding [28], as a powerful and flexible tool. An advantage of the Bayesian MCMC framework is that it is relatively easy to implement measures of model quality such as posterior predictive ability (PPA) checks.
In this paper, four different candidate models for the number of spots were compared. Poisson and ZIP models were considered, with the log of the Poisson parameter of each of the models regressed on environmental and genetic effects. The two models were extended further to include a random residual in the regression, aimed to capture overdispersion other than that due to extra zeros. Two of the models were selected and fitted to a sample of Corriedale sheep to obtain estimates of population parameters.

Simulation
Four different scenarios (H1-H4) were simulated as described in Table I. The rationale underlying the models is that the observed number of spots in each animal follows a Poisson distribution with the logarithm of its parameter expressed as a linear model. The Poisson distribution does not accommodate well the overdispersion caused by excess zeros, so a ZIP model is a reasonable competitor. Furthermore, the parameter of the Poisson distribution represents the expected propensity of spots, so an additional error (residual) term at the regression level allows modelling individual differences in propensity. The two models (Poisson and ZIP), each with or without residuals, give the four models (P, Z, Pe and Ze) studied.
Data were generated from either ZIP (H1, H2) or Poisson (H3, H4) distributions; the log of the Poisson parameter contained (H2, H4) or did not contain (H1, H3) a random residual. In all four models, the ram effects were assumed to follow independent normal distributions with null mean and variance r 2 ram ; the residual was independent and identically distributed as e i;j;k $ N ð0; r 2 e Þ (H2, H4).

382
In each scenario, 100 datasets (replicates) were randomly generated, with 1000 observations each. For each animal, the covariate age was randomly sampled, resembling the distribution of the age in the observed data. Forty rams (sires) were sampled in each dataset and randomly assigned to observations. In each scenario, the true parameters were selected to resemble the observed distribution of spots.

Models fitted in the simulation
Four models were fitted to the simulated data (Z, Ze, P and Pe), each matching a specific scenario, as shown in Table I. Models are connected as illustrated in Figure 2. A path between two models involves fixing or adding a single parameter. Preliminary analysis indicated that flock and year effects (and their interaction) had minor importance, so these factors were not included in the simulations. However, when models were fitted to the real data, the regression models included flock and year effects.

Bayesian computation
Parameter inference was done using the OpenBUGS software [31]. Vague priors were assigned to represent initial uncertainty. A normal distribution centred at zero with precision 0.01 was used for location parameters, while a Gamma (0.01, 0.01) distribution was assumed for each of the two variance parameters. Several different hyper-parameter values were assigned in pilot runs, with the only observable difference being the time needed to attain convergence. For each scenario and model, the burn-in period was determined from preliminary runs, based on four chains, starting at different points. Final runs were performed with two chains each. The burn-in period was of 10 000 iterations, and samples were obtained from the following 10 000 iterations, without thinning. An exception was model Pe in scenario H4, where the required burn-in period was 30 000 iterations.
The b's are unknown regressions.
Poisson and ZIP models for spots in sheep

End points for model comparison
Models were contrasted first through simulated data (comparing true and estimated parameter values) and by using the deviance information criterion (DIC), estimates of marginal likelihoods with the method of Newton and Raftery [22] and via PPA. The DIC [30] was obtained directly from OpenBUGS. PPA was patterned after Sorensen and Waagepetersen [29]. Suppose that for a model M, h  as discrepancy statistic, where k k in the numerator is the mean, and ffiffiffiffi ffi k k p in the denominator is the standard deviation; for the ZIP model, the mean and standard deviation were replaced by their corresponding values. For simplicity, dependence between observations from the same animal was ignored, so that the only source of correlation considered was that resulting from a half-sib family structure. Clearly, the limited dataset precludes precise estimation of genetic parameters, but this was not an objective of this study.

Simulations
For each scenario simulated, the results are presented for the ''true'' model and for the other three models. Values of the DIC (highlighting pD, the ''effective number of parameters'') and of the difference statistic used for PPA are shown in Table II.

DIC
In scenario H1, where Z is the true model, Pe performed better (lower DIC) than the true model, in spite of the penalty resulting from a larger pD (number of parameters). The value of nearly 400 effective parameters in 1000 observations indicates that very few observations clustered under the same Poisson distribution. Models with residuals had a higher pD but lower deviance. Except for the P model, the other specifications had similar DIC, at least in the light of the between replicates standard deviation. Clearly the Poisson model was the worst under the ''true ZIP'' scenario.
In scenario H2 (Ze is the true model), Pe was, again, better than the true model and the picture with respect to pD was as in the H1 scenario, although differences between models with residuals were smaller. Models without residuals had the poorest performance; P had the worse DIC.

Poisson and ZIP models for spots in sheep
Under H3 (P is the true model) all DIC values were similar. Models with residuals had smaller pD than in scenarios H1, H2, probably due to the simpler nature of this simulation scenario. Finally, in scenario H4, the true model (Pe) was best under the DIC, followed by Ze.
The global picture is clearer when the number of times (in 100 simulations) in which each model had the smallest DIC was considered (Tab. III). Model Pe outperformed other models except under H3. Notably, DIC selected the right model only in 172 out of 400 comparisons (43% of the time).

PPA
Values of the PPA difference statistic close to zero indicate essentially no differences between observed and predicted responses. In regard to the PPA results, the true model always predicted best, and this was essentially true for all scenarios (Tab. II). The pure Poisson model (P) performed badly in ZIP scenarios (H1 and H2), while Ze did reasonably well in all four scenarios. In H3, PPA was similar for all models. The problem with this criterion seems to be its low discriminative power, relative to its high standard deviations over replications. Alternatives to PPA are cross-validation techniques, but these were not considered due to computational expense.

Marginal likelihood
It was impossible to calculate the Bayes factor for several pairs of models, given the huge differences in marginal likelihoods. For this reason, only estimates of marginal log-likelihood are presented for each model and scenario (Tab. IV). On the basis of this criterion, Pe was the best model in all scenarios.

Parameter inference
As expected, parameter estimates were in agreement with their ''true'' values when a model matched its corresponding scenario (Tab. V). However, when models pertained to a different scenario, their performances were markedly different. Regressions on age were well inferred, but estimated intercepts b 0 were severely understated when Poisson models were applied to ZIP scenarios. Model Ze estimates were always in agreement with the ''true'' values, regardless of the scenario. Pe model estimates of intercept and of the residual variance were strongly biased in ZIP scenarios. Finally, models with residuals estimated the sire (ram) variance well.
The ability of different models to predict breeding values is of interest. Given that the ''true'' values of rams were known, their Spearman rank correlation with the predicted values (posterior mean) was calculated for each combination of scenarios and models. Ze was best in all scenarios, but differences were small (Fig. 3). The ZIP model performed reasonably well in all scenarios, while pure Poisson models did well only in their own scenarios, with median correlations between 0.41 and 0.53.

Field data
Given the simulation results, data were analysed using the two ''best'' models, Pe and Ze, including flock, year, flock · year and age effect as covariates, and ram as a genetic effect. The effects were removed successively from each main model and DIC was computed; full models were found to be the best. Since the effects of flock, year and flock-year interaction were negligible, these are not reported in Table VI. DIC (pD) was 1103 (200.6) and 1070 (201.8) for the Ze and Pe models, respectively. The number of spots increased with age, and there was no strong evidence of inflation at zero (the estimate of the probability of the perfect state, h, was 0.05). It seems that most of the overdispersion is due to unaccounted for between-individual variability.  In standard animal breeding theory, in a sire model, heritability is defined as: Estimates of heritability (in a log-scale) under the two models considered are displayed in Table VI. Posterior medians were 0.25 and 0.17 for the Ze and Pe model, respectively, but the distribution was very skewed due to the few rams (19) used in the study.

DISCUSSION
While easier to measure than the number of dark fibres per animal, modelling the number of spots poses several challenges in regard to standard methodology of animal breeding. It is very difficult to obtain a good fit of the data with simple linear models with normal distributions for the random effects. Frequently used Box-Cox transformations (e.g., log or reciprocal) cannot be used, given the number of zeros. Additionally, there is the issue of an excess of zeros relative to Poisson sampling. One attractive model for dealing with this is the ZIP. One can think of a fraction h of ''perfect'' animals that will never develop spots, whereas others only will develop spots at random, following a Poisson distribution with parameter k. Moreover, variation in the k's can be accounted for by a model including environmental and genetic factors, such as age, flock, year or ram.
We first compared the performance of four models (i.e., Poisson and ZIP with and without an error term in the regression) using a simulation that resembled the field data structure. Based on the end points considered, two ''competitive''  [5]. This provides a very flexible structure, which explains at least partially, the good performance of these two models in the simulation results.
Using the DIC, the Pe model was best in most scenarios, in spite of a larger pD, a term that penalises the likelihood and that is associated with the effective number of parameters of the model. The results of the marginal likelihoods also supported this interpretation (Tab. IV), while the difference statistic for the posterior predictive ability (DPPA) displayed low discriminative power. However, the Pe model did not produce good estimates of the intercept (number of spots at age = 0, i.e. birth) and of the residual variance, when data were generated from ZIP distributions (Tab. V). This suggests that even the results of DIC should be viewed with some caution. Furthermore, as pointed out by several discussants in Spiegelhalter et al. [30], the DIC may underpenalise model complexity. Several alternative versions of DIC were proposed by Celeux et al. [4] to address models with missing data or mixtures of distributions. However, despite important differences in performance, each alternative proposed has its own drawbacks and no single solution emerges as unanimously appropriate.
On the contrary, Ze was robust across all situations, since it estimated the true parameters well. A ''stable model'' is appealing under practical conditions. In animal breeding, a stable model with good predictive ability is desired. All models produced a good agreement between ''true'' and predicted breeding values, especially Ze, which maintained its ability across scenarios.
Based on the simulation, Pe and Ze were chosen to analyse the field dataset. As expected, under the DIC Pe outperformed Ze. However, parameter estimates were similar (Tab. VI). This may be explained by the fact that the estimate of h in Ze was low, pointing to a relatively small effect of ''perfect'' individuals on inference when the Poisson model includes a residual.
A series of environmental and genetic factors may be related to the number of spots. Simple observation (even in humans) suggests that this number increases with age and environmental stress factors (e.g., solar irradiation can be invoked as causative agents [7,8,14,23]). Variability in the underlying genetic mechanisms responsible for the spots is likely, at least in different races, as well as in susceptibility to environmental stress factors.
The age of the animals was the main environmental factor to consider, consistently. However, it is not known if this relationship arises from an intrinsic ageing process independent of environmental factors, or if environmental stressors such as sun irradiation drive the process. Anyhow, it is possible to envisage management measures aiming to reduce incidence of dark fibres. If an intrinsic ageing process is the main factor, reducing the age at shearing could be a practice to take into consideration. This requires additional research.
In an animal breeding context, genetic and environmental variances are extremely important since they define heritability, a key parameter used to select among breeding strategies. The meaning of heritability in non-linear hierarchical models, such as Pe or Ze, is not straightforward. However, the magnitude of heritability suggests scope for genetic selection. Different simulations indicated that predicted breeding values (for log k) were in good agreement with ''true'' values, so these models are probably useful for selection purposes.
In summary, hierarchical models for count data were studied with the aim of defining strategies for reducing incidence of dark fibres in wool from Corriedale sheep. ZIP and Poisson models with random residuals performed better than their counterparts without residuals. Ageing related processes seem to drive the number of dark spots in sheep, and further research should be done to address this underlying phenomenon.