A Graphical Tool for Assessing the Suitability of a Count Regression Model

Whilst many numeric methods, such as AIC and deviance, exist for assessing or comparing model ﬁt, diagrammatic methods are few. We present here a diagnostic plot, which we refer to as a ‘Quantile Band plot’, that may be used to visually assess the suitability of a given count data model. In the case of diagnosed model inadequacy, the plot has the unique feature of conveying precise information on the character of the violation, hence pointing the data analyst towards a potentially better model choice.


Introduction
Consider univariate count data Y = {Y 1 , . . ., Y n }, possibly accompanied by associated covariate vectors x i ∈ R p , i = 1. . . ., n.Any attempt at drawing inferential conclusions from these observations will require, at first place, a choice of the response distribution, which determines the likelihood function and hence impacts on all further analysis, whether this is done in a Bayesian or frequentist framework.Typical choices for count data distributions would be for instance the Poisson distribution (which describes the number of independent events occurring at constant rate within a certain time interval), the Negative Binomial distribution (which can be considered as an overdispersed extension of the Poisson distribution), the Poisson Inverse Gaussian distribution, the Zero-Inflated Poisson distribution, and many more.
Despite the availability of this plethora of possible count data models, applied users will generally strive for using the most simple model wherever possible, which, in this framework, is the Poisson distribution.The immediate question is then whether a particular choice of distribution is appropriate; which can be phrased as the question of whether the observed data are 'plausible' given the properties of the count data distribution being postulated.One can disassemble this question as asking whether the observed number of zeros is plausible given the postulated count data model (in other words, whether the number of zeros is 'consistent' with the number of zeros predicted by this model), whether the number of ones is plausible given the postulated model, and so on.Clearly, due to the inherent randomness of the system, there will generally be more than one 'number of counts k' which is plausible, and henceforth the task will be to provide an appropriate range of plausible numbers of counts, as a function of k.The methodology presented in this manuscript will do exactly this, and it will summarize the information, for all count values k = 0, 1, 2, 3 . .., in a novel diagrammatic tool.
The resulting diagnostic plot may be used to visually assess the suitability of a given count data model.If it is determined that the model in question is not suitable, the plot has the unique feature of conveying precise information on the character of the violation.For instance, this tool will, at a glance, allow the user to detect not only the presence of zeroinflation (more zeros than could be plausibly expected under the given count data model), but also yet poorly explored features such as for instance '3-deflation' (less observations of the count 3 than would be expected under the count data model) or '10-inflation'.This turns out to be highly relevant for assessing the presence of digit preference in a given count data set.The insights gained through this diagnostic plot should also help the data analyst in making a more informed model choice; for instance it is clear that in the case of detected zero-inflation in a Poisson model, a zero-inflated Poisson (ZIP) model should be more adequate.
Denote the postulated response distribution by G(θ i ), where θ i ∈ R d is a multivariate parameter vector, where each of the d components may or may not depend on covariates.It is helpful to think of θ i as being composed of a mean component, µ i = E(Y i |x i ), with x i ∈ R p , and the remaining (shape, scale, dispersion,...) parameters bundled in δ i , so that θ i = (µ i , δ i ), though it is not strictly necessary that any of the d components actually represents the mean.The δ i may depend on the same or other covariates as µ i , but will often be assumed not to depend on covariates at all.In the special case that G is the Poisson or geometric distribution, δ i is empty.
The question of interest is whether the data Y are plausible given the distributional assumption G; that is, whether it can be plausibly assumed that Y have in fact been generated from G. At some occasion, the parameters θ i may be fixed and known (even if they depend on i), but more often they will be unknown and need to be estimated from the data.We will, initially, not distinguish between these two cases, that is we assume that routines to obtain estimates θi are readily available.
We would like to make the point that we do not consider the developed graphical tool as a technique to test for the adequacy of the predictor specification; for instance, of the mean component µ i = E(Y i |x i ).Obviously, this specification is an important part of the model, but the overall distribution of the number of zeros, 1's, etc, will generally not depend strongly on it.Our concern is the suitability of the distribution G, given a certain choice of θ i .
To fix the notation more precisely, denote p i (k) = P (k|θ i ) the probability of observing the count k under covariate x i and model G, which can be estimated by pi (k) = P (k| θi ) from the fitted model.For instance, in the special case that G(θ i ) corresponds to Pois(µ i ), one has pi This scenario is discussed in Wilson and Einbeck (2018) with focus on the case k = 0.This manuscript generalizes those ideas to general k and G and proposes a generic diagrammatic tool.
Denote by N (k) the 'counts of counts', that is the number of occurrences of count k among the data in Y, where k≥0 N (k) = n.It is clear that, for fixed k, N (k) can be described by a sum of n Bernoulli trials with success probabilities p i (k), i = 1, . . ., n.The resulting distribution, of which one can think as a Binomial distribution with unequal success probabilities, is known as a Poisson-Binomial distribution (Chen and Liu 1997), some properties of which we summarize in Appendix A. Hence, for any choice of k and G, a range of plausible values of N (k) can be obtained from this distribution, using fitted success probabilities pi (k) as model parameters.By doing this for a range of values of k, one can draw diagrams which give envelopes for plausible values of N (k) which can then be compared to the true values.For reasons that will become clear in later sections, we refer to such diagrams as Quantile Band plots.
The remainder of this exposition is organized as follows.The graphical tool will be presented and explained in systematic form in Section 2, using an example involving digit preference for illustrative purposes.Computational details of the methodology as well as the problem of parameter estimation are deferred to Section 3. Further examples, including real data examples, follow in Section 4, before the paper is concluded in Section 5. Some complementary technicalities and definitions are included in the Appendices.

