Econometric analysis of non-standard count data

(1)

by

Ryan T. Godwin BA, Queen’s University, 2006 MA, University of Saskatchewan, 2008 A Thesis Submitted in Partial Fulfillment

of the Requirements for the Degree of DOCTOR OF PHILOSOPHY in the Department of Economics

 Ryan T. Godwin, 2012 University of Victoria

(2)

Supervisory Committee

Econometric Analysis of Non-Standard Count Data by

Ryan T. Godwin

BAH, Queen’s University, 2006 MA, University of Saskatchewan, 2008

Supervisory Committee

Dr. David E.A. Giles, Department of Economics Supervisor

Dr. Judith A. Clarke, Department of Economics Departmental Member

Dr. Kenneth G. Stewart, Department of Economics Departmental Member

Dr. Min Tsao, Department of Mathematics and Statistics Outside Member

(3)

Abstract

Supervisory Committee

Dr. David E.A. Giles, Department of Economics Supervisor

Dr. Judith A. Clarke, Department of Economics Departmental Member

Dr. Kenneth G. Stewart, Department of Economics Departmental Member

Dr. Min Tsao, Department of Mathematics and Statistics Outside Member

This thesis discusses various issues in the estimation of models for count data. In the first part of the thesis, we derive an analytic expression for the bias of the maximum

likelihood estimator (MLE) of the parameter in a doubly-truncated Poisson distribution, which proves highly effective as a means of bias correction. We explore the

circumstances under which bias is likely to be problematic, and provide some indication of the statistical significance of the bias. Over a range of sample sizes, our method outperforms the alternative of bias correction via the parametric bootstrap. We show that MLEs obtained from sample sizes which elicit appreciable bias also have sampling distributions which are unsuited to be approximated by large-sample asymptotics, and bootstrapping confidence intervals around our bias-adjusted estimator is preferred, as two tiers of bootstrapping may incur a heavy computational burden.

Modelling count data where the counts are strictly positive is often accomplished using a positive Poisson distribution. Inspection of the data sometimes reveals an excess of ones, analogous to zero-inflation in a regular Poisson model. The latter situation has well developed methods for modelling and testing, such as the zero-inflated Poisson (ZIP) model, and a score test for zero-inflation in a ZIP model. The issue of count inflation in a positive Poisson distribution does not seem to have been considered in a similar way. In the second part of the thesis, we propose a one-inflated positive Poisson (OIPP) model, and develop a score test to determine whether there are “too many” ones for a positive Poisson model to fit well. We explore the performance of our score test, and compare it

(4)

to a likelihood ratio test, via Monte Carlo simulation. We find that the score test performs well, and that the OIPP model may be useful in many cases.

The third part of the thesis considers the possibility of one-inflation in zero-truncated data, when overdispersion is present. We propose a new model to deal with such a phenomenon, the one-inflated zero-truncated negative binomial (OIZTNB) model. The finite sample properties of the maximum likelihood estimators for the parameters of such a model are discussed. This Chapter considers likelihood ratio tests which assist in specifying the OIZTNB model, and investigates the finite sample properties of such tests. The OIZTNB model is illustrated using the medpar data set, which describes the hospital length of stay for a set of patients in Arizona. This is a data set that is widely used to highlight the merits of the zero-truncated negative binomial (ZTNB) model. We find that our OIZTNB model fits the data better than does the ZTNB model, and this leads us to conclude that the data are generated by a one-inflated process.

(5)

List of Tables

Table 1. %Bias, %MSE for , , and , and Pitman’s Nearness for and ... 31

Table 2. Mean confidence intervals and associated coverage probabilities for asymptotic confidence intervals around , and for bootstrapped confidence intervals around and . c = 0, d = 3, = 5, n = 25. ... 33

Table 3. c = 2, d = 8, = 20. Measures of the statistical significance of the bias correction. ... 35

Table 4. Number of diseased coronary vessels found at angiography (reproduced from Nademanee et al., Fig. 2). ... 36

Table 5. Estimates of using , , and , with associated bootstrap confidence intervals for and . ... 37

Table 6. Fitted frequencies, and p-values for Pearson’s goodness of fit test, using , , and . ... 38

Table 7. Simulation results for n = 25, λ = 0.25, compared to selected reproduced results from Mathews and Appleton (1993). ... 40

Table 8. Empirical Prob(type I error) for the score and likelihood ratio tests, when the DGP is consistent with the null hypothesis of no one-inflation, in the positive Poisson model. ... 60

Table 9. Estimation results for a positive Poisson regression on claims, from 1963 to 1973, subcategory amusement devices. ... 71

Table 10. Estimation results for a one-inflated positive Poisson regression on claims, from 1963 to 1973, subcategory amusement devices. ... 72

Table 11. Estimation results for a positive Poisson regression on claims, from 1963 to 1981, subcategory organic compounds. ... 74

Table 12. Estimation results for a one-inflated positive Poisson regression on claims, from 1963 to 1981, subcategory organic compounds. ... 74

Table 13. %Bias and %MSE for , , and , nrep = 10,000. ... 87

Table 14. Empirical Prob(type I error) and empirical type I critical values, for likelihood ratio tests where H0: OIPP model and HA: OIZTNB model. ... 93

Table 15. Empirical Prob(type I error) and empirical type I critical values, for likelihood ratio tests where H0: ZTNB model and HA: OIZTNB model. ... 95

Table 16. Empirical power of the likelihood ratio test where H0: OIPP model and HA: OIZTNB model. λ = 2, ω = 0.1. ... 98

Table 17. Empirical power of the likelihood ratio test where H0: ZTNB model and HA: OIZTNB model. λ = 6, α = 0.5. ... 100

Table 18. Estimation results for the medpar data set, under OIPP, ZTNB, and OIZTNB models. ... 104

Table B.1.1. %Bias, %MSE for , , and , and Pitman’s Nearness for and ... 122

Table B.1.2. ... 123

(8)

Table B.2.1. Mean confidence intervals and associated coverage probabilities for asymptotic confidence intervals around , and for bootstrapped confidence intervals

around and . c = 0, d = 3, = 0.5 ... 125 Table B.2.2. c = 0, d = 3, = 2 ... 126 Table B.2.3. c = 0, d = 3, = 5 ... 127 Table B.2.4. c = 2, d = 8, = 1.3 ... 128 Table B.2.5. c = 2, d = 8, = 5.3 ... 129 Table B.2.6. c = 2, d = 8, = 20 ... 130 Table B.2.7. c = 1, d = ∞, = 0.75 ... 131

Table B.3.1. c = 0, d = 3, = 0.5. % of lower C.I.s which are negative, using asymptotic standard errors around . ... 132

Table B.3.2. c = 0, d = 3, = 2 ... 132 Table B.3.3. c = 0, d = 3, = 5 ... 132 Table B.3.4. c = 2, d = 8, = 1.3 ... 132 Table B.3.5. c = 2, d = 8, = 5.3 ... 132 Table B.3.6. c = 2, d = 8, = 20 ... 132 Table B.3.7. c = 1, d = ∞, = 0.75 ... 132

Table B.4.1. c = 0, d = 3, = 0.5. Measures of fit. ... 133

