5
$\begingroup$

I have large data set of $\approx 10^6$ points where each point contains the information of a

(year, count)

of a particular event. There are many counts for a each year. Looking at the histograms of counts for a given year, I've conjectured that they follow a Poisson distribution. By generating the Poisson distribution with the maximum likelihood parameter $P(\lambda_{MLE})$, $\lambda_{MLE}=\frac{1}{n}\sum{k_i}$. I think what I need to do next is compare the values of this expected distribution to my normalized histogram and compute a $\chi^2$ statistic. Wiki however warns that:

In cases where the expected value, E, is found to be small (indicating a small underlying population probability, and/or a small number of observations), the normal approximation of the multinomial distribution can fail, and in such cases it is found to be more appropriate to use the G-test, a likelihood ratio-based test statistic.

This will certainly be the case for larger values of $\lambda$ at the more extreme counts. How do I calculate how "good" my data will fit to the Poisson distribution? Ideally what I would like is to show a plot of ($\lambda$ vs year) with an error bar related to this "goodness".

$\endgroup$
2
$\begingroup$

The direct answer to the question - how to deal with small expected counts - is that one might either

(a) combine ranges of cells at the end (a very common approach),

(b) use a different (and perhaps better) goodness of fit test, or

(c) consider dropping the chi-square approximation, and see if one can deal with the discrete distribution of the test statistic more directly, perhaps by simulation.

Approach (a) can be found in many texts. There are many ways you can go about combining cells, but many people simply work from one end or the other, combining cells either to the left, or to the right, until the expected counts are sufficiently high for their purpose.

(However, the most commonly cited rule of thumb for the expected number - that it should be at least 5 for the chi-square approximation to hold - is unnecessarily strict for the sort of approximation most people would require. Many subsequent papers have suggested less stringent rules.)

The other answer by user36381 suggests that with such a large sample size, goodness of fit tests are almost certain to reject; this is true. However, I'm not so sure comparing to other reference distributions will help, since they, too, would almost certainly be rejected by a decent goodness of fit test.

(Why are you testing whether it's Poisson? If you have around a million data points, the sample itself contains a lot of information about distributional shape - do you actually need a name for the distribution?)

$\endgroup$
  • $\begingroup$ Per your last point, I don't really need to prove if the distribution is Poisson, but for all my samples the mean and the variance were very close suggesting that it would be a good fit. Since the Poisson distribution is one-parameter, it makes for a nice plot versus time since it conveys both the mean and the variance (or at least close to it). Since I'm only plotting $\lambda$ vs time, I needed some way of estimating how accurate this was and that a single $\lambda$ was essentially the same as showing each distribution. $\endgroup$ – Hooked Dec 18 '13 at 20:21
2
$\begingroup$

This doesn't entirely answer your question regarding how to implement the goodness of fit test. But with a sample size that large, your data will almost certainly fail any test comparing it to a Poisson distribution even if your distribution closely resembles a Poisson distribution (unless, of course, you simulated Poisson data). The power of your test to detect minute deviations from the Poisson distribution would be very high.

It may be better to compare it to other reference distributions to see which distribution it fits best.

$\endgroup$
  • $\begingroup$ See the following link (particularly Warning 2) regarding the use of distribution goodness of fit tests for large samples. It has some good examples: blog.minitab.com/blog/statistics-and-quality-data-analysis/… $\endgroup$ – Underminer Dec 18 '13 at 17:58
  • $\begingroup$ So to get at my question, would the deviation of the target distribution to my normalized dataset (e.g. the raw score of $\chi^2$) be a good error bar to put on the chart? $\endgroup$ – Hooked Dec 18 '13 at 18:35
  • $\begingroup$ I'm not exactly sure what you want the error bars to represent. The original post says "goodness" (of fit), but goodness of fit does not necessarily represent error for the value of the plot, but rather "error" for the assumption you are making about the data used to construct the bar. You could give error bars that give relative fit to one another, but would not make sense to be plotted on the graph. $\endgroup$ – Underminer Dec 18 '13 at 21:34
  • $\begingroup$ I would suggest visual inspection of the QQ plot (for each year) to determine appropriateness of the Poisson assumption, and then report $\lambda$ values. $\endgroup$ – Underminer Dec 18 '13 at 21:44
1
$\begingroup$

One way to go is by using the "deviance residual" from glms. This is a likelihood ratio for each count. You can easily get these from the glm () function in R. The glm you would fit is one with no intercept, and dummy variables for each year. However with a million observstions this may run into memory problems.

You can calculate the deviance residuals for the poisson model as $$ d_{ti}^2=2 (\hat{\lambda}_t-n_{ti})+2n_{ti}\log\left (\frac {n_{ti}}{\hat {\lambda}_t}\right) $$ I have use "t" to denote the year ("time"), "i" to denote the observations, $ n_{ti} $ is the "count" variable, and $\hat {\lambda}_t $ is the MLE for year "t". This is approximately equal to the usual "pearson chi-square" residual when lambda is large, but has more "robust" behaviour for small lambda due to the logarithm. The standard pearson residual is $r_{ti}=\frac {(\hat{\lambda}_t-n_{ti})^2}{\hat{\lambda}_t} $ which diverges like $\frac {1}{\hat{\lambda}_t}$ - much faster than logarithm.

Now if your Poisson model is adequate you should roughly have $\sum_{ti} d_{ti}^2\approx df $ where "df" is the number of observed counts minus the number of lambdas you fitted.

Another way you can test the poisson model is to fit a negative binomial model, and see how close the "concentration" parameter is to infinity.

$\endgroup$
  • $\begingroup$ Thanks! There is a lot to look over here but it looks very helpful. Do you have a reference or a link for the equation of the "deviance residuals for the poisson model" or is it standard knowledge in a stat textbook? $\endgroup$ – Hooked Dec 19 '13 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.