Quantile band plots
Based on the principles outlined above, we propose here a diagnostic plot to visually assess the suitability of a given model for the data.We firstly present the algorithm for the construction of the plot.

Algorithm for plot construction
The first aspect to decide on is the range K = [k min , k max ] of count values that is to be assessed.Typical choices would be k min = 0 and k max = max(Y) (and this is what will be used by the default in the graphical tool).Other choices may be preferable in specific circumstances.
A diagnostic plot may be constructed as follows.(The items labelled by a * symbol are to be understood as optional.While the Quantile Band plot as advocated in this work includes the execution of these optional items, there may be certain situation when the data analyst might prefer to omit them; for instance if the quantitative information on the count frequencies is to be conveyed through the plot.)Specification Determine the model G(θ i ) for the data Y.Obtain estimates θi , i = 1, . . ., n where required.

Create graph
(i) Plot the functions b α/2 (k) and b 1−α/2 (k) versus k.Then add to the plot the observed shifted counts, A(k), of the observed data Y. [If item (iii) above has not been carried out, replace b and A by q and N , respectively, and in this case one may optionally add m(k) to the plot.]*(ii) Rotate the plot by 90 degrees, so that k is orientated along the vertical axis.
If the data is consistent with the distribution fitted, the curve A(k) [N (k), respectively] should (largely) stay within the bands b α/2 (k) and b 1−α/2 (k) [q α/2 (k) and q 1−α/2 (k), respectively].If the data is not consistent with the distribution fitted then A(k) [N (k)] is likely not to stay within these bands.Informally we will refer to the line representing the b(k)'s as the upper and lower quantile bands, and the line representing the A(k)'s as the count-line.
Note that there are several possible choices of how exactly to compute quantiles for a discrete distribution.For our purposes, the quantiles employed are the mid-quantiles, which are discussed in Section 3.1.