Table B.5.1. c = 0, d = 3, = 0.5. Statistical significance of O(n-1) bias correction, relative to both and ’s sampling distribution... 135

Table C.1.1. Empirical Prob(type I error) for the score and likelihood ratio tests of no 1-inflation in the positive Poisson model, for nominal significances = 0.01, 0.05, 0.1. The case of no covariates. ... 137

Table C.1.2. Covariates according to design 1A. ... 138

Table C.1.3. Covariates according to design 1B. ... 139

Table C.1.4. Covariates according to design 2. ... 140

Table F.1. Empirical power of the likelihood ratio test of H0: OIPP model versus HA: OIZTNB model. λ = 2, ω = 0.1. ... 191

(9)

Table F.5. Empirical power of the likelihood ratio test of H0: ZTNB model versus HA:

OIZTNB model. λ = 2, α = 0.5. ... 195 Table F.6. Empirical power of the likelihood ratio test of H0: ZTNB model versus HA:

OIZTNB model. λ = 2, α = 2. ... 196 Table F.7. Empirical power of the likelihood ratio test of H0: ZTNB model versus HA:

OIZTNB model. λ = 6, α = 0.5. ... 197 Table F.8. Empirical power of the likelihood ratio test of H0: ZTNB model versus HA:

(10)

List of Figures

Figure 1. Venn diagram conveying how various count data issues are considered in the

following Chapters. ... 3

Figure 2. Empirical Prob(type I error) for the score and likelihood ratio tests, when the DGP is consistent with the null hypothesis of no one-inflation, in the positive Poisson model. ... 61

Figure 3. Empirical power of the score and likelihood ratio tests of no 1-inflation in the positive Poisson model, for various values of and , for the case of no covariates, for n = 250, and with nominal significance = 0.05. ... 63

Figure 4. Empirical power of van den Broek’s score test, for various and , and in the case of no covariates. n = 50. ... 68

Figure 5. Histogram of claims for observations from 1963 to 1973, subcategory amusement devices. ... 70

Figure 6. Histogram of claims for observations from 1963 to 1981, subcategory organic compounds. ... 73

Figure 7. The nesting of various zero-truncated count models. ... 90

Figure 8. Empirical Prob(type I error) for likelihood ratio tests where H0: OIPP model and HA: OIZTNB model. ... 94

Figure 9. Empirical Prob(type I error) for likelihood ratio tests where H0: ZTNB model and HA: OIZTNB model. ... 96

Figure 10. Empirical power of the likelihood ratio test where H0: OIPP model and HA: OIZTNB model, using nominal significance of 0.05. λ = 2, ω = 0.1 and 0.2. ... 99

Figure 11. Empirical power of the likelihood ratio test where H0: ZTNB model and HA: OIZTNB model, using nominal significance of 0.05. λ = 6, α = 0.5 and 2. ... 101

Figure 12. Actual proportion, and predicted probabilities of los, from the estimates from OIPP, ZTNB, and OIZTNB models. ... 107

Figure 13. Venn diagram conveying the issues addressed by the OIPP and OIZTNB models, compared to the issues addressed by existing count data models. ... 111

Figure A.1. Asymptotic O(n-1) bias correction plots ... 119

Figure A.2 ... 120

Figure A.3 ... 121

Figure C.2.1. Graphical representation of empirical Prob(type I error) for the score and likelihood ratio tests of no 1-inflation in the positive Poisson model, for nominal significances = 0.01, 0.05, 0.1. (Graphical representation of table in C.1.1). The case of no covariates. ... 141

Figure C.2.2. Covariates according to design 1A. ... 142

Figure C.2.3. Covariates according to design 1B. ... 143

Figure C.2.4. Covariates according to design 2. ... 144

Figure D.1.1. Empirical power of the score and likelihood ratio tests of no one-inflation in the positive Poisson model, for various values of n, , and , for various designs of , and with nominal significance = 0.05. The case of no covariates. n = 25. ... 145

(11)

Figure D.1.3. The case of no covariates. n = 100. ... 147

Figure D.2.1. Covariates according to Design 1A. n = 25. ... 151

Figure D.3.1. Covariates according to Design 1B. n = 25. ... 157

Figure D.4.1. Covariates according to Design 2. n = 25. ... 163

Figure E.1.1. Adjusted empirical power of the score and likelihood ratio tests of no 1-inflation in the positive Poisson model, for various values of , , and , for various designs of , and with nominal significance = 0.05. The case of no covariates. n = 100. ... 169

Figure E.1.2. The case of no covariates. n = 250. ... 170

Figure E.2.1. Covariates according to Design 1A. n = 25. ... 173

Figure E.3.1. Covariates according to Design 1B. n = 25. ... 179

Figure E.4.1. Covariates according to Design 2. n = 25. ... 185

(12)

(13)

Acknowledgments

I am deeply indebted to my supervisor, Dr. David Giles, not only for guidance on this thesis, but also as a teacher of econometrics, and as a mentor in many things. I am thankful to Dr. Judith Clarke for many helpful suggestions on earlier drafts, and also for my econometrics education. I am thankful for the comments received while presenting job market papers at the University of Manitoba and the University of Saskatchewan, for they greatly improved the direction of the final work of my thesis. I am grateful for the support received by many professors at the University of Victoria, and especially the Department of Economics staff. I am also thankful for helpful discussions with Jacob Schwartz. The patience and support of friends and family has been greatly appreciated, especially that of my partner, Ashley.

(14)

Dedication

(15)

Chapter One: General Introduction

There are several themes which are common to the following three Chapters of this thesis, the most prevalent being that the statistical models discussed pertain to count data. Data are typically considered to be “counts” when they are comprised of non-negative integers which are cardinal in nature. Count data arise in many areas, including medicine, demography, and various fields of economics. For example, the number of decayed, missing and filled teeth is a count which is an important indicator of the dental status of a person (Böhning et al., 1999). In demography, individual fertility is often measured by a count of children (Winklemann, 2008, p. 257). In health economics, count data may occur when describing the frequency of a health problem, or when observing health care utilization (Winklemann, 2008, p. 254).

When analyzing count data, a model should be specified that recognizes the limited support and discrete nature of the distribution for the variable being modelled. For example, the linear regression model ignores the non-negativity of count data, and can lead to “significant deficiencies” (Cameron and Trivedi, 1998, p. 2) if misapplied to count data. Furthermore, by ignoring the discrete nature of count data, the linear

regression model is unable to provide estimates for certain probabilities of events, which is typically the goal of the modelling exercise (Greene, 2008, p. 906).

The benchmark model for the analysis of count data is the Poisson model (Greene, 2008, p. 542). The Poisson distribution is typically used to model the number of events which

(16)

occur over a specific time period, with a classic example being the number of calls arriving at a switchboard over a certain interval. There are a number of stochastic

processes which can underlie a Poisson process. In one case, a Poisson process arises as a binomial limit. In another case, a Poisson process arises when inter-arrival times are independently exponentially distributed (Winklemann, 2008, p. 17).

