Distributions You Can Count On …But What’s the Point? ^†

Abstract

The Poisson regression model remains an important tool in the econometric analysis of count data. In a pioneering contribution to the econometric analysis of such models, Lung-Fei Lee presented a specification test for a Poisson model against a broad class of discrete distributions sometimes called the Katz family. Two members of this alternative class are the binomial and negative binomial distributions, which are commonly used with count data to allow for under- and over-dispersion, respectively. In this paper we explore the structure of other distributions within the class and their suitability as alternatives to the Poisson model. Potential difficulties with the Katz likelihood leads us to investigate a class of point optimal tests of the Poisson assumption against the alternative of over-dispersion in both the regression and intercept only cases. In a simulation study, we compare score tests of ‘Poisson-ness’ with various point optimal tests, based on the Katz family, and conclude that it is possible to choose a point optimal test which is better in the intercept only case, although the nuisance parameters arising in the regression case are problematic. One possible cause is poor choice of the point at which to optimize. Consequently, we explore the use of Hellinger distance to aid this choice. Ultimately we conclude that score tests remain the most practical approach to testing for over-dispersion in this context.

1. Introduction

The well-known Pearson family of continuous distributions, originally explored by Pearson (1895), is comprised of any solution to a particular differential equation. In his PhD thesis, Katz (1945) explored a family of discrete distributions that are solutions to a difference equation analogous to the Pearson differential equation.1 The Pearson family is a collection of four-parameter distributions and specializations thereof. Katz (1965, p. 175) observes that certain specializations ‘produce simpler and more manageable classes’ and restricts attention to a set of one- and two-parameter distributions. In particular, his restrictions result in a family of distributions that nest the two-parameter binomial and negative binomial (or Pascal) distributions, together with the one-parameter Poisson distribution.2 A defining characteristic of these distributions is that they arise when certain parameters, or parameter ratios, take integer values and so represent a set of measure zero in respect of the set of family members, which are defined in terms of real-valued parameters. The Katz family of distributions has proved important in the analysis of count data. It provides a framework within which practitioners can extend simple Poisson models to models that allow for individual heterogeneity, using the Poisson regression model (PRM). The PRM can, in turn, be extended to models that allow for either over-dispersion, using the negative binomial regression model (NBRM), or under-dispersion, using the binomial regression model. We shall, for the most part, defer consideration of under-dispersion to another time.

The problems of modelling and testing for over-dispersion have proved important in the count data literature. Essentially concurrently, papers by Cameron and Trivedi (1986); Lee (1986) and Lawless (1987a, 1987b) made substantial contributions to the literature on inference in the PRM, the NBRM, and testing for over-dispersion, with both Cameron and Trivedi (1986) and Lee (1986), in particular, couching substantial parts of their analysis within the context of the Katz family of distributions. This class of distributions is interesting because the binomial and negative binomial distributions are alternative specifications to the Poisson that allow under- and over-dispersion, respectively. Subsequent contributions to this literature include Dean (1992); Dean and Lawless (1989); Qu et al. (1990) and Fang (2003). Collectively they have explored likelihood ratio (LR), Lagrange multiplier (LM), and Wald tests, together with tests based on generalised method of moments (GMM) estimators, for over-dispersion in the PRM. Fang (2003) concludes that his preferred GMM test is that based on the fewest over-identifying assumptions offering essentially the same power as tests based on more over-identifying restrictions but having the greatest ease of calculation.3 Interestingly, this preferred test is that originally proposed by Katz (1965), on an ad hoc basis.

In this paper, we investigate a new family of tests for over-dispersion in the PRM by exploring point optimal tests where the alternative hypothesis lies in the Katz family of distributions. An analysis of the Katz likelihood reveals that maximum likelihood estimation may be problematic in the over-dispersed case, suggesting that the use of point optimal tests may have value. For overviews of the use of point optimal tests in econometrics see King (1987) and King and Sriananthakumar (2015). To the best of our knowledge they have not previously been used in the context of testing for over-dispersion.

This paper can be thought of as being comprised of three main parts. The first part provides a very brief description of the family of distributions introduced by Katz (1965), the second explores the role that these distributions can play in extending the PRM to allow for over-dispersion and, finally, we introduce a new class of point-optimal tests for over-dispersion. Specifically, in Section 2 we explore the Katz family of distributions, although most of the analysis is relegated to the Appendix A while Section 3 explores the PRM and NBRM. In particular, we highlight that the typical treatments of the NBRM really have little to do with what might be thought of as the canonical negative binomial distribution. Section 4 then focuses on the problem of testing for over-dispersion and the structure of the Katz likelihood. It is here that we introduce our family of point optimal tests and explore their small sample characteristics relative to some existing tests via a simulation study. We find that it is possible to choose a point optimal test which is better in the intercept only case, although the regression case proves problematic. One possible source of weakness in our point optimal tests is the choice of ‘point’ at which to optimize. In Section 5 we explore the use of Hellinger distance as a device to assist in choice of point. Although exact calculation of the Hellinger distance in this context is not tractable, it is straight-forward to obtain bounds on the distance. Using the upper bound, we find that the implied optimal points are extremely close to zero, implying that use of the score test is close to the optimal strategy in this context and so our advice to practitioners is to continue to use score tests to test for over-dispersion in this context. Section 6 concludes.

2. The Katz Family of Distributions

Among his many and varied interests, Karl Pearson was concerned with the problem of modelling (possibly) asymmetric empirical distributions. To this end, he developed a four-parameter family of skewed continuous distributions as solutions to a particular differential equation (Pearson 1895). The idea being that the distributions might be fitted to any data set using the method of moments approach that he had developed earlier (Pearson 1894). Perhaps surprisingly, the motivation for the choice of differential equation came from a difference equation that could be used to generate the hypergeometric distribution (Pearson 1895, pp.360–361); that is, from a discrete distribution. If we let $p\left(y\right)\equiv P[Y=y]$ denote the probability that the discrete random variable Y takes a value $y\in \mathbb{Y}$, where $\mathbb{Y}$ denotes the support of the distribution of Y, then the form of this difference equation was

$$\frac{p\left(y\right)-p(y-1)}{p(y-1)}=\frac{a-y}{b_0+b_{1}y+b_2{y}^2},$$

(1)

with a, $b_0$, $b_1$, and $b_2$ denoting the various parameters of the distribution. We note that this expression is of the form

$$p\left(y\right)/p(y-1)=P\left(y\right)/Q\left(y\right),$$

where P and Q are polynomials in y, and remark in passing that the sequence of probabilities so-defined are hypergeometric in that the ratio of adjacent terms in the sequence can be expressed as a ratio of polynomials in the index y.