Illustration via digit preference
Digit Preference refers to the mis-reporting of some numbers in favour of "preferred numbers".An early example in the literature is Myers (2002) who reports a tendency in the the 1910, 1920, and 1930 U.S. censuses to report ages of 20 as 21 and ages of 31 as 30.Camarda, Eilers, and Gampe (2002) provide an extensive list of references to literature concerning digit preference, and report that maybe the most common form of the phenomenon is "heaping" of data at multiples of 5.
We present here a hypothetical example that simulates a situation where the number of, say, annual theatre visits follows a Poisson distribution with parameter 7, however in practice the following occurs: • visit counts of 0, 1 and 2 are correctly reported; • visit counts that are a multiple of 5 are correctly reported; • visit counts that are within 2 of a multiple of 5 are reported correctly with probability φ and as that multiple of 5 with probability 1 − φ.The question of interest is whether the data could plausibly have been generated from a Poisson distribution.Therefore, the above procedure is applied to the data of Table 1, choosing K = [0, 16] as the range of counts of interest, and omitting the two 'starred' items for now.
The values of N (k), q 0.025 (k) and q 0.975 (k), for k = 0, 1, . . ., 7 are given in Table 2; the full table is given in Table 9 of Appendix D. From the left-hand part of this table one can construct the diagnostic plot displayed in Figure 1.Clearly N (5) > q 0.975 (5) and N (10) > q 0.975 (10) as is illustrated by the red line which represents the N (k) lying above the green line representing the q 0.975 (k) for k = 5 and k = 10.Similarly the red line representing the N (k) lies beneath the green line representing the q 0.025 (k) for some other values of k, indicating the unsuitability of a Poisson model here.
Figure 1: Diagnostic plot for the data of Table 1.We refer to plots in this form as raw Quantile Band plots.
Whilst the plot of Figure 1 contains complete information, the range of N (k) leads to the bands described by N (k), q 0.025 (k) and q 0.975 (k) being somewhat close.For data where the range of N (k) is large, it frequently occurs that these bands become extremely close and difficult to distinguish (examples of such plots are given in Figure 9 in Appendix E).
A superior plot may be obtained by including the 'starred items' from Section 2.1.After subtracting the median from all other quantities we arrive at the information displayed in the right part of Table 2, and hence the Quantile Band plot in the top panel of Figure 2. Similar to the display of boxplots, which can be presented either vertically or horizontally, one may display the Quantile Band plot in either orientation.We prefer the vertically rotated version, as shown in the bottom panel of Figure 2, and will use this version throughout the remainder of this paper.We acknowledge however that others may prefer the horizontal version, and also that there are occasions when practical considerations such as availability of space in a publication render the horizontal version preferable.One immediately concludes from this plot that there is evidence of inflation of multiples of five, and also some evidence of deflation of the counts four and nine (which is arguably an artifact of the former).Overall, this gives evidence that the Poisson assumption is not adequate.We further discuss this example in Section 4.1.
Enhanced interpretations on other diagnostic aspects can often be drawn from the specific nature of the pattern; these will be discussed in the examples of Section 4. We will turn now to the question of how exactly the quantiles of the distribution of N (k) are obtained.

Quantiles and mid-quantiles
We consider now the question of how exactly to compute the quantiles of the distribution of N (k).Recall that for each (fixed) value of k, N (k) is described as a Poisson-Binomial distribution with parameters pi (k), i = 1, . . ., n.We denote by F the distribution function of The γ-quantile of this distribution is traditionally defined as where t is a non-negative real value.For the purposes of Section 2 one will have either γ = α/2 or γ = 1 − α/2.For the Poisson-Binomial distribution, these quantiles can be obtained straightforwardly from the package poibin (Hong 2013).However, as Ma, Genton, and Parzen (2011) have pointed out, there are "many drawbacks" with using traditional quantiles for discrete distributions.Firstly, the quantile function is discontinuous, leading to unfavourable theoretical properties.Secondly, it causes interpretational problems: For instance, for all observable values t ∈ N = {0, . . .n} of N (k), the sum of the (right-handed) p-value and the left-handed quantile turns out to be larger than 1, since the probability mass at t is 'double-counted'.
The suggested improvement by Ma et al. (2011) is to split the probability mass at discrete values accordingly, which gives rise to the definition of the mid-distribution function F mid (t) = F (t) − 0.5p(t), where p(t) is the corresponding point mass at t (which is equal to 0 for all t ∈ N ).
This concept then leads to the definition of mid-quantiles, and is also naturally related to the notion of mid-p-values (Franck 1986).The mid-quantile-function, which we denote by Q mid (γ) henceforth, is constructed so that, at the points of observable values t , one has Q mid (F mid (t )) = t , and a piecewise linear interpolation in between (Ma et al, 2011).Since the exact mathematical expression of mid-quantiles in itself is rather complicated, we have deferred this alongside with a tutorial to Appendix B.
In all computations leading to Quantile Band plots in this paper, we have used for step (ii) in Section 2 the convention q γ (k) = Q mid (γ), with the distribution F which gives rise to this quantile being the distribution of N (k) corresponding to an underlying count distribution G.
Of course, also the medians m(k) are computed as 'mid-medians' in this manner.Wilson and Einbeck (2018) constructed mid-quantiles in the special case k = 0 and G = Pois, and referred to the resulting intervals [Q mid (α/2), Q mid (1 − α/2)) as mid-quantile intervals (MQI).In this spirit, we will also refer to tables which provide information as in the left half of Table 2 as mid-quantile tables.