The Poisson model is likely so important due to its relation to many other count models, rather than for its practical usefulness. The Poisson model is a specific case of many other more general models. That is, the Poisson model is nested by the Poisson-log-normal, negative binomial, generalized event count model, hurdle Poisson, zero-inflated Poisson, and doubly-truncated Poisson models (Winklemann, 2008, p.272). The elegance of the Poisson distribution contributes to its continued use as a building block for more complicated, yet empirically more relevant models. All of the models in this thesis, including the ones introduced, may attribute the Poisson model to their genesis.

While many models that build upon the Poisson distribution have subtly different characteristics (such as the negative binomial model), a characteristic which is common to the following three Chapters is almost immediately apparent from the data. Non-standard count data is a more accurate term for the type of data considered here, a title which implies immediate departure from the Poisson distribution. All three chapters in this thesis deal with count data where some of the integers on the non-negative real line, below or above a specific terminal, are impossible to observe. In Chapter 2, an arbitrary general case is considered, where the lower and upper terminals may take on any value.

(17)

In ze h a th o d is re C F fo If pr n Chapter 3 a ero truncatio owever. In C certain coun han count tru f this inferen eveloped an s the presenc eason for dis Chapter 3, an igure 1. Ven ollowing Ch f the non-sta rocess, then and Chapter on is conside Chapter 3, w nt must be ad uncation, and nce is the ma d discussed. ce of overdis scarding the nd Chapter 4 nn diagram c apters. andard featur their dismis r 4, only the ered. Trunca we consider th dded to the P d must be ca ain purpose o A further n spersion in th basic Poisso address the conveying ho

res noted abo ssal would ha special (and ation is not th he occurrenc Poisson proc arefully infer of Chapter 3 on-standard he data, whic on model. Fi issues, and c ow various c

ove are inde ave importan

d most freque he only

non-ce of count i cess. Count i rred from the 3, where test feature whi ch is likely t igure 1 below combination count data is ed relevant t nt implicatio ently encoun -standard fea inflation, wh inflation is m e data. The a ting procedu ch is addres to be the mo w illustrates n of issues, d

ssues are con

to the data g ons for the p

ntered) case ature allowed here the exce much more s accomplishm ures are sed in Chapt ost common how Chapte discussed ab nsidered in th generating properties of of d for, ess of subtle ment ter 4 er 2, ove. he

(18)

maximum likelihood estimation. Maximum likelihood estimation is also a central theme to this thesis, as it is the only estimation methodology considered. The method of

maximum likelihood estimation is primarily due to the work of R. A. Fisher, between the years 1912 and 1922, with some of the work being anticipated by Y. S. Edgeworth in 1908 (Aldrich, 1997). Among the properties of a maximum likelihood estimator, perhaps the most important are its (at least weak) consistency, and its asymptotic efficiency. Given that certain regularity conditions hold (White, 1982, p. 2), a maximum likelihood estimator is both weakly consistent and asymptotically efficient. Maximum likelihood estimation requires an explicit assumption of the parametric family responsible for the data generating process. Should the parametric family be chosen incorrectly, a

misspecification occurs, with the likely result being the loss of the desirable large-sample properties of the MLE. Dismissing any of the non-standard features of count data

discussed above would entail such a misspecification. While maximum likelihood

estimators are consistent (and asymptotically unbiased), they may or may not be biased in finite samples. This property is investigated in Chapter 2, in the context of the maximum likelihood estimator for the doubly-truncated Poisson distribution. The finite-sample bias of the estimator is explored, as are two different methods of bias correction. Bias is also explored in Chapter 4.

The Monte Carlo method provides a laboratory in which the properties of estimators and tests can be explored. This method is used extensively throughout this thesis. Monte Carlo experiments are used in Chapter 2 to assess the bias of three competing estimators; in Chapter 3 to compare two competing test statistics; and in Chapter 4 to determine the

(19)

viability of a new statistical model for a particular type of count data. The Monte Carlo method is associated with repetitive calculations performed by a computer, although the method itself is older than the computer. Monte Carlo methods are merely assisted by the fast generation of random numbers and quick computation enabled by computers. The Monte Carlo method has been used in many fields of research, including physics, engineering, statistics, and economics.

The Monte Carlo method was used as early as 1933 by Enrico Fermi, and likely

contributed to the work that won him the Nobel Prize in 1938 (Anderson, 1986, p. 99). It was not until 1947 that the Monte Carlo method was thusly named by Stanislaw Ulam, Nicolas Metropolis, and John von Neumann, after Stanislaw Ulam’s gambling uncle (Metropolis, 1987). It was a secret moniker for a method which would assist in the calculations necessary for nuclear fusion. The spirit of the Monte Carlo method is well captured in the sentiments of Stanislaw Ulam, as he recalls his first thoughts and attempts at practising the method. He was trying to determine the chances that a hand of solitaire would come out successfully. He wondered if the most practical method would be to deal one-hundred hands, and simply observe the outcome (Eckhardt, 1987). The use of

random generation and repetitive calculation are the two central tenets to Monte Carlo experimentation, a method which has flourished since the first electronic computer was built in 1945, and a method which has had a profound impact on mathematics and statistics since. “At long last, mathematics achieved a certain parity - the twofold aspect of experiment and theory - that all other sciences enjoy” (Metropolis, 1987, p.130).

(20)

This thesis proceeds as follows. In Chapter 2, we discuss the finite sample properties of a bias-corrected estimator for the parameter in the doubly-truncated Poisson distribution. In Chapter 3, we derive and examine a test for detecting inflation at the one count in the positive Poisson distribution (a special case of the doubly-truncated Poisson distribution). In Chapter 4, we propose an alternative model to the one discussed in Chapter 3, and discuss how to determine its appropriateness through specification testing. Chapter 5 provides some general concluding remarks.

(21)

Chapter Two: Bias Reduction for the Maximum

Likelihood Estimator of the Doubly-Truncated Poisson

Distribution

1. Introduction

This Chapter is concerned primarily with the performance of the maximum likelihood estimator for the shape parameter of a doubly-truncated Poisson distribution, with results extending to the more frequently applied positive Poisson distribution. Although the MLEs cannot be expressed in a closed form solution, we examine bias and MSE through Monte Carlo simulation, and obtain an analytic O(n-1) bias expression. We propose a simple bias-corrected estimator, discuss under what circumstances bias is likely to be problematic, and argue, for several reasons, that employing the analytic bias correction is preferred to the alternative of bias correction using the parametric bootstrap.

According to Firth (1993), in a model with a p-dimensional parameter vector , the asymptotic bias of the MLE may be written as

⋯ ,

where is a single bias term in an infinite sequence. The bias of an MLE is usually ignored on the grounds that it will be small compared with the standard errors (Cordeiro and McCullagh, 1991). In small enough sample sizes, however, bias may be significant. There are several techniques for determining bias, where focus is usually placed on the

(22)