Pearson did not pursue a discrete analogue to his family of distributions. Indeed, apart from some incidental investigations along these lines, Katz (1945) provided the first detailed analysis of the family of distributions arising from (1) although, apart from some abstracts (Katz 1946, 1948), it was not until Katz (1965) that this material was published. In the event, Katz (1965) focussed on a two-parameter special case of (1),4 which he expressed in the form

$$\frac{p(y+1)}{p\left(y\right)}=\frac{\lambda +y\gamma }{y+1},\lambda \in \mathsf{\Lambda }\subseteq \mathbb{R},\gamma \in \mathsf{\Gamma }\subseteq \mathbb{R},$$

(2)

with $y\in \mathbb{Y}\subseteq {\mathbb{Z}}_0$, where ${\mathbb{Z}}_0$ denotes the set of non-negative integers and subject to the usual axiomatic properties of probability:

$$\begin{array}{c}0\le p\left(y\right)\le 1\end{array}$$

(3)

$$\begin{array}{c}\sum _{y\in \mathbb{Y}}p\left(y\right)=1.\end{array}$$

(4)

As we demonstrate in the Appendix A, there are circumstances where both $\mathsf{\Lambda }$ and $\mathsf{\Gamma }$ may include values that are positive, negative, or zero.

Although Katz himself included zero in $\mathcal{Y}$, subsequent literature has not always done so, choosing instead to focus on the difference equation (2), whilst still referring to the resulting distributions as members of the Katz family; see, for example, Sundt and Jewell (1981); Willmot (1988) and Miller (1998). We too shall proceed in this latter manner, focussing on the what we call left-truncated Katz distributions, that include the original definition of Katz (1965) as the special case of no left-truncation. We relegate the technical analysis of these distributions to the Appendix A, which also gathers a number of other properties of this family of distributions.

Two members of this family that will be of particular interest to us are the Poisson and negative binomial distributions, which are commonly encountered in the modelling of counts and the possibility of over-dispersion. The probability mass functions (pmfs) of these distributions are the form:

$$p\left(y\right)=\left\{\begin{array}{cc}{\displaystyle \frac{{\lambda }^y{e}^{-\lambda }}{y!},}\hfill & \gamma =0,\hfill \\ {\displaystyle \frac{{(1-\gamma )}^{\lambda /\gamma }{\gamma }^y{\left(\frac{\lambda }{\gamma }\right)}_y}{y!},}\hfill & 0<\gamma <1,\hfill \end{array}\right.y\in \{0,1,2,\dots \},$$

(5)

respectively.5 Evidently, when $\gamma =0$, $p\left(y\right)$ is the pmf of the Poisson distribution. When $\gamma >0$, if $\lambda /\gamma =r$ is integer then $p\left(y\right)$ yields a standard representation of the negative binomial pmf, where the probability of success in any given trial is $\pi =1-\gamma$. Even if $\lambda /\gamma =\tau$ is not integer, $p\left(y\right)$ is still the pmf of a negative binomial distribution — see, for example, (9) — although the interpretation of $\lambda /\gamma$ differs between the two cases.6 We shall, hereafter, denote the Poisson distribution with parameter $\lambda$, $\mathcal{P}\left(\lambda \right)$, and the negative binomial with parameters $\tau$ and $\pi$, $\mathcal{NB}(\tau ,\pi )$.

Before moving on, let us consider the well-known Poisson approximation to the negative binomial. A common statement of this result is $\mathcal{NB}(\tau ,\pi )\to \mathcal{P}\left(\lambda \right)$ as $\tau \to \infty$ provided $\lambda =\tau (1-\pi )$ remains fixed. That is, $\pi \to 1$ at the same rate as $\tau$ diverges. One advantage of the parameterization adopted in (5) is that the somewhat convoluted requirement on how the parameters evolve in the approximation readily reduces to $\gamma \to 0^{+}$ for fixed $\lambda$.7

3. The Poisson, Negative Binomial, and Katz Regression Models

3.1. The Poisson Regression Model

The PRM extends the Poisson distribution to allow for individual heterogeneity. It has played an important role in the analysis of count data in both econometrics and statistics — early references include Gart (1964); Jorgenson (1961), and Haight (1967, chp. 5) — and is readily available in standard software such as MATLABmathsizesmall, Stata, and R. The use of the PRM in econometrics became increasingly widespread following the significant contributions of Gilbert (1979, 1982) and Hausman et al. (1984). Recent summaries can be found in Greene (2007); Winkelmann (2008) and Cameron and Trivedi (2013).

The PRM is obtained from the Poisson distribution by replacing the fixed parameter $\lambda$ with a function, denoted ${\lambda }_i$ say, of the k-vector of characteristics $x_i$ that can vary across individuals. Specifically, in the language of generalized linear models (GLIMs), we have the link function

$$ln{\lambda }_i=x_i^{\top }\beta ,$$

(6)

with regression coefficients $\beta$. The work of Nelder and Wedderburn (1972) and Frome et al. (1973) shows how iterated least squares methods can be used to obtain maximum likelihood estimates of $\beta$; see also McCullagh and Nelder (1989).

One shortcoming of the PRM is the implied equality of mean and variance that is characteristic of the Poisson distribution. Specifically, on replacing $\lambda$ with ${\lambda }_i$ in (Appendix A.2.2), we obtain8

$$\mathrm{E}\left[Y_i\mid {\lambda }_i\right]=\mathrm{V}\left[Y_i\mid {\lambda }_i\right]={\lambda }_i.$$

(7)

This is at odds with the observation that variability typically exceeds location in real world data, a feature known as over-dispersion. A common response to concerns about over-dispersion has been to explore extensions to the Poisson model that allow for different means and variances. To the extent that the Poisson regression model can be nested in such generalizations, this approach provides a framework within which one might test for either over-dispersion or underdispersion, although we will not explore this latter case here.

The fundamental characteristic of the PRM is that it is a function of the linear index, $x_i^{\top }\beta$, only through the ‘parameter’ ${\lambda }_i$, as per (7). In the next two sub-sections we will consider different extensions to this model, the first being the classical NBRM and the second being what we dub the Katz regression model (KRM). Both models extend the PRM by nesting it within a richer model with an additional ‘parameter’. An important distinction between the NBRM and the KRM is the role of ${\lambda }_i$. In the case of the NBRM, ${\lambda }_i$ remains the conditional mean of the count $Y_i$, whereas this is not the case in the KRM. A second distinction between the models is that the additional ‘parameter’ is typically treated as being a function of the linear index in the NBRM whereas in our treatment of the KRM it is not, it is a genuine parameter, although it is easy to envisage extensions where that requirement is relaxed.

3.2. The Classical Negative Binomial Regression Model

There are numerous paths leading to what might reasonably be called a negative binomial regression model.9 This is due, at least in part, to the variety of ways in which one might generate a negative binomial distribution. For example, Boswell and Patil (1970) provide 15 different derivations and, of course, there is a variety of parameterizations of the negative binomial distribution that can also lead to differences. Below we explore a fairly commonly adopted approach and consider some of its implications.

Our starting point is the following observation, originally due to Greenwood and Yule (1920). Suppose that $Y\mid \theta \sim \mathcal{P}\left(\theta \right)$, where $\theta$ is a random variable whose distribution is gamma with shape ($\tau$) and rate ($\eta$) parameters, written $\theta \sim \mathcal{G}(\eta ,\tau )$, so that the corresponding density function is,10

$$g(\theta ;\eta ,\tau )={\eta }^{\tau }{\theta }^{\tau -1}exp\{-\theta \eta \}/\mathsf{\Gamma }\left(\tau \right),\tau >0,\eta >0,\theta >0,$$

(8)

with $\mathrm{E}\left[\theta \right]=\tau /\eta$ and $\mathrm{V}\left[\theta \right]=\tau /{\eta }^2$.11 Then, we obtain an unconditional distribution for Y on averaging with respect to $\theta >0$, so that

$$\begin{array}{cc}\hfill \mathrm{Prob}\left(Y=y\mid \eta ,\tau \right)& ={\int }_{\theta >0}f(y\mid \theta )g(\theta ;\eta ,\tau )\mathrm{d}\theta \hfill \\ \hfill & =\frac{\mathsf{\Gamma }(y+\tau )}{y!\mathsf{\Gamma }\left(\tau \right)}{\left(\frac{1}{1+\eta }\right)}^y{\left(\frac{\eta }{1+\eta }\right)}^{\tau }.\hfill \end{array}$$

(9)

If one imposes the restriction $\tau \equiv r\in {\mathbb{N}}^{+}$, where ${\mathbb{N}}^{+}$ denotes the set of positive integers (or the natural numbers), then this is simply a form of the negative binomial (Pascal) pmf, with $\pi =\eta /(1+\eta )$, see (A5). Note that

$$\begin{array}{cc}\hfill \mathrm{E}\left[Y\right]& =\tau /\eta =\lambda \left(\mathrm{say}\right),\hfill \end{array}$$

the same as for the gamma distribution (8), and that

$$\begin{array}{cc}\hfill \mathrm{V}\left[Y\right]& =\frac{\tau }{\eta }+\frac{\tau }{{\eta }^2}=\mathrm{E}\left[Y\right]+{\tau }^{-1}{\left(\mathrm{E}\left[Y\right]\right)}^2=\lambda +{\tau }^{-1}{\lambda }^2.\hfill \end{array}$$

One posible path to a NBRM is to extend the analysis of Greenwood and Yule (1920) to allow for individual heterogeneity; we follow the treatment of Cameron and Trivedi (1986, p. 32). Specifically, we replace $\theta$ by ${\theta }_i$, where

$$ln{\theta }_i=x_i^{\top }\beta +{ϵ}_i,$$

(10)

with ${ϵ}_i$ a disturbance term reflecting unobservables. Cameron and Trivedi (1986) then assume that either ${ϵ}_i$, or ‘equivalently’ ${\theta }_i$, have a gamma distribution, conditional on the regressors. Their analysis then proceeds under the latter assumption, which is completely analogous to the developments of Greenwood and Yule (1920). Specifically, letting ${\theta }_i\mid x_i\sim \mathcal{G}({\eta }_i,{\tau }_i)$ yields

$$\mathrm{Prob}\left(Y_i=y\mid x_i;{\eta }_i,{\tau }_i\right)=\frac{\mathsf{\Gamma }(y+{\tau }_i)}{y!\mathsf{\Gamma }\left({\tau }_i\right)}{\left(\frac{1}{1+{\eta }_i}\right)}^y{\left(\frac{{\eta }_i}{1+{\eta }_i}\right)}^{{\tau }_i}.$$

(11)

Moreover,

$$\begin{array}{c}\mathrm{E}\left[Y_i\mid x_i;{\eta }_i,{\tau }_i\right]=\frac{{\tau }_i}{{\eta }_i}={\lambda }_i\left(\mathrm{say}\right),\end{array}$$

(12)

and

$$\begin{array}{c}\mathrm{V}\left[Y_i\mid x_i;{\eta }_i,{\tau }_i\right]=\frac{{\tau }_i}{{\eta }_i}+\frac{{\tau }_i}{{\eta }_i^2}={\lambda }_i+{\tau }_i^{-1}{\lambda }_i^2.\end{array}$$

(13)

It is immediately obvious that, in the final analysis, the functional form of (10) is a complete irrelevance, with only the parameters of the mixing Gamma distribution of any importance and we have made no assumptions about them beyond allowing the possibility of varying at the individual level. From here, Cameron and Trivedi (1986) argue that a variety of models are available on defining

$${\tau }_i={\alpha }^{-1}{(\mathrm{E}\left[Y_i\mid x_i\right];{\eta }_i,{\tau }_i)}^k$$

(14)

for $\alpha >0$ and arbitrary constant k, so that

$$\mathrm{V}\left[Y_i\mid x_i;{\eta }_i,{\tau }_i\right]=\mathrm{E}\left[Y_i\mid x_i;{\eta }_i,{\tau }_i\right]+\alpha {\left(\mathrm{E}\left[Y_i\mid x_i;{\eta }_i,{\tau }_i\right]\right)}^{2-k}.$$

Special cases of importance are then the Negbin I model (obtained when $k=1$) and the Negbin II model ($k=0$),12 of which the latter is probably the more popular in the literature. This model nests the PRM as a limiting case where $\alpha \to 0$ from above, a testable proposition, which is equivalent to ${\tau }_i$ diverging to ∞ for all i.

The specification (10) becomes more relevant if, instead, we assume that $h_i=e^{{ϵ}_i}\mid x_i\sim \mathcal{G}({\eta }_i,{\tau }_i)$ rather than ${\theta }_i$. Define ${\delta }_i=exp\left\{x_i^{\top }\beta \right\}$, so that ${\theta }_i={\delta }_i{h}_i$. Thus, conditional on $x_i$, ${\theta }_i$ is a scaled Gamma random variate which is, itself, a Gamma random variate. From the properties of the Gamma distribution we have immediately that ${\theta }_i|x_i\sim \mathcal{G}({\eta }_i/{\delta }_i,{\tau }_i)$. Moreover, analogs of results (11)–(13) are immediately available in this case on replacing ${\eta }_i$ by ${\eta }_i/{\delta }_i$. In short, the differing distributional assumptions are ‘equivalent’ in that the structure of the results is the same in both cases, however, it is only in this latter case that (10) has any relevance, through the presence of ${\delta }_i$ in the various expressions.

The attraction of the formulation (11)–(13) of Cameron and Trivedi (1986) is its close resemblance to a GLIM, which simplifies estimation.13 As noted above, the null of a PRM obtains as ${\tau }_i\to \infty$, however, results in a relatively odd PRM with a potentially unbounded mean, unless ${\eta }_i$ is diverging to infinity at the same rate as is ${\tau }_i$. Moreover, there is a Davies-type problem relating to the separate identification of both $\tau$ and $\eta$ when the null is true. Greene (2008) camouflages this difference by imposing the restriction $\eta =\tau$.14 He refers to this restriction as being mean preserving, by which is meant that when $\eta =\tau$, $\mathrm{E}\left[Y_i\mid x_i;\eta ,\tau \right]={\delta }_i$, as would be the case in the PRM. We should note that, in order to generate the same class of models as do Cameron and Trivedi (1986), Greene (2008) also allows $\tau$ to be replaced by ${\tau }_i$, as defined by (14), but this means that $\eta$ must be replaced by ${\tau }_i$ too.15 Of course, the restriction that ${\tau }_i={\eta }_i$ reduces the two parameter mixing gamma distribution to a single parameter distribution with the loss of modelling flexibility that implies. However, without this restriction, the conditional mean of $Y_i$ is other than ${\delta }_i$.

3.3. The Katz Regression Model

The fundamental difference between the NBRM, as described in the above, and what we refer to as the Katz regression model (KRM) lies in the generation of the underlying distribution. Specifically, the Katz family of distributions is not generated via a mixing argument and so, in contrast to the NBRM, the probabilistic quantities of interest (pmfs and moments) are not functions of the parameters of the mixing distribution; see (11)–(13). In this sense, the parameterization of the Katz family is more natural than that of the NBRM. Directly analogously with the PRM, the KRM can be generated from (A6) simply by replacing $\lambda$ by ${\lambda }_i$, as per (6), which is analogous to our earlier development of the PRM.16 Equally, one might explore models that see $\gamma$ replaced by functions of regressors, ${\gamma }_i$ say, although we will not. Note that the conditional mean and variance of this distribution are given by (Appendix A.2.2), with $\lambda$ replaced by ${\lambda }_i$. Contrast this structure with that for the NBRM described above. There we saw that the conditional mean of the dependent variable was not varying with $\gamma$, being a function of the linear index $x_i^{\top }\beta$ alone. Similarly, by construction, the variance exceeded the mean of the dependent variable, but the reduction to a Poisson model requires the shape parameter $\tau$ of the mixing gamma distribution to be unbounded, which yields a degenerate distribution for given rate parameter $\eta$.

It is clear that it is not necessarily desirable to preserve the mean, in Greene’s sense of equating ${\tau }_i$ and ${\eta }_i$(Greene 2008), because, as $\gamma$ increases, the mean for both the NBRM and the KRM should be decreasing relative to that of the PRM.

We note in passing that this is the model that underlies the generalized event count (GEC) model of King (1989); this model was also considered by Ghahfarokhi et al. (2008). That they obtained more complicated models than that proposed here, resulting in the models being less popular than the NBRM in practice, stems from the fact that they did not have (A6) as the pmf implied by (2), which in part is due to working with (2) rather than (A1).

4. Testing for Over-dispersion in Poisson Regression model

There is a vast literature addressing the problem of over-dispersion and how to test for it. We will not attempt to provide a comprehensive survey of this literature, focussing instead on a few key contributions, although it should be noted that most of the references cited so far will have some discussion of the problem. We shall break our discussion into two parts. First, we shall restrict attention to the case where the only regressor in the model is an intercept, so that ${\lambda }_i=\lambda$ is a constant. Then we will extend the analysis to allow for additional regressors. In each case, the null hypothesis will be that the data have been generated by the PRM. We investigate the performance of tests whose preferred alternative is, variously, that the data have come from one of the Negbin I, Negbin II, or Katz regression models.17

4.1. The Katz Likelihood

For any positive real number $n=\lambda /\gamma$, with $\lambda >0$ and $0<\gamma <1$, the Katz pmf is given by

$$\mathrm{Prob}\left(Y=k\right)=n(n+1)(n+2)\dots (n+k-1)\frac{{\gamma }^k{(1-\gamma )}^n}{k!}$$

and hence

$$\begin{array}{cc}\hfill L\left(y_1,\dots ,y_N\mid \lambda ,\gamma \right)& =\prod _{i=1}^N\left[n(n+1)(n+2)\dots (n+y_i-1)\frac{{\gamma }^{y_i}{(1-\gamma )}^n}{y_i!}\right]\hfill \\ \hfill & =\left[\prod _{i=1}^N\prod _{s=0}^{y_i-1}(n+s)\right]{(1-\gamma )}^{Nn}{\gamma }^{{\sum }_{i=1}^N{y}_i}{\left[\prod _{i=1}^N{y}_i!\right]}^{-1},\hfill \end{array}$$

(15)

where products of the form ${\prod }_{s=0}^{y_i-1}=1$ for $y_i=0$. In textbook cases, where n is known, the first and last terms in (15) are functions of the data only and ${\sum }_{i=1}^N{y}_i$ is sufficient by the factorization theorem. But in the current context, when n is not known, there is no reduction to a fixed dimensional sufficient statistic. Even if $\lambda$ is known there is no sufficiency reduction; the ratio n is required. Only the entire sample (or the order statistics) are sufficient but even they are are not complete so that different parameter configurations may give rise to the same data. We may surmise this from the likelihood (15) since any combination of $\lambda$ and $\gamma$ that preserves n gives the same likelihood. Adopting the convention that ${\sum }_{s=0}^{y_i-1}log(n+s)=0$ when $y_i=0$, the log likelihood is

$$\begin{array}{cc}\hfill logL& =\sum _{i=1}^N\sum _{s=0}^{y_i-1}log\left(n+s\right)+\sum _{i=1}^{N}log\left(\frac{{\gamma }^{y_i}}{y_i!}\right)+\sum _{i=1}^{N}log\left({\left(1-\gamma \right)}^n\right)\hfill \\ \hfill & =\sum _{i=1}^N\sum _{s=0}^{y_i-1}log\left(\lambda +s\gamma \right)+Nlog{\left(1-\gamma \right)}^{\frac{\lambda }{\gamma }}-\sum _{i=1}^{N}log\left(y_i!\right),\hfill \end{array}$$

(16)

where the second line follows by substituting for n and simplifying.

In fact, (nonlinear) maximum likelihood estimators for n do not exist when the sample variance is less than the mean, that is, $s_y^2\le \overline{y}$ (see Al-Khasawneh (2010) and the references therein). Note that, even when drawing from a negative binomial, which by definition is over-dispersed, many individual samples will exhibit under-dispersion. We can see the difficulty explicitly by looking at simple moment estimates for $\lambda$ and $\gamma$, that is, solve (Appendix A.2.2) to get

$$\begin{array}{cc}\hfill \widehat{\lambda }& =\overline{y}\left(1-\widehat{\gamma }\right)\hfill \\ \hfill \widehat{\gamma }& =1-\left(\overline{y}/s_y^2\right).\hfill \end{array}$$

Hence, problems arise when $s_y^2\le \overline{y}$ since then $\widehat{\gamma }\le 0$, which is illegtimate when investigating over-dispersion. Even if $s_y^2\ge \overline{y}$, convergence issues arise if the difference is not great. This suggests that test procedures that use maximum likelihood estimates, such as Wald or Likelihood Ratio tests, can be problematic and that there may be a role for point optimal approaches.

4.2. Point Optimal Tests

Point optimal tests have had a long and varied career in econometrics; see King (1987) and King and Sriananthakumar (2015) for an overview. These tests optimize power at a particular parameter value under the alternative, the idea being to have good power at a point where incorrectly accepting the null really matters. This is in contrast to, say, a score test that is locally best, in that it has the steepest power function local to the null hypothesis. Although not an undesireable property in any way, the practical difference between a null model and some other model local to the null is often, although not always, vanishingly small. So, optimizing the ability to distiguish between such null model and another local to it is not necessarily all that desireable a property. Moreover, there is implicit in such an approach the notion that the power function will be monotonically increasing, which ideally it should be, and that it will remain near the power envelope as the data generating process diverges from the null. In many cases this is indeed what happens, although we know that power functions are likely to cross, as otherwise the test would be uniformly most powerful, which is a very rare property indeed. The divergence between the power function of a score test and the power envelope is then something that requires exploration on a case by case basis and we will explore this below.

The log likelihood ratio of the $Katz-NB$ alternative to the Poisson (P) null is written

$$\begin{array}{cc}\hfill LLR\left({\lambda }_1,{\lambda }_0,\gamma \right)& =\sum _{i=1}^N\sum _{s=0}^{y_i-1}log\left({\lambda }_1+s\gamma \right)+Nlog{\left(1-\gamma \right)}^{\frac{{\lambda }_1}{\gamma }}\hfill \\ \hfill & -\left(log{\lambda }_0\sum _{i=1}^N{y}_i-N{\lambda }_0\right).\hfill \end{array}$$

Assuming that both distributions are fully specified (with $\lambda ={\lambda }_0={\lambda }_1$), the Neyman-Pearson Lemma states that the $UMP$ test of $\gamma =0$ versus $\gamma ={\gamma }_1$ is given by $LLR\left(\lambda ,\lambda ,{\gamma }_1\right)$. Hence assuming $\lambda$ known, the so-called power envelope is determined by computing $LLR\left(\lambda ,\lambda ,\gamma \right)$ over a range of values of $\gamma \in \left(0,1\right)$. A $PO$ test is constructed by choosing a fixed $\gamma ={\gamma }_{PO}$ to be a ‘representative’ value under the alternative $Katz-NB$ distribution, giving $LLR\left(\lambda ,\lambda ,{\gamma }_{PO}\right)$. It is desirable that ${\gamma }_{PO}$ be chosen so that the power of the test $LLR\left(\lambda ,\lambda ,{\gamma }_{PO}\right)$ is as close as possible to the power of the family of tests $LLR\left(\lambda ,\lambda ,\gamma \right)$, $\gamma \in \left(0,1\right)$. That is, ideally, ${\gamma }_{PO}$ is chosen so that the power function of the resulting test is as close to the power envelope as possible.

4.3. Score Test

A common alternative to likelihood ratio approaches, which does not require maximum likelihood estimation of the parameter of interest, is to construct optimal tests local to the null $\gamma =0$. The so-called eficient score tests, or simply score test, are derived by differentiating the log likelihood with respect to $\gamma$ and then setting $\gamma =0$. Such score tests are easily found for the Katz family using (16); see, for example, Katz (1965) and Lee (1986). Specifically, the score test is

$$S\left(\lambda \right)=\frac{1}{2\lambda }\sum _{i=1}^N\left[y_i\left(y_i-1\right)-{\lambda }^2\right].$$

(17)

This test was originally proposed by Katz (1965) on heuristic grounds and by Lee (1986) as a formal score test.18^,19

4.4. Simulation Experiments

4.4.1. The Unconditional Model

In this section we simulate the powers of the point optimal and score tests and compare them to the benchmark power envelope. First we give details for the power envelope and this is followed by a description of the operational tests. We present results for a sample size of $N=50$ throughout.

In Figure 1, Figure 2 and Figure 3 below, we consider a range of values for $\lambda \in (0,8)$ and $\gamma \in (0,1)$ and look at the relative performance of the $PO$ and S tests and the power envelope. For each $(\lambda ,\gamma )$ pair we generate samples from the Katz family, $K(\lambda ,\gamma )$. This is efficiently accomplished using $p\left(0\right)={(1-\gamma )}^{\frac{\lambda }{\gamma }}$ along with the defining recurrence $p(y+1)=(\lambda +\gamma y)/(y+1)p\left(y\right)$ and then sampling from the inverse cumulative distribution function. Setting $\gamma =ϵ$, for some very small positive $ϵ$, effectively simulates from the Poisson null. The null critical values are computed by simulation to avoid asymptotic approximations. This means that the sizes of tests are accurate (up to simulation error) and hence that the power comparisons are meaningful in smaller sample sizes.

For the power envelope, we simulate from the null $P\left(\lambda \right)=K(\lambda ,ϵ)$ distribution, compute 10,000 values of $LLR(\lambda ,\lambda ,\gamma )$ and extract the $95\%$ quantile as a critical value, $cv$. The $DGP$ is the Katz family $K(\lambda ,\gamma )$ and we simulate $10,000$ replicates from the DGP and count the percentage of times the $PO$ statistic, $LLR(\lambda ,\lambda ,\gamma )$, exceeded the $cv$ to calculate the power envelope.

The operational $PO$ test estimates ${\lambda }_0$, ${\lambda }_1$ and fixes $\gamma$ at ${\gamma }_{PO}$. To do this, we use simple moment based estimators, that is, compute ${\widehat{\lambda }}_{NB}={\widehat{\lambda }}_P(1-{\widehat{\gamma }}_{NB})$, where ${\widehat{\lambda }}_P=\overline{y}$ and ${\widehat{\gamma }}_{NB}=1-\left(\overline{y}/s_y^2\right)$ using the mean and variance of the data at hand. Should ${\widehat{\gamma }}_{NB}$ stray negative, we truncate and set ${\widehat{\gamma }}_{NB}=ϵ$. The $PO$ test is computed as $LLR\left({\widehat{\lambda }}_P,{\widehat{\lambda }}_{NB},{\gamma }_{PO}\right)$ while the score test, using (17), is $S\left({\widehat{\lambda }}_P\right)$. The null is Poisson $P\left(\lambda \right)=K(\lambda ,ϵ)$ and we simulate from the null, computing $LLR\left({\widehat{\lambda }}_P,{\widehat{\lambda }}_{NB},{\gamma }_{PO}\right)$ and $S\left({\widehat{\lambda }}_P\right)$ for each realization, to get $5\%$ $cv$’s. We calculate the power by simulating from the DGP $K(\lambda ,\gamma )$.

It is helpful to view Figure 1, Figure 2 and Figure 3 in the light of the cross sections displayed in Figure 4. It is clear that, for small values of $\gamma$, no test can be expected to perform well close to the null. Equally, for large values of $\gamma$, with high degrees of over-dispersion, all reasonable tests can be expected to be powerful. Thus, for small and large degrees of over-dispersion we expect to see little difference in the performance of the envelope and the $PO$ test as the left panel of Figure 1 attests. For moderate values of $\gamma$, the power of the $PO$ test can be significantly smaller than the envelope as the coloured scale suggests. In the left panel the difference in power between the $PO$ and score tests is plotted for ${\gamma }_{PO}=0.1$. The differences are not large but the $PO$ test uniformly dominates as the scale indicates. Figure 2 plots the same surfaces for ${\gamma }_{PO}=0.3$ and the interesting feature is that both tests perform similarly but none dominates the other. In Figure 3, ${\gamma }_{PO}=0.5$ and the score test dominates by a small margin except where the $DGP$ corresponds to ${\gamma }_{PO}$.

Since the shapes of the surfaces are quite smooth over $\lambda$, we plot a cross-section at $\lambda =5$ to look at absolute performance. There are three panels in Figure 4 each corresponding to a value of ${\gamma }_{PO}=0.1$, $0.3$, $0.5$. The power envelope is shown in red, with those of the PO test in blue and the Score test in orange. Also shown (vertically) is the PO point ${\gamma }_{PO}$.

The power envelope reaches unity at around $\gamma =0.2$. This corresponds to a degree of over-dispersion, in the $Katz-NB$ distribution, of ${\sigma }^2/\mu =1/\left(1-\gamma \right)=1.25$. For the tests to reach equivalent power requires $\gamma =0.6$ with ${\sigma }^2/\mu =2.5$, roughly, and $\gamma =0.9$ is required at ${\sigma }^2/\mu =10$. So, neither test can match the envelope unless the degree of over-dispersion is quite large. For ${\gamma }_{PO}$ less than $0.3$ the $PO$ test performs better uniformly, at $0.3$ they perform equally well and for ${\gamma }_{PO}>0.3$ the score test is better. Thus, a choice of ${\gamma }_{PO}$ which is small will uniformly dominate.

4.4.2. The Katz Regression

In practice, the analysis of over-dispersion often takes place when covariates need to be taken into account. As explained in Section 3 there are many ways in which this may be approached. We work directly from the definition of the Katz family rather than mix over a kernel Poisson distribution. The log of the likelihood takes the form

$$LL=\sum _{i=1}^N\sum _{s=0}^{y_i-1}log({\lambda }_i+s\gamma )+\sum _{i=1}^{N}log{(1-\gamma )}^{\frac{{\lambda }_i}{\gamma }},$$

with ${\lambda }_i=exp\left({\beta }_1{x}_i\right)$ varying and $\gamma$ fixed. This gives $\mathrm{E}\left[Y_i|x_i\right]={\mu }_i={\lambda }_i/(1-\gamma )$ and $\mathrm{V}\left[Y_i|x_i\right]={\lambda }_i/{(1-\gamma )}^2$.

The $PO$ test, using ${\gamma }_{PO}$, is based on the log likelihood ratio,

$$\begin{array}{cc}\hfill LLR({\lambda }_0,{\lambda }_1,{\gamma }_{PO})& =\sum _{i=1}^N\sum _{s=0}^{y_i-1}log({\lambda }_{1,i}+s{\gamma }_{PO})+\sum _{i=1}^N\frac{{\lambda }_{1,i}}{{\gamma }_{PO}}log(1-{\gamma }_{PO})\hfill \\ \hfill & -\left(\sum _{i=1}^N{y}_{i}log{\lambda }_{0i}-\sum _{i=1}^N{\lambda }_{0i}\right)\hfill \end{array}$$

and the $PO$ test needs to estimate the parameters ${\lambda }_0$ and ${\lambda }_1$. Estimating ${\lambda }_{1,i}=exp\left({\beta }_1{x}_i\right)$ may be problematic as trying to fit a $NB$ regression when the data is Poisson can lead to identification/convergence problems, exacerbated by the fact that fitting these types of regressions requires nonlinear maximum likelihood estimation. We avoided this issue in the last sub-section by using Katz moment estimators. Here, we use a regression version of the same idea. First, estimate a Poisson regression $P\left({\mu }_i^{*}\left(x_i\right)\right)$ (including a constant) which will return the mean estimate $exp\left(b_0+b_1{x}_i\right)$. Noting that we can write ${\mu }_i=exp\left({\beta }_1{x}_i\right)/\left(1-\gamma \right)=exp\left({\beta }_0+{\beta }_1{x}_i\right)$, where ${\beta }_0=-log\left(1-\gamma \right)$, we can set ${\widehat{\beta }}_1=b_1$ and hence ${\widehat{\lambda }}_{1i}=exp\left({\widehat{\beta }}_1{x}_i\right)$ to give the vector ${\widehat{\lambda }}_1$. To get ${\widehat{\lambda }}_{0,i}$ we fit the Poisson without the constant term which returns the estimate $exp\left(b_2{x}_i\right)$, which gives ${\widehat{\lambda }}_{0,i}=exp\left(b_2{x}_i\right)$, and hence $LLR\left({\widehat{\lambda }}_0,{\widehat{\lambda }}_1,{\gamma }_{PO}\right)$.

As a comparator to the $PO$ statistic, in the regression setting, we use the score test of Dean and Lawless (1989), which avoids the potential difficulties associated with maximum likelihood estimation. Thus, we estimate the Poisson regression $P\left(\mu \left(x_i\right)\right)$ (with a constant) to get the vector of predictors ${\widehat{\mu }}_i$ and, using $z_i={\left(y_i-{\widehat{\mu }}_i\right)}^2-y_i$, the test $S\left(\widehat{\mu }\right)$ is computed as the t-statistic in the regression of $z_i$ on ${\widehat{\mu }}_i$.20 Again critical values are computed by simulation. The null is generated as $P\left({\lambda }_i\right)$, with ${\lambda }_i=exp\left({\beta }_1{x}_i\right)$ used to keep the means of the counts low. We used $x_i{\sim }^{i}log\left(P\left(2\right)+1\right)$ and the $x_i$ are kept fixed under replication. We generate simulated critical values, based on $10,000$ replications, for the tests $S\left(\widehat{\mu }\right)$ and $LLR\left({\widehat{\lambda }}_0,{\widehat{\lambda }}_1,{\gamma }_{PO}\right)$. To compute powers, the DGP $K\left({\lambda }_i,\gamma \right)=NB(n_i,p)$ is used, where $n_i={\lambda }_i/\gamma$ and $\gamma$ takes a selection of values in $\left(0,1\right)$. As usual, $\pi =1-\gamma$. The results are presented in Figure 5.

The $PO$ test performs badly for very high degrees of over-dispersion when ${\gamma }_{PO}$ is less than $0.5$ approximately and is dominated by the score test for ${\gamma }_{PO}$ greater than $0.5$. However, the choice ${\gamma }_{PO}=0.5$ does lead to superior $PO$ performance albeit by not a great margin.

4.4.3. Summary

Our experimental results are mixed. In the unconditional model, the point optimal tests appeared to work best when $\gamma$ was small, with their performance deteriorating relative to the score test as $\gamma$ increased. In the regression model, the score test outperformed the point optimal tests suggesting that the null distribution of the score tests was more robust to the presence of nuisance parameters than was that of the point optimal tests, with none of the test statistics being pivotal. However, these rankings were also sensitive to the choice of point. This begs the question as to whether or not we are choosing the ‘point’ for the point optimal tests in a sensible way. It is to this question that we turn in the next section.

5. Hellinger Distance

The reasoning behind the use of point optimal tests is to put power where it is of greatest practical use. The immediate problem facing the use of point optimal tests is where to place the ‘point’. Sometimes the testing problem suggests a solution. Other times the choice is less clear and is often based an the outcome from a simulation study ‘run-off’, making the results somewhat ad hoc. The attraction of point optimal test in the context of testing for over-dispersion is that the parameter space of interest, namely that of $\gamma$, is bounded and so there is some hope of finding an appropriate point. One way of defining appropriate in this context is where the distribution under the alternative starts to depart from that under the null in some substantial way. The questions then reduces to one of how we might measure such a departure. In this section we explore the use of Hellinger distance $\left(\mathcal{H}\right)$ for this purpose. We think that this is a novel use of such a distance measure and is of independent interest. We do not, however, assert that Hellinger distance is the only choice or even the best choice in this context, but it does yield some interesting results.

To begin, various definitions of Hellinger distance are available.21 Originally proposed in an integral form by Hellinger (1909), we will work with the following discrete variant:

Definition 1

(Hellinger Distance for Discrete Random Variables). The squared Hellinger distance between these two discrete distributions P and Q is

$${\mathcal{H}}^2=\frac{1}{2}\sum _{j=1}^k{\left(\sqrt{p_j}-\sqrt{q_j}\right)}^2=1-\sum _{j=1}^k\sqrt{p_j{q}_j},$$

(18)

where $P=(p_1,\dots ,p_k)$ and $Q=(q_1,\dots ,q_k)$.

We note in passing that the Hellinger distance is bounded, $0\le \mathcal{H}\le 1\Rightarrow 0\le {\mathcal{H}}^2\le 1$. $\mathcal{H}=1$ iff P assigns zero probability to anywhere that Q assigns positive probability and $\mathcal{H}=0$ iff $P=Q$.

5.1. The Poisson Distribution

By way of example, to illustrate the basic idea and to help calibrate the procedure, suppose that we choose as our base case a Poisson distribution with parameter ${\lambda }_0$ so that the implied standard deviation is $\sqrt{{\lambda }_0}$. Writing $\mathrm{Prob}\left(X=x\mid \lambda \right)\equiv \mathcal{P}\left(\lambda \right)$, we are going to explore the behaviour of $\mathcal{H}$ as we compare $\mathcal{P}\left({\lambda }_0\right)$ with $\mathcal{P}\left({\lambda }_1\right)$ for various $({\lambda }_0,{\lambda }_1)$. When comparing Poisson distributions, the squared Hellinger distance is readily shown to

$${\mathcal{H}}^2=1-exp\left\{-\frac{1}{2}{\left(\sqrt{{\lambda }_1}-\sqrt{{\lambda }_0}\right)}^2\right\}.$$

(19)

Figure 6 provides some insight into the sensitivity of Poisson pmfs to changes in parameter values when the parameters are small and includes examples that are variously skewed to the right, (roughly) symmetric, and skewed to the left. Observe that, here we have used $\lambda =1$ as the base case and that, as $\lambda$ increases, it is by one standard deviation each time and so these changes are quite dramatic.

In Figure 7 we present values for $\mathcal{H}$ for various ${\lambda }_0$ and ${\lambda }_1$. The dashed and dotted lines correspond to Hellinger distances of 0.1 and 0.05, respectively.

We observe that $\mathcal{H}$ is asymmetric in $\lambda$ for all ${\lambda }_0$ considered, which reflects the skewed nature of Poisson distributions. Note that, as ${\lambda }_0$ increases, so too does standard deviation of the base distribution. As this happens a given value of $\mathcal{H}$ will admit great differences between ${\lambda }_0$ and ${\lambda }_1$. For example, when ${\lambda }_0=1$ a Hellinger distance of $0.1$ or greater is achieved for any $0.8\approx L\le {\lambda }_1\le U\approx 1.3$. In contrast, when ${\lambda }_0=9$, $\mathcal{H}\le 0.1$ for all $8.2\approx L\le {\lambda }_1\le U\approx 9.9$, which is a much wider interval than the previous case. Given that we are seeking to construct point optimal tests that compete with locally best tests, these results suggest that we need to be looking at points for which the Hellinger distance is quite small.

5.2. The Katz Distribution

We will take the Poisson ($\gamma =0$) as our base model. Moreover, as Poisson-ness, or otherwise, is completely determined by the value of $\gamma$, we will hold $\lambda$ fixed across models. Here the support under both null (equi-dispersion) and alternative (over-dispersion) is $y\in \mathbb{Y}=\{0,1,2,\dots \}$ and so

$$\begin{array}{cc}\hfill {\mathcal{H}}^2& =1-\sum _{y=0}^{\infty }\sqrt{\frac{e^{-\lambda }{\lambda }^y}{y!}\frac{{\left(\frac{\lambda }{\gamma }\right)}_y{\gamma }^y{(1-\gamma )}^{\lambda /\gamma }}{y!}}\hfill \\ \hfill & =1-{\left[e^{-1}{(1-\gamma )}^{1/\gamma }\right]}^{\lambda /2}\sum _{y=0}^{\infty }\frac{1}{y!}\sqrt{{\left(\frac{\lambda }{\gamma }\right)}_y{\left(\lambda \gamma \right)}^y}.\hfill \end{array}$$

Although not amenable to direct solution we notice that

$${\left(\frac{\lambda }{\gamma }\right)}_y{\gamma }^y={\gamma }^y\prod _{j=0}^{y-1}\left(\frac{\lambda }{\gamma }+j\right)=\prod _{j=0}^{y-1}\left(\lambda +j\gamma \right)>{\lambda }^y,\mathrm{for}\mathrm{all}0<\gamma <1.$$

(20)

Therefore,

$$\begin{array}{cc}\hfill {\mathcal{H}}^2& >1-{\left[e^{-1}{(1-\gamma )}^{1/\gamma }\right]}^{\lambda /2}\sum _{y=0}^{\infty }\frac{1}{y!}\sqrt{{\left[{\left(\frac{\lambda }{\gamma }\right)}_y{\gamma }^y\right]}^2}\hfill \\ \hfill & =1-e^{-\lambda /2}{(1-\gamma )}^{-\lambda /\left(2\gamma \right)}=h_L^2,\left(\mathrm{say}\right).\hfill \end{array}$$

(21)

We can solve this non-linear equation for ${\gamma }_L$ numerically for given $\lambda$ and $h_L$. Some results are reported in

Table 1. Values for ${\gamma }_L$ Obtained From (21) For Given $h_L$ and $\lambda$ (scaled by a factor of ${10}^{12}$).

${\mathit{h}}_{\mathit{L}}$∖$\mathit{\lambda }$	1	2	3	4	5
0.02	0.0348	0.0493	0.1261	0.1556	0.0147
0.04	0.0019	0.0121	0.0131	0.0250	0.0148
0.06	0.0032	0.0014	0.0115	0.0109	0.0136
0.08	0.0008	0.0033	0.0018	0.0061	0.0140
0.1	0.0027	0.0006	0.0022	0.0057	0.0072

.

We see that all values of ${\gamma }_L$ are positive, albeit extremely to zero. Alternatively, from (20) we also have the result

$$\begin{array}{cc}\hfill {\mathcal{H}}^2& <1-{\left[e^{-1}{(1-\gamma )}^{1/\gamma }\right]}^{\lambda /2}\sum _{y=0}^{\infty }\frac{{\lambda }^y}{y!}=1-e^{\lambda /2}{(1-\gamma )}^{\lambda /\left(2\gamma \right)}=h_U^2,\left(\mathrm{say}\right).\hfill \end{array}$$

(22)

Solutions to (22) for various $\lambda$ and $h_U$ are given in

Table 2. Values for ${\gamma }_U$ Obtained From (22) For Given $h_U$ and $\lambda$.

${\mathit{h}}_{\mathit{U}}\setminus \mathit{\lambda }$	1	2	3	4	5
0.02	0.0016	0.0008	0.0005	0.0004	0.0003
0.04	0.0064	0.0032	0.0021	0.0016	0.0013
0.06	0.0143	0.0072	0.0048	0.0036	0.0029
0.08	0.0252	0.0127	0.0085	0.0064	0.0051
0.1	0.0391	0.0198	0.0133	0.0100	0.0080

.

We see that ${\gamma }_U$ is monotonically increasing in $h_U$ but monotonically decreasing in $\lambda$. That is, once $\lambda$ becomes sufficiently large, even small departures of the hellinger distance from zero are consistent with $\gamma >0$.

All of the above said, however, the over-riding conclusion is that the optimal ‘points’ are going to be sufficiently close to zero that it is not clear that there is much benefit over just using the the score test, which is essentially point optimal at $\gamma =0$. The main reason for such a conclusion is our earlier results indicating that, in the regression context, the score test is much less subject to the influence of the regression coefficients, which are nuisance parameters in this testing problem.

6. Conclusions

At a fundamental level, this paper explores the use of point optimal tests in the problem of testing for over dispersion. Our basis of comparison is the score test of Lee (1986), which is the same as the earlier method of moments test proposed by Katz (1965). Our findings are somewhat disappointing and we are unable to recommend that practitioners change their current practices as the performance of the point optimal tests is, at best, mixed. It may be possible to improve the performance of the point optimal tests by a more refined analysis of (i) the problem of nuisance parameters and (ii) the construction of p-values, along the lines suggested by King and Sriananthakumar (2015). This we leave for further work.

Along the way, the paper has made two other contributions. First, in the Appendix A we have provided a reasonably exhaustive treatment of the family of distributions consistent with the difference equation of Katz (1965). To the best of our knowledge this treatment extends all known earlier results by allowing for arbitrary points of left truncation. This expands the class of distributions originally considered by Katz (1965), which can be characterized as including zero in the support of the count variable. The treatment is closest to that of Willmot (1988), although there are differences in the mode of analysis and he restricts attention to extensions where only zero is omitted from the support of the count variable. We note in passing that right truncation is a much easier problem to deal with as it neither expands nor contracts the members in the family, in the way that left-truncation does. Its only consequence is the introduction of a scale factor equal to $1-R$, where R denotes the upper tail probability that has been truncated.

The other contribution that we have made is to introduce the use of Hellinger distance as a metric by which one might settle on the ‘points’ characterizing point optimal tests. This is novel and allows a more systematic treatment than the grid searches that have characterized such choices in the past.