Parameter estimation
The view taken in our approach is that the production of the Quantile Band plot succeeds the parameter estimation.In other words, it is assumed that, once a user has specified a certain count data model G(θ i ), a routine to estimate θ i is readily available.The distribution G and the estimated θi then serve as input to the production of the Quantile Band plot.In this sense, the estimation of the θ i is not considered as an intrinsic part of the methodology as such.Some words on parameter estimation still appear in order.
Usually, the θi will be estimated through Maximum Likelihood (ML).A frequent scenario is where θ i = (µ i , δ i ), where δ i is some dispersion or shape parameter.This would be the case, for instance, for the Negative Binomial (Type I or II), Neyman Type A, the Pòlya-Aeppli, or the Poisson-Inverse Gauss distribution.For, say, the latter case one may consider a model with constant dispersion index δ i ≡ δ, and a log-link for the mean regression parameter, i.e. log(µ , for some vector of predictors x i ∈ R p .The actual vector of parameters to estimate would then be (β, δ) ∈ R p+1 .But of course, it is also possible to describe the dispersion or other parameters through appropriate linear predictor terms, which may involve the same or different covariates as the mean function.
Having obtained the underlying estimates ( β, δ), one can immediately obtain the estimates θi via the pre-specified predictor configurations, which then allows production of the pi (k) = P (k| θi ), and hence execution of the machinery outlined in Section 2.
The actual calculation of the involved MLE's will usually be carried out through software.For instance, in the case that G corresponds to a Poisson or Binomial distribution, in which cases d = 1, estimation of θ i reduces to fitting a generalized linear model (McCullagh and Nelder 1989) and hence can be carried out using the glm function in R. A wide range of further multi-parameter count data distributions can be fitted using functionalities provided by the R packages VGAM (Yee 2010) and gamlss (Rigby and Stasinopoulos 2005).Estimation routines for some further count distributions, including the ones mentioned earlier in this subsection, which are relevant for specific applications for instance in dosimetry, are provided in form of R code in the supplementary material of Oliveira, Einbeck, Higueras, Ainsbury, Puig, and Rothkamm (2016).
However, it needs to be pointed out at this occasion that Maximum Likelihood is not the only way to estimate model parameters.For instance, it is long known that the ML estimate of the Poisson mean, that is, the whole sample mean, can perform very poorly if the data are zero-inflated or zero-deflated (or, to condense these two terms, 'zero-modified', see also da Silva, Ribeiro, Conceiçáo, Andrade, and Louzada ( 2018)).In response to such problems, Plackett (1953), Irwin (1959) and Ridout and C.B. (1992) present formulae for estimating the Poisson parameter from the mean of the positive data values; that is from the zero-truncated data.These estimators reduce the bias of the ML estimator of the mean parameter, but turn out to be less precise than the maximum likelihood estimator.Wilson and Einbeck (2018) went one step further and suggested (for scenarios in which zero-modification is expected or suspected) to balance bias and variance of the Poisson mean estimate through a weighted mean of the whole sample estimator and the zero-truncated estimator, with a weight of 2/3 for the whole sample mean performing favorably in simulation studies.In the case of actual zero-modification, this 'hybrid' estimator will considerably reduce the bias, but in its absence it will not behave notably differently than the whole sample estimator.It is emphasized that the problem being solved through such measures is intrinsic to zero counts.For inflation or deflation of higher counts k, the impact on the estimates of µ i is much less severe since the effects tend to cancel out over neighbouring counts.For multi-parameter distributions, the second parameter can absorb the overdispersion created through the excess zeros to some extent.Hence, in line with this reasoning, we use the hybrid estimation technique in the applications in Section 4 only for one-parameter distributions, i.e. the Poisson and the geometric distribution.The example in Section 4.3 will illustrate the impact of the different mean estimators on the Quantile Band plot explicitly.

Multiple testing issues
One may argue that due to the consideration of a sequence of mid-quantile intervals for N (0), N (1), . . .one has to account for multiple testing issues.It is certainly true that the count-line being outside of the mid-quantile bands at N (0) is equivalent to determining that there is zero-modification relative to the null model, as in the test of Wilson and Einbeck (2018), similar statements being possible for other N (k)'s, and that if several such tests were performed simultaneously then multiple testing issues would arise.The pragmatic view is that our proposed plot should not be considered as a testing procedure, but as a simple diagrammatic tool which supports the data analyst in identifying potential model inadequacies, similar in spirit to a QQ plot.

Digit preference revisited
We saw in Section 2.2 that a Poisson model is unsuitable for the data of Table 1.Whilst in the introduction to that example the data generating mechanism is informally described, even if this had not been the case the Quantile Band plots of Figure 2 clearly illustrate that the observed numbers of multiples of five are considerably greater than is to be expected under a Poisson model, and most other values somewhat less than would be expected, which would lead the analyst to suspect that digit preference for multiples of five is a feature of the data.The mid-quantile table corresponding to the Poisson model modified for digit preference is given in Table 10, Appendix D. Agresti (2002) discusses the fitting of a variety of models to data concerning the responses to the question: "Within the last twelve months, how many people have you known personally that were victims of homicide?"The respondents were classified as Black or White according to their race.The data is summarised here in Table 3.  4. (For the Geometric and Poisson models the mean is modelled by race using a log-link, for the PIG and negative binomial models both parameters are modelled by race using log-links, and for for the ZIP model, the mean and zero-inflation parameters are both modelled by race using log and logit-links respectively.)Quantile Band plots corresponding to the various models are illustrated in Figure 4, with the corresponding mid-quantile tables provided in Appendix D (Tables 11 to 16).It is apparent that the one parameter geometric and Poisson models are inadequate; we do see however from the plots for the Poisson and geometric models that whilst both underestimate the numbers of zeros, and overestimate the number of 1's, the over-estimation of 1's is not as severe for the geometric model as it is with the Poisson, but the over estimation of larger values is more severe under the geometric model than under the Poisson.It is most interesting to compare the two-parameter models.In all cases the observed counts all fall within their respective mid-quantile intervals, so it may be argued that all four are suitable.Reflecting its position as the poorest of the two parameter models considered here by both AIC and BIC criteria, it is noticeable that the PIG model considerably overestimates the number of 1's, and considerably underestimates the number of 2's and 3's.The Quantile Band plots for the two types of negative binomial are identical; it is notable that the fits of the negative binomial models and the zero-inflated Poisson models are similar under both the AIC and BIC criteria, it is apparent however that the ZIP model slightly underestimates the numbers of 1's, and overestimates the numbers of 2's compared to the observed data, the reverse being the case for the negative binomial models.

Frequency of homicides
We include in Figure 9, Appendix E the mid-quantile band plots obtained when Poisson and PIG models are fitted to the data of Table 3.These plots illustrate the superiority of the median-adjusted plots of Figure 4.

Choosing between Poisson and ZIP models
In Section 3.2 we discussed the estimation of model parameters, and mentioned three estimators of the parameter of a Poisson model that have been proposed in the literature.We present here an example that further explores this issue, and illustrates the usefulness of Quantile Band plots as a tool in determining the appropriateness of a given model, under consideration of different possible estimates of the model parameters.
Consider the data of Table 5, which are generated by concatenating 200 zeros to a sample of 800 data drawn from a Pois(0.5)distribution, thus the data are zero-inflated by construction.The mean of these data is 0.435 and the variance 0.484, indicating at first glance that a Poisson model may not be unreasonable.
Table 5: Simulated data with zero-inflation 0 1 2 3 4 5 663 256 67 12 1 1 This is now a situation as touched upon in Section 3.2, where the Poisson mean estimate is possibly affected by the presence of zero-inflation.Hence, we consider in this example the application of three different estimators of the Poisson mean, which are the whole sample mean (μ W ), zero-truncated mean (μ T ), and hybrid mean (μ H ), respectively.Table 6 gives the corresponding three estimates along with the log-likelihoods of the Poisson models when using those estimates.It is apparent that, under the log-likelihood criterion, the Pois(μ W ) model has the best fit (indeed, as μW is the maximum likelihood estimator this must be the case), the fit of the Pois(μ H ) model is only slightly poorer, and that of the Pois(μ T ) model more so.The log-likelihood statistics on their own however do not say anything about the suitability of the models: it is possible that all or none are compatible with the observed data.
Table 6: Model fits using whole sample mean (μ W ), zero-truncated mean (μ T ), and hybrid mean (μ H ), respectively.μ log-likelihood μW = 0.435 −873.0 μT = 0.534 −882.9 μH = 0.468 −874.2 Therefore, we proceed with the production of Quantile Band plots, which are depicted for the three Poisson models as well as the ZIP model in Figure 5.The two plots in the top row show, interestingly, that both a Poisson model with mean parameter of 0.435 and a ZIP model are compatible with the data.What has happened here is that the presence of extra zeros has reduced the value of μW (the whole sample mean) to a value which renders the numbers of observed zeros and 1's compatible with a Pois(μ W ) model.There is no contradiction here, and both statements are correct: a zero-inflated model with small mean parameter is simply hard to distinguish from a Poisson model with even smaller mean parameter (note that the Poisson mean estimate under the ZIP model is 0.533 and the estimated zero-inflation parameter is 0.184).
The bottom plots show that Poisson models with larger mean parameters, as obtained through the truncated and hybrid estimators, are not compatible with the data.Also these are true statements: The plots simply assert that Poisson models with means of 0.468 and 0.534, respectively, are incompatible with the data.However, arising from this is the interesting question which of the three Poisson-based plots tells the most meaningful story in terms of the data generating mechanism.From this angle, the hybrid and the truncated estimators do the more useful job, by producing Poisson mean estimates which are closer to the true value (of 0.5), allowing for the indication of the presence of inflation of zeros relative to the Poisson model.
Finally, it is worth noting that, while the count-line of the top-left diagram of Figure 5 corresponding to the Pois(μ W ) model remains within the quantile bands b 0.025 and b 0.975 , it does hint that somewhat more zeros and considerably less 1's than expected under the Poisson model are present in the observed data.This pattern drawn by the count-line is a characteristic of zero-inflation, which may lead the analyst to speculate that the data is in fact Table 7: Frequency of dicentric chromosomes after acute homogeneous in vitro exposure to doses between 0 and 4.5Gy of Cobalt-60 γ-rays.(This corresponds to data set A1 in the notation of Oliveira et al. (2016), where also the reference for the data source is provided.) Frequency of counts dose 0 1 2 3 4 5 0.00 2591 1 0 0 0 0 0.25 2185 8 0 0 0 0 0.75 2550 44 1 0 0 0 1.00 2231 54 2 0 0 0 1.50 1712 96 3 0 0 0 2.50 1196 123 7 1 0 0 3.00 1070 320 41 6 1 0 4.50 895 360 110 25 5 1 zero-inflated, and hence to proceed with fitting a ZIP model which then delivers the nicely behaved Quantile Band plot as in the top right panel.

Biodosimetry data
We consider data consisting of n = 14430 chromosome aberration counts previously studied by Oliveira et al. (2016).The covariate dose, with values between 0 and 4.5Gy, gives the radiation dose applied in vitro to blood sample cells, causing DNA damage in form of doublestrand breaks.When incorrectly repaired by the cellular DNA-damage response mechanism, this can lead to dicentric chromosomes which can be counted under a microscope.That is, each examined blood sample cell contributes, for known covariate dose, exactly one count observation.For this data set, the counts take values in the range from 0 to 5. Data of this type have been fitted traditionally through Poisson regression models, though deviation from the Poisson property, and specifically the presence of excess zero counts, has been regularly reported in the literature, see e.g.Puig and Barquineiro (2011).
Table 7 displays the data under investigation, and Figure 6 contains the Quantile Band plots obtained when Poisson and zero-inflated Poisson models, using a log-link and quadratic polynomial for dose, but constant zero-inflation parameter in the ZIP case, are fitted to these data.The left hand plot clearly indicates the unsuitability of the Poisson model, whereas the right hand plot indicates that ZIP is suitable.Oliveira et al. (2016) carried out an extensive analysis of this data set, applying several statistical tests and model selection criteria in order to decide for an adequate modelling strategy.Specifically, they found that a Negative Binomial type 2 (hereafter NB2) model returned the lowest AIC (7489.1),closely followed by a ZIP model (AIC=7490.4).Other models considered included the Poisson as reference model (AIC=7504.7),and a Poisson Inverse Gaussian (PIG; AIC=7495.2).
The two plots in Figure 7 corresponding to the NB2 and PIG models, respectively, illustrate cases where the adjusted observed data line, A(k), remains close to the centre line.For the NB2, all observations lie between the 43% and 57% quantiles of their respective Poisson-Binomial distribution.Hence, there is less random variation amongst observed counts than would be expected under the NB2 model, most likely indicating that the variance of the fitted model is inflated in order to accommodate the number of observed zeros.A similar effect is observed for the PIG model.In summary, these plots suggest that the ZIP model is the most adequate model for these data, deviating from what would be concluded by looking at a single-number model selection criterion such as AIC.3).For ease of comparisons, the horizontal axes for the one parameter Poisson and geometric models are drawn to the same scale, and the horizontal axes for the four two-parameter models are drawn to the same scale.q q q q q q −50 0 50 0 1 2 3 4 5