first O(n-1) term, as the higher order terms disappear quickly as n grows. Bartlett (1953) was perhaps the first to provide a first order bias term, for a maximum likelihood estimator in a model involving one parameter. Once bias is determined, a bias-corrected estimator may be obtained by estimating the bias term at the MLE, and then subtracting this estimated bias from the MLE itself. Many developments have since been made; see Shenton and Bowman (1963, 1977) for example. Multi-parameter cases have been studied by Shenton and Wallington (1962) and by Cox and Snell (1968), for example. While there is a substantial literature extending these works, of particular import to this Chapter is the contribution of Cordeiro and Klein (1994), where analytical determination of the O(n-1) bias expression is somewhat simplified. Applications of Cordeiro and Klein’s (1994) procedure for finding the O(n-1) bias expression may be found in Giles (2010), Giles (2012), Giles et al. (2011), and Schwartz et al. (2011).

An alternative to the above “corrective” bias reduction method is a “preventive” approach proposed by Firth (1993), where the score vector is altered before the MLE is derived. However, this is not considered here. Other approaches include using resampling methods such as the jackknife or bootstrap. The jackknife (Quenouille, 1956) and the bootstrap (Efron, 1979) allow for O(n-1) bias corrections which can be implemented by brute force computing. The wide applicability of the bootstrap is an advantage over the more rigorous procedure unique to each distribution in the Cox and Snell/Cordeiro and Klein method. Once the analytic work has been undertaken, however, and the O(n-1_{) bias}

term has been determined, preference of one method over another is a matter of finite sample comparisons. This Chapter will provide first order bias terms for the MLEs for

(23)

the doubly-truncated and positive Poisson models, and compare analytic bias corrections to those enacted by a parametric bootstrap. We will also study various statistics which are constructed from these estimators.

The Poisson distribution is a discrete probability distribution expressing the probability of a number of events occurring over some specified interval, and has been a benchmark for the analysis of discrete count models. There have been several modifications made to the Poisson model in order to allow for the data to have various idiosyncrasies. When there is a fundamental difference in the data-generating process at the zero outcome, a preponderance of zeros perhaps, a hurdle model (Mullahy, 1986) or zero-inflated Poisson model might be employed. Examples of the need to accommodate excess zeros might arise in modelling manufacturing defects (Lambert, 1992), or in the demand for recreational boating trips (Gurmu and Trivedi, 1996).

Another distortion to the Poisson model occurs when the counts above or below a specific terminal conglomerate, requiring the use of a censored Poisson model. Survey data, for example, are typical of this tendency, when respondents may choose “X or more.” Terza (1985) considers such an application, where respondents could answer “three or more” at most, when asked how often they frequent a shopping area. Brännäs (1992) considers an employment questionnaire where the number of unemployment spells over a three year period could be at most “four or more.” Censoring may also occur when exact counts above a terminal value are infeasible, as may be the case in counting the number of ticks on sheep.

(24)

Yet another possible distinction from the regular Poisson model, and the model under consideration in this Chapter, is when the data generating process excludes the

observation of count outcomes outside of a specific interval. For example, observations on the number of items purchased when there is a “limit per customer” would exclude the zero class (since customers who do not buy the item will not be identified), and would exclude observations above the limit. Similarly, the number of serious criminal offences committed over a certain time period by felons under California’s three strikes law, is characteristic of a double truncation. Data of this type might be considered to follow a doubly-truncated Poisson distribution. Problems of truncation in a Poisson distribution, of one form or another, had early consideration by Bliss (1948), David and Johnson (1952), Moore (1952, 1954), Rider (1953), and Plackett (1953). However, it was not until Cohen (1954) that various forms of truncation were concisely generalized into the doubly-truncated Poisson model.

The doubly-truncated Poisson model encompasses singly-left or right truncation, as well as the regular Poisson model. However, the most frequent truncation phenomenon to arise is where the Poisson distribution otherwise appropriately describes the data generating process but for the omission of the zero-class. This is a special case of the doubly-truncated Poisson distribution, and is referred to as the zero-deleted or positive Poisson distribution. The prevalence of this form of truncation may be attributed to the tendency of an observational apparatus to become active only when at least one event occurs (Johnson and Kotz, 1969). For example, the size of groups observed in public places (Coleman and James, 1961), the number of accidents per worker in a factory

(25)

(David and Johnson, 1952), or the number of visits to a hospital, could be described by a positive Poisson model.

In Section 2 we present some background results for the doubly-truncated Poisson model, and in Section 3 we discuss maximum likelihood estimation of the model. In Section 4 we derive an analytic bias-corrected maximum likelihood estimator. In Section 5 we present the results from Monte Carlo experiments which assess the performance of our estimator. Section 6 presents some illustrative examples, and Section 7 provides some concluding remarks.

2. Background

2.1 Doubly-Truncated Poisson Model

From Cohen (1954, p. 160) the probability mass function for a doubly truncated Poisson random variable, y, is:

1

! ; , ,

(1)

where c and d correspond to the left and right truncations respectively, and

! ,

is the cumulative probability that c or more occurrences will be observed in the full Poisson distribution. Cohen’s notation is appropriate to the primary intent of his paper; namely, facilitating the non-linear computation of the MLE in a doubly-truncated Poisson model. We find it convenient to use the result that:

(26)

1

! , allowing (1) to be rewritten as:

! ! . (2)

The first two moments of this distribution are readily shown to be:

∑ _{1 !} ∑ _! , (3) ∑ _{1 !} ∑ _! ,

giving the variance as:

∑ _{1 !}

∑ _!

∑ _{1 !}

∑ _!

. (4)

The regular Poisson model assumes an equality between mean and variance, that is, the model assumes equidispersion. In the doubly truncated Poisson model, the mean and variance are characteristic of underdispersion. The difference in the relationship between the mean and variance for the two models provides an intuitive reason as to why fitting a regular Poisson model to truncated data is a fundamental misspecification. A model which assumes equidispersion should not be used to describe data which arises from a process characteristic of underdispersion.

(27)

2.2 Positive Poisson Model

As noted above, a special case of the doubly-truncated Poisson distribution is the positive Poisson distribution. As this special case encompasses the bulk of empirical applications, in this Chapter the positive Poisson distribution has been considered somewhat

separately. There are two advantages to this approach. First, the bias correction formula derived for the doubly-truncated case (in section 4.4) can be somewhat refined, and the coding of various functions associated with the MLE can be greatly simplified. Second, it can be shown that the Hessian of the log-likelihood is unambiguously negative, hence the MLE is unique. In the case of double-truncation, we do not have as tidy a proof regarding uniqueness, only a strong argument.

The mass function for a random variable, y, which follows a positive Poisson distribution is

1 ! . (5)

The first two moments are:

1 , (6)

1 1 , giving the variance as:

(28)

As can be seen from (6) and (7), the positive Poisson model exhibits underdispersion. The rightmost term on the right side of (7) may be considered an “adjustment factor,” and must be acknowledged in tests of overdispersion (Gurmu, 1991).

3. Maximum Likelihood Estimation

3.1 Maximum Likelihood Estimation of the Doubly-Truncated Poisson

Parameter

The sole estimation methodology under consideration in this Chapter (and entire thesis) is maximum likelihood estimation, as it is germane to the O(n-1) bias correction to follow, and it has been the dominant methodology for Poisson and modified Poisson models in the literature.

Under independent sampling from (1), Cohen writes the log-likelihood as

ln 1 ln ln ! ,

which may be rewritten as:

ln ln

! .

The maximum likelihood estimator, , is found by solving the following first-order condition:

∑ _{1 !}

∑ _!

(29)

Here we note that there is no closed form solution for , hence estimation must be accomplished numerically.

3.2 Uniqueness of the Doubly-Truncated Poisson MLE

To ensure the log-likelihood has been maximized, it would be sufficient to determine that the following scalar Hessian quantity is negative:

∑ _{2 !}

∑ _!

∑ _{1 !}

∑ _!

. (8)

However, we cannot sign the above term, until having first observed the sample

average, . This may warrant a certain amount of concern when estimation is undertaken. We cannot preclude multiple maxima, for example. Furthermore, the very possibility of a positive sign casts doubt on the appropriateness of a normal approximation for the

sampling distribution of , and on the validity of asymptotic standard errors and confidence intervals (Pratt, 1981, p.103).

Fortunately, examining the expectation of the Hessian permits insight suggesting the log-likelihood will be concave, provided we are comfortable with the sophistication of the numerical method used to obtain . Using (3) we can write:

∑ _{1 !} ∑ ! ∑ _{2 !} ∑ ! ∑ _{1 !} ∑ ! ,

(30)

where henceforth it is understood that summations occur over j, from c to d. We wish to show that the RHS of the above expression is negative, or:

∑ _{1 !} ∑ ! ∑ _{2 !} ∑ ! ∑ _{1 !} ∑ ! .

Equivalently, we need to show that:

! 1 ! 2 ! 1 ! , or: 1 ! 1 ! 1 1 ! 1 ! , 1 1 ! 1 ! 1 1 ! 1 ! , 1 1 ! 1 ! 1 ! .

Now consider _! and _!. Then by the Cauchy-Schwarz inequality we have:

,

and the expected value of the second derivative of the log-likelihood is negative. One should wonder how this is of help, since the Hessian of the log-likelihood evaluated at the MLE should still require inspection of to determine concavity. By taking the

(31)

substitution of ∑ !

∑ _! for is unambiguously negative. As was noted by David and

Johnson (1952), solving the first-order condition amounts to equating the sample mean value of to its expectation, and so is found by the solution of

∑ _{1 !}

∑ ! .

Insofar as our numerical methods are accurate, successful convergence will imply that, when evaluated at the MLE, the above first-order condition will be very nearly satisfied. Hence, proving negativity of the expectation of the Hessian is very near to proving concavity of the log-likelihood. In 1,750,000 simulations over seven different parameter configurations, the maximum deviation between the expected Hessian and the actual Hessian was 0.504%, which only signifies that asymptotic standard errors are unlikely to be affected by the choice of Hessian. In the simulations, the Hessian was always found to be negative.

Although the evidence above is quite compelling in terms of establishing the negativity of the Hessian, an alternative approach could be considered. Namely, the probability of the scalar Hessian being negative could be evaluated. However, the complexity of the expression in (8) suggests that this probability would have to be approximated by a suitable expansion, and it is unclear what would be considered a sufficiently small probability to indicate a “clear” result. Accordingly, this possibility is not pursued here, but remains a matter for future investigation.

(32)

3.3 Maximum Likelihood Estimation of the Positive Poisson Parameter

Under independent sampling, from (5), the log-likelihood function is

ln ln 1 ! .

The maximum likelihood estimator, , is found by solving

1 0 . (9)

As noted above, solving the first order condition has the interpretation of finding a value for which equates the mean of the distribution to the mean of the data, so the maximum likelihood estimator is also the method of moments estimator. Also, without exception, (9) may be rewritten as:

1 ∑

,

where it may be recalled that the LHS of the above first order condition is equal to (6), the mean of the distribution.

3.4 Uniqueness of the Positive Poisson MLE

The Hessian of the log-likelihood is given by

1 .

A sufficient condition for uniqueness of the MLE is that the above expression should be negative. The Hessian will become more negative the higher is the sample average, . Since it is impossible to observe a 0 in a positive Poisson model, has a lower bound of 1. A condition for uniqueness may then be written as:

(33)

1

. (10)

Rewriting (10) as:

1 , (11)

and using a Taylor series expansion on both sides of (11) we get: 1

2 2! 2 3! 2 ⋯ 1 2! 3! 4! ⋯ . (12)

Subtracting the LHS of (12) from (12) we can write:

0 2 1

1 ! 2 . (13)

By examining the term inside the numerator of the summation in (13), it can be seen that 2 1, for all 1. Hence, the inequality in (10) holds, and the MLE for the shape parameter in a positive Poisson model is unique.

4. Bias of the Maximum Likelihood Estimators

4.1 Background

Let l(θ) be the log-likelihood function based on a sample of n observations, with p-dimensional parameter vector, θ. l(θ) is assumed to be regular with respect to all

derivatives up to and including the third order. The joint cumulants of l(θ) are defined as follows:

; , 1,2, … . , (14)

(34)

, ; , , 1,2, … . , . (16)

The derivatives of the cumulants are defined as:

; , , 1,2, … . , . (17)

All of the expressions in (14) - (17) are assumed to be O(n), and (14) is the , element of Fisher’s expected information matrix, . Cox and Snell (1968) showed that when the sample data are independent, the bias of the sth element of is:

0.5 _, ; 1,2, … . , , (18)

where is the , element of the inverse of . Extending the work of Cox and Snell (1968), Cordeiro and Klein (1994) show that (18) can be rewritten, where the data need not be independent provided that all of the k terms are O(n):

0.5 ; 1,2, … . , . (19)

Note that the Cordeiro and Klein version does not involve the terms as defined in (16). Now, let 0.5 , for , , 1,2, … . , ; and define the following matrices:

; , , 1,2, … . , (20)

… … . | . (21)

Equations (20) and (21) allow (19) to be rewritten as:

. (22)

Finally, we can define the “bias-corrected” MLE for θ as:

(35)

where | and | . Note how (22) and (23) exhibit the merit of the Cox and Snell/Cordeiro and Klein procedure; they can be evaluated even when the MLE cannot be expressed in a closed form solution, as (22) and (23) are achieved solely through manipulation of the cumulants of the likelihood function. They do not require an explicit algebraic formula for .

4.2 On Treating Doubly-Truncated and Positive Poisson Separately

A practitioner working with singly-left-truncated data ( ∞) may easily implement the bias correction formula provided for doubly-truncated data, by selecting a sufficiently large value for d. However, when ∞ is implicitly considered, the bias correction formula becomes much more elegant, and significant coding hurdles are eliminated. For example, in the doubly-truncated case, terms with a negative factorial in the denominator sometimes arise. This is dealt with by simply setting those terms equal to zero, using the convention that the factorial of a negative number is equal to 1 0⁄ . Hence, a term with a negative factorial in the numerator is undefined, but a term with a negative factorial in the denominator becomes zero. As computer algorithms do not take account of this

convention, such terms with negative factorials must be sought out “manually.” Great care must be taken in coding not only the bias correction formula provided in this

Chapter, but also the gradient of the log-likelihood, the Hessian and expected Hessian of the log-likelihood, and the mean and variance of the distribution implied by the MLE. In short, one who requires a zero-deleted Poisson distribution would be ill-advised to use results obtained for doubly-truncated Poisson.

(36)

To exhibit the Cox-Snell/Cordeiro-Klein bias correction formula in a more agreeable setting to begin with, the procedure is applied to the positive Poisson distribution in section 4.3. Section 4.4 provides straightforward results for the more general doubly-truncated Poisson distribution.

4.3 Bias of the MLE in the Positive Poisson Distribution

In order to obtain Cordeiro and Klein’s bias expression (22), we need the first three derivatives of the log-likelihood with respect to the shape parameter, :

∑ 1 ∑ 1 2 ∑ 1 1 .

Using the result that / 1 , we can then obtain the following expressions: 1 1 , and 2 1 1 1 . So, the expected information measure is:

1

1 ,

(37)

1 1 1 1 , and 0.5 2 2 2 1 .

The final step is to pre-multiply A by and post-multiply by , which in the single parameter case, amounts to dividing A by . Thus, Cordeiro and Klein’s (1994) bias expression, to O(n-1_{) is:}

1 2 2

2 1 .

A sufficient condition for the above O(n-1) bias term to be negative is that:

2 2 0 . (24)

Since (24) may be rewritten as 1

! ; 0 ,

the O(n-1) bias for the MLE in a positive Poisson model is strictly negative.

4.4 Bias of the MLE in the Doubly-Truncated Poisson Distribution

In addition to the second derivative of the log-likelihood (8), we require:

2 ∑ 3 ! ∑ ! 3 ∑ _{1 !}∑ _{2 !} ∑ ! 2 ∑ _{1 !} ∑ ! . (25)

(38)

∑ _{1 !} ∑ ! ∑ _{2 !} ∑ ! ∑ _{1 !} ∑ ! , 2 ∑ _{1 !} ∑ ! ∑ _{3 !} ∑ ! 3 ∑ _{1 !}∑ _{2 !} ∑ ! 2 ∑ _{1 !} ∑ ! .

The expected information measure is . Additionally, we have:

∑ _{1 !} ∑ ! ∑ _{1 !}∑ 2_! ∑ ! ∑ _{3 !} ∑ ! 3 ∑ _{1 !}∑ _{2 !} ∑ ! 2 ∑ _{1 !} ∑ ! , and ∑ _{1 !} ∑ _{1 !} ∑ _{1 !}∑ 2_! ∑ ! ∑ _{1 !} ∑ ! ∑ _{3 !} 2 ∑ ! 3 ∑ _{1 !}∑ _{2 !} 2 ∑ ! ∑ _{1 !} ∑ ! ,

giving the O(n-1) bias expression of the MLE from a doubly-truncated Poisson model as:

(39)

1 ∑ jjλ 1 ! ∑ λj! ∑ _j λ_{1 ! ∑} j 2 λ_j! ∑ λj! ∑ _j λ_{1 !} ∑ λj! ∑ _jλ _{3 !} 2∑ λj! 3∑ _jλ _{1 ! ∑} _jλ _{2 !} 2 ∑ λj! ∑ _jλ _{1 !} ∑ λj! ∑ _j λ_{1 !} ∑ λj! 2∑ _j λ_{1 ! ∑} _jλ _{2 !} ∑ λj! ∑λj! 2∑ _j λ_{1 ! ∑} _jλ _{1 !} ∑ λj! ∑λj! ∑ _jλ _{2 !} ∑ λj! 2∑ _jλ _{2 ! ∑} _jλ _{1 !} ∑ λj! ∑ _jλ _{1 !} ∑ λj! . (26)

Plots of some numerical evaluations of this expression may be found in Appendix A. Simulations reveal that (26) will be negative (positive) when is such that the mean of the distribution is close to the left (right) truncation. These simulations also suggests that the main factor (other than sample size) in determining the severity of the O(n-1_{) bias, is}

the nearness of the mean of the distribution to a left or right truncation value. From Figures A.1 and A.2 in Appendix A, we see that a low mean is not associated with bias when there is no left truncation (c = 0), but as the mean moves towards the right truncation, the positive bias becomes worse. In Figure A.3, we see that nearness of to c is associated with negative bias, while the bias becomes positive and is at its worst when nears d.

(40)

5. Monte Carlo Experiments

5.1 Experimental Design

The following Monte Carlo experiments compare the performance of various statistics based on , , and ; the MLE, analytic bias-adjusted MLE, and (parametric) bootstrap bias-adjusted MLE, respectively. The three estimators themselves are compared in terms of percent bias, percent mean-square error (MSE), and Pitman’s measure of nearness (PN)1, while other statistics are compared by examining means and coverage

probabilities.

The Monte Carlo experiments were undertaken using the R statistical environment (R, 2010). While randomly generated zero-deleted Poisson variates are available in the VGAM (Yee, 2011) package, doubly-truncated variates are not. Doubly-truncated Poisson random variates were generated by using the cumulative mass function to transform a randomly generated 0,1 variate. This inversion method was verified by checking the empirical sample moments against the moments of the theoretical

distribution.

Maximizing the log-likelihoods in R was achieved using the optimize algorithm in the

stats (R, 2010) package, a routine which searches over an interval for a maximum. While

several methods were attempted, this simple line search method proved faster and as accurate as using the first-order condition in a Newton-Raphson algorithm. The lower

(41)

and upper bounds of the search interval were set at 0.01 and 99, respectively. There was no indication that a corner solution had occurred in any of the experiments.

The parametric bootstrap was used as an alternative bias correction method against which the analytic bias correction could be compared. The simple, non-parametric bootstrap, involves resampling from the initial data. Alternatively, the parametric bootstrap differs in that resampling occurs from the fitted model. Hall (1998, p.930) and Davison and Hinkley (1997, p.15) describe the differences between the simple bootstrap and the parametric bootstrap. When a fully specified parametric model is known to be the data generating process, as in the case of maximum likelihood estimation, it is generally thought that the parametric bootstrap is more appropriate than the nonparametric bootstrap since more information is included in the resampling plan.

The parametric bootstrap bias-adjusted estimator was obtained in the following way. First, the MLE was obtained. The fitted model, obtained by substituting into the cumulative density function of the doubly truncated Poisson model, is used to generate

NB = 1000 bootstrap samples. The parametric bootstrap bias-adjusted estimator, , is then

obtained as 2 1⁄ ∑ , where is the MLE of obtained from the

jth of the NB bootstrap samples. Ideally, NB should be as large as computing and time

constraints allow. Our choice of NB = 1000 took time constraints into consideration,

however, preliminary experimentation suggested that increasing NB beyond 1000 made

very little difference to . While the parametric bootstrap performed better than the non-parametric bootstrap in this setting, further refinements may be possible. It should

(42)

therefore be noted that the Monte Carlo results may favour the analytic adjustment, when comparing the analytic adjustment to the bootstrap adjustment. We have elected to employ a bootstrap technique which is straightforward and likely to be used in practice. One possible improvement on the parametric bootstrap could be made by instead considering the double-bootstrap (Beran, 1988). As the bootstrapping technique employed in this Chapter assumes that the bias is a constant, another possible

improvement is to use a linear bootstrap correction, proposed by MacKinnon and Smith (1998). The linear bootstrap has been shown to be extremely effective (Giles et al., 2011), but is not considered here.

We now describe the essential structure of a Monte Carlo experiment, by describing one repetition of the computer program used to enact an experiment. Having set the true parameter and sample size n for the desired data generating process, a random sample,

y, is generated. The estimates , , and are then obtained. Asymptotic confidence

intervals are constructed around . Parametric bootstrapped confidence intervals are constructed around and , using 1000 bootstrap samples. In the case of , bootstrapped confidence intervals are accomplished using the same bootstrapped statistics used to determine . In the case of , confidence intervals are achieved by generating 1000 random samples as if were the true data generating process, calculating the MLE and bias-correcting using (26), for each sample. In both bootstrapped confidence intervals, we simply replicate the calculation of the estimate in the same way it was determined in the first place. It is not feasible to bootstrap confidence intervals around , due to time constraints, since replicating involves two layers of bootstrapping. Using the

(43)

constructed confidence intervals, we record whether or not the true parameter lies within the interval, so that we can examine coverage probability. For the asymptotic confidence interval, we note whether or not the left confidence interval is negative. We also record how extreme and are relative to each other. This is accomplished by placing and in the other’s bootstrapped statistics, which were obtained when constructing the

bootstrap confidence intervals. Finally, the three estimates , , and are used to generate predicted probabilities at each count, which are compared to the actual

proportions at each count, in a goodness of fit test. The above outlined process is repeated 50,000 times for each experiment.

The reader may be somewhat surprised at the rather large values of the Poisson parameter that were used in the Monte Carlo experiments. When a truncation from above occurs, the implication of a higher is that it brings the mean of the distribution closer to the right truncation, d. As → ∞, → . As noted earlier, it appears that the closer the mean of the distribution is to d, the larger is the O(n-1) bias.

5.2 Monte Carlo Results

Table 1 reports some selected results for percent-bias and percent-MSE, while the results of experiments for seven different parameter configurations may be found in Appendix B. Table 1 presents a situation where the analytic adjustment performs much better than the parametric bootstrap, when there is an appreciable amount of bias for n = 10. In fact, the parametric bootstrap makes the bias worse, although it reduces the %MSE at n = 10. The fact that one of the bias correction techniques makes bias worse in small samples is

(44)

not unprecedented. Both bias correction techniques eliminate an O(n-1) bias component; in small samples the remaining components may be of the opposite sign, and together larger in magnitude. As n increases, the parametric bootstrap and analytic

bias-adjustment become increasingly similar. This is because the O(n-1) bias component accounts for a larger portion of the bias as the sample size grows.

Tables B.1.1 and B.1.2 in Appendix B exhibit the tendency of the bias to become worse as the mean of the distribution moves toward the right truncation, d. When the mean of the distribution is not close to the right truncation d, (for example when c = 0, d = 3, = 0.5 and = 2, and when c = 2, d = 8, = 1.3 and = 5.3), the bootstrap and analytic bias-adjustment perform very similarly, and there does not appear to be a clear advantage to one method over another. In Table B.1.3, where is chosen such that the mean of the positive Poisson distribution is close to the (only) left-truncation, it can also be seen that the bootstrap and analytic bias-adjustments are very similar. In the cases mentioned above, the percent-bias is hardly of concern. It is not until the bias becomes substantial that the analytic method dominates in small samples. Pitman’s measure of nearness also supports the above observations.

Pitman’s measure of nearness can be used to compare two competing estimators. It expresses the probability that one estimator is “nearer” (in absolute terms) to the true parameter than the competing estimator. In five out of the seven parameter configurations examined, Pitman’s measure of nearness either favours or heavily favours , when n = 10

(45)

and n = 15. Out of 35 calculations of Pitman’s nearness, is favoured 21 times, and is never heavily favoured.

Table 1. %Bias, %MSE for , , and , and Pitman’s Nearness for and

c = 0, d = 3, = 5, = 2.4 n [%MSE] %Bias [%MSE] %Bias [%MSE] %Bias

Pitman’s Nearness2 . 10 20.883 [59.676] -9.528 [9.666] -29.937 [12.271] 0.872 15 12.918 [28.406] -2.300 [9.938] -9.726 [7.493] 0.595 25 _[10.830]6.733 _[6.909]-0.490 _[6.272]-2.174 0.473 50 _[3.878]3.002 _[3.290]-0.156 _[3.234]-0.408 0.484 100 _[1.748]1.415 _[1.620]-0.081 _[1.616]-0.136 0.510 While the choice of bias-adjusted estimator may be clear when the mean of the distribution is close to the right truncation, we have noted that there is no clear choice otherwise, in terms of percent-bias and percent-MSE. There may be computational reasons for preferring the use of the analytic adjustment, however. Namely, when the practitioner wishes to bootstrap standard errors, confidence intervals, or p-values, in

addition to bias-adjusting the MLE.

The bootstrap builds up an empirical distribution for a statistic by calculating many statistics in the same manner in which the original statistic was calculated. This empirical distribution may then be used to correct bias, or to determine a p-value, for example. If the original statistic is (the bootstrap bias-adjusted MLE), for example, constructing the

(46)

empirical distribution for becomes computationally burdensome. If one were to choose 1000 bootstraps, constructing the empirical distribution for would entail the solving of 1 million non-linear equations. If the analytic bias-adjustment is used instead, the practitioner does not face the same restraints should they wish to bootstrap a p-value.

Tables B.2.1 to B.2.7, in Appendix B, show the mean (over 50,000 Monte Carlo

repetitions) confidence intervals and associated coverage probabilities for , using both asymptotic standard errors3 and bootstrapped confidence intervals. Also reported are the mean bootstrapped confidence intervals and associated coverage probability for . Recall that bootstrapping confidence intervals around is not at all feasible in a Monte Carlo setting. As can be seen, the bootstrapped confidence intervals can be quite different from the ones obtained through asymptotic theory. In most cases, confidence intervals

bootstrapped around are much narrower, while maintaining nominal coverage probability.

Table 2 reproduces a portion of Table B.2.3 from Appendix B. The results from Table 1,

n = 25, do not clearly indicate that the analytic bias-adjustment is preferred. In the

situation depicted in Table 2, however, it appears that both bias-adjustment and the bootstrapping of confidence intervals or p-values would be required. Bootstrapped confidence intervals around are markedly different from the ones constructed via asymptotics, suggesting that at n = 25 the asymptotic normal distribution may be a poor approximation to the finite sample distribution of .

(47)

Table 2. Mean confidence intervals and associated coverage probabilities for asymptotic confidence intervals around , and for bootstrapped confidence intervals around and . c = 0, d = 3, = 5, n = 25.

Mean confidence interval:

99% C.I. 95% C.I. 90% C.I. estimate Asy.se.( ) 1.669 2.546 2.994 5.337 7.679 8.128 9.005 Boot. ( ) 2.957 3.380 3.622 9.152 10.435 14.021 Boot. ( ) 2.676 3.062 3.271 4.975 7.334 8.083 9.832 C.P. (%) at nominal: 90% 95% 99% Asy.se.( ) 90.55 93.74 97.65 Boot. ( ) 88.83 93.75 98.03 Boot. ( ) 89.14 94.48 99.06

Tables B.3.1 to B.3.7 provide some further evidence towards the need to bootstrap confidence intervals, showing the percentage of asymptotic confidence intervals which extend into the negative domain. This could be problematic for hypothesis testing when is near its boundary of zero, and is evidence that the asymptotic normal distribution is a poor approximation for the sampling distribution of when is small. The occurrence of a negative confidence interval can be quite frequent, especially when is small to begin with. Together with Tables B.2.1-B.2.7, these results suggest that bootstrapping

confidence intervals is more appropriate than relying on asymptotics, at least for sample sizes under 25, in the present context.

Tables B.4.1 – B.4.7 report the squared differences of the fitted probabilities from the theoretical ones, summed across each count, and averaged over Monte Carlo replications.

(48)

These squared differences are determined for the three different estimators, suggesting that using a bias-adjusted estimator improves the fit of the model marginally, although there is no clear difference between using or . To further explore the performance of fitted probabilities using , , and , we also report coverage probabilities for Pearson’s goodness-of-fit test. The goodness-of-fit test has null hypothesis : the fitted model generated the data, and alternative hypothesis : not . The goodness-of-fit test

statistic is calculated by ∑ , where are the observed frequencies and are the expected frequencies determined by the fitted model, and k are the number of unique counts. is asymptotically chi-square distributed with degrees of freedom equal to k – 1, the number of unique counts minus the number of estimated parameters. For example, when c = 0 and d = 3, there are only four unique counts, and the degrees of freedom of are 3. When d = ∞, k must be observed from the sample, y. The Monte Carlo program calculates the test statistic in each repetition and notes whether the null hypothesis is rejected or not. Over 50,000 replications, this information is used to

determine coverage probability. The reported coverage probabilities are very close to nominal when c = 0, d = 3, and when c = 2, d = 8, = 5.3. Otherwise, coverage

probability is generally quite poor, even when n = 100. These findings may suggest that an alternative method of goodness-of-fit testing is warranted, in the present context.

As noted, a criticism which arises in the bias correction literature is that the bias of an estimator is likely to be small relative to the estimator’s standard error. To address this criticism, we provide two different bootstrap p-values, which are computed by examining where and fall in each other’s empirical distribution, as determined by 1000 bootstrap

(49)

statistics. These two p-values are denoted Pr , Pr and Pr , Pr , where Pr denotes that the empirical distribution is determined from the calculation of many ’s. The use of the min operator ensures that the reported p-value is always below 0.5. These p-values are averaged over 50,000 Monte Carlo repetitions, and provide a direct measure of the average statistical significance of the analytic O(n-1_{) bias correction. Table 3 below reports p-values for one parameter}

configuration, and the remaining results may be found in Tables B.5.1 – B.5.7.

Table 3. c = 2, d = 8, = 20. Measures of the statistical significance of the bias correction. n Pr , Pr Pr , Pr 10 0.238 0.107 15 0.288 0.183 25 0.335 0.261 50 0.382 0.331 100 0.417 0.381

The analytically-bias-corrected estimator can be significantly different from the uncorrected one, especially in small samples. Comparing the columns of Table 3 also exhibits the relatively smaller variance of .

(50)

6. Illustrative Examples

6.1 Number of Diseased Coronary Vessels

The following illustrative example uses data from Figure 2 of Nademanee et al. (1987, p.6), on the number of diseased coronary vessels found at angiography. Data is

reproduced below, in Table 4.

Table 4. Number of diseased coronary vessels found at angiography (reproduced from Nademanee et al., Fig. 2).

# diseased

vessels All groups Group A* Group B* Group C*

0 3 3 0 0 1 9 6 1 2 2 14 4 5 5 3 13 1 2 10 *Group A = no transient myocardial ischemia; Group B = transient myocardial ischemia

< 60 minutes/24 h; Group C = transient myocardial ischemia > 60 minutes/24 h

The data describes the number of diseased coronary vessels (diseased arteries which supply blood to the heart) which are found at angiography, divided into groups based on varying degrees of self-reported transient myocardial ischemia. Angiography is an

invasive procedure which allows visualization of the vessel, in order to determine if there is obstruction of blood flow to the heart muscle. Transient myocardial ischemia (angina) is essentially pain in the heart due to lack of oxygen. Estimates of for the different groups, along with confidence intervals, are provided in Table 5, using the three different estimators discussed.

(51)

Except for Group C, at first glance there does not seem to be an appreciable difference between the estimates. Upon closer examination, however, it appears that some of the confidence intervals can be quite different. In particular, the upper (bootstrapped) critical values for Group B experience quite a shift. Under maximum likelihood estimation, rather than our bias-adjusted maximum likelihood estimation, we would be much less likely to conclude that Groups B and C are statistically different. We should also note that Group C, the group with mean closest to the right truncation, experiences the most drastic change out of all the groups.

Table 5. Estimates of using , , and , with associated bootstrap confidence intervals for and . Group 99% C.I. 95% C.I. 90% C.I. All 2.570 2.793 2.957 3.906 5.272 5.755 7.019 A 0.659 0.738 0.819 1.372 2.135 2.298 2.891 B 1.244 1.846 2.096 3.644 8.594 12.71 24.85 C 3.135 3.610 3.892 6.131 10.87 13.48 17.79 All 2.524 2.664 2.816 3.798 5.068 5.509 6.028 A 0.576 0.728 0.807 1.336 2.047 2.196 2.534 B 1.193 1.536 1.731 3.183 4.950 5.827 6.470 C 2.329 2.976 3.176 5.488 8.746 10.16 11.92 All 3.798 A 1.367 B 2.968 C 5.154

In Table 6, we show how the different estimators affect the fit of the model. There are minor changes in the fitted frequencies, with the models appearing to do remarkably well.

Econometric analysis of non-standard count data

Supervisory Committee

Abstract

Table of Contents

List of Tables

List of Figures

Acknowledgments

Dedication

Chapter One: General Introduction

Chapter Two: Bias Reduction for the Maximum

Likelihood Estimator of the Doubly-Truncated Poisson

Distribution

1. Introduction

2. Background

2.1 Doubly-Truncated Poisson Model

2.2 Positive Poisson Model

3. Maximum Likelihood Estimation

3.1 Maximum Likelihood Estimation of the Doubly-Truncated Poisson

Parameter

3.2 Uniqueness of the Doubly-Truncated Poisson MLE

3.3 Maximum Likelihood Estimation of the Positive Poisson Parameter

3.4 Uniqueness of the Positive Poisson MLE

4. Bias of the Maximum Likelihood Estimators

4.1 Background

4.2 On Treating Doubly-Truncated and Positive Poisson Separately

4.3 Bias of the MLE in the Positive Poisson Distribution

4.4 Bias of the MLE in the Doubly-Truncated Poisson Distribution

5. Monte Carlo Experiments

5.1 Experimental Design

5.2 Monte Carlo Results

6. Illustrative Examples

6.1 Number of Diseased Coronary Vessels