Poisson Model Truncated Sample Estimator
A(k)=N(k)−m(k) RESPONSE q q q q q q −50 0 50

Conclusion
Whilst the principal purpose of the Quantile Band plots presented here is, as indicated in the title of this paper, to assess the suitability of a count regression model to a given data set, not to determine the 'best' model, the plots are an extremely useful tool to help answer the question: 'what is the best model?'.It is sometimes overlooked that the purpose of fitting a model to a given sample of data is not to find the model that 'best fits' the sample, but to attempt to discover a model that is, to paraphrase George Box, probably incorrect but of use as a model for the data from which the sample was taken.Whilst likelihood-based methods such as log-likelihood, AIC, BIC and various other 'information criteria' have an enormous role to play in determining and estimating the coefficients of 'the most suitable model', the fact that they frequently disagree on the nature of that model illustrates that none of them are perfect arbitrators.The graphical tool presented in this paper is an alternative approach to assessing model fit, which may be used on its own, or in conjunction with other methods.Its unique feature is that it enables the user to determine whether the observed frequency of a given count in the data is compatible with that to be expected under a given distribution; as stated in Section 2 if the majority of the frequencies of observed counts are indicated as being compatible with the model under consideration (i.e. they lie within the mid-quantile bands), then the model is likely to be appropriate for the data.Of course, under this procedure several models may be deemed 'appropriate'.This is well illustrated by Figure 4 which indicates that all four of the two parameter models considered are appropriate.A strength of the plots is that if they indicate non-suitability of a model, they also indicate the nature of the unsuitability, for example the top two plots of Figure 4 show that Poisson or geometric models are unsuitable as they severely underestimate the amounts of zeros, and overestimate the amount of 1's.Exactly the same information is contained in the horizontal and vertical forms of the plot.Raw Quantile Band plots (for example Figure 1) actually contain more information than the median-adjusted forms, but frequently are impractical (see Figure 9).
Quantile Band plots contain more information however than whether the observed counts are compatible with the stated model.As discussed in the various examples of Section 4 they indicate which values of k are over-or underestimated by the model; the example of Section 4.3 shows how the plots may be used to help determine the most suitable parameter for a model, which may differ from the maximum likelihood estimate.In the example of Section 4.4 the count-lines of the Quantile Band plots for the negative binomial and PIG models of Figure 7 indicate very little variation of the observed counts about the median values, possibly indicating that the counts exhibit less random variation than expected under these models, leading the researcher to speculate as to the reason for the possibly larger than necessary estimates of the dispersion parameter.
Whilst the authors advise the use of mid-quantiles as outlined in Section 3.1 for the construction of the quantile bands, other forms of quantiles may also be adopted.An attractive alternative may also be the use of expectiles, especially with view to their uniqueness property for discrete data (Eilers 2013).At present the plots are only proposed for use with univariate count models, but could be extended to continuous models by binning data, or to multivariate models by increasing the dimension of the plots.
R Code for the production of Quantile Band plots (in all three versions; that is with or without the starred items from Section 2.1) will be made available as supplementary material accompanying this publication.
partly supported by CRoNoS COST Action IC1408.The authors are grateful to a referee and the Editor for insightful comments.quantile functions of Bin(7,0.35)p q q q q q q q q q q q q q q q q q qq q q q q q q q traditional quantiles mid−quantiles R code for the computation of mid-quantiles from any probability mass function will be provided by the authors as supplementary material.

Density of Poisson model modified for digit preference
This probability mass function was used for generating data under digit preference in Section 2.2.

otherwise
Mid-quantile tables

Figure 2 :
Figure 2: Construction of Quantile Band plots for simulated covariate-free data.Subtracting the blue coloured median curve m(k) in the plot of Figure 1 from all other curves gives the horizontal version of the plot (top), which after rotation leads to the vertical Quantile Band plot (bottom).
When the role of the distribution G is taken by a Poisson model that incorporates digit preference (φ = 0.8), that is a modified Poisson model with probability density function g as given in Appendix C, one obtains the Quantile Band plot displayed in Figure3.We note that for k > 3 the quantile bands are wider at points corresponding to multiples of five, and narrower at other points, reflecting the modified densities of the Poisson model which incorporates digit preference.The count-line is now interior to the quantile bands at all points, indicating that b 0.025 (k) < A(k) < b 0.975 (k), for all values of k, and hence the suitability of the modified Poisson model.

Figure 3 :
Figure 3: Quantile Band plot for modified Poisson model (incorporating digit preference) The relative fits of Geometric, Poisson, Poisson Inverse Gaussian (PIG), Negative Binomial (Type I), Negative Binomial (Type II) and Zero-Inflated Poisson (ZIP) models are shown Table 3: Number of victims of murder known in past year,

Figure 4 :
Figure 4: Quantile Band Plots for the homicide data (Table3).For ease of comparisons, the horizontal axes for the one parameter Poisson and geometric models are drawn to the same scale, and the horizontal axes for the four two-parameter models are drawn to the same scale.

Table 1 :
Reported numbers of theatre visits

Table 2 :
Data of Table1with upper and lower quantiles for N (k) and A(k) (first 8 rows).

Table 4 :
Values of information criteria for several count regression models fitted to homicide data.

Table 8 :
Elements of the computation of mid-quantiles, for the binomial